OpenMHD: Godunov-type code formagnetic reconnection and MHD problems

Seiji ZENITANI
Kobe University[1] Zenitani, Phys. Plasmas 22, 032114 (2015)[2] Zenitani, Astrophysics Source Code Library, ascl:1604.001 (2016)

MHD shocks in reconnection
vx
5

764
128z

3

2016

19C

A

914

15

E

13

9

1

10
B

Petschekslow shocks

87

18

114aD

104b

175 6123
16x

2
1. Petschek slow shock (Petschek 1964)
2. outer shell = slow shock (Ugai 1995)
3. intermediate shock (Abe & Hoshino 2001) or slow shock (Saito et al. 1995)
4a fast shock (Forbes & Priest 1983)
4b oblique shock & Mach disk (Takasao et al. 2015)
5. looptop front (Ugai 1987)
6. tangential discontinuity
7. post-plasmoid vertical slow shock (Zenitani et al. 2010)
8. outer vertical slow shock (Zenitani & Miyoshi 2011)
9. fast-mode wave front (Saito et al. 1995)
10.overexpanded shock-diamonds (Zenitani et al. 2010)
A. reconnection inflow
11.underexpanded shock-diamonds (Zenitani 2015)
B. outflow jet
12.slow expansion wavefront (Zenitani 2015)C. post-plasmoid reverse flow
13.contact discontinuity (Zenitani & Miyoshi 2011, 2015)
D. internal flow
14.contact discontinuity (Zenitani 2015)
E. flapping jet (KH instability)
15.vortex-driven shocklets (Miura 1992, 1995)
16.contact discontinuity (Zenitani 2015)17.rotational discontinuity [in guide-field reconnection] (Petschek & Thorne 1967)
18.conduction front [with heat conduction] (Yokoyama & Shibata 1997)
19.forward head shock [in asymmetric reconnection] (Nitta et al. 2016)

Zenitani & Miyoshi 2011
Zenitani 2015

OpenMHD - public code
• MHD code by Zenitani &Miyoshi (2011, 2015) withrecent improvements• Fortran 90, 2D (Internally 3D),
HLLD, MUSCL, MPI+OpenMP,MPI-IO, GPU (CUDA)• Visualization code in IDL andPython3 (matplotlib)

https://sci.nao.ac.jp/MEMBER/zenitani/openmhd-e.html

MHD equations• Resistive MHD equations• Ohm’s law• Hyperbolic divergence cleaning (Dedner+ 2002)• Second-order TVD RK + Strang splitting

Riemann solver
UUUx

HLLD solver• Proposed by Miyoshi & Kusano 2005
• An advanced variant of HLL solver

UUA problem in the HLLD scheme• In the parallel direction, CA = Cfast for β<1.2.
• RD sometimes outruns the Riemann fan.
• We offer the problem parameters.
β=0.3

β=2.0

B
BCfastCA**CfastCA! left
vx0 = -3.550886791794632E-002
vy0 = 0.202198834859736
vz0 = 0.d0
pr0 = 0.573857561003025
ro0 = 1.46274860814133
bx0 = 1.d0
by0 = 0.535958263440366E-004
bz0 = 0.d0! rightvx0 = -3.563873070012032E-002
vy0 = 0.201451798427819
vz0 = 0.d0
pr0 = 0.575034711174729
ro0 = 2.46444099174925
bx0 = 1.d0
by0 = -4.588121811862009E-004
bz0 = 0.d0

Solutions• 1. HLLC=HLLD hybrid solver•

Positivity•

If clauses - not ideal for optimization• 2. Faster Cfast (Miyoshi 2015 private comm.)• 3. Adjust wavespeeds: SL*=max(SL,SL*), SR*=min(SR*,SR)•Good for optimization•

Theoreticians are unhappy

Cfast

CfastCACA*

δv/CA~10-3-4*

OpenMHD/2019 (HLLD) - Profiler output
Orszag-Tang vortex

Magnetic reconnection

1%1%
2%
3%
4%

1%
2%
2%
2%
3%
3%

19%

14%

70%

ﬂux̲solver
limiterset̲dtu2v
tvd̲rk22
tvd̲rk21
ﬂux̲glmothergfortran 4.8

71%

ﬂux̲solver
limiterset̲dttvd̲rk22
ﬂux̲resistive
tvd̲rk21
u2v
ﬂux̲glmotherPerformance tuning (1/2)• Loop indexes• Specifying intent attributes in subroutines•

In particular, intent(in) and intent(out) work great.
subroutine bc(U,ix,jx)
implicit none
include 'param.h'
integer :: ix, jx

subroutine bc(U,ix,jx)
implicit noneinclude 'param.h'
integer, intent(in) :: ix, jx• Reducing if-statements•

This does not improve the performance on Intel computers,
but it works well on FX computers•

minmod limiter:if( gA*gB .le. 0 ) thengrad = 0.d0
else
if( gA .gt. 0 ) thengrad = 0.5d0 * min(gA,gB)
else
grad = 0.5d0 * max(gA,gB)
endif
endif•

grad = (sign(0.25d0,gA)+sign(0.25d0,gB))*min(abs(gA),abs(gB))

MC limiter:
grad = (sign(0.5d0,gA)+sign(0.5d0,gB))*min(abs(gA),abs(gB),0.25d0*abs(gA+gB))

Performance tuning (2/2)• Reducing square roots and divisions•[before] Fastest fast-mode speed
vfL = sqrt( ( (f1+B2) + sqrt(max( (f1+B2)**2-f2, 0.d0 ))) / ( 2*VL(i,j,ro) ))vfR = sqrt( ( (f1+B2) + sqrt(max( (f1+B2)**2-f2, 0.d0 ))) / ( 2*VR(i,j,ro) ))f1 = max( vfL, vfR )•[after]aL = ( (f1+B2) + sqrt(max( (f1+B2)**2-f2, 0.d0 ))) / ( 2*VL(i,j,ro) )aR = ( (f1+B2) + sqrt(max( (f1+B2)**2-f2, 0.d0 ))) / ( 2*VR(i,j,ro) )f1 = sqrt( max( aL, aR ) )• Using register variables• Using collapse option in OpenMP directives!$omp parallel do private(i,j,k) collapse(2)do k=1,var1
do j=2,jx-1
do i=2,ix-1• Using good compilers•

Intel Fortran compiler is usually better than gfortran•

Gfortran 5.4 or 6.x is substantially better than gfortran 4.8

Improving performance - slope limiters
minmodMCfast
500000(a+b) - min(0.d0,a,b) - max(0.d0,a,b)

400000(sign(0.5d0,a)+sign(0.5d0,b))* min( abs(a),abs(b) )300000sign(1.d0,a) * max( 0.d0,min(abs(a),sign(1.d0,a)*b) )200000100000If statements

0.95 0.96 0.97 0.98 0.99
Intel (gfortran -O2)

fast

Intel (ifort)

1

1.01 1.02 1.03 1.04

SPARC (Fujitsu fortran)

0FX

If statements
minmod(2a,minmod(2b,c))(sign(0.5d0,a)+sign(0.5d0,b)) * min( abs(2a),abs(2b), abs(c) )Improving performance - slope limiters
If statements

min, sign etc.

300000300000270000270000240000240000210000210000180000180000150000150000FX
minmod

MC

fast

FX
Koren

MPI nonblocking communication

1time• Other operation while
sending/receiving data(communication hiding)2(nonblocking)(blocking)call mpi_sendrecv(.. )call mpi_sendrecv(.. )12call mpi_irecv(.. )call mpi_isend(.. )call mpi_irecv(.. )call mpi_isend(.. )

time
1
1
2
2

call mpi_waitall(.. )1call mpi_waitall(.. )2

Parallel I/O (MPI-IO)• Multiple MPI processes concurrently access data in a single file.
• The way to access the shared file is defined by MPI_File_set_view.

MPI
process

fileMPI
process

MPI
process

MPI
process

Parallel efficiency / strong scaling• Very nice scaling up to 5000 cores; moderate scaling at 10000 cores• 94% efficiency at 3000 cores with flat MPI• 98% at 5000 physical cores with hyperthreading (@Xeon Phi KNL)
timeefficiency

1000

Ideal
2D: N x 1

750

2D: N x 10
2D: N x 50

500

1.2

2omp+2HT

10.8
0.6
0.4

250

0.200300

3000

phys. cores3003000phys. cores

OpenMHD/2019

Performance: years of improvement
KDK-A (Xeon Phi KNL)NAOJ XC30 (Xeon)

JAXA FX100 (SPARC)

50343037.5

25.5

22.525171512.5

8.5

fast 7.502018 2016 2015 2014 2011

HLL/HLLD

02018 2016 2015 2014 2011

HLL/HLLD
Optimized for various CPU architectures

02018 2016 2015 2014 2011

HLL/HLLD

GPU version• OpenMHD also runs on NVIDIA GPUs
• Partially written by CUDA FortranO-T vortex
800 x 800

1000x 40750
500
2500CPUCPU x4 RTX 2080 GV10030MagneticReconnection
2000 x 500

231580

CPU x4 RTX 2080 GP100

GV100fastWeb site• Search “OpenMHD” with Google.
• Online documents [in Japanese]. English translation in progress.
• Gateway to a mailing list (41 users).

https://sci.nao.ac.jp/MEMBER/zenitani/openmhd-e.html

Built-in problems• Riemann problem• Solar-wind problem• KH & secondary instabilities• Orszag-Tang vortex• Magnetic reconnection

Source code management• Our repository is hosted by GitHub.
• It is open to public -- pull requests are always welcome!

Research papers based on OpenMHD• Asymmetric reconnection• Plasmoid-dominated reconnection

Shimizu et al. 2017 PoPHosseinpour et al. 2018 PoP
Nitta et al. 2016, 2019 ApJ•

Used in college & graduate school•

Kobe University and Ehime UniversityReferences• [1] S. Zenitani & T. Miyoshi, Phys. Plasmas 18, 022105 (2011)• [2] S. Zenitani, Phys. Plasmas 22, 032114 (2015)• [3] A. Dedner et al., J. Comput. Phys. 175, 645 (2002)• [4] T. Miyoshi & K. Kusano, J. Comput. Phys. 208, 315 (2005)• [5] 銭谷誠司, 生存圏研究 13, 27 (2017)

URLs• https://ascl.net/1604.001• https://github.com/zenitani/OpenMHD• https://sci.nao.ac.jp/MEMBER/zenitani/openmhd-e.html