MHD constrained transport calculations create asymmetric floating-point errors

[felker]

Bug report

Summary of issue:

This has been a known issue since April 2018, but I am filing a bug report now in order to document it and help track down the source of the problem soon.

The constrained transport scheme introduces small asymmetries in problems that should be perfectly symmetric: first in the B field, then in the Hydro variables. The high-order solver quickly amplifies these differences into visible asymmetries. I have confirmed that these differences occur even with the second-order scheme with time/xorder=2 PLM and time/integrator=vl2; solving the Orszag-Tang problem at high-resolution will produce visible asymmetries even at second-order accuracy.

Partial fixes have already been applied to various MHD algorithmic components:

• HLLD: If you consider the Orszag-Tang problem, but set B=0 to formulate a purely hydrodynamics problem, it still should maintain 90-degree rotational symmetry. Originally HLLE was the only MHD Riemann solver which produced solutions that satisfied this. PR Improve symmetry-preservation of HLLD floating point operations #120 fixed a minor bias in HLLD on 2018-04-28.
• Averaging field components on longitudinal faces: PR Use simple average when computing cell-centered B  #122 was merged on 2018-05-03 which, for example, fixed the initial condition in the MHD Rotor problem.

Nevertheless, symmetry is still imperfect in the Orszag-Tang solutions, and others have reported asymmetries with a Magnetized Noh implosion problem variant solved using the second order scheme with PPM reconstruction:

The issue is almost certainly due to the left-to-right associativity of floating-point arithmetic in the calculations of the electric field gradients and the upwinded average at the corner:

athena/src/field/calculate_corner_e.cpp

Lines 92 to 112 in b60e890

	// integrate E3 to corner using SG07
	for (int k=ks; k<=ke; ++k) {
	for (int j=js; j<=je+1; ++j) {
	#pragma omp simd
	for (int i=is; i<=ie+1; ++i) {
	Real de3_l2 = (1.0-w_x1f(k,j-1,i))*(e3_x2f(k,j,i ) - cc_e_(k,j-1,i )) +
	( w_x1f(k,j-1,i))*(e3_x2f(k,j,i-1) - cc_e_(k,j-1,i-1));
	Real de3_r2 = (1.0-w_x1f(k,j ,i))*(e3_x2f(k,j,i ) - cc_e_(k,j ,i )) +
	( w_x1f(k,j ,i))*(e3_x2f(k,j,i-1) - cc_e_(k,j ,i-1));
	Real de3_l1 = (1.0-w_x2f(k,j,i-1))*(e3_x1f(k,j ,i) - cc_e_(k,j ,i-1)) +
	( w_x2f(k,j,i-1))*(e3_x1f(k,j-1,i) - cc_e_(k,j-1,i-1));
	Real de3_r1 = (1.0-w_x2f(k,j,i ))*(e3_x1f(k,j ,i) - cc_e_(k,j ,i )) +
	( w_x2f(k,j,i ))*(e3_x1f(k,j-1,i) - cc_e_(k,j-1,i ));
	e3(k,j,i) = 0.25*(de3_l1 + de3_r1 + de3_l2 + de3_r2 + e3_x2f(k,j,i-1) +
	e3_x2f(k,j,i) + e3_x1f(k,j-1,i) + e3_x1f(k,j,i));
	}
	}}

However, the CT operations contain many more terms than the analogous sources of asymmetry fixed in #144 (viscosity), #98 (PPM), #112 (x1f calculations), #114 (AMR/SMR prolongation and restriction operators) since those are all 1D stencils. So, rearranging the terms to guarantee that there is no directional bias to the floating-point round-off error will be more difficult.

Proposed steps:

• Start debugging with the simpler mirror symmetry of the MHD R-T problem before attempting to fix the rotational symmetry of Orszag-Tang