Comparison of algorithm in Aerosol and Aghora for compressible flows

21 Pages·2017·1.45 MB·English

Comparison of algorithm in Aerosol and Aghora for
compressible flows
Dragan Mbengoue, Damien Genet, Cedric Lachat, Emeric Martin, Maxime
Mogé, Vincent Perrier, Florent Renac, François Rue, Mario Ricchiuto
To cite this version:
Dragan Mbengoue, Damien Genet, Cedric Lachat, Emeric Martin, Maxime Mogé, et al.. Comparison
of algorithm in Aerosol and Aghora for compressible flows. [Research Report] RR-8200, 2013. ￿hal-
00773531v1￿
HAL Id: hal-00773531
https://hal.inria.fr/hal-00773531v1
Submitted on 14 Jan 2013 (v1), last revised 17 Jan 2013 (v2)
HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est
archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
Comparison of algorithm
in Aerosol and Aghora
for compressible flows
D.A. Mbengoue, D. Genet, C. Lachat, E. Martin, M. Mogé, V. Perrier,
F. Renac, F. Rué, M. Ricchiuto
G
N
E
+
R
F
0--
0
2
8
R--
R
A/
RI
N
I
N
RESEARCH R
S
I
REPORT
9
9
N° 8200 3
6
9-
4
January2013 2
0
N
Project-TeamsBACCHUSand SS
I
CAGIRE

Comparison of algorithm in Aerosol and Aghora for
compressible flows
D.A.Mbengoue∗,D.Genet†,C.Lachat‡,E.Martin§,M.Mogé¶,V.
Perrier(cid:107),F.Renac∗∗,F.Ru醆,M.Ricchiuto‡‡
Project-TeamsBACCHUSandCAGIRE
ResearchReport n°8200—January2013—16pages
Abstract: This article summarizes the work done within the COLARGOL project during CEMRACS
2012. The aim of this project is to compare the implementations of high order finite elements methods
forcompressibleflows thathavebeendevelopedat ONERAandatINRIAfor aboutoneyear, withinthe
AGHORAandAEROSOLlibraries.
Key-words: Finiteelements,parallelperformance,compressibleflows
∗INRIABordeauxSud-Ouest,EquipeBACCHUS
†INRIABordeauxSud-Ouest,EquipeBACCHUS
‡INRIABordeauxSud-Ouest,EquipeBACCHUS
§ONERA
¶LMAP,UniversitédePauetdesPaysdel’Adour
(cid:107)LMAP,UniversitédePauetdesPaysdel’Adour,INRIABordeauxSud-Ouest,EquipeCAGIRE
∗∗ONERA
††INRIA,serviceSED
‡‡INRIABordeauxSud-Ouest,EquipeBACCHUS
RESEARCHCENTRE
BORDEAUX–SUD-OUEST
351,CoursdelaLibération
BâtimentA29
33405TalenceCedex
Comparison of algorithm in Aerosol and Aghora for
compressible flows
Résumé: CetarticlerésumeletravaileffectuédurantleprojetCOLARGOLpendantle
CEMRACS2012.Lebutdeceprojetestdecomparerlesimplémentationsdeméthodes
d’élémentsfinisd’ordreélevépourlesfluidescompressiblesdevéloppéesàl’ONERA
etàl’INRIAdepuisenvironunandanslesbibliothèquesAGHORAetAEROSOL.
Mots-clés: Finiteelements,parallelperformance,compressibleflows
COLARGOL 3
Contents
1 DiscontinuousGalerkinmethodsandlibraries 4
1.1 TheEulermodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Runge-KuttadiscontinuousGalerkinformulation . . . . . . . . . . . 5
1.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4 Remarksonquadratureformulas . . . . . . . . . . . . . . . . . . . . 7
2 TheAGHORAandAEROSOLlibraries 7
2.1 TheAGHORAlibrary . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 TheAEROSOLlibrary . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.1 Underthebonnet . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2.2 ThePAMPAlibrary . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Comparisonofimplementations . . . . . . . . . . . . . . . . . . . . 11
2.3.1 Howtotakeintoaccountthegeometry . . . . . . . . . . . . 11
2.3.2 Computationofbasis . . . . . . . . . . . . . . . . . . . . . . 11
2.3.3 Communicationhandling . . . . . . . . . . . . . . . . . . . . 12
3 Comparisons 12
3.1 Numericaltests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2 Yeevortex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Introduction
Ontheneedforhighorderinaerodynamics
Overthelastthreedecades,mostofthesimulationsfordesigningaircrafthavebeenre-
ducedtosecondorderfinitevolumesmethods,mainlywithRANS(ReynoldsAveraged
Navier-Stokes)modelsforturbulence. Theusualtrickforhavinganaccuratesolution
withaRANSsimulationconsistsinperformingasequenceofcomputationsonmore
and more refined meshes. However, the RANS approach can be inefficient for simu-
latingsomeproblemsthatareintrisicallytimedependent,becauseitisatimeaveraged
model,andthusstationary. Fortimedependentsimulations,themostaccuratemethod
is called DNS (Direct Numerical Simulation) and consists in meshing the domain at
theturbulentscale,withoutanyturbulencemodel. However,duetotherequiredmesh
fineness,DNSisoftenconsideredonlyforlowReynoldsnumbers(sincethenumberof
cellsincreasesasRe9/4)and/orforacademicconfigurations,forwhichefficientfinite
differencesmethodscanbeused. Anintermediatesolution,betweenDNSandRANS
approachconsistsinapplyingaspatialfiltertotheNavier-Stokesequations.Theresult-
ingmodelistimedependent, butshallbeclosedfortakingintoaccountsmallscales.
ThismethodiscalledtheLargeEddySimulation(LES)approach. ForbothDNSand
LES, the accuracy yielded by adaptive mesh refinement is harder to implement be-
causetheproblemistimedependent. Itwouldthereforeimplytheoreticallytoperform
ameshadaptationateachtimestep,whichisverycostly.Italsoraisesissuesregarding
dynamicloadbalancingasfarasparallelcomputingisconcerned.
Howtoderiveahighordermethod
Mathematicallyspeaking,asufficientconditionforhavingahighorderapproximation
of a function is to find a piecewise polynomial approximation of this function, for
RRn°8200
4 D.A.Mbengoueetal
examplebyinterpolationorL2-projection. Oncethetypeofapproximationischosen,
anumericalschememustbederivedfromthisapproximation. Thefollowingmethods
canbeconsidered.
• Highorderfinitevolumemethods
In finite volume approximations, the solution is considered as piecewise con-
stant. For getting a high order method, a high order polynomial representation
can be computed, by interpolating the values on the neighboring cells. For ex-
ample, in one dimension, a second order polynomial can be computed on one
cellbyinterpolatingthevaluesontheleftandrightcells. Thiscanbegeneral-
ized to higher dimensions and to higher order polynomial approximations, and
thepricetopayforhavingahighorderrepresentationistovisitmoreandmore
neighbors. Moreover, in a hyperbolic framework, the aim is to rebuild a non
oscillatoryinterpolation[?],whichcomplicatestheproblem. Thismethodisnot
compact,andthereforeisnotwellsuitedtoparallelcomputing.
• ContinuousGalerkinmethod
Inthismethod,thesolutionisapproximatedbypiecewisecontinuouspolynomial
functions. ThenumericalschemeisthenobtainedbywritingtheL2orthogonal-
itybetweentheapproximationbasisandtheequationprojectedontheapproxi-
mation space. It is well adapted to some problems such as purely parabolic or
elliptic problems. For Euler equations, and more generally for hyperbolic sys-
tems (even linear), this method is known to be unstable. This method can be
stabilized, for example with the SUPG method [?]. Nevertheless, this method
dependsonparametersthatcanbehardtotune.
• ResidualDistributionschemes
ThedevelopmentoffluctuationschemesbeganwiththeearlyworkofPhilRoe[?].
Residual distribution is a weighted residual approach in which local discrete
equations are derived as a sum of elemental contributions proportional to ele-
mentintegralsoftheequations(thecellresiduals)viamatrixweights.Thebasics
ofthemethodarethoroughlydiscussedin[?]. Asincontinuousanddiscontin-
uous finite element methods, the key toward high accuracy is the use of a high
orderpolynomialrepresentationoftheunknownsinthecomputationofthecell
residuals[?]. Althoughthisapproachhasshowngreatpotentialinsteadyappli-
cations [?, ?], its current state of the art [?] shows that further work is needed
tobringthemethodtothelevelofmaturityofmorepopulartechniques,suchas
DiscontinuousGalerkin.
• DiscontinuousGalerkinmethod
ThedevelopmentofdiscontinuousGalerkinfornonlinearhyperbolicequations
began in [?]. Its stabilization for flows with shocks was developed in the 90’s,
mainlybyCockburnandShu(see[?]forareview). Inthesametime,asolution
fordealingwithNavier-Stokesequationswasproposedin[?]. Inspiteofitscost
(ithasmuchmoredegreesoffreedomthanclassicalcontinuousfiniteelements
methods), it is attractive because it is naturally L2 stable for linear problems,
because a cell entropy inequality can be proven [?], and also because it has a
compactstencilsothatithasagoodbehaviorinparallelenvironment[?]. The
success of these methods lies in their flexibility thanks to their high degree of
locality. ThesepropertiesmaketheDGmethodwellsuitedtoparallelcomput-
ing,hp-refinement,hp-multigrid,unstructuredmeshes,theweakapplicationof
boundaryconditions,etc.
Inria
COLARGOL 5
The development of different high order methods for aerospace application was the
topicoftheEuropeanprojectADIGMA,andwereferto[?]forrecentdevelopmentson
thistopic.
1 Discontinuous Galerkin methods and libraries
1.1 TheEulermodel
Let Ω ⊂ Rd be a bounded domain where d is the space dimension and consider the
followingproblem
∂ u+∇·f(u)=0, inΩ×(0,∞) , (1)
t
with initial condition u(·,0) = u (·) in Ω and appropriate boundary conditions pre-
0
scribed on ∂Ω. The vector u = (ρ,ρv,ρE)(cid:62) represents the conservative variables
with ρ the density, v ∈ Rd the velocity vector and E = ε(p,ρ)+v2/2 where ε is
the specific internal energy. Here, we suppose that the fluid follows the perfect gas
equation
p
ε(p,ρ)= ,
(γ−1)ρ
wherepdenotesthepressureandγistheratioofspecificheats. Thenonlinearconvec-
tivefluxesin(1)aredefinedby
 
ρv(cid:62)
 
f(u)= ρvv(cid:62)+pI  . (2)
 
 
(ρE+p)v(cid:62)
Theproblem(1)ishyperbolicprovidedtheconservativevariablevectortakesvaluesin
thesetofadmissiblestates
v2
Ψ ={u∈Rd+2 : ρ>0, E− >0} . (3)
2
1.2 Runge-KuttadiscontinuousGalerkinformulation
ThediscontinuousGalerkinmethodisafiniteelementmethodinwhichtheweakfor-
mulationisprojectedonaspaceofdiscontinuouspiecewisepolynomialsoftheprob-
lem (1). The domain Ω is partitioned into a shape-regular mesh Ω consisting of
h
nonoverlapping and nonempty elements κ of characteristic size h := min{h ,κ ∈
κ
Ω }whereh isad-dimensionalmeasureofκ.WedefinethesetsE andE ofinterior
h κ i b
andboundaryfacesinΩ ,respectively.
h
Welookforapproximatesolutionsinthefunctionspaceofdiscontinuouspolyno-
mialsVp = {φ ∈ L2(Ω ) : φ| ◦F−1 ∈ Pp(κˆ), ∀κ ∈ Ω },wherePp(κˆ)denotes
h h κ κ h
the polynomial space associated to the reference element κˆ corresponding to the ele-
mentκ. Eachphysicalelementκistheimageofoneofthefollowingreferenceshapes
κˆ through the mapping F : simplex (line, triangle or tetrahedron), tensor elements
κ
(quadrangle, hexaedron), prism or pyramid. The space Pp(κˆ) might be composed of
the set of polynomials with degree lower or equal to p in the case of simplex, or of
tensorproductofonedimensionalpolynomialswithdegreelowerorequaltopinthe
caseoftensorelements, orofatensorofonedimensionalbasisandtwodimensional
RRn°8200
6 D.A.Mbengoueetal
e
v - v +
h h
κ+
κ-
n

Figure1:Innerandexteriorelementsκ+andκ−anddefinitionoftracesv±ontheinterfaceeandofthe
h
unitoutwardnormalvectorn.
basis on a triangle if κˆ is a prism, or a conical product of one dimensional basis and
two dimensional basis on a quadrangle if κˆ is a pyramid. The numerical solution of
equation(1)issoughtundertheform
(cid:88)Nκ
u (x,t)= φl(x)Ul(t), ∀x∈κ, κ∈Ω , ∀t≥0 , (4)
h κ κ h
l=1
where (Ul) are the degrees of freedom in the element κ. The semi-discrete
κ 1≤l≤Nκ
formofequation(1)reads: findu in[Vp]d+2suchthatforallv inVpwehave
h h h h
(cid:90)
v ∂ u dx+B (u ,v )=0 . (5)
h t h h h h
Ωh
Hereafter,wewillusethenotation[[φ]]=φ+−φ−whichdenotesthejumpoperator
definedforagiveninterfacee ∈ E . Here,φ+ andφ− arethetracesofanyquantityφ
i
ontheinterfaceetakenfromwithintheinterioroftheelementκ+ andtheinteriorof
theneighboringelementκ−,respectively(seeFigure1).
Thespacediscretizationinequation(5)canthenbewrittenas
(cid:18)un+1−un(cid:19)
M h h +B (un,v )=0 . (6)
κ δt h h h
fortheexplicitEulertimestepping,theextensiontootherexplicittimesteppingbeing
straightforward. In(6)M denotesthelocalmassmatrixdefinedas
κ
(cid:90)
∀(i,j)∈[0;N ]×[0;N ] M(i,j) = φi(x)φj(x)dx ,
κ κ κ κ κ
κ
andthespacediscretizationoperatorB isdefinedby
h
(cid:90)
B (u ,v ) = − f(u )·∇v dV
h h h h h
Ωh
(cid:90)
+ [[v ]]˜f(u+,u−,n)dS
h h h
Ei
(cid:90)
+ v+f(u (u+,n))·ndS , (7)
h b h
Eb
Inria
COLARGOL 7
wherendenotestheunitoutwardnormalvectortoanelementκ+ (seeFigure1)and
u isanappropriateoperatorwhichallowstoimposeboundaryconditionsonE . The
b b
numerical flux ˜f may be chosen to be any monotonic Lipschitz function satisfying
consistency,conservativityandentropydissipativityproperties(see[?]forinstance).
Throughoutthisstudy,thesemi-discreteequation(5)isadvancedintimebymeans
ofanexplicittreatmentwiththird-orderstrongstabilitypreservingRunge-Kuttameth-
ods[?,?].
1.3 Algorithm
AsanexplicitRunge-Kuttatimesteppingwaschosen,thealgorithmconsistsincom-
putingthespatialresidual,andtheninupdatingtheintermediate(orthefinal)stepsof
theRunge-Kuttamethodbyinvertingthemassmatrix. InthediscontinuousGalerkin
method,themassmatrixisblockdiagonal. Itcanactuallybediagonalprovidedanor-
thogonalbasisisusedoneachelement. Wepointoutthatthisstepofthemethoddoes
notrequireanycommunicationinadistributedmemoryenvironmentifallelementsare
onthesameprocess(whichisalwaysthecase). Wenowdwellonthedifferentloops
thatmustbemadeforcomputingthelocalspatialresidualdefinedon(7). Wedetailthe
computationsneededfortheelementsloop,theotherloopsbeingsimilar.
Weareinterestedincomputing
(cid:90)
f(u )·∇v dV
h h
κ
forallbasisfunctionv ofκ. UsingthedefinitionofthebasisgiveninSection1.2
h
(cid:90) (cid:90)
f(u )·∇φidV = f((cid:80)Nκ φl(x)Ul)·∇φidV .
h κ l=1 κ κ κ
κ κ
AsinSection1.2,wedenotebyF themapfromκtoκˆ. Intheintegral,wechangethe
κ
variablexbyx=F (xˆ),sothat
κ
(cid:90) (cid:16) (cid:17) (cid:90) (cid:16) (cid:17)
f (cid:80)Nκ φl(x)Ul ·∇φidV = f (cid:80)Nκ φl ◦F (xˆ)Ul ·DF−1∇φˆi |detDF | dxˆ
l=1 κ κ κ l=1 κ κ κ κ κˆ κ
κ (cid:90)κˆ (cid:16) (cid:17)
= f (cid:80)Nκ φˆl(xˆ)Ul ·DF−1∇φˆi |detDF | dxˆ .
l=1 κˆ κ κ κ κ
κˆ
An approximated quadrature formula is used for computing this integral. We denote
byxˆ thequadraturepointsandbyω theweights.
α α
(cid:90) f(cid:32)(cid:88)Nκ φl(x)Ul(cid:33)·∇φidV ≈(cid:88)ω f(cid:32)(cid:88)Nκ φˆl(xˆ )Ul(cid:33)·DF−1(xˆ )∇φˆi(xˆ )|detDF (xˆ )|
κ κ κ α κˆ α κ κ α κ α κ α
κ l=1 α l=1
1.4 Remarksonquadratureformulas
Forlinearproblems,quadratureformulascanbechosenforbeingexact. Fornonlinear
problems,[?]suggests,foranapproximationofdegreep,totakea(2p)thorderformula
forcells,anda(2p+1)thorderforfaces.
For hypercube shapes, the optimal set of points (i.e. the one ensuring the highest
degree with a given number of points) is the Gauss set of points. For other shapes,
the optimal set of points is often unknown. A systematic way of deriving quadrature
RRn°8200