High Order Methods with Adaptive Mesh Refinement for the Solution of - - PowerPoint PPT Presentation

High Order Methods with Adaptive Mesh Refinement for the Solution of the Relativistic MHD equations

Collaborator: Michael Dumbser, University of Trento

Olindo Zanotti

University of Trento 14th International Conference on Hyperbolic problems: Theory, Numerics, Applications. June 25th 2012.

Outline

• Part I: Astrophysical motivations for considering

Resistive Relativistic MHD

• Part II: The mathematics of Resistive Relativistic

MHD: a hyperbolic system with stiff source terms

• Part III: A numerical scheme for RRMHD
• High order reconstruction
• Using local space-time Discontinuous Galerkin for treating

stiffness in the source terms

• Construction of 1-step time evolution schemes
• Adaptive Mesh Refinement
• Part IV: High Lundquist number relativistic

magnetic reconnection

Why Resistive MHD?

In most cases, the ideal MHD approximation of infinite conductivity is correct: magnetic Reynolds and Lundquist number

ratio between diffusion timescale and advection timescale

In some cases, however, such an approximation is completely wrong and resistivity must be taken into account. This is particularly the case when magnetic reconnection takes place

1 1,       L v S L v R

A m

(Lyubarski 2005)

B J B B v B

t

                   ) (

2



Where relativistic magnetic reconnection?

 Giant flares in gamma-ray repeaters (Lyutikov 2003, 2006)  Current sheets at the Y point in pulsars magnetospheres  Dissipation of alternating fields at the termination shock of a relativistic pulsar wind (Petri & Lyubarski 2007)  Gamma-ray burst jets, where particle acceleration induced by magnetic reconnection is supposed to take place (Drenkhahn & Spruit 2002 , McKinney & Uzdensky 2010) Accretion disc coronae of AGN where magnetic loops emerge from the disc via buoyancy instability (di Matteo 1998, Jaroschek 2004)

The Relativistic plasma

• The energy-momentum tensor of a relativistic plasma is:
• The Faraday EM tensor satisfies Maxwell’s equations
• Unlike ideal MHD, the current is not a derived quantity, and it is

provided after specifying Ohm’s law

• In ideal MHD one simply has

2 2 2 2

4 / ) ( ) ( dz dy dx dt dx dx g g F F F F pg u u h T T T

em m

          

              



    

    

F J F

] ) ( [ v v E B v E v J

c

                

t B c E B J c t E c B E                              1 , 4 1 , 4  

 

2 /

  

 F F 



B v E      

assume 3+1 formalism

Ohm’s law

] ) ( )[ ( ] ) ( [

2

v v B E v B B E v v E B v E v q J                                 

The electric conductivity is actually a tensor

time collision , 1 ,

2

                     m e m e

) (

          

      b u b b g e u J

e

    

RRMHD ideal RMHD

) (                                             

i i t c i i t i i t i i k j ijk i t i k j ijk i t i i t i j i j t i i t

J q E B J B e E E e B S U Z S Dv D   

 

i j j i j i j i i j k j ijk i i

E B p E E B B v v h Z E B p h U B E v h S D       2 / 2 / 2 / 2 /

2 2 2 2 2 2 2

                

• where

                            

i i t k j ijk i i k j ijk i t i i t i j i j t i i t

B B v e E E e B S U Z S Dv D ) (

Hyperbolic Divergence Cleaning approach of Dedner et al. JCP, 175, 645, 2002

Of course all of this holds also in curved spacetime

) ( 2 1

2 2 2 2

B E p W h U B E v W h S W D               

+ evolution equations for the electromagnetic field

RRMHD: a system with potential stiff source terms!

] ) ( [ v v E B v E v q J J B E

t

                             

In the limit of the electric field evolves with a hyperbolic equation, not a parabolic one like in classical MHD

 

In the limit of the source becomes stiff

  

Traditional approaches:

• 1. Strang-splitting, i.e. operator splitting, (Komissarov 2007,

Zenitani et al. 2009)

• 2. Implicit-Explicit Runge Kutta (Palenzuela et al 2009)

1 ) (

2

        B B v B

t

     

d x t

u s u f u    

a

) ( 1 ) (     

Timescale of dissipative process Timescale of advection

Computational strategy

• 1. Apply a reconstruction operator to the Discontinuous

Galerkin scheme at the beginning of each time step

• 2. Provide

the time evolution

• f

the reconstructed polynomials by solving the local space-time Discontinuous Galerkin predictor step

• 3. Solve the Riemann problem by some flux formula
• 4. Perform a one-step time update from time level n to n+1

with quadrature (or quadrature free) formulation

Reconstruction: the scheme (I)

         ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (

) ( 1 ) ( 4 ) ( 1 ) ( 3 ) ( 1 ) ( 2 ) ( 1 ) ( 1 ) ( 4 ) ( 1 ) ( 3 ) ( 1 ) ( 2 ) ( 1 ) ( 1 ) ( 4 ) ( 1 ) ( 3 ) ( 1 ) ( 2 ) ( 1 m m m m m m m m m m m m m m m m m m m m m

Z Z Z Z Z Z Z z Y Y Y Y Y Y Y y X X X X X X X x                     

) ( 1 m N m

Q Q

E



  

) , (

) (

   

m

Q x x 

Set up a triangulation of space:

Physical coordinates reference coordinates

z

x

y

1 2 3 4 4 3 1 2   

M NP

P

] 1 , [ , ,    

M NP

P

) ( ˆ ) ( ) , (

1 n n l L l l n h

t u t u

N

   





  ) ( ˆ ) ( ) , (

1 n n l L l l n h

t w t w

M

   





 

At time t^n the vector of conserved quantities is represented by piecewise polynomials of degree N

l



l



and coincide up to order N, while

N M 

n m n m d

i

Q n m

    



, , ,  

From it, we reconstruct over polynomials of degree M :

) ( ) 1 ( ! 1 d M M d LM     

Choose

l



among Legendre polynomials

Number of degrees of freedom

Reconstruction: the scheme (II)

Traditional Galerkin representation

M NP

P

The reconstruction is obtained through L2 projection of the unknown polynomials Weak form of the identity of the reconstructed solution of degree M and the original numerical solution of degree N:

) ( m i Q h k Q h k

S Q d u d w

i i

    

 

   

Computed with Gaussian quadrature Computed analytically

As many equations as the element of the stencil. For stability reasons, the number of elements in the stencil is chosen larger than Solved using the constrained least squares technique Choose a stencil for the reconstruction:



e

n k k j m

T S

1 )) ( ( ) ( 



The number of cells in the stencil is taken as

Reconstruction: the scheme (III)

n

N m l m l

L l u w    1 for

) ( ) (

  ) / CEILING( /

N M N M

L L n L L n  

To recap:

Reconstruction is applied to polynomials of degree N and generates polynomials of degree M

class hybrid Galerkin classical FV classical       N M M N N

Dumbser et al (2008) JCP, 227, 8209

Time evolution of the reconstructed polynomial

) ( ) ( W S W F t W      

 

      S W F W ) (





x J W tS S J W tF F

T

        

 

 ), ( , ) (

 

      S W F W

k h h k h k

, ) ( , ,    



) ˆ ( ˆ ) ˆ ( ˆ ˆ

l l l l l l h l l l l l l h l l l h

W S S S W F F F W W

     

    

         

Use the local space-time Discontinuous Galerkin approach

      d d g f g f

Q Q

 ) , ( ) , ( ,

1

 



Multiply by the space-time test function and integrate over the space-time reference control volume

) , (     

k k 

 

    



 d g f g f

Q Q 

 ) , ( ) , ( ,

….after integration by parts in time:

   

 

        S W F W w W

k h h k k h h k h k

, ) ( , , , ,

1

     



This polynomial is obtained by the high order reconstruction. At this point we use the expansions over the polynomial basis

M NP

P

   

  

              

l l k l l k n l l k l k l l k

S F w W ˆ , ˆ , ˆ , ˆ , ,

1

                               



1

K

F



K

M

i l n l i l i l

F K K w F K S M K W

, 1 1 1 1 1 , 1 1 1

ˆ ) ( ˆ ) ( ˆ ) ( ˆ

      

  



The unknown

….the source term is taken implicitly as:

In practice, the relevant source term will involve the derivative of the electric current with respect to the primitive variables!

 

i l i l i l i l

W W W S t S S ˆ ˆ ˆ ˆ

1 , 1 ,

     

    1 

             V W V S W S

We solve the algebraic equation by first looking for a first guess at first order:

n l

w W S W ˆ ) (  



1 1

1

   M K K



The resulting cell average is used as initial guess in the full equation

W

i l n l i l i l

F K K w F K S M K W

, 1 1 1 1 1 , 1 1 1

ˆ ) ( ˆ ) ( ˆ ) ( ˆ

      

  



Once we have

l l l h

W W ˆ ) , ( ) , (



     

we update in time through

 

i i i n i n i

S t f f x t W W       

   2 / 1 2 / 1 1

with

  

 

  1 1 1 1 2 / 1

)) , ( ( )) , ( ), , 1 ( (        d d W S S d W W f f

i i i i h i



 

 

2 / 1 2 / 1

) , ( 1

i i

x x n n i

dx t x W x W

Computed with Gaussian quadrature One step time update!

Corrector: one step time update

Adaptive Mesh Refinement: I

We followed the standard approach by Berger-Oliger-Colella, adopting a set of hierarchical RECTANGULAR grids. CALL ComputeFluxes(reflev) CALL MPIExchange_duh CALL Update(reflev) ! CALL EvolveVirtual(reflev) CALL Average(reflev) ! IF(reflev.EQ.0) THEN CALL ComputeTimeStep(amax) DO ireflev = 1, CurrentRefLevel CALL TimeEvolution(ireflev) CALL Project(ireflev) CALL MPIExchange_qh ENDDO ENDIF

Mignone et al. ApJ (2012)

Adaptive Mesh Refinement: II

Relativistic magnetic reconnection

Initial conditions (relativistic Harris model)

) 2 tanh( ) 2 ( cosh 1 ) 2 ( cosh 1 2 1

2 2 2

         

y x y m m m

v v x B B x p p x B p       

Anomalous resistivity

Zanotti & Dumbser MNRAS (2012)

Dynamical evolution

Two blobs of matter are ejected along the y direction while the magnetic field dissipates in a X type topology

Zanotti & Dumbser MNRAS (2012)

Plasma acceleration

Plasma acceleration

At high magnetizations the dependence is in agreement with the theoretical prediction:

2 / 1 m

  

Zanotti & Dumbser MNRAS (2012)

High Lundquist number magnetic reconnection

What happens when the background resistivity is decreased? (the resistivity jump between the background and the anomalous spot is increased?)

Zanotti & Dumbser MNRAS (2012)

High Lundquist number magnetic reconnection

What happens when the background resistivity is decreased? (the resistivity jump between the background and the anomalous spot is increased?) A tearing instability is produced…..

Zanotti & Dumbser MNRAS (2012)

High Lundquist number magnetic reconnection

Newtonian regime (Samtaney et al 2009): Relativistic regime:

8

10 ~

c

S S 

4

10 ~

c

S S 

Zanotti & Dumbser MNRAS (2012)

Conclusions

Future directions: inclusion of more physical Ohm’s laws and generalization to curved spacetimes in general relativity.

• This is a promising period for numerical relativistic MHD. Having

included resistivity allows for a larger class of problems to be studied.

• Disc accretion onto Black Holes
• Modelling of current sheets around rotating neutron stars
• References:
• Dumbser, M., Balsara, D., Toro, E., Munz, C., JCP, 227, 8209, (2008)
• Dumbser, M., Enaux, C., Toro, E., JCP, 227, 3971, (2008)
• Palenzuela, C., et al., MNRAS, 339, 1727, (2009)
• Dumbser, M., Zanotti, O., JCP, 228, 6991, (2009)
• Zenitani, Hesse, Klimas, ApJ, 696, 1385, (2009)
• Zanotti, O. Dumbser, M., (2011) arXiv:1103.5924