Transcript A Godunov method for ideal MHD using Constrained Transport

Studies of the MRI with a New
Godunov Scheme for MHD:
Jim Stone & Tom Gardiner
Princeton University
Recent collaborators:
John Hawley (UVa)
Peter Teuben (UMd)
Outline of Talk
A. Motivation
B. Basic Elements of the Algorithm
C. Some Tests & Applications
D. Evolution of vortices in disks
E. MRI, with comparison to ZEUS
ZEUS-like algorithms have been used successfully to
study the MRI in both 3D global and local simulations
Local simulation
Global simulation
Hawley, Gammie, & Balbus 1995; 1996; Stone et al. 1996;
Armitage 1998; etc.
These methods have been extended to study the radiation dominated
inner regions of disks around compact objects. Equations of radiation
MHD:
(Stone, Mihalas, & Norman 1992)
Use ZEUS with fluxlimited diffusion
module (Turner &
Stone 2001)
e.g. photon bubble instability in accretion disk atmospheres
QuickTime™ and a
TIFF (LZW) decompressor
are needed to see this picture.
Gammie (1998) and Blaes &
Socrates (2001) have shown
magnetosonic waves are linearly
unstable in radiation dominated
atmospheres
Turner et al. (2004) have shown they evolve into shocks in nonlinear regime:
Limitation is FLD
So why develop a new MHD code?
Global model of geometrically thin (H/R << 1) disk covering 10H
in R, 10H in Z, and 2p in azimuth with resolution of shearing box
(128 grid points/H) would be easier with nested grids.
Nested (and adaptive) grids work best with single-step Eulerian
methods based on the conservative form
MHD algorithms in ZEUS are 15+ years old - a new code could
take advantage of developments in numerical MHD since then.
Our Choice: higher-order Godunov methods combined with
Constrained Transport (CT)
B. Basic elements of the algorithm
A variety of authors have combined CT with Godunov
schemes previously • Ryu, Miniati, Jones, & Frank 1998
• Dai & Woodward 1998
• Balsara & Spicer 1999
• Toth 2000
• Londrillo & Del Zanna 2000
• Pen, Arras, & Wong 2003
However, scheme developed here differs in:
1. method by which EMFs are computed at corners.
2. calculation of PPM interface states in MHD
3. extension of unsplit integrator to MHD
Gardiner & Stone 2005
The CT Algorithm
Ez,i
1 2, j 1
Ez,i
Ez,i , j
1 2, j 1 2
Ez,i
12
E z,i
1 2, j
1, j 1 2
• Finite Volume / Godunov
algorithm gives E-field at
face centers.
• “CT Algorithm” defines
E-field at grid cell
corners.
• Arithmetic averaging: 2D
plane-parallel flow does
not reduce to equivalent
1D problem
• Algorithms which
reconstruct E-field at
corner are superior
Gardiner & Stone 2005
Simple advection tests demonstrate
differences
Field Loop Advection (b = 106): MUSCL - Hancock
QuickTime™ and a
GIF decompressor
are needed to see this picture.
Arithmetic average
QuickTime™ and a
GIF decompressor
are needed to see this picture.
Gardiner & Stone 2005
Directionally unsplit integration
CornerTransportUpwind [Colella 1991] (12 R-solves in 3D)
• Optimally Stable, CFL < 1
• Complex & Expensive for MHD...
CTU (6 R-solves)
• Stable for CFL < 1/2
• Relatively Simple...
MUSCL-Hancock
• Stable for CFL < 1/2
• Very Simple, but diffusive...
C. Some Tests & Applications
1.
2.
3.
4.
5.
6.
7.
8.
Convergence Rate of Linear Waves
Nonlinear Circularly Polarized Alfven Wave
RJ Riemann problems rotated to grid
Advection of a Field Loop (Current Cylinder)
Current Sheet
Hydro and MHD Blast Waves
Hydro implosion
RT and KH instability
See http://www.astro.princeton.edu/~jstone/tests for more tests.
Linear Wave Convergence
(2N x N x N) Grid
3D RT Instability (200x200x300 grid)
Goal: measure growth rate of fingers, structure in 3D
A test
Multi-mode in hydrodynamics
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
An application
Multi-mode in with strong B
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Codes are publicly available
Code, documentation, and tests posted on web.
Download a copy from
www.astro.princeton.edu/~jstone/athena.html
Current status:
• 1D version publicly available (C and F95)
• 2D version publicly available (C and F95)
• 3D version being tested on applications (C only)
Expect to release parallelized 3D Cartesian grid code in ~1yr
D. Evolution of vorticity in
hydrodynamical shearing disk
• HBS 1996 showed nonlinear random motions decay in
hydrodynamical shearing box using PPM
• Recent work has renewed interest in evolution of vortices in disks,
especially transit amplification of leading->trailing waves
2D incompressible
Umurhan & Regev 2004; Yecko 2004
2D compressible
Johnson & Gammie 2005
3D anelastic in stratified disks
Barranco & Marcus 2005
Our goal: Study evolution of vortices in 3D compressible disks at
the highest resolutions possible, closely following JG2005
A test: shear amplification of a single, leading,
incompressible vortex
Initial vorticity
(kx/ky)0 = 4
KE is amplified by ~ (kx/ky)2 : can code reproduce large amplification
factors?
Johnson &
Gammie 2005
Dimensionless time:
t = 1.5W t + kx0/ky
There is no aliasing at shearing-sheet boundaries
using this single-step unsplit integrator
Trailing wave never seeds
new leading waves
y
x
Evolution of specific vorticity
Initially W(k) ~ k-5/3
2D 5122 simulation
Box size: 4Hx4H
QuickTime™ and a
Evolved
to t=100
Video decompressor
are needed to see this picture.
Comparison of evolution in 2D and 3D
Wz initialized identically in 2D and 3D
Random Vz added in 3D
2D grid: 2562 (4Hx4H) 3D grid: 2562x64 (4Hx4Hx1H)
2D simulation
T = 0.5W t
T = 2W t
T = 5W t
T = 10W t
T = 20W t
In 2D, long-lived vortices emerge
In 3D, vorticity and turbulence decays much more rapidly
Evolution of KE: 2D versus 3D
NB: Initial evolution identical
Evolution of Stress: 2D versus 3D
In 2D, long-lived vortices emerge
In 3D, vorticity (and KE and stress) decays much more rapidly
Further 3D runs in both stratified and unstratified boxes warranted
E. MHD studies of the MRI
Start from a vertical field with zero net flux: Bz=B0sin(2px)
Sustained turbulence not possible in 2D – dissipation rate
after saturation is sensitive to numerical dissipation
2D MRI
Animation of angular velocity fluctuations: dVy=Vy+1.5W0x
shows saturation of MRI and decay in 2D
3rd order reconstruction,
2562 Grid
bmin=4000, orbits 2-10
QuickTime™ and a
GIF decompressor
are needed to see this picture.
Magnetic Energy Evolution
ZEUS vs. Athena
Numerical dissipation is ~1.5 times smaller with
CTU & 3rd order reconstruction than ZEUS.
3D MRI
Animation of angular velocity fluctuations: dVy=Vy+1.5W0x
Initial Field Geometry is Uniform By
CTU with 3rd order
reconstruction,
128 x 256 x 128 Grid
bmin= 100, orbits 4-20
QuickTime™ and a
Video decompressor
are needed to see this picture.
Goal: Since Athena is
strictly conservative, can
measure spectrum of T
fluctuations from
dissipation of turbulence
Requires adding opticallythin radiative cooling
Stress & Energy for <Bz>  0
No qualitative difference
with ZEUS results (Hawley,
Gammie, & Balbus 1995)
Energy Conservation
Saturation implies:
Previously observed for ReM = 1 (Sano et al. 2001).
Dependence of saturated state on cooling
Red line: no cooling; Green line: tcool = Q
Internal energy
Reynolds stress
Magnetic energy
Maxwell stress
Cooling has almost no effect except on internal energy
Density Fluctuations
• Compressible fluctuations lead to spiral waves & M ~ 1.5 shocks
• Temperature fluctuations are dominated by compressive waves.
• Will this be a generic result in global disks?
Going beyond MHD
If mean-free-path of particles is long compared to gyroradius,
anisotropic transport coefficients can be important (Braginskii 1965)
Generically, astrophysical plasmas are diffuse, and kinetic effects
may be important in some circumstances
These effects can produce qualitative changes to the dynamics:
e.g. with anisotropic heat conduction, the convective
instability criterion becomes dT/dz < 0 (Balbus 2000)
Evolution of field lines in atmosphere with dS/dz > 0 and dT/dz < 0
including anisotropic heat conduction,
Stable layer
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
Unstable
z
Stable layer
Details: 2D 1282 , initial b=100, Bx proportional to sin(z), fixed T at
vertical boundaries, isotropic conduction near boundaries
Summary
Developed a new Godunov scheme for MHD
Code now being used for various 3D applications
• Evolution of vortices in disks
• Energy dissipation in MRI turbulence
Further algorithm development is planned
• Curvilinear grids
• Nested and adaptive grids
Future applications
• 3D global models of thin disks
• Fragmentation and collapse in self-gravitating MHD turbulence
Beyond ideal MHD
• Kinetic effects in diffuse plasmas
• Particle acceleration (CRs are important!)
• Dust dynamics in protoplanetary disks
A. Astrophysical motivation
• Accretion disks
• Star formation
• Protoplanetary disks
• X-ray clusters
Numerical challenges
• 3D MHD
• Non-ideal MHD
• Radiation hydrodynamics (Prad >> Pgas)
• Nested and adaptive grids
• Kinetic MHD effects