Transcript Document

Athena中算法与GLM九个方程
刘静静
4月7号
主要内容





阶段任务
背景(CESE+HLL)
Athena中的不同算法
GLM-MHD九个方程
参考文献
阶段任务
对程序CESE（in the near sun)and HLL(off-sun)的一些修改，从以下三个方
面：
 在off-sun部分使用参考Athena code中的不同算法试验
 两部分使用不同的参数无量纲化，在off-sun部分使用6Rs处的太阳分参
数无量纲化
 考虑将MHD方程组改成GLM九个方程，模拟太阳风
背景(CESE+HLL）




从太阳到地球的计算区域，分成了近太阳和远太阳两部分区域
在近太阳区域，使用阴阳网格和CESE算法，在远太阳区域使用自适应（AMR）
网格和HLL算法
阴阳网格可以避免奇点和在极区网格集中问题，更好地在太阳表面达到高精度
的解；AMR可以自动的捕捉等离子体流的特征（far-field)，例如日球层电流片
和激波，并且可以节省计算资源
控制方程：
 t     u  0


1


 t u      uu   p  B12  B1  B 0   B1B1  B1B 0  B 0 B1  
2




j0  B 0   g Ω (Ω r )  2 Ω  u  B(  B )


1

 t e     e  p  B12  B1  B 0 u  u  B1 B1  B 0  
2



B
 B1  1  E  j0  u  g Ω(Ω r )  Qe  u  B(  B )
t
 t B1    uB  BU    t B 0  u  B     B 
其中 E  u  B, j0    B0 ,磁场分成背景磁场和扰动磁场
B  B0  B1
Athena
Athena：



The equations: ideal MHD
The numerical algorithms in Athna are based on directionally unsplit,higher order
Godunov methods
Discretization 1)mass,momentum,energy:finite volume





t  n1 2
t  n1 2
t  n1 2
U in, j ,1k  U in, j ,k 
Fi 1 2, j ,k  Fi n112,2j ,k 
Gi , j 1 2,k  Gin, j 1122,k 
H i , j ,k 1 2  H in, j ,1k 21 2
x
y
z
n
z k 1 2 y j 1 2 xi 1 2 
1
    U x, y, z, t n dxdydz
Volume-averaged
with U i , j ,k 
xyz zk 1 2 y j1 2 xi1 2
 n 1 2
t n1 z k 1 2 y j 1 2 
1
Fi 1 2, j ,k 

F xi 1 2 , y, z, t dydzdt
tyz t n zk 1 2 y j1 2
time- and area n 1 2
t n1 z k 1 2 xi 1 2 
1
Gi , j 1 2,k 

G x, y j 1 2 , z , t dxdzdt averaged fluxes
txz t n zk 1 2 xi1 2
 n 1 2
t n1 y j 1 2 xi 1 2 
1
H i , j ,k 1 2 
  n   H x, y, zk 1 2 , t dxdydt
txy t y j1 2 xi1 2
2)magnetic field:based on area rather than volumes averages








Athena
 n1
n
t n1 2
t
Bx ,i 1 2, j ,k  Bx ,i , j ,k 
 z ,i 1 2, j 1 2,k   zn,i 1122, j 1 2,k   yn,i1122, j ,k 1 2   yn,i1122, j ,k 1 2
y
z

with


Bxn,i 1 2, j ,k 
n 1 2
x ,i , j 1 2 , k 1 2



1 zk 1 2 y j1 2
Bx xi 1 2 , y, z, t n dydz


z
y
yz k 1 2 j1 2


area-averaged
1 t n1 xi1 2
electromotive force averaged along




x
,
y
,
z
,
t
dxdt
x
j

1
2
k

1
2
n
the appropriate line element
xt t xi1 2
advantages:1) ideal for use AMR
2)Superior for shock capturing and evolving the contact and rotational
discontinuties
Athna
The chart for the steps in the 2D algorithm in Athena
Athena
The algorithm for computing MHD interface states:piecewise contant(first-order)
reconstruction,piecewise linear(second-order)resconstruction,piecewise parabolic(thirdorder) reconstruction
 The algorithm for computing fluxes:HLL solvers,Roe’s method
Remarks：1)the reconstruction used in Athena require characteristic variables and a
characte-ristic evolution of the linearized systerm
2)The Godunov methods do not require expensive solvers based on complex characteristic
decompositions

Data reconstruction

Piecewise constant reconstruction: assume the primitive variables are piecewise
constant within each cell
u L,i 1 2  ui 1, u R,i 1 2  ui

Piecewise linear reconstruction: assume the primitive variables vary linearly within
each cell
x  xi
ui  x   ui 
where
ui is
ui
x
a limited slope for the cell,two types of limiter as follow:
1


ui  min mod ui 1  ui , ui 1  ui 1 , ui  ui 1 ,  1,2
2


ui 

2ui 1  ui ui  ui 1   
ui 1  ui   ui  ui 1 
2
2
1
ui1  ui 1 
 2
WENO reconstruction: can achieve higher than second order
Basic ideal:several cells can formulate a module S kr (r denotes the number of cells
formulated the module,k denotes the total number of modules,different modules have
different interpolation polynomials Pk x 
Data reconstruction

the total polynomial R(x) of reconstruction is a convex combination of the
above polynomials Pj(x),
R x    w P  x 
where w j is weight cofficient, w  0,  w  1
k 1
j
j 0
j
k 1
j
j 0
j
The WENO reconstruction can have (2k-1) order, and is non-oscillatory,but the
computation is complex
Godunov Fluxes


First proposed by Godunov S.K. in 1959
The basic idea: at t n ,in each cell the primitive variables are constant. At the
interface bewteen the neighbor cells xi 1, xi and xi , xi 1  , there is a initial discontinuity
uin1 2 , x  xi
ui  x    n
ui 1 2 , x  xi
then formulate a local Riemann problem bewteen the neighbour cells.
 Godunov methods do require expensive solvers based on complex characteristic
decompositions and capture high quality shock
 HLL-family solvers:
 f L , if 1  0

F  f k , if k k 1  0
 f , if   0
N
 R
Roe’s method



An useful linearization for the MHD equations
Include all the characteristics of the systerm,and less diffusive and more accurate
for intermediate waves
Jacobian is evaluated using an average state(Roe average)
  L R
    
H   H   H      
B   B   B      
v

L vL  R v R
L
L
L
L
R
R
where H  E  P   is the enthalpy
the Roe fluxes are simply:
L
R
R
R
L
L
R
R
1

F Roe   f Lf R   k L k  U L  U R R k 
2
k

Disadvantage: may return negative densities or pressures
HLL

assuming an average intermediate state between the fastest and slowest
waves
are the minimum signal speed
and the maximum signal speed

intermediate state

S RU R  S LU L  FR  FL
U 
SR  SL
the HLL fluxes
SL , SR

FHLL
 FL if SL  0

 F if SL  0  SR
F if S  0
R
 R
F 
S R FL  S L FR  S R S L U R  U L 
SR  SL
HLL
Remarks:
 must be estimated appropriately
Davis
S R  max l U L , l U R , S L  min m U L , m U R 
Einfeldt et al

 
 min  U ,  U 
S R  max l U L , l U Roe
SL
Roe
m
m
R
 The solver is fast and do not need the characteristic decomposition
too diffusive and cannot resolve isolated contact discontinuities very well
HLLE



Using a singal constant intermediate state computed from a conservative average
Do not require a characteristic decomposition of MHD equations
The HLLE flux at the interface:
HLLE
i 1 2
F



b  f L,i 1 2  b f R,i 1 2
b  b


b b 
ui  ui 1 
 
b  b
where f L,i 1 2  f u L,i 1 2 , f R,i 1 2  f u R,i 1 2 are the fluxes evaluated using the left and
right states of the conserved variables,and



 

b   max max M , vR  c f R ,0 , b   min min l , vL  c f L ,0
If both M  0 and vR  c f
R
 0 (or L  0 and vL  c f L  0
),the HLLE flux will be f L,i 1 2 or f R,i 1 2 
the HLLE can guarantee the pressure and density is positive,but in the
multiple dimensions,it does not necessarily guarantee.What’s more,the HLLE
neglects the Alfven,slow magnetosonic,and contact waves.

HLLC

the intermediate states in the Riemann fan are separated into two intermediate
states by a contact discontinuity can resolve isolated contact discontinuities
exactly
S M be evaluated from HLLaverage
SM 


vn   S R  vnR  RvnR  S L  vnL  LvnL
S R  vnR  R  S L  vnL  L

 pT1R  pT1L
S R  vnR  R  S L  vnL  L
the numerical flux of HLLC
FHLLC
 FL if SL  0
 F  if S  0  S

L
M
  L
 FR if SM  0  S R

 FR if SR  0

 S U
FL  FL  S L U L  U L
FR  FR
L

R
U R


HLLC

Positively conservative
S L  uL 
a:the soud speed

p
 1
 1
aL , S R  u R 
aR
2
2

HLLC can dramatically improve the results of the HLL solver, and has much less
computational time than the HLLEM
HLLD


Five-wave Riemann solver for
MHD,HLL-Discontinuities solver
Composed of four intermediate
states
indicate the speeds of the
fast magnetosonic waves, Alfvén
waves, and entropy wave

 min v
S , S , S M


S R  max vnR  c f nR , vnL  c f nL
SL
S R  S M 
nR
Bn
 R
 c f nR , vnL  c f nL
, S L  S M 
Bn
 L
SM


vn 



S R  vnR  R vnR  S L  vnL  LvnL
S R  vnR  R  S L  vnL  L
 pT1R  pT1L
S R  vnR  R  S L  vnL  L
HLLD

The numerical flux vector of the HLLD Riemann solver for MHD equations
FHLLC
 FnL if
 F  if
 nL


 FnL if
  
 FnR if
 FnR if


 FnR if
SL  0
SL  0  S L
SL  0  S M
SM  0  S R
SR  0  S R
SR  0
 Positively conservative
S L  uL 
 1
 1
c f , S R  uR 
cf
2
2
c f L , c f R :the fast magnetosonic speed
L
R
两部分无量纲化

在off-sun部分使用6Rs处太阳风参数做无量纲化Cartesian部
分，利用初始时的流场和磁场参数Parker解在6Rs处的值
Call parker（r6rs,ur,gamma0,T0)

进一步改进：在给定的网格上，指定半径处的太阳风参数，
分两部分无量纲化
subroutine nearpoint(r)
subroutine
nearpoint(r)
给定的半径
距离最小的点
利用Parker解求出所得到点的太阳风参数，进行无量纲化
GLM-MHD



GLM(generalized Lagrange multiplier)
The form of equations
Solver for GLM-MHD
GLM



Coupling the divergence constraint by introducing a generalized Lagrange
multiplier
the divergence errors are transported to the domain boundaries with the
maximal admissible speed and are damped at the same time
Magnetic induction equations are replaced:





 t B    uBT  Bu T    0

D     B  0
Different choices for the linear operator D
D   0
Elliptic correction:
Parabolic correction: D   1  , c  0, 
p
c 2p
Hyperbolic correction: D   1   , c  0,  
t
h
2
ch
Combination of parabolic and hyperbolic ansatz:
D  
1
1
  2
2 t
ch
cp
方程组形式

MHD方程组变为

 t    x ( u x )  0
  T 
1  2   T 

 t u      u u   p  B   BB   0
2 








 t B    u B T  Bu T   0



1 2     
 t e   e  p  B u  B u  B   0
2 



ch2
2
 t  ch   B   2 
cp
方程形式


The MHD equations can be symmetrized by adding some hyperbollic terms on the
right-hand side
The changed eqations:

 t    x ( u x )  0
 
  
1    

 t u      u u T   p  B 2   BB T     B B
2 






 
 T
T
 t B    u B  Bu     B u


  

1      
 t e   e  p  B 2 u  Bu  B     B u  B   
2 


2

c

2
 t  ch   B  u     h2 
cp
Remarks : 1)call the equations the extended GLM(EGLM) formulation of MHD
equations
2) c p significantly depends on the grid size and the scheme used, c p is a
function of ch
方程的特征值


The eigenvalues of the GLM-MHD coinside with the ordinary MHD waves plus
two additional modes  ch ,for a total of nine characteristic waves
For one dimensional ,x direction
1,9  ch , 2,8  vx  c f , 3,7  vx  ca , 4,6  vx  cs , 5  vx

ch  max vx  c f , x , v y  c f , y , vz  c f , z

remarks：1)show that the system is hyperbolic
2)only the waves traveling with speeds
or 
1，9 can
carry a change in
Bx
The solver of GLM-MHD


Solver for the GLM-MHD without additional source
Treat the 2 2 linear system given by the B and  from the other ordinary 7wave MHD equations in an operator-split fashion
U n1  S t 2 At S t 2 U n
where S and A are the advection and source step operators separately
1)Advection step :based on the corner transport upwind(CTU) method, second order
accurate discretization
 Fin1122  Fin1122 G nj 11 22  G nj 11 22 H nk 11 22  H nk 11 22 
U n 1  U n  t n 



x
y
z


Where F,G,H are the numerical fluxes computed by solving a Riemann problem
between suitable time-centered left and right states





Fin1122  R Vin,1 2 , Vin1,12 , G in1122  R Vjn,1 2 , Vjn11,2 , Hin1122  R Vkn,1 2 , Vkn11,2
R(·, ·) denotes the flux obtained by means of a Riemann solver, Vkn,1 2 ,Vkn11,2 are
computed via a Taylor expansion in the direction normal to a given interface

The solver for GLM-MHD
2)Source step: solver the initial value problem without the
B
term
c

 
t
c
2
h
2
p
can be integrated exactly for a time increment

 t    0  exp  
n

ch 
c
, with   h h , h  minx, y, z 
h t 
cp
Remarks:
numerical experiments indicate that the divergence errors are mininized when the  lies
in the range [0,1]

is an unphysical variable,the initial condition
step

 0 given
by the output of the most recent
condition for  :assume that the behavior of  and  at the boundary is
identical,use a homogeneous Dirichlet condition ,nonreflecting boundary condition
Boundary
参考文献

Xueshang Feng,Shaohua Zhang,Changqing Xiang,Liping Yang,Chaowei Jiang,”A
Hybrid Solar Wind Model of CESE+HLL Method with Yin-Yang Overset Grid and
AMR Grid”

Takahiro Miyoshi,Naoki Terada,”The HLLD Approximate Riemann Solver for
Magnetospheric Simulation”
Takahiro Miyoshi,Kanya Kusano,”A multi-state HLL approximate solver for ideal
magnetohydrodynamics”
A.Mignone,G.Bodo,”PLUTO:A NUMERICAL CODE FOR COMPUTATIONAL
ASTROPHYSICS”
Shengtai Li,”An HLLC Riemann solver for magneto-hydrodynamics”
James M.Stone,Thomas A.Gardiner,”ATHENA：A NEW CODE FOR
ASTROPHYSICAL MHD”
A.Dedner,F.Kemm,”Hyperbolic Divergence Cleaning for the MHD Equations”
Andrea Mignone,Petros Tzeferacos,Gianluigi Bodo,”High-order conservative finite
difference GLM-MHD schemes for cell-centered MHD






参考文献
Andrea Mignone,Petros Tzeferacos,”A Second-Order Unsplit Godunov Scheme for
Cell-Centeres MHD:CTU-GLM scheme”
 Shengtai Li,Hui Li,”A Modern Code for Solving Magneto-hydrodynamics or Hydrodynamics Equations
 Dinshaw S.Balsara,”Multidimensional HLLE Riemann Solver;Application to Euler and
Magnetohydrodynamic Flows”
