MHD reproducibility

[tomidakn]

This commit improves reproducibility of parallel MHD simulations with a new EMF correction method, taking a maximum value among coinciding edges between MeshBlocks. See issue #645.

Prerequisite checklist

• My code follows the Athena++ Style Guide
• My change requires a change to the documentation.
• I have updated the documentation in the Wiki accordingly.
• I have added tests to cover my changes.
• All new and existing tests passed.

Description

Previously we simply used average, but such operations are not reproducible because they depend on the ordering of arrival of MPI communications. This implementation is simpler than the previous one, as it does not need to count the number of coinciding edges.

Note 1: this improves the reproducibility between runs, but it does NOT reproduce calculations with different MeshBlock configuration (e.g. 1 MeshBlock vs multiple MeshBlocks).

Note 2: the operation can be minimum instead of maximum, or max/min in absolute values. I even tried switching max/min at every cycle, but I could not find any meaningful difference between them. All we need is to use exactly same EMF for coinciding edges.

Note 3: I now feel this implementation is better than the previous one, but I'm not 100% sure yet.

This also contains minor leftover fixes.

[tomidakn]

I will add a note on Wiki as the code becomes different from what is written on the paper once this is merged.

[felker]

  Your Athena++ distribution has now been configured with the following options: 
  Problem generator:            jgg
  Coordinate system:            cartesian
  Equation of state:            isothermal
  Riemann solver:               hlld
  Magnetic fields:              ON
  Number of scalars:            0
  Number of chemical species:   0
  Special relativity:           OFF
  General relativity:           OFF
  Radiative Transfer:           OFF
  Implicit Radiation:           OFF
  Cosmic Ray Transport:         OFF
  Cosmic Ray Diffusion:         OFF
  Frame transformations:        OFF
  Self-Gravity:                 OFF
  Super-Time-Stepping:          OFF
  Chemistry:                    OFF
  KIDA rates:                   OFF
  ChemRadiation:                OFF
  chem_ode_solver:              OFF
  Debug flags:                  OFF
  Code coverage flags:          OFF
  Linker flags:                   -lfftw3
  Floating-point precision:     double
  Number of ghost cells:        4
  MPI parallelism:              OFF
  OpenMP parallelism:           OFF
  FFT:                          ON
  HDF5 output:                  OFF
  Compiler:                     g++
  Compilation command:          g++  -O3 -std=c++11

Setup complete, entering main loop...


Terminating on time limit
time=3 cycle=653
tlim=3 nlim=1000

zone-cycles = 5349376
cpu time used  = 11.7239
zone-cycles/cpu_second = 456281

Setup complete, entering main loop...


Terminating on time limit
time=3 cycle=646
tlim=3 nlim=1000

zone-cycles = 5292032
cpu time used  = 13.8666
zone-cycles/cpu_second = 381638
[MHD_SHWAVE]: L1-Norm Errors
[MHD_SHWAVE]: w/o Orbital Advection
[MHD_SHWAVE]: epsilon_c = 0.9165251299739616, xi_dBx = 1.7225109800993983
[MHD_SHWAVE]: dByc Error is more than 1%.
[MHD_SHWAVE]: xi_dBx Error is more than 1%.
[MHD_SHWAVE]: w/  Orbital Advection
[MHD_SHWAVE]: epsilon_c = 0.010621588626771223, xi_dBx = 0.02034164709843645
[MHD_SHWAVE]: dByc Error is more than 1%.
[MHD_SHWAVE]: xi_dBx Error is more than 1%.
shearingbox.mhd_shwave test: prepare(), run(), analyze() finished
  Your Athena++ distribution has now been configured with the following options: 
  Problem generator:            hb3
  Coordinate system:            cartesian
  Equation of state:            isothermal
  Riemann solver:               hlld
  Magnetic fields:              ON
  Number of scalars:            0
  Number of chemical species:   0
  Special relativity:           OFF
  General relativity:           OFF
  Radiative Transfer:           OFF
  Implicit Radiation:           OFF
  Cosmic Ray Transport:         OFF
  Cosmic Ray Diffusion:         OFF
  Frame transformations:        OFF
  Self-Gravity:                 OFF
  Super-Time-Stepping:          OFF
  Chemistry:                    OFF
  KIDA rates:                   OFF
  ChemRadiation:                OFF
  chem_ode_solver:              OFF
  Debug flags:                  OFF
  Code coverage flags:          OFF
  Linker flags:                  
  Floating-point precision:     double
  Number of ghost cells:        2
  MPI parallelism:              OFF
  OpenMP parallelism:           OFF
  FFT:                          OFF
  HDF5 output:                  OFF
  Compiler:                     g++
  Compilation command:          g++  -O3 -std=c++11
iso_cs = 0.00408
d0     = 1
p0     = 1.66464e-05
B0     = 9.12316e-05
ipert  = 1
ifield = 1
beta   = 4000
[hb3.cpp]: [Lx,Lz,Ly] = [1,1,1]

Setup complete, entering main loop...


Terminating on time limit
time=50265.5 cycle=35171
tlim=50265.5 nlim=500000

zone-cycles = 144060416
cpu time used  = 73.6399
zone-cycles/cpu_second = 1.95628e+06
[MRI-2D]: Ref(stress,ME,ratio) = 4.308659080942425e-05 0.00015249213554582367 3.539201702457957
[MRI-2D]: New(stress,ME,ratio) = 4.398143733780276e-05 0.00015356316377865482 3.4915449124411584
shearingbox.mri2d test: prepare(), run(), analyze() finished
  Your Athena++ distribution has now been configured with the following options: 
  Problem generator:            ssheet
  Coordinate system:            cartesian
  Equation of state:            isothermal
  Riemann solver:               hlle
  Magnetic fields:              OFF
  Number of scalars:            0
  Number of chemical species:   0
  Special relativity:           OFF
  General relativity:           OFF
  Radiative Transfer:           OFF
  Implicit Radiation:           OFF
  Cosmic Ray Transport:         OFF
  Cosmic Ray Diffusion:         OFF
  Frame transformations:        OFF
  Self-Gravity:                 OFF
  Super-Time-Stepping:          OFF
  Chemistry:                    OFF
  KIDA rates:                   OFF
  ChemRadiation:                OFF
  chem_ode_solver:              OFF
  Debug flags:                  OFF
  Code coverage flags:          OFF
  Linker flags:                  
  Floating-point precision:     double
  Number of ghost cells:        2
  MPI parallelism:              OFF
  OpenMP parallelism:           OFF
  FFT:                          OFF
  HDF5 output:                  OFF
  Compiler:                     g++
  Compilation command:          g++  -O3 -std=c++11
iso_cs = 0.001
d0 = 1
p0 = 1e-06
ipert  = 3
[ssheet.cpp]: [Lx,Ly,Lz] = [4,4,1]

Setup complete, entering main loop...


Terminating on time limit
time=8000 cycle=1266
tlim=8000 nlim=2000

zone-cycles = 5185536
cpu time used  = 1.12162
zone-cycles/cpu_second = 4.62327e+06
iso_cs = 0.001
d0 = 1
p0 = 1e-06
ipert  = 3
[ssheet.cpp]: [Lx,Ly,Lz] = [4,4,1]

Setup complete, entering main loop...


Terminating on time limit
time=8000 cycle=321
tlim=8000 nlim=2000

zone-cycles = 1314816
cpu time used  = 0.357559
zone-cycles/cpu_second = 3.6772e+06
[SSHEET_SHWAVE]: L1-Norm Errors
[SSHEET_SHWAVE]: w/o Orbital Advection
[SSHEET_SHWAVE]: epsilon_c = 0.107596209104392, epsilon_s = 0.002195062813586482
[SSHEET_SHWAVE]: w/  Orbital Advection
[SSHEET_SHWAVE]: epsilon_c = 0.091714869167537, epsilon_s = 0.0002732013850596203
shearingbox.ssheet test: prepare(), run(), analyze() finished

Results:
    shearingbox.mhd_shwave: failed; time elapsed: 39.1 s
    shearingbox.mri2d: passed; time elapsed: 87.2 s
    shearingbox.ssheet: passed; time elapsed: 15.1 s

Summary: 2 out of 3 tests passed

Traceback (most recent call last):
  File "./run_tests.py", line 331, in <module>
    main(**vars(args))
  File "./run_tests.py", line 198, in main
    raise TestError()
__main__.TestError

real	2m22.003s
user	6m20.525s
sys	0m28.447s

[tomidakn]

I'll look into it.