François De Vleeschouwer
email: fdevleeschouwer@ulg.ac.be

Etienne Juvigné

Virginie Renson

José Antonio Naranjo

MINERAL CHEMISTRY OF THE LLAIMA PUMICE, SOUTHERN CHILE: EVIDENCE FOR MAGMA MIXING

Keywords: Chile, Lake District, Holocene, volcano, pumice, olivine

is paper details the mineralogy and the chemistry of the Llaima pumice deposit from in southern Chile. is andesitic to dacitic pumice was erupted from Llaima Volcano. is eruption took place in the early Holocene through at least two pulses that are well recorded in deposits around the Icalma-Galletué area. In 2001, the Llaima pumice was sampled from six outcrops. e mineralogy has been studied under polarizing microscope and electron microprobe (EPMA). e Llaima pumice contains phenocrystals of orthopyroxene, clinopyroxene, olivine and titano-magnetite. e olivine is characterized by two chemically distinguishable populations (i.e. Fo 32 to Fo 35 and Fo 76 to Fo 82 ). e respective chemical compositions of these two distinct olivines emphazise that the eruption was triggered by magma mixing.

RÉSUMÉ. Chimie des minéraux de la ponce Llaima, Sud du Chili : Evidence d'un mélange magmatique. Ce travail détaille les compositions minéralogique et chimique de la ponce Llaima. Celle-ci est andésitique à dacitique, et résulte d'une éruption du volcan Llaima [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF]. Cette éruption s'est eff ectuée en deux temps durant la première moitié de l'Holocène. En 2001, la ponce Llaima a été échantillonnée dans six coupes. Elle est étudiée au microscope polarisant et à la microsonde électronique. La ponce Llaima contient des phénocristaux d'orthopyroxène, de clinopyroxène, d'olivine et de titano-magnétite. Les deux populations d'olivine, chimiquement distinguables (i.e. Fo 32 à Fo 35 et Fo 76 à Fo 82 ) attestent que l'éruption fut déclenchée par un mélange de magma.

Mots-clés : Chili, Région des Lacs, Holocène, volcan, ponce, olivine in other environment (i.e. "dry" outcrops). Both lake catchments are covered by 4 -6 m of Holocene sandy to silty deposits in which two major pumice layers are interbedded (Pl. 1A). ese two pumice layers have been previously described. [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF] demonstrated that the lower pumice layer (Llaima pumice) was generated by Llaima Volcano (38° 42'S -71°44'W) at about 8830 ± 80 a BP (9900 ± 300 a cal BP) during a plinian eruption of its second eruptive phase called phase Llaima II. During this eruptive phase, a column of 4 km 3 of pumice reached around 40 km high and was dispered by wind towards the south-east. is explosive stage ended with the collapse of the plinian column generating a 6-8 m thick base surge deposit. e Upper pumice layer (Alpehué pumice) took place during the last vigorous plinian eruption (column of 44 km high, 7.5 km 3 of airfall pumice and 0.4 km 3 of pyrclastic fl ow) of Caldeira de Sollipulli (38° 58'S -71°31'W) causing the formation of crater Alpehué at 2900 a BP (3020 ± 60 a cal BP) [START_REF] Gilbert | Non-explosive, constructional evolution of the ice fi lled caldera at Volcán Sollipulli, Chile[END_REF][START_REF] Naranjo | Volcanismo explosivo reciente en la caldera del volcán Sollipulli, Andes del Sur (39°S)[END_REF]. In this short contribution, we document mineralogical data for

Introduction

e studied area is located in the central part of the Southern Volcanic Zone (SVZ) (López-Escóbar et al., 1977;orpe & Francis, 1979;López-Escóbar, 1984). e SVZ rocks are dominated by basalts with subsidiary basaltic andesites and dacites rich in orthopyroxene, clinopyroxene and Fe-rich olivine (López-Escóbar, op. cit.;[START_REF] Futa | Sr and Nd isotopic and trace elements compositions of Quaternary volcanic centres of the Southern Andes[END_REF]. e occurrence of Fe-rich olivine is sometimes interpreted as resulting from partially assimilated inclusions of basaltic material [START_REF] Gerlach | Recent volcanism in the Puyehue-Cordon Caulle Region, Southern Andes, Chile (40.5°S): Petrogenesis of evolved lavas[END_REF]. Lakes Icalma and Galletué are two glacial lakes (~4 km long) located in la Región de la Araucaña (Chile, 37°0' -39°30' S; Fig. 1). ey are surrounded westwards by a volcanic fi eld that has been active throughout the Quaternary (e.g. [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF][START_REF] Naranjo | Volcanismo explosivo reciente en la caldera del volcán Sollipulli, Andes del Sur (39°S)[END_REF][START_REF] Gonzalezferran | Volcanes de Chile[END_REF]. Hence a lot of tephra falls were transported by dominating Westerlies winds.

e eruptions that emplaced these tephra falls are well recorded in peat bogs or lake sediments which are good sedimentary traps, but only the thickest ones are recorded the Llaima pumice. is provides some new insights into the mechanism of the eruption.

Material and method

Sample collection

During the austral summer of 2002, a fi eld campaign was conducted in the Icalma and Galletué watersheds. Five outcrops (e.g. Pl. 1.A) around the two lakes were sampled from the diff erent watersheds (Fig. 1C). An additional outcrop was sampled westwards of the Icalma catchment (Fig. 1C, location CO2). Samples (around 1 dm 3 of bulk material) in each tephra bed were taken in respect of signifi cant macroscopical changes of texture, grain size or colour occurs. us, the number and position of samples vary in the diff erent studied outcrops. Because of the dry summer climate, land reclamation (particularly for rearing), and high drainage (enhanced by the permeable nature of the volcanic layers), only the two lacustrine infi lling peat bogs that subsist downstream of Lake Galletué (Fig. 1) were cored. Coring was carried out with a Russian type steel corer [START_REF] Belokopytkov | Giktorf 's peat borers[END_REF][START_REF] Jowsey | An improved peat sampler[END_REF]. ese peat bogs are relatively heterogeneous and display alternating humid and dry zones. us, each peat bog was studied and precored in order to determine the deepest site that could record the longest peat deposit.

Sample preparation and optical analyses

Outcrop samples containing particles coarser than 420 µm were roughly crushed and separated in four equal parts. An aliquot, made up of a small part of each quarter, was then re-crushed and sieved to isolate the grain size fraction between 105 and 420 µm. Samples containing only particles smaller than 420 µm were sieved without preliminary crushing to isolate the grain size fraction between 105 and 420 µm. In the cores from peat-bogs, the Alpehué pumice was visible to the naked eye (Pl. 1B). e Llaima pumice was not present, because these bogs are located outside the dispersion lobe. Each tephra bed was washed and treated in the same way as outcrop samples. In all cases, smear slides were mounted as follows:

(1) A small part of each 105 -420 µm fraction was mounted on a slide in Canada balsam and all the grains were determined in terms of their mineralogical nature. (2) Heavy minerals were separated by centrifuging 0.5 cm 3 of the 105-420 µm fraction in test tubes containing 10 ml bromoform (density = 2800 kg.m -3 ). Each tube's bottom was frozen in liquid nitrogen, allowing top the liquid containing the light fraction to be entirely discarded. e bromoform ice, containing the dense fraction was thawed, fi ltered and dried. Dense minerals were mounted on a slide. About 300 grains on each slides were identifi ed and counted under polarising microscope using the 'Van Harten rubber technique' [START_REF] Van Harten | On the estimation of relative grain frequencies in heavy mineral slides[END_REF].

Chemical analyses

Monomineral spots from both pumice layers were isolated by hand-picking mineral grains in the 105-420 µm fraction. All samples were mounted in PETREPOXY ® resin.

Major elements analyses were performed with a Cameca Camebax SX 50 electron microprobe from the C.I.A. of Louvain-la-Neuve University (Belgium). Analytical conditions are 15 kV, 10 nA current (analytical time: 10 seconds; bulk standard analyses used for calibration available from the author).

14 C age dating

One sample of a charcoal-rich humus layer outcropping 15 cm above the Llaima pumice at location CO2 was collected for 14 C-dating. Dating was performed by A.M.S. (1.5 SDH-Pelletron Model "Compact Carbon AMS") at the Poznań Radiocarbon Laboratory (Poland). e sample was acid-base digested to remove any humic acid and fulvic acid, or any secondary carbonate deposited in the sediment's organic fraction. It was then burned to obtain purifi ed CO 2 for introduction into the spectrometer's carbon cathode. Calibrations were made according to [START_REF] Stuiver | INTCAL98 Radiocarbon Age Calibration, 24,000-0 cal BP[END_REF].

Results

14 C age dating

e charcoal-rich humus was dated at 9,030 ± 45 C yr BP (Poz-252), corresponding to bracketed calibrated ages of 10,180 ± 70 yr cal BP. Since our sample was taken just above the tephra, that result should be considered as a minimal age.

Llaima pumice

In the Icalma and Galletué catchments, the Llaima pumice is up to 100 cm thick. It is a yellowish to brownish coarse grained pumice layer with variable amounts of lithic fragments (1-39% in CO2, 3-17% in OC1, 0-57% in OC2; De Vleeschouwer, 2002). e pumice mineralogy is characterised by colourless glass shards, xenomorphic plagioclase, automorphic orthopyroxene and clinopyroxene, apatite, opaque minerals and in addition, a high amount of hypautomorphic to automorphic olivine (Pl. 1B-C). All these crystals are phenocrysts and their concentrations vary from 2 to 13% in respect to the diff erent outcrops and samples. e olivine population can be divided into two groups based on colour: olivine phenocrysts are either colourless or honey-yellow. Glass shards are optically yellowish to brown, refl ecting a relatively advanced degradation stage. ese colours are darker at the deposit's base and top. is could be explained by a stagnation of percolating water at the deposit base, degrading the pumice more intensively. 0,01 0,00 0,01 0,00 0,02 0,00 0,01 0,00

Al 3+ (IV) 0,03 0,03 0,01 0,02 0,12 0,04 0,03 0,01

Al 3+ (VI) 0,01 0,01 0,01 0,01 0,00 0,00 0,01 0,01 Fe 3+ 0,01 0,02 0,01 0,01 0,09 0,03 0,03 0,02 Fe 2+ 0,70 0,24 0,78 0,05 0,19 0,08 0,35 0,01 Mn 2+ 0,01 0,01 0,01 0,01

Mg 2+

1,20 0,20 1,14 0,05 0,82 0,05 0,75 0,01 Ca 2+ 0,06 0,01 0,05 0,00 0,85 0,04 0,83 0,01 Na + 0,00 0,00 0,00 0,00 0,02 0,00 0,02 0,00 K + 0,00 0,00 0,00 0,00 0,00 0,00 0,00 0,00 Cr 3+ 0,00 0,00 0,00 0,00 Table 1. Mean values and cationic proportions of chemical analyses for glass shards and minerals from the Llaima pumice (the Alpehué pumice is given for comparison).

e degradation of the deposit top can be accounted for meteoritic weathering of the tephra. A fresh outbreak of lithics unaccompanied by any erosion structure is observed in the pumice at three locations. In fact, lithics concentration increases from 1 to 32%, 115 cm above the base of the pumice at CO2; from 3 to 11%, 10 cm above the base of the pumice at OC1, and from 0 to 57%, 25 cm above the base of the pumice at OC2 (Fig. 2). Since the increase in lithic content is abrupt, it may refl ect a structural change in the source vent as the explosions blowed off older rock fragments from the partially collapsed crater and upper conduit. In such context, Llaima pumice can be considered as deposited in two pulses during a single plinian eruption phase. Llaima pumice has an andesitic to dacitic composition (60-65% SiO 2 ), as already shown by [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF]. Glass shards are rhyodacitic (70-74% SiO 2 ; Fig. 3), transparent minerals are made of andesine (Fig. 4A) and fl uorapatite. Neither K-feldspar nor quartz was detected. Orthopyroxene is hypersthene and clinopyroxenes are made of Ca-augite and diopside (Fig. 4B). Opaque minerals are titano-magnetites (Fig. 4C). e transparent olivine is made of Fo 76 to Fo 82 whereas the honey yellow one is made of Fo 32 to Fo 35 . Mean chemical analyses and cationic proportions are reported in Table 1.

Discussion

Age of the eruption

In spite of the fact that the age obtained in this work (10,180 ± 70 yr cal BP) is a minimal age, it is consistent with the age obtained by [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF]: 9900 ± 300 yr cal BP for a charcoal taken above the pumice layer. Both results place the tephra-fall around 10,000 yr cal BP. However, [START_REF] Mardones | Gestion de los recursos hidricos de la cuenca del Rio Bio-Bio y del area marina costera adyacente[END_REF] have dated a silty to clayey sample that was taken just below a pumice tephra layer interpreted as the equivalent of Llaima pumice in a core of a nearby peat bog (Guallelmapu, S 38°72' -W 71° 13', pumice depth between 95 and 105 cm, sample depth: 125 cm). It has given an age of 4580 ± 70 yr BP (5225 ± 275 yr cal BP). Provided the tephra correlation is accepted, such a diff erence between [START_REF] Mardones | Gestion de los recursos hidricos de la cuenca del Rio Bio-Bio y del area marina costera adyacente[END_REF] and our dating cannot be reconciled. Mineralogical, chemical and age consistencies between Naranjo & Moreno work (1991) and ours implies that the Mardones et al. (op. cit.) date is not useable, probably because of a young carbon contamination by water percolation or presence of rootlets through the pumice deposit.

Eruption process

e key founding of the present study is the occurrence of two diff erent olivine populations in the same tephra as well as olivine in rhyodacitic liquids. Olivine had already been optically identifi ed in the Llaima pumice by [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF]. Our investigation shows that these olivine populations are optically but also chemically diff erent. e Fe/Mg ratio for both olivine populations and for the Llaima pumice's glass shards is reported on Fig. 4D using the terminology [START_REF] Roeder | Olivine-Liquid Equilibrium[END_REF].

e F o 32-35 fi eld is overlapped by the glass shard's one. e Fo 76-82 is isolated. Both observations demonstrate that Fo 32-35 crystallised in equilibrium with the magma and that Fo 76-82 is in disequilibrium. is is consistent with the fact that even if olivine is more common in basic magmas, Fe-rich olivine is known to be a common mineral phase in the dacitic liquid from this region [START_REF] Futa | Sr and Nd isotopic and trace elements compositions of Quaternary volcanic centres of the Southern Andes[END_REF]López-Escóbar, 1984). Moreover, olivine with bimodal distributions (Fo 24-26 and Fo 74-87 ) have already been observed in andesites from Puyehue volcano [START_REF] Gerlach | Recent volcanism in the Puyehue-Cordon Caulle Region, Southern Andes, Chile (40.5°S): Petrogenesis of evolved lavas[END_REF]. ese andesites from Puyehue volcano are interpreted to be disequilibrium phenocrysts resulted from a mixing between a basaltic magma and an evolved magma to produce supercooled mixtures. e structural study of the Llaima pumice had already suggested a mixing with a more basic (i.e. hotter) magma at the bottom of the Llaima magmatic chamber [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF].

e presence of Fo 76-82 is consistent with this hypothesis. e eruption of the Llaima pumice was caused by mixing between a resting acid magma with a replenishment batch of more basic magma containing Fo 76-82 .

Summary and conclusion

Two major tephras blanketed the catchments of lakes Icalma and Galletué (southern Chile) during the Holocene. e lower one, the Llaima pumice [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF], is an olivine-bearing andesite to dacite with pyroxene and titano-magnetite. It erupted in two pulses of a single eruption phase, from Volcano Llaima during the early Holocene (10,180 ± 70 yr cal BP -this work, 9900 ± 300 yr cal BP - [START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF].

e eruption was triggered by magma mixing as suggested in the previous study of Naranjo & Moreno (op. cit.) and by the two populations of completely diff erent compositions of olivine, one of which did not crystallise in equilibrium with the magmatic liquid.

Figure 1 .

 1 Figure 1. Schematic map of studied area. Black dots: sampled sites, OC1, OC2, OC3, OC4, CO1, CO2: outcrops; PB1 and PB2: peat bogs. e Llaima volcano is 40 km westwards.

Figure 2 .

 2 Figure 2. Detailed logs and components for the Llaima Pumice sections CO2, OC1 and OC2. A fresh lithics outbreak is observed above the middle of the pumice deposit. Lithics are mainly composed of fragments of poorly diff erenciated magma (i.e. basaltic scoria s.l.) along with a small fraction of schists.

Figure 4 .

 4 Figure 4. Chemical analyses of minerals found in the Llaima Pumice reported in corresponding ternary diagrams. A. feldspars diagram. B. pyroxenes diagram according to Morimoto (1988). C. oxides diagram. D. Fe/Mg ratios diagram according to[START_REF] Roeder | Olivine-Liquid Equilibrium[END_REF]. Note that the Fo 76-82 fi eld is strongly isolated from the Llaima Pumice glass shards' fi eld.

Figure 3 .

 3 Figure 3. Llaima Pumice glass shards: chemical analyses reported in the TAS vs. SiO 2 diagram (Le Bas et al., 1986); Black diamonds: whole rock analyses from the Llaima pumice[START_REF] Naranjo | Actividad explosiva postglacial en el Volcán Llaima, Andes del Sur (38°45'S)[END_REF]

56 Cationic proportions calculated on a fi ve cations basis Cationic proportions calculated on a three cations basis

	Llaïma glass shards Mean (n=20) Llaïma Mean (n=2) Std dev Orthopyroxene Alpehué glass shards Std dev Alpehué Mean (n=18) Mean (n=5) Std dev Std dev	Clinopyroxene Std dev Llaïma Mean (n=26) Mean (n=5) Alpehué Std dev
	SiO 2 SiO 2 TiO 2 TiO 2 Al 2 O 3 Al 2 O 3 FeOt FeO	72,14 51,81 0,41 0,35 15,49 0,92 3,42 22,19	1,24 1,53 0,13 0,02 0,83 0,49 0,56 6,00	74,88 52,28 0,21 0,18 14,12 0,54 1,50 24,74	0,27 0,47 0,08 0,14 0,09 0,23 0,12 1,44	50,26 0,64 2,76 8,62	0,80 0,13 0,80 2,03	52,06 0,26 0,85 11,86		0,22 0,06 0,16 0,23
	MnO Fe 2 O 3 MgO MgO	0,60 0,29 0,45 21,32	0,10 0,40 0,16 4,45	0,26 0,24 20,04	0,02 0,02 1,18	0,28 14,68	0,17 0,91	13,28		0,23
	CaO CaO	2,31 1,50	0,34 0,40	1,10 1,34	0,05 0,06	21,05	1,04	20,58		0,24
	Na 2 O Na 2 O K 2 O K 2 O Total Cr 2 O 3	3,28 0,00 1,84 0,00 99,94 0,00	1,01 0,00 0,21 0,00 0,00	3,26 0,00 4,24 0,00 99,81	0,77 0,00 0,10 0,00	0,29 0,00 0,08	0,04 0,01 0,15	0,33 0,00		0,03 0,00
	Total	98,37	1,27	99,11	0,66	98,67	0,51	99,21		0,21
	Si 4+	1,97	Plagioclase Cationic proportions calculated on a four cations basis 0,03 1,99 0,02 1,89 0,03 Olivine (Fo 32-35 )	1,97	Olivine (Fo 76-82 ) 0,01
	Ti 4+	Llaïma Mean (n=12) Std dev	Alpehué Mean (n=10) Std dev		Llaïma Mean (n=19) Mean (n=24) Llaïma Std dev
	SiO 2	58,11	4,14	58,44	1,77	SiO 2	32,22	38,80	0,79
	TiO 2	0,00	0,00	0,00	0,00	TiO 2	0,00	0,01	0,03
	Al 2 O 3	25,32	2,50	25,07	1,05	Al 2 O 3	0,00	0,01	0,02
	FeO	0,32	0,13	0,29	0,07	FeO	50,61	18,99	2,93
	MnO	----	MnO	1,00	0,10	0,17
	MgO	0,02	0,02	0,01	0,01	MgO	14,41	41,09	2,22
	CaO	8,19	3,07	7,79	1,20	CaO	0,20	0,21	0,03
	Na 2 O	6,74	1,75	6,84	0,64	Na 2 O	0,00	0,00	0,01
	K 2 O	0,17	0,06	0,42	0,12	K 2 O	0,00	0,00	0,01
	Cr 2 O 3	----	Cr 2 O 3	0,00	0,00	0,00
	Total 0,Si 4+ 98,86 0,45 98,95 0,30 Total 98,45 99,21 2,63 0,16 2,65 0,07 Si 1,00 1,00 0,01
	Ti 4+	0,00	0,00	0,00	0,00	Ti	0,00	0,00	0,00
	Al 3+	1,37	0,16	1,35	0,07	Al 3+	0,00	0,00	0,00
	Fe 3+	0,01	0,01	0,01	0,00	Fe 3+	0,00	0,00	0,00
	Mn 2+	----	Fe 2+	1,31	0,41	0,07
	Mg 2+	0,00	0,00	0,00	0,00	Mn 2+	0,03	0,00	0,00
	Ca 2+	0,40	0,16	0,38	0,06	Mg 2+	0,66	1,58	0,07
	Na +	0,59	0,15	0,60	0,05	Ca 2+	0,01	0,01	0,00
	K +	0,00	0,00	0,02	0,01	Na +	0,00	0,00	0,00
	Cr 3+	----	K +	0,00	0,00	0,00
						Cr 3+	0,00	0,00	0,00
	Ab	59,16	15,19	59,86	5,47	Fo	32,42	79,29	3,52
	An	39,88	15,51	37,72	5,88
	Or	0,96	0,37	2,42	0,72

Acknowledgements. is work was funded as O.S.T.C. project EV/10/12B (coord. M. De Batist, R.C.M.G., University of Ghent, Belgium) and the Lefranc Foundation. Fondecyt Project No. 1960186 is also acknowledged. We thank N. Fagel and S. Bertrand (U.R.A.P. University of Liège, Belgium) for helping to improve the paper. e paper has greatly benefi ted from the reviews of G. Ernst and M. Pino.


Salli Moustafa

Ivan Dutka-Malen

Laurent Plagne

Angélique Ponçot

Pierre Ramet

Shared Memory Parallelism for 3D Cartesian Discrete Ordinates Solver

Keywords: Deterministic transport, SN, multicore processor, wavefront, TBB, SIMD

This paper describes the performance of DOMINO, a 3D Cartesian SN Solver that implements two nested levels of parallelism (multicore+SIMD) on shared memory computation nodes. DOMINO solves a 3D full core PWR k ef f problem involving 26 groups, 288 angular directions (S16), 46×10 6 spatial cells and 1×10 12 DoFs within 12 hours on a 32 core SMP node. This represents a sustained performance of 184 GFlops and 31% of the SMP node peak perfomance.

I. The DOMINO COCAGNE SN Solver

COCAGNE is the name of the new nuclear core simulation system developed at EDF R&D. Recently, a shared-memory parallel solver named DOMINO based on a discrete ordinates method has been introduced into the COCAGNE system.

DOMINO Numerical Schemes

DOMINO implements the Discrete Ordinates Method for the stationary neutron transport in a multi-dimensional Cartesian geometry. The quadrature used is the level symmetric one. In three-dimensions, this quadrature leads to N(N + 2) discrete directions where N stands for the method's order. The space discretization corresponds to the Diamond Differencing scheme DD0. The DD0 element has 1 moment and 3 mesh-edge incoming fluxes per cell. The spatial problem is solved through the well-known sweep algorithm. We use a classical Diffusion Synthetic Acceleration (DSA) method. The sweep operation is the most computationally intensive portion of DOMINO. For each incoming direction of the angular quadrature, the angular and volumic flux Degrees of Freedom (DoFs) of each spatial mesh cell must be updated. For simplicity's sake, we illustrate the parallel sweep on a 2D example with a DD0 spatial discretization scheme. Let {c i j } be the n x × n y cells of a 2D Cartesian spatial mesh. Let us consider Ω d an angular direction incoming from the bottom left corner. Let { d ψ x i j , d ψ y i j , φ i j } be the corresponding incoming angular and volumic flux DoFs. In the sweep operation one has to treat each cell c i j by updating the volumic φ i j and the angular outcoming d ψ x i+1 j and d ψ y i j+1 DoFs that depend on the three d ψ x i j , d ψ y i j and φ i j input values. This implies an ordering constraint for the sweep operation: a cell c i j can only be treated if the following two conditions are fulfilled:

Shared Memory Parallel Implementation of DOMINO

c i-1, j is already treated or i = 0, c i, j-1 is already treated or j = 0.
Obviously c 00 is the first cell that can be treated but the second can be either c 10 or c 01 . . . or both can be treated in parallel. We rely on the Intel TBB primitive parallel_do [START_REF] Robison | Optimization via reflection on work stealing in TBB[END_REF] that enables a dynamic scheduling of the parallel tasks to implement the sweep. This parallel function allows a set of threads to execute the tasks from a task list which is dynamically updated. In the beginning, the task list contains the cell c 00 . One of the running threads treats this cell and updates the task list to {c 10 , c 01 } and so on. In order to reduce the overhead due to the thread scheduling, the cells are not treated individually but they are packed into groups of cells, called macrocells in the DOMINO implementation. All the cores of the supercomputing nodes are used in this parallel wavefront algorithm. Figure 1 illustrates this parallel sweep strategy.

During each sweep, the angular directions that belong to the same octant are computed via SIMD instructions allowing to perform simultaneously two double precision or four simple precision floating point operations.

In this paper the computations are made with one of the com-puting nodes of the EDF R&D Ivanoé supercomputer [START_REF] Robison | Optimization via reflection on work stealing in TBB[END_REF] . This 1024 GB node is a 4-socket Intel X7560 processor (2.26GHz), 8 cores per socket, and 4-way SIMD units totalizing 128 parallel FPUs that can be used concurrently. As a consequence, one can expect a 128 fold speed-up for an ideal parallel implementation of DOMINO compared to a sequential non-SIMD run time.

II. Scalability and Efficiency

We analyze the performance of the sweep operation for a simple 2-group k ef f computation with 480 3 spatial cells with different angular quadrature orders. Figure 2 and Table 1 sumarize the performance of the 2-level parallel implementation (mul-ticore+SIMD) of DOMINO. Since the arithmetic intensity of the computation increases with the direction number, the sweep performance improves with the angular quadrature order. S2 and S4 timings are exactly the same because we use SIMD units to treat 4 angular directions that belong to the same octant at the same time. GFlops value is estimated by dividing the floating point operation number by the completion time:

GFlops = 25 × N cells × N dir Time in nanosecond ,
where the factor 25 accounts for the 20 add/mult operations and 1 division per cell (we count 5 flops for the division).

In the S16 case (N dir = 288), the sweep performance reaches 216 GFlops which corresponds to 37% of the peak performance given by (Single Precision): Peak = ncores × frequency × SIMD width × SIMD units = 32 × 2.26 × 4 × 2 = 578 GFlops.

III. DOMINO/MCNP Comparison

The benchmark used for calculations is described in reference. [START_REF] Courau | Hybrid Parallel Code Acceleration Methods in Full-Core Reactor Physics Calculations[END_REF] It corresponds to a simplified 3D PWR first core loaded with 3 different types of fuel assemblies characterized by a specific Uranium-235 enrichment (low, medium and highly enriched uranium). The calculation mesh is based on the cell pin mesh (289×289) and along the z-axis, the height is subdived into 70 meshes for the 26-group calculation. For spatially converged results, meshes are refined by 2 in directions x, y and z. The angular discretization corresponds to a S16 level symmetric angular quadrature, which results in 288 directions in 3D. To simplify the comparisons to MCNP results, the 3D 26-group fluxes are integrated over energy and space. Practically, a group collapsing is performed from 26 to 2 energy groups with a boundary fixed at 0.625 eV, and a spatial integration over each pin-cell along the z-axis is made. Comparisons between DOMINO and MCNP presented in Table 2 are computed in the form:

∆k ef f = k DOMINO ef f -k MCNP ef f k MCNP ef f × 10 5 in pcm, | δφ g |= max {p i } | φ DOMINO g (p i ) -φ MCNP g (p i ) φ MCNP g (p i ) |,
where p i stands for the pin cell index. This SN k ef f calculation involved 1.05 × 10 12 DoFs obtained within 126 power iterations for a 10 -5 stopping criterium. The DOMINO implementation reached 184.3 GFlops which corresponds to 31% of the SMP node peak performance.

Encouraged by these node performances, we are currently working on extending DOMINO to distributed memory machines.

Figure 1 :

 1 Figure 1: Wave front algorithm applied to the sweep across the mesh macrocells. The incoming fluxes come from the bottom left corner. The treated macrocells are pink colored while the white macrocells are still waiting. Only blue dashed border macrocells are ready for treatment and their respective indexes are stored in a task list from which a set of threads concurrently pick their work item. This list is dynamically updated each time a new macrocell has been treated.

Figure 2 :

 2 Figure 2: Time per sweep over all angular directions as a function of the core number (mesh: 480 × 480 × 480).

Table 1 :

 1 Time per sweep over all angular directions as a function of the core number (mesh: 480 × 480 × 480).

	1 Ranked 149 in Nov 2012 TOP 500 List

Table 2 :

 2 Comparison to MCNP/DOMINO: S16, 26-group, 578 × 578 × 140 spatial cells. 2.26GHz SMP node with 32 cores


Salli Moustafa

Ivan Dutka-Malen

Laurent Plagne

Angélique Ponçot

Pierre Ramet

Shared Memory Parallelism for 3D Cartesian Discrete Ordinates Solver

Keywords: Deterministic transport, SN, multicore processor, wavefront, TBB, SIMD

This paper describes the design and the performance of DOMINO, a 3D Cartesian SN solver that implements two nested levels of parallelism (multicore+SIMD) on shared memory computation nodes. DOMINO is written in C++,a multi-paradigm programming language that enables the use of powerful and generic parallel programming tools such as Intel TBB and Eigen. These two libraries allow us to combine multi-thread parallelism with vector operations in an efficient and yet portable way. As a result, DOMINO can exploit the full power of modern multi-core processors and is able to tackle very large simulations, that usually require large HPC clusters, using a single computing node. For example, DOMINO solves a 3D full core PWR eigenvalue problem involving 26 energy groups, 288 angular directions (S16), 46 ⇥ 10 6 spatial cells and 1 ⇥ 10 12 DoFs within 11 hours on a single 32-core SMP node. This represents a sustained performance of 235 GFlops and 40.74% of the SMP node peak perfomance for the DOMINO sweep implementation. The very high Flops/Watt ratio of DOMINO makes it a very interesting building block for a future many-nodes nuclear simulation tool.

I. INTRODUCTION Industrial Context

As part of its activity, EDF R&D is developing a new nuclear core simulation code named COCAGNE. This code relies on DIABOLO, a Simplified PN (SPN) method to compute the neutron flux inside the core for eigenvalue calculations. [START_REF] Kirschenmann | Parallel SPN on Multi-Core CPUS and Many-Core GPUS[END_REF] In order to assess the accuracy of SPN results, a 3D Cartesian model of PWR nuclear cores has been designed and a reference neutron flux inside this core has been computed with the MCNP5 Monte Carlo transport code. [START_REF] Courau | Hybrid Parallel Code Acceleration Methods in Full-Core Reactor Physics Calculations[END_REF] This kind of 3D whole core probabilistic evaluation of the flux is computationally very demanding. An efficient deterministic approach is therefore required to reduce the computation effort dedicated to reference simulations. In this paper we introduce DOMINO, a new shared memory parallel 3D Cartesian SN solver specialized for PWR core reactivity computations which is fully integrated in the COCAGNE system.

A First Step Toward Efficient Use of the HPC Clusters

Our aim is to build a massively parallel deterministic neutron transport solver such as Denovo [START_REF] Davidson | Massively Parallel, Three-Dimensional Transport Solutions for the k-Eigenvalue Problem[END_REF] which is under active development at ORNL 1 , allowing simulations on large distributed memory supercomputers. To achieve this goal, we first focus on designing an efficient implementation for one supercomputing node. Modern supercomputers are generally built upon hierarchical architectures that can be characterized by their number n nodes of computing nodes, the number n sockets of processors in each node, the number n cores of cores in each processor, and [START_REF] Kirschenmann | Parallel SPN on Multi-Core CPUS and Many-Core GPUS[END_REF] Oak Ridge National Laboratory finally the number n simd of Single Instruction on Multiple Data (SIMD) units in each core. Although we plan to use these four levels of parallelism in the future, the DOMINO implementation is, at the present time, limited to the inner three levels (multi-socket computer nodes) which can be addressed with shared memory parallel programming paradigms. Obviously, an efficient implementation for supercomputing nodes is required to attain an efficient use of many-node HPC facilities. This hierarchical path for building massively parallel solvers on modern supercomputers is well illustrated in ref [START_REF] Petrini | Multicore surprises: Lessons learned from optimizing Sweep3D on the Cell Broadband Engine[END_REF] where the authors give algorithmic and implementation strategies for optimizing the Sweep3D [START_REF] Koch | Solution of the firstorder form of the 3-D discrete ordinates equation on a massively parallel processor[END_REF] kernel on the IBM Cell BE processor.

Proper Tools For Multi-Core Processors Programming DOMINO is written in C++, a multi-paradigm programming language that enables the use powerful and generic parallel programming tools such as Intel TBB [START_REF]Intel TBB: a C++ template library for task parallelism[END_REF] and Eigen. [START_REF]Eigen: a C++ template library for linear algebra[END_REF] These two libraries allow us to combine multi-thread parallelism with vector operations in an efficient and yet portable way. As a result, DOMINO can exploit the full power of modern multicore processors and is able to tackle very large simulations, that usually require large HPC clusters, using a single computing node. For example, DOMINO solves a 3D full core PWR k eff problem involving 26 energy groups, 288 angular directions (S 16 ), 46 ⇥ 10 6 spatial cells and 1 ⇥ 10 12 DoFs within 11 hours on a single 32-core SMP node. This represents a sustained performance of 235 GFlops and 40.74% of the SMP node peak performance for the DOMINO sweep implementation. The very high Flops/Watt ratio of DOMINO makes it a very interesting building block for a future distributed memory nuclear simulation tool.

Paper Outline

In this paper, we describe the design of DOMINO and assess its performance. Section II briefly introduces the discrete ordinates numerical scheme and algorithms used in DOMINO. Section III describes the two nested levels of parallelism (mul-ticore+SIMD) of DOMINO. Then section IV analyses the parallel scalability of the sweep kernel that represents the computationally most demanding part of code. In section V we present 3D full core PWR k eff computations carried out with DOMINO. Both k eff and flux accuracies are assessed by comparison with Monte-Carlo MCNP computations. Then we compare the performances of DOMINO with those of PENTRAN [START_REF] Courau | 3D Neutron Transport and HPC: A PWR Full Core Calculation Using PENTRAN SN Code and IBM BLUEGENE/P Computers[END_REF] and Denovo, [START_REF] Evans | Threedimensional full core power calculations for pressurized water reactors[END_REF] two massively parallel deterministic transport solvers. Section VI concludes the paper by describing our ongoing work devoted to extending DOMINO to distributed memory machines.

II. THE DOMINO COCAGNE SN SOLVER

COCAGNE is the name of the new nuclear core simulation system developed at EDF R&D. The COCAGNE system mostly relies on DIABOLO, an approximate solver based on the Simplified PN equations (SPN) [START_REF] Kirschenmann | Parallel SPN on Multi-Core CPUS and Many-Core GPUS[END_REF][START_REF] Courau | DOMINO: A Fast 3D Cartesian Discrete Ordiantes Solver for Reference PWR Simulations and SPN Validations[END_REF] for its industrial simulations. Recently, an additional reference solver named DOMINO and based on a discrete ordinates method has been introduced into the COCAGNE system. These two solvers share a common interface and a similar design within the COCAGNE system. In particular, both DOMINO and DIABOLO are built upon Legolas++,aC ++ library dedicated to solving multilevel blocked linear algebra systems. [START_REF] Plagne | Generic Programming for Deterministic Neutron Transport Codes[END_REF] This design similarity facilitates the interoperability between the two solvers. This section starts with a description of the numerical schemes used in DOMINO and then gives an overview of its algorithmic structure.

DOMINO Numerical Schemes

DOMINO (Discrete Ordinates Method In NeutrOnics) implements the Discrete Ordinates Method for the stationary neutron transport equation (1) in multi-dimensional Cartesian geometries. This equation depends on energy (E), the angular directions ( Ω) and the spatial position (r) of the particles, defining respectively the velocity, propagation direction and the spatial localization of particles.

Ω • -! r (r, E, Ω) +Σ t (r, E) (r, E, Ω) = Z 1 0 dE 0 Z S 2 d Ω0 Σ s (r, E 0 ! E, Ω0 • Ω) (r, E 0 , Ω0 ) + 1 k χ(E) 4⇡ Z 1 0 dE 0 Z S 2 d Ω0 ⌫Σ f (r, E 0 ) (r, E 0 , Ω0 ). ( 1
)
This is an eigenvalue problem with k representing the multiplication of neutrons in successive fission generations; it can be rewritten as:

H (r, E, Ω) = 1 k F (r, E, Ω), (2)
where H and F represent transport and fission operators. The problem (2) is solved using an inverse power algorithm, which leads to the computation of the neutron flux and k, by iterating on the fission term:

H n+1 = 1 k n F n , k n+1 = k n < F n+1 , F n+1 > < F n+1 , F n > .
Each power iteration solves a multi-group problem by using the Gauss-Seidel algorithm, which is the same as solving G one-group space-angle problems:

H gg g (r, Ω) = - X g 0 ,g H gg 0 g 0 (r, Ω) + S g (r), (3)
where S g (r) represents the fission sources for the group g.W e refer to these one-group problems as monogroup equations.

In the S N method, each of these equations is discretized on a finite number of angular directions defined by the quadrature formula used. DOMINO supports both Level Symmetric and Gauss-Legendre quadrature formulas, which lead respectively to N(N + 2) and m ⇥ n angular directions; N stands for the Level Symmetric quadrature formula order, m and n represent respectively the number of azimutal and polar directions used in the Gauss-Legendre quadrature formula. For both quadrature formulas, each angular direction is associated to a weight w i for integral calculation on the unit sphere S 2 . Hence, we have the following transport equations coupled by the scattering term:

L n+1 = B -Rφ n ,φ n+1 (r) = ndir X j=1 w j n+1 (r, Ω j ),
where L is the spatial stream matrix, composed of the diagonal of H, R = H -L is the spatial scattering matrix, and B = S the source term; ndir is the total number of angular directions depending on the quadrature formula used.

The general structure of the DOMINO solver is summarized by Algorithm 1. The first level of iterations solves the eigenvalue problem by an inverse power algorithm with a Chebychev acceleration. Then a Gauss-Seidel algorithm is used to solve the multigroup problem coming from the discretization of the energetic variable. Inside each multigroup iteration, scattering iterations with a Diffusion Synthetic Acceleration, [START_REF] Larsen | Unconditionally stable diffusion synthetic acceleration methods for the slab geometry discrete ordinates equations[END_REF] solve the spatial problem by the well-known sweep algorithm.

The space discretization corresponds to the Diamond Differencing scheme (DD). [START_REF] Martin | A three-dimensional high-order diamond differencing discretization with a consistent acceleration scheme[END_REF] Up to now, only the order 0 has been implemented. In three dimension, the DD0 element has 1 moment and 3 mesh-edge incoming fluxes per cell as indicated in Figure 1. The aim of this algorithm is to solve the space-angle problem, by inverting the streaming operator of the monogroup transport equation: lines 0-0 of Algorithm 1. The volumic flux computation inside the cell c (line 0) needs to know incoming data for this cell: incoming angular fluxes, L , B , F , the total cross section Σ t and the source term S . Outgoing angular flux are then updated on lines 0, 0 and 0. At this point, we can observe that once a cell has been swept for a given octant, incoming

Q ext = vS [g] + X g 0 ,g Σ g 0 !g s .φ g 0 ; B Scattering iterations while Non convergence do Q = Q ext +Σ g!g s φ g ; 0 Ωk . r k +Σ k = Q 8 Ωk 2 S N ; 0 φ g = X Ω k ! k k ; C = X g ⌫Σ f,
B c = (i, j, k) 0 forall the d 2 Directions[o] do B d = (⌫, ⌘, ⇠) 0 ✏ x = 2⌫ ∆x ; ✏ y = 2⌘ ∆y ; ✏ z = 2⇠ ∆z ; 0 [o][c][d] = ✏ x L +✏ y B +✏ z F +S ✏ x +✏ y +✏ z +Σ t ; 0 R [o][c][d] = 2 [o][c][d] -L [o][c][d]; 0 T [o][c][d] = 2 [o][c][d] -B [o][c][d]; 0 BF [o][c][d] = 2 [o][c][d] -F [o][c][d]; 0 φ[k][ j][i] = φ[k][ j][i] + [o][c][d] ⇤ ![d];
angular flux are not longer used. This property allows an optimization on the memory footprint of the solver: incoming and outgoing angular fluxes are stored on the same memory location, which dramatically increases the arithmetic intensity of the code. This feature is the key factor that allows obtaining a code that uses the full potential of modern multi-core architectures (see section 3). Finally we add the contribution of the volumic flux to the scalar flux on line 0, using the weight associated to the direction

d, ![d].
Let us count the total number of arithmetic operations per angular direction and per spatial cell in this sweep algorithm. We have 20 add/mult operations (the quantities [START_REF] Courau | Hybrid Parallel Code Acceleration Methods in Full-Core Reactor Physics Calculations[END_REF] ∆u , u = x, y, z can be computed once per spatial cell) each corresponding to 1 (flop) and 1 floating point division (line 0). The question how many flops to count for one division operation is a tricky one as the answer depends on the target architecture. In Ref, [START_REF] Vladimirov | Arithmetics on Intel's Sandy Bridge and Westmere CPUs: not all FLOPS are created equal[END_REF] authors show that the Sandy Bridge microarchitecture gives no performance gain for the division compared to the Nehalem microarchitecture, that is to say it costs exactly the same time to perform 8 packed single precision floating point divisions using AVX as 4 packed single precision floating point divisions using SSE. For our studies, we choose to count 5 flops per floating point division; this leads to a total of 25 flops for the sweep operation, determining the arithmetic intensity.

The Critical Arithmetic Intensity Issue

During the last decade, successive supercomputer node generations brought a regular and impressive improvement of their peak performance. Since their operating frequency remained almost unchanged, the multi-socket processors peak performance resulted mainly from their parallelism. In the context of transport simulations, the parallel computing power of a processor is proportional to the number n FPU of its parallel floating points units (FPU):

n FPU = n sockets ⇥ n cores ⇥ n simd .
Surprisingly enough, the computer bandwidth that measures the maximal data flow between the computer RAM and the FPUs did not increase as fast as the peak performance. The consequence of the broadening gap between the node bandwidth and its peak performance is to put a dramatic emphasis on the arithmetic intensity of the algorithms. If the arithmetic intensity i of a given algorithm, defined by: i = Number of floating points operations Number of RAM access (Read+Write) , is lower than a critical value i c , then its performance does not depend on the computational power of a target processor but mainly on the system memory bandwidth: the algorithm is then said to be memory bound. The critical arithmetic intensity of processors increases with each generation causing an ever larger fraction of algorithms to be memory bound. As a consequence of this processor evolution, the whole design of DOMINO is aimed at maximizing the arithmetic intensity in order to exploit the full parallel power of multicore processors.

III. TWO LEVELS OF PARALLELISM INSIDE THE DOMINO SOLVER

A shared memory supercomputing node can be seen as a multicore processor. To achieve a full utilization of its computing resources, we need to use both the multicore parallelism and the vector level parallelism also known as the Single Instruction Multiple Data model (SIMD). This latter level is a real bottleneck for performance issues, and its implementation constraints, involving hardware specifications of the chip, give guidelines to the algorithmic design of an application.

The Multicore Parallelism

The sweep operation is the most computationally intensive portion of DOMINO. For each incoming direction of the angular quadrature, the angular and volumic flux Degrees of Freedom (DoFs) of each spatial mesh cell must be updated. For simplicity's sake, we base the following sweep parallel description on a 2D example with a DD0 spatial discretization scheme. Let {c ij } be the n x ⇥ n y cells of a 2D Cartesian spatial mesh. Let us consider Ω d an angular direction coming from the bottom left corner. Let { d x ij , d y ij ,φ ij } be the corresponding incoming angular and volumic flux DoFs. In the sweep operation one has to process each cell c ij by updating the volumic φ ij and the angular outgoing d x i+1 j and d y ij+1 DoFs that depend on the three d x ij , d y ij and φ ij input values. This implies an ordering constraint for the sweep operation: a cell c ij can only be processed if the following two conditions are fulfilled: c i-1, j has been already processed or i = 0, c i, j-1 has been already processed or j = 0.

Obviously c 00 is the first cell that can be processed but the second can be either c 10 or c 01 . . . or both can be processed in parallel. We rely on the Intel TBB primitive parallel_do [START_REF] Robison | Optimization via reflection on work stealing in TBB[END_REF] that enables a dynamic scheduling of the parallel tasks to implement the sweep. This parallel function allows a pool of threads to execute the tasks from a task list which is dynamically updated. In the beginning, the task list contains the cell c 00 . One of the running threads processes this cell and updates the task list to {c 10 , c 01 } and so on. In order to reduce the overhead due to the thread scheduling, the cells are not processed individually but they are packed into groups of cells, called macrocells in the DOMINO implementation. The choice of the macrocell size depends of the available cache memory. In order to benefit from data reuse, data must be kept as long as possible in cache. If the macrocell size is too small, overheads due to thread scheduling deteriorate performances; the same consequence occurs when it is too big because of cache misses. For our studies, we perform an heuristic for the target architecture and determine the best size by trying several sizes. Figure 3 illustrates this parallel sweep strategy. All the n sockets ⇥ n cores cores of the supercomputing node are used in this parallel wavefront algorithm.

When we are using the vacuum condition, with no external neutron source, we can extract more parallelism from the sweep operation. In fact, we can start the sweep for all 8 octants in parallel as indicated on Figure 4, provided the scalar flux reduction is performed by using mutexes, in order to avoid race conditions which can lead to wrong results. This corresponds to the parallelization of the most external forall loop of the Algorithm 2. By so doing, the potential parallelism is multiplied by 8 which can significantly improve the speed-up especially when a large number of cores are used to sweep small geometries. It is an ongoing work to extend this strategy to symmetric boundary conditions. Basically you start the sweep of an octant as soon as its first macrocell has been swept for the previous octant: it follows a pipeline of the sweep over the octants.

SIMD Parallelism

During each sweep, the angular directions that belong to the same octant are computed via SIMD instructions allowing to process simultaneously up to 8 simple precision floating point operations when using Intel AVX (Advanced Vector Extensions) instructions. In this section we recall the basic of the SIMD programming model with emphasis on Intel SSE and AVX; then we present how we handle these instructions in DOMINO.

Figure 4: Parallelization of the sweep operation: all 4 octants are swept simultaneously. Pink colored macrocells are concurrently processed using mutexes: the first thread that starts the sweep of a given macrocell activates a lock that prevents other threads from starting to process the same macrocell and at the end, the lock is released so that another thread can start processing this macrocell.

SIMD Programming Model Vectorization is a computing model which consists in applying a single instruction to multiple data (SIMD). Most modern processors provide dedicated functional units that implement such instructions. Exploitation of these units can be done by relying on automatic vectorization by the compiler. However, while automatic vectorization is possible for simple loops, it is generally not the case for complex kernels. Meanwhile, by explicitly using assembly instructions or compiler intrinsics corresponding to the target architecture, one can make use of these vector units with the expense of some hardware constraints specific to vector instructions. In fact, vector instructions operate on packed data loaded inside specialized registers of fixed size. To perform fast load and store operations, data items need to be well aligned on cache boundary: 16 bytes for Intel SSE and 32 bytes for Intel AVX. Sometimes we have to resort to padding to satisfy this requirement. As an example, when loading 256 bits packet data with Intel AVX in a contiguous memory region of size 256 + 32 ⇥ 3 = 352 bits, we need to extend this region with 512 -352 = 160 bits at the memory allocation stage. Attention should be ported on padding as it increases global memory consumption and useless computations.

Angular Vectorization For the octant currently being swept, inside each spatial cell, we compute simultaneously several angular directions belonging to the octant (the forall loop on line 0 of Algorithm 2). SIMD instructions operate on packs of directions. The pack size, which depends on the precision and on the target architecture, is 8 in single precision and 4 in double precision on AVX enabled processors; on SSE enabled processors this pack size is divided by 2. Handling any angular quadrature order requires us to set up a padding system: for example, S 16 Level Symmetric angular quadrature formula gives 16(16 + 2) = 288 angular directions or 288/8 = 36 directions per octant; when using single precision Intel AVX, as 36 is not a multiple of 8, we perform 40 angular directions processing per spatial cell corresponding to an efficiency of 36/40 = 0.9.

On the other hand, product quadrature formulas, such as Gauss-Legendre, give more flexibility to overcome this limitation as you can choose a combination of the number of azimutal and polar directions that give a total number of directions divisible by, say, 8.

C++ programming language is well suited to deal with this vectorization procedure. Instead of hard-coding compiler intrinsics in the code, the arithmetic operators needed are overloaded by their corresponding intrinsics; we rely on Eigen2 for doing this job. Hence we are able to benefit from new vector instructions of the next generation of processors without having to modify the code, improving maintainability and readability of the code. Listing 1 shows a snapshot of the SIMD implementation of the sweep algorithm, which features some constructs of then Eigen library.

IV. Sweep Scalability and Efficiency

We analyze the performance of the sweep operation for a simple 2-group k ef f computation with 480 ⇥ 480 ⇥ 480 spatial cells with different angular quadratures on the A and B computing nodes described in Table 1 Figure 5 summarizes the performance of the 2-level parallel implementation (multicore+SIMD) of DOMINO for different Level-Symmetric angular quadratures on Node A. Since the arithmetic intensity of the computation increases with the number of direction, the sweep performance improves with the angular quadrature order. Because of their relatively low arithmetic intensities, the parallel speed-ups of the S 2 and S 4 computations saturate above 16 cores while it is not the case for the S 8-16 quadratures. Note that S 2 and S 4 timings are exactly the same due to our padding strategy (see section 2). Here we use SIMD units that process at least 4 angular directions that belong to the same octant simultaneously. Since a given S N Level-Symmetric defines d o = N(N + 2)/8 angular directions per octant, we obtain d o = 1 (resp. d o = 3) for the S 2 (resp. S 4 ) quadrature leading to 4 -1 = 3 (resp 4 -3 = 1) padded directions. In order to assess the sweep perfomance we compare the corresponding GFLOPS (Giga FLOating Point operations per Second) to the theoretical Peak performance of the computing Nodes. These two metrics are defined as:

• GFLOPS value is estimated by dividing the number of floating point operations by the completion time:

GFLOPS = 25 ⇥ N cells ⇥ N dir Time in nanoseconds , (4)
where the factor 25 is the total number of the floating point operations in the sweep operation, counted in section 2.

• The peak performance of a supercomputing node is evaluated by taking into account all the operations that all the available floating point arithmetic units can complete per clock cycle:

Peak = ncores ⇥ frequency ⇥ SIMD width ⇥ SIMD units. (5)
Table 2 summarizes the 32-core performance corresponding to Figure 5. While the S 2 and S 4 timings are identical (27.9s), the corresponding GFLOPS figures (11.2 and 33.8) differ since we only count useful angular directions (1 and 3) in the GFLOPS evaluation (Eq 4). In the S 16 case (N dir = 288), the sweep performance reaches 248.9 GFLOPS which corresponds to 43% of the Node A peak performance in single precision.

Impact of Parallelizing Over the Octants

Figure 6 shows the performance behavior that one can observe when allowing to compute all the octants in parallel as explained in the previous section. For the large spatial mesh, the parallelism over the octant is negligible. On the contrary, for the smaller mesh the additional amount of task-parallelism helps to maintain a speed-up almost linear up to 32 cores. The same experiments were carried out using Node B. Here we use Gauss-Legendre quadratures. The corresponding result are summarized in Figure 7. This nodes accepts AVX instructions but is still compatible with SSE4 version of the sweep. Scalar, SSE and AVX versions of the sweep have been evaluated on Node B using 16 threads. The results are summarized in Table 3. One can see that, in this case, the performance does not increase much when switching from SSE to AVX. We are currently investigating this issue. The performance of the sweep operation reaches 237.7 GFLOPS which corresponds to 35.71% of the peak performance of the supercomputing node B in single precision. As a general trend, the 16-core Sandy Bridge node leads to sweep performances that are on par with the 32-core Nehalem node.

In the next section we show that the good performance of the sweep kernel allows DOMINO to solve 3D PWR reactivity problems efficiently.

V. 3D PWR Full Core Performance

The benchmark used for calculations is described in Ref. [START_REF] Courau | Hybrid Parallel Code Acceleration Methods in Full-Core Reactor Physics Calculations[END_REF] It corresponds to a simplified 3D PWR first core loaded with 3 different types of fuel assemblies characterized by a specific 235 U enrichment (low, medium and highly enriched uranium). No inserted control device is considered in this core model. Along the z-axis, the 360 cm assembly is axially reflected with 30 cm of water which results in a total core height of 420 cm. The 3 types of fuel assemblies appear in Figure 8 where the central assembly corresponds to the lowest enrichment, while the last row of fuel assemblies have the highest enrichment to flatten the neutron flux.

Each fuel assembly is a 17 ⇥ 17 array with a lattice pitch of 1.26 cm that contains 264 fuel pins and 25 water holes. The boundary condition associated with this benchmark problem is a pure leakage without any incoming angular flux. The associated nuclear data, an 8-group and a 26-group libraries, were obtained from a fuel assembly heterogeneous transport calculation performed with the cell code DRAGON. [START_REF] Marleau | A User's Guide for DRAGON 3.05[END_REF] A 26-group computation has been carried out with DOMINO using a S 16 Level-Symmetric angular quadrature. Table 4 summarizes the results of this k eff computation. We assess the accuracy of DOMINO by comparing obtained neutron flux and k eff to a reference obtained with MCNP5. To make easier the comparisons to MCNP results, the 3D 26-group fluxes are integrated over energy and space. Practically, a group collapsing is performed from 26 to 2 energy groups with a boundary fixed 4 are computed in the form:

∆k ef f = k DOMINO ef f -k MCNP ef f k MCNP ef f ⇥ 10 5 in pcm, |δφ g | = max {p i } | φ DOMINO g (p i ) -φ MCNP g (p i ) φ MCNP g (p i ) |,
where p i stands for the pin cell index. The 26-group k ef f computation completes in less than 11h, and corresponds to a sweep performance of 235 GFLOPS which is equivalent to 40.74% of the peak performance of the supercomputing node A. In order to evaluate the performance of the complete DOMINO k ef f computation we divide the sweep operation number (Eq.4) by the complete wall-clock time [START_REF] Davidson | Massively Parallel, Three-Dimensional Transport Solutions for the k-Eigenvalue Problem[END_REF] . It leads to an average performance of 111 GFLOPS.

26-group

N
Based on our 3D PWR model, the comparison with MCNP and the performance evaluation in terms of GFLOPS show that DOMINO is both accurate and computationally efficient. However it is important to compare DOMINO to other deterministic solvers in order to check if all the counted floating point operations correspond to a truly efficient algorithm.

DOMINO Compared With Pentran And Denovo

There exists several deterministic neutron transport solvers based on discrete ordinates method; we propose here a comparison of DOMINO with two solvers that can address 3D Cartesian meshes.

Pentran PENTRAN [START_REF] Courau | 3D Neutron Transport and HPC: A PWR Full Core Calculation Using PENTRAN SN Code and IBM BLUEGENE/P Computers[END_REF] is a 3D discrete ordinates code that can automatically distribute the problem phase space among the angular (A), energy (E), and spatial variables (S) using a parallel memory and task allocation across a virtual processor array (VPA). One has the possibility to focus on one of the axes of the VPA by assigning a specific weight for this axis; the parallelism is implemented using a pure MPI programming model to exploit both shared and distributed memory supercomputers.

We compare S 8 DOMINO and PENTRAN k eff computations for the 8-group version of our PWR Benchmark. [START_REF] Courau | 3D Neutron Transport and HPC: A PWR Full Core Calculation Using PENTRAN SN Code and IBM BLUEGENE/P Computers[END_REF] The results, obtained on Ivanoe, the EDF Intel Xeon based supercomputer [START_REF] Petrini | Multicore surprises: Lessons learned from optimizing Sweep3D on the Cell Broadband Engine[END_REF] , are summarized in Table 5. Note that different values for the stopping criterion associated to the flux (✏ φ ) has been used for the different computations. The PENTRAN results have been extracted from ref. [START_REF] Courau | Hybrid Parallel Code Acceleration Methods in Full-Core Reactor Physics Calculations[END_REF] Using a single node DOMINO is faster (67 min) than PEN-TRAN running on 289 nodes (4752 min). Actually, PENTRAN is a general-purpose transport code and is not specially adapted to 3D PWR simulations. In particular, PENTRAN lacks an efficient acceleration scheme that is required to cope efficiently with highly diffusive media. In order to compare the performances of the parallel implementations of the codes, a nonaccelerated DOMINO computation has been carried out (third column in Table 5). This non-accelerated DOMINO computation remains faster (573 min) than PENTRAN.

Several factors can explain this large performance difference. The PENTRAN code is built upon a single programming paradigm and does not explicitly address shared memory and vector issues that tend to become more and more important on modern processors. In addition, the spatial discretization scheme of DOMINO (DD0) is simpler than the DTW scheme used in PENTRAN and should consume less CPU cycles. At last it is difficult to evaluate the efficiency penalty arising from using a distributed architecture (289 nodes) compared to our shared memory approach. We are about to complete a distributed version of DOMINO that would allow us to conclude on this point.

PENTRAN DOMINO DOMINO

No Accel PWR 8g [START_REF] Courau | 3D Neutron Transport and HPC: A PWR Full Core Calculation Using PENTRAN SN Code and IBM BLUEGENE/P Computers[END_REF] PWR

Denovo

Denovo [START_REF] Evans | Threedimensional full core power calculations for pressurized water reactors[END_REF] is another deterministic neutron transport solver based on the discrete ordinates method, for radiation shielding and reactor physics applications under active development at ORNL. It implements a multilevel parallel decomposition on the phase space. This decomposition allows concurrency over energy in addition to space-angle parallelism. The spatial parallelism uses the KBA-based parallel decomposition; eigenvalue solvers implemented in Denovo are power iteration, as in DOMINO, and an Arnoldi solver. Extracted from ref, [START_REF] Evans | Threedimensional full core power calculations for pressurized water reactors[END_REF] we compare S 12 calculations performed with the Denovo solver on a 2-group version of the PWR-900 benchmark [START_REF] Courau | Specifications of a 3D PWR Core Benchmark for Neutron Transport[END_REF] with DOMINO. The Denovo computation was run on the Jaguar XT5 supercomputer (18688 compute nodes, each with dual 2.6 GHz AMD 6-core Istanbul processor). Please note that this machine differs from Ivanoe and that the performance comparison is not very precise. In addition, the axial meshes used in DOMINO and Denovo are slightly different (756 vs 700). Nonetheless, we hope that this comparison provides a correct trend on how DOMINO compares to Denovo. The results are summarized in Table 6.

Denovo is able to use efficiently a large number of nodes and thus solves the 2-group problem much faster than DOMINO. This illustrates our need to address many-node computations with DOMINO. As a positive point, we confirm that DOMINO makes an efficient use of the available 32 cores. Although the number of solved DoF per minute and per core (DoF/min/core) of DENOVO (1.88 ⇥ 10 6 ) is lower than the one of DOMINO (64.7 ⇥ 10 6 ), one should keep in mind that Amdahl's law imposes this metric to decrease with the number of cores. Note that this number (64.7 ⇥ 10 6 ) is significantly higher than what we measure for the 8-group computation (16.7 ⇥ 10 6 ). Indeed one needs more inner-group iterations when the group number DENOVO DOMINO PWR 2g [START_REF] Davidson | Massively Parallel, Three-Dimensional Transport Solutions for the k-Eigenvalue Problem[END_REF] PWR 2g increases.

VI. Conclusion

The two-level (multi-core+SIMD) parallel implementation of the sweep algorithm in DOMINO has been described. It leads to a very efficient deterministic Cartesian transport solver that runs on shared-memory HPC nodes. Compared to Monte-Carlo solution (MCNP), DOMINO exhibits very accurate results both on 8 and 26 energy groups 3D PWR Benchmark. On comparable 3D nuclear core reactivity computations, DOMINO has shown to be very fast compared to PENTRAN and not very far from Denovo, a solver that addresses much larger HPC devices. We believe that the very high Flops/Watt ratio of DOMINO makes it a very promising building block for a future distributed memory nuclear simulation tool.

Figure 1 :

 1 Figure 1: Degrees of freedom for the DD0 element

Figure 2 :

 2 Figure 2: Data dependencies over one spatial cell in 2D 2. The Sweep Algorithm Before going into details of the parallel implementation, Algorithm 2 describes the sweep operation which represents the most time-consuming operation inside the solver.The aim of this algorithm is to solve the space-angle problem, by inverting the streaming operator of the monogroup transport equation: lines 0-0 of Algorithm 1. The volumic flux computation inside the cell c (line 0) needs to know incoming data for this cell: incoming angular fluxes, L , B , F , the total cross section Σ t and the source term S . Outgoing angular flux are then updated on lines 0, 0 and 0. At this point, we can observe that once a cell has been swept for a given octant, incoming

Algorithm 2 :

 2 The Sweep Operation forall the o 2 Octants do forall the c 2 Cells do

Figure 3 :

 3 Figure3: Wave front algorithm applied to sweep across the mesh macrocells. The incoming fluxes come from the bottom left corner. The processed macrocells are pink colored while the white macrocells are still waiting. Only blue dashed border macrocells are ready for processing and their respective indexes are stored in a task list from which a pool of threads concurrently pick their work item. This list is dynamically updated each time a new macrocell has been processed.

Figure 5 :

 5 Figure 5: Time per sweep over all angular directions as a function of the core number on Node A (mesh: 480 ⇥ 480 ⇥ 480).

Figure 6 :

 6 Figure 6: Performance impact of the octant parallelization (mesh: 120 ⇥ 120 ⇥ 120 and 480 ⇥ 480 ⇥ 480, using S 8 Level Symmetric quadrature on Node A)

Figure 7 :

 7 Figure 7: Time per sweep over all angular directions, using Gauss-Legendre quadrature formula, as a function of the core number on Node B (mesh: 480 ⇥ 480 ⇥ 480)

Figure 8 :

 8 Figure 8: 2D-View of the PWR Core Model

Table 1 :

 1 . Characteristics of two computing nodes. The theoretical peak performance figures are computed from Eq. (5).

	Node A	Node B

Table 2 :

 2 Time per sweep over all angular directions using 1 and 32 cores on Node A (mesh: 480 ⇥ 480 ⇥ 480).

		N dir	Seq.	32-Core Speed GFLOPS % Peak
			Time (s) Time (s)	Up	(Eq 4)	Perf.
	S2	8 ⇥ 1	27.9	1.96	14.2	11.2	1.94%
	S4	8 ⇥ 3	27.9	1.96	14.2	33.8	5.8%
	S8 8 ⇥ 10	41.5	1.84	22.6	120.4	20.8%
	S12 8 ⇥ 21	65.0	2.4	27.0	190.6	32.9%
	S16 8 ⇥ 36	89.7	3.2	28.0	248.9	43.0%

Table 3 :

 3 Impact of SIMD parallelism for a 32 ⇥ 16 Gauss-Legendre quadrature sweep performance (mesh: 480 ⇥ 480 ⇥ 480).

Table 4 :

 4 3D PWR k eff computation with S 16 angular quadrature on Node A (see Table1). The k eff tolerance is set to 1 ⇥ 10 -5 . The flux and k eff obtained with DOMINO are compared to a MCNP5 reference.

Table 5 :

 5 DOMINO-PENTRAN Comparison for a S 8 8-group 3D PWR k eff computation.

Table 6 :

 6 DOMINO-Denovo Comparison for a S 12 2-group 3D PWR k eff computation.

	Computer	Jaguar	Ivanoe
	XY Mesh	578 ⇥ 578 578 ⇥ 578
	Z Mesh	700	756
	N dir	168	168
	N group	2	2
	k eff tol N DoF (⇥10 9 )	1 ⇥ 10 -3 78.6	1 ⇥ 10 -5 84.9
	N cores	20400	32
	Wall-clock
	time (min)	2.05	41
	DoF/min (⇥10 6 )	38330	2069
	DoF/min/core (⇥10 6 )	1.88	64.7

Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.

Here we neglect all the other flops of the rest of the code.

Ivanoe was ranked 213 in the June 2013 Top

list.

ACKNOWLEDGMENTS

The authors wish to thank Stanislas Odinot from Intel Corporation for providing access to different Intel devices. We also thank our reviewers who encouraged us to compare DOMINO with other transport codes.


G Vermot Des Roches

-E Balmes
email: balmes@sdtools.com

-H Ben Dhia

- R Lemaire
email: remi.lemaire@fr.bosch.com

- T Pasquet

Compatibility measure and penalized contact resolution for incompatible interfaces

Keywords: Penalized contact, Regularization, (1-ǫ)-Compatibility, Non-conforming mesh MOTS-CLÉS : Contact pénalisé, Regularisation, (1-ǫ)-Compatiblité, Maillages non compatibles

Handling of large industrial mechanical assemblies implies structure interactions commonly modeled with contact formulations. In cases where component interfaces are discretized using non conforming meshes, classical contact solutions have difficulties producing correct contact pressure fields. The method presented in this paper gives a relevant measure of interface compatibility and shows how it can be exploited to obtain regular contact pressures or limit over-integration in the contact formulation.

RÉSUMÉ. La gestion d'assemblages mécaniques complexes nécessite souvent des modèles d'interaction de structure par une formulation en contact. Lorsque la discrétisation élément finis des interfaces est non compatible, les résolutions classiques peuvent êtres non régulières. La méthode présentée dans ce papier donne une mesure pertinente de la compatibilité des interfaces et montre comment l'exploiter pour régulariser les pressions de contact ou limiter la sur-intégration dans la formulation du contact.

Introduction

The constant improvement of computational power allows the implementation of numerical prototyping at every stage of the design cycle, thus reducing costs and development time. In particular, simulation of complex mechanical assemblies with the finite element method is now widespread in industry. Such models however imply the ability to handle components interactions of discretized interfaces at the system level.

Interaction between moving parts typically requires contact models which can be split in two main categories, using either Lagrange or penalized formulations. Lagrange formulations enforce strict non-penetration, while penalization authorizes a relative level of interface interpenetration through a pressure-gap relationship. Both models are an idealization of contact, indeed, the reality of contacting surfaces implies interpenetration seen at the microscopic level as the compression of contacting asperities. Non linear penalization formulations, which is the modeling context of this study, can then be seen as a characterization at the system level of the local contact behavior.

External modeling and model generation constraints usually found for industrial models (automatic mesh generation at the component level, mesh provided by contractors,... ) often lead to non-conforming meshes for the contacting interfaces. The definition of master/slave surfaces is also a question that is often with no clear answer, in particular with partially covering surfaces due to large displacements or the presence of holes.

The problem of implementing bilateral or contact coupling for non-conforming meshes has been widely studied in the finite element [START_REF] Babuska | The finite element method with Lagrange multiplier[END_REF][START_REF] Bernadi | Domain decomposition by the mortar element method[END_REF][START_REF] Kim | Interface Element Method (IEM) for a partitioned system with non-matching interfaces[END_REF] and component mode synthesis [START_REF] Balmes | Use of generalized interface degrees of freedom in component mode synthesis[END_REF][START_REF] De Klerk | General Framework for Dynamic Substructuring : History, Review and Classification of Techniques[END_REF] literature. The strategy presented here is based on the notions of pseudo compatibility introduced in (Ben [START_REF] Dhia | Mesure de compatibilité et application aux problèmes de sousstructuration[END_REF] extended to tridimensional meshed volumes and the need to use specific integration rules for interfaces [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF].

Section 2 illustrates the problematic and proposes a measure of non conforming interface compatibility which can be taken as a quality indicator. The measure is directly exploited in section 3 to regularize contact solutions by describing interface fields as a linear combination of quasi-compatible shapes, through either primal or dual formulations. Section 4 eventually provides applications to an academic model showing large mesh refinement differences and to an industrial pad/disc model from Bosch, showing large non compatible interfaces.

A measure of interface compatibility

A compatibility measure is an indicator of the interface ability to transmit a field from one surface to the other. For a general non conforming interface, not all field distributions can be equally represented by both surfaces. A measure of compatibility thus needs to evaluate differences between the representations of a given field on both surfaces.

Illustration of the problem

The usual finite element formulation of penalized contact can yield poor results in the case of non conforming interface meshes, especially when mesh sizes differ in contacting surfaces and partial element coverage is present -e.g. due to holes or global motion. The academic illustration presented in figure 1 is a contact between two cubes in vertical translation. The lower cube is clamped at his base, drilled and finely meshed. The upper cube is plain, coarsely meshed and a pressure is applied to its top surface.

Figure 1. Example of incompatible mesh not covering. Reference solution is plotted in the middle, against the 1 (right) and 4 (far right) point rule integration solution using the upper cube as master Figure 1 presents displacement results obtained for different formulations. Using the coarse upper cube as the master means that a lot of slave nodes are not matched (or not seen by any master contact point). The richer shapes of Γ 2 cannot be represented on Γ 1 and are therefore not seen nor constrained by the contact formulation on Γ 1 . This pattern is close to hourglass modes observed for under-integrated elements.

Introducing very rich integration avoids oscillations but can significantly augment computational costs and leads to excessive stiffness, or locking, when zero gap is enforced strictly [START_REF] Dhia | Mesure de compatibilité et application aux problèmes de sousstructuration[END_REF]. The following sections will introduce an energy function for the gap and show how this can be used to quantify and possibly relax compatibility of displacements on the interfaces.

Definition of gap energy

Contact between two solids Ω 1 and Ω 2 , is defined between two surfaces Γ 1 and Γ 2 . The compatibility is computed on the subset of effective contact Γ of Γ 1 × Γ 2 , as shown figure 2. Any contact formulation needs to evaluate the gap field between two surfaces, which is defined for a displacement u(q 1 ) of Γ 1 and u(q 2 ) of Γ 2 , where q 1 (resp. q 2 ) is the discretized displacement on Γ 1 (resp. Γ 2 ) as

{g(u(q1), u(q 2 ))} = {u(q 1 ) -u(q 2 )} T {n} [1]
The first step is to define a scalar product on Γ 1 × Γ 2 that is a norm for the gap.

The strain energy of a penalized contact with a uniform contact stiffness density k is chosen here. For a displacement u, this energy is given by

E p = Γ kg (u(q 1 ), u(q 2 )) 2 dS [2]
This scalar product must be approximated as precisely as possible [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF]. The computation strategy chosen here is an automated Delaunay triangulation over the subset Γ of Γ 1 × Γ 2 in effective contact. The scalar product is computed by numerical quadrature from the mesh of Γ

E p = xi kw i J(x i )g i (u i (q 1 ), u i (q 2 )) 2 [3]
In practice, the gap is computed on N Γ integration points (Gauss) and linearly linked to the normal displacement observation matrices C 1

N NΓ×NΓ 1 (resp. C 2 N NΓ×NΓ 2 ) of Γ 1 (resp. Γ 2 ) on Γ.
This yields the discretized gap formulation

{g} i = C 1 Ni {q 1 } -C 2 Ni {q 2 } [4]
The contact strain energy is then of the form

E p = q 1 q 2 T C 1 Ni -C 2 Ni T \ kw i J(x i ) \ C 1 Ni -C 2 Ni q 1 q 2 [5]
which defines the scalar product on Γ. The scalar product matrix is noted [A] and is written by blocks

E p = {q} T [A] {q} = q 1 q 2 T A 11 A T 21 A 21 A 22 q 1 q 2 [6]

Robust compatibility computation

In [START_REF] Dhia | Mesure de compatibilité et application aux problèmes de sousstructuration[END_REF] the so-called (1 -ǫ)-compatibility is computed by considering the norm difference between a displacement q 1 of an interface Γ 1 and its projection to the facing interface π 1 2 {q 1 } ∈ Γ 2 . The projection, defined as the vector of Γ 2 that minimizes the gap energy, verifies

[A 21 ] {q 1 } -[A 22 ] π 1 2 {q 1 } = 0 [7]
This suggests a formulation of the (1 -ǫ)-compatibility as a Rayleigh quotient, as

C 1 2 ({q 1 }) 2 = π 1 2 {q 1 } A22 {q 1 } A11 = {q 1 } T [A 21 ] T [A 22 ] -1 [A 21 ] {q 1 } {q 1 } T [A 11 ] {q 1 } [8]
Such formulation raises robustness issues as partial element covering tends to yield ill conditioning of [A 11 ] and [A 22 ]. The formulation presented here computes the singular value decomposition (SVD) of [A] (rather than

[A 21 ] suggested in (Ben Dhia et al., 2003)). As [A] is symmetric, [A] = . . . u 1 i u 2 i . . . \ σ i\ {U } T [9]
The SVD generates a displacement basis {u i } of the coupled interface with the level of gap energy generated σ i . Low singular values thus characterize compatible interface displacements. The compatible displacement fields obtained are then known to be represented on both surfaces well enough -which will be quantified by a tolerance ǫ.

A new definition of interface compatibility is suggested as

C 1 2 (ǫ) = card i/ σi maxi σi ≤ ǫ min i=1,2 (card(Γ i ))
[10]

It is useful to note that the gap observation (4) can be performed with the full contact set. Unmatched contact DOF will generate zero terms on the diagonal of [A]. No conditioning issue affects the SVD computation which would output in that case a null (to numerical precision) singular value associated with shapes showing unmatched points displacements only. The description of the unmatched point displacements is then a full recombination of the finite element basis. The refinement of Γ is nevertheless necessary for a good description of the non spurious contact points.

Using (1-ǫ) vector pairs to solve contact problems

The compatibility measure through an SVD generates couples of quasi-compatible displacements on both sides of an interface. The following sections discuss their use for contact problems.

Primal and dual formulation

The primal formulation constrains the interface displacement to the subspace of quasi-compatible displacements. The model DOF are sorted to segregate the interface DOF (q i1 and q i2 ) from the rest (q c ), thus allowing the projection to the generalized quasi-compatible interface DOF (q ǫ1 and q ǫ2 ). Regular contact force distributions are then a consequence of the imposed displacement regularity.

   q c q i1 q i2    =   I 0 0 0 U 1 ǫ 0 0 0 U 2 ǫ      q c q ǫ1 q ǫ2    = [T ǫ ]    q c q ǫ1 q ǫ2    [11]
The gap observation is consequently projected and can be written

{g} = C 1 N i U 1 ǫ {q ǫ1 } -C 2 N i U 2 ǫ {q ǫ2 } [12]
Exact compatibility (zero gap on Γ) is known to generate locking issues, as discussed for example in [START_REF] Balmes | Use of generalized interface degrees of freedom in component mode synthesis[END_REF][START_REF] Dhia | Mesure de compatibilité et application aux problèmes de sousstructuration[END_REF]. Introducing quasicompatible displacements relaxes the constraint while avoiding large relative displacements.

In a dual formulation, the projection on the (1-ǫ)-compatible displacement pairs of the contact forces is assumed to be zero. This clearly avoids local stress concentrations associated with locking but also leads to larger gaps. Dual formulations also have the advantage of preserving the use of physical displacement DOFs.

Regularized gap function for non-linear contact laws

In general penalized contact approaches, the load is a non linear function of the gap. Exponential functions are in particular used for brake squeal studies (Vermot des [START_REF] Vermot Des Roches | Time simulation of squeal phenomena in realistic brake models[END_REF]. In a dual formulation, locking is avoided but the gap is no longer regular. It is thus proposed to use a regularized gap observation on Γ to compute contact forces.

For a surface Γ n the regularized displacement is a linear combination of the trace on Γ n of the compatible vectors,

[U n ǫ ] {q ǫ n } NΓ n ×Nǫ = [U n ǫ ] NΓ n ×Nǫ {α n } Nǫ [13]
where {α n } is a solution of the minimization problem

{α n } = min {α} [U n ǫ ] {α} -{q n } K [14]
This defines a pseudo-inverse, as the solution searched is the closest to {q n } possible, relatively to a norm K. The choice of K is fully open. It cannot be a restriction of [A] due to conditioning problems, but considering a norm in displacement or strain on Γ n seems meaningful. The identity was chosen here so that the modified gap observation is

{g} i = C 1 N i U 1 ǫ U 1 ǫ + {q 1 } -C 2 N i U 2 ǫ U 2 ǫ + {q 2 } [15]
Projection of both sides of the interface can be considered. In master/slave configurations, a slave surface only projection is sufficient to suppress the observation of non quasi-compatible displacements.

Opportunities for contact under-integration

Obtaining proper quasi-compatible modes requires a good representation of the scalar product presented in section 2.2. Relatively rich integration rules on both surfaces are therefore needed for these computation steps. However, the selection of quasi-compatible modes typically restricts the number of needed interface degrees of freedom significantly. The number of integration points needed for contact can consequently be decreased as information redundancy occurs. This is particularly interesting in vibration studies for which the contact surface is assumed constant during long time simulations, so that coarser rules will decrease computation time.

The proposed strategy is to use a fine integration rule to compute [A] and U i ǫ defined in (8) and ( 9), then to switch to a coarser rule for the evaluation of the regularized gap (15) and resulting contact loads. Such under-integration could be optimized by the use of the Orthogonal Maximum Sequence method [START_REF] Balmes | Orthogonal Maximum Sequence Sensor Placements Algorithms for modal tests, expansion and visibility[END_REF] to create an optimal set of contact points associated to the quasi-compatible shapes to be observed.

Illustrations

The concepts of sections 2 and 3 are applied to the academic model which features great mesh refinement differences, and to an industrial model featuring a relatively large interface. Regularization patterns are presented on the first case while under-integration strategies in the scope of an improvement of dynamical behavior are discussed on the latter model.

A drilled cube example

This section illustrates the resolution method presented to the drilled cube example, using the upper cube as contact master, which yields poor results when no regularization is considered. The first step to the resolution is the creation of Γ, and the SVD computation of [A], as plotted in figure 3. The gap is computed on the Gauss integration points of Γ for which a Delaunay triangulation is performed for visualization purposes. The compatible shapes shown in figure 3a or 3b correspond to very small gaps. The incompatible shapes in figure 3c and 3d show a high singular value (8% and 12%) and correspond to a large gap visible on Γ. In this problem, 25 quasi-compatible shapes at most exist since this is the number of DOF of the coarse interface. A clear jump can be seen in figure 4 after this threshold.

To illustrate the regularized gap approach, the projected observation

C 1 N U 1 ǫ U 2 ǫ
+ is plotted and compared to C 1 N in figure 5. The observation is the dual of the contact pressure as they are directly related by the pressure-gap contact law. The regularization method clearly distributes the gap observation to the contact point vicinity, instead of having a point-to-matched-point observation. Basic solutions provide poor contact fields, using a richer integration rule for that case is however an improvement. Indeed the approached contact area seen by the covered contact points increases by 25% when switching from the center point rule to the 2x2 points rule.

The primal resolution directly constraints displacements but can raise locking issues. No such pattern is observed for the cube and direct displacement restriction yields in this case better contact pressure fields. It is however costlier to implement as the interface projection would typically require DOF reordering prior to solve.

Application to an industrial brake model

State-of-the-art brake models used in this section are provided by Bosch. The system is meshed automatically by component using 10 node tetrahedron elements, resulting in non conforming interfaces. Figure 7 shows the components retained, the disc and a pad with its backplate and lining. Pressure is applied using the trace of the caliper fingers on the backplate.

Static responses are properly handled here using a rich contact integration rule (even over-integration). For applications in dynamics, one seeks to compute small oscillations around a bilateral contact state. The measure of incompatibility introduced in the paper shows that although the mesh sizes are similar, compatibility is not very good. The combination of incompatible mesh and rich integration, then leads to locking issues which will be illustrated at the end of this section. Applying the concepts presented in sections 2.2 leads to the compatible singular values and associated compatibility measure shown in in figure 8. The first 10% of compatible shapes are free movements for the uncovered nodes of the partially covered elements. A threshold ǫ value is seen before more quasi-compatible shapes are found. More than 70% of shapes are above a ǫ = 10 -4 threshold. Thus although mesh sizes are similar, the compatibility is low. Enforcing strict displacement equality is thus expected to induce locking. The resolution of a contact problem with an exponential law is here satisfying using a rich integration rule (6 points for a 6 nodes triangle), but requires 1824 contact points. A center point integration rule uses only 394 points but results in strong contact irregularities, as illustrated in figure 10. As well documented in the literature, over-integration is a good solution to solve static contact problems. In dynamics, enforcing permanent contact exactly for incompatible meshes leads to locking. As an illustration, figure 11 shows a transfer function from the vertical displacements of a disc point to a pad point. Modes are computed enforcing perfect bilateral contact (no gap) of the contact interface through either the basic contact normal displacement observation or the modified one (equation ( 15)), using ǫ = 10 -3 . Left: Dynamic behavior differences between the use of a 6 point integration rule contact integration law, and the regularization strategy based on a 1 point integration rule for ǫ = 10 -3 (42% of quasi-compatible shapes kept). Right: Frequency differences observed

The over-integrated solution shows positive frequency shifts from 6 to 30%, corresponding to resonance shifts visible in the figure even at low frequencies and a clear stiffening, or locking, of the structure. The application of interest being brake squeal simulations (Vermot des [START_REF] Vermot Des Roches | Time simulation of squeal phenomena in realistic brake models[END_REF], the typical frequency range of interest goes up to 16kHz and the great differences observed in this range must clearly be addressed.

In a penalty based contact, the existence of a tangent stiffness may alleviate some of the problems. Some of the locking may disappear for a soft contact. The nonphysical nature of such contact is however not very satisfactory.

Conclusion

The article presented a general methodology to deal with field transmission in non conforming interface meshes. A scalar product on the interface is used to define compatible motion. Using a very rich integration for the evaluation of the gap energy, avoids stress concentrations associated with integration points in coarse meshes and is applicable to partially covered interfaces and cases where no obvious slave/master strategy exists.

Using rich integration the contact must be verified at a high number of points so that only strictly compatible deformations of the interface are allowed. This leads to an over stiffening of the model, known as locking, which is particularly visible for applications in dynamics where the loads transmitted in the interface are not as smooth as in contact problems.

A singular value decomposition of the gap energy operator generates pairs of (1-ǫ)compatible deformations which properly handles difficult cases of very partial overlap. The relation between total number of DOF and number of (1-ǫ)-compatible vectors is a direct numerical indicator of mesh compatibility. The industrial application shown illustrates the fact that meshes with comparable element sizes can have relatively low compatibility.

(1-ǫ)-compatible vector pairs can be used to avoid inappropriate gaps in a primal method that only allows compatible motion or to compute a regular pressure field by introducing a regularized gap function (dual method). An academic example illustrated how this is applicable to pathological cases such as gross under-integration which can be associated with the selection of a coarse master mesh. Another use might be the introduction of a coarse contact rule for applications in long transient simulations.

The need to relax continuity to avoid locking was finally illustrated for a real example. This leads to correct coupled behavior while avoiding pathologies of underintegration. It clearly appears that locking is particularly important for contact between solids. And that soft contact law, by opposition to Langrangian formulations or penalty approaches with stiff contact, will be less sensitive to the problem. The trade-off between the local nature of gap evaluations and the global nature of (1-ǫ)compatible deformations still needs further discussion.

Figure 2 .

 2 Figure 2. Definition of a domain Γ for the compatibility computation

Figure 3 .

 3 Figure 3. Shapes resulting from the SVD of [A] plotted on Γ 1 and Γ 2 . The gap is computed and plotted on the Gauss points of Γ

Figure 4 .

 4 Figure 4. Singular values of [A] in increasing incompatibility

Figure 5 .

 5 Figure 5. Initial and regularized observation for a given Gauss integration point (o) of the master surface

  (a) Center rule (b) Center Dual Reg. (c) Center Primal Reg.

Figure 6 .

 6 Figure 6. Gap observations obtained using the center integration rule for contact, observed on the rich 3x3 rule Gauss points Delaunay triangulation, ǫ = 10 -5

Figure 7 .

 7 Figure 7. Pad/Disc model and superimposed interface discretization

Figure 8 .

 8 Figure 8. Left: Singular values of [A] -Right: (1-ǫ)-compatibility as function of ǫ Sample shapes resulting from the SVD are plotted in figure 9. Shapes are plotted on Γ 1 and Γ 2 which are the upper and lower surfaces. The gap between both surfaces displacement is plotted on Γ shown at the intermediate height.

Figure 9 .

 9 Figure 9. Sample quasi-compatible shapes. Γ is represented as a Delaunay triangulation of its Gauss points

Figure 10 .

 10 Figure 10. Basic contact pressure solution [MPa] for rich and poor integration rulesFigure 10c clearly illustrates that the regularized gap function avoids large oscillations and yields smoother results when a small ǫ (10 -5 ) is considered. The pressure does not however exactly converge towards the rich 6 point per element rule. Further

  Figure11. Left: Dynamic behavior differences between the use of a 6 point integration rule contact integration law, and the regularization strategy based on a 1 point integration rule for ǫ = 10 -3 (42% of quasi-compatible shapes kept). Right: Frequency differences observed