Bug: SUPERBEE reverses decreasing slopes

[aliencaocao]

In

MFC/src/simulation/m_muscl.fpp

Lines 201 to 212 in 511cda5

	else if (muscl_lim == muscl_lim_superbee) then ! SUPERBEE
	if (slopeL*slopeR > muscl_eps) then
	slope = -1._wp*min(-min(2._wp*abs(slopeL), abs(slopeR)), -min(abs(slopeL), &
	& 2._wp*abs(slopeR)))
	end if
	end if
	! reconstruct from left side
	vL_rs_vf_x(j, k, l, i) = v_rs_ws_muscl(${SF('')}$, i) - (5.e-1_wp*slope)
	! reconstruct from the right side
	vR_rs_vf_x(j, k, l, i) = v_rs_ws_muscl(${SF('')}$, i) + (5.e-1_wp*slope)

, SUPERBEE computes the slope magnitude without restoring its sign.

Cell values (3, 2, 1) give slope +1 instead of -1, reversing the reconstructed interface states. This affects second-order SUPERBEE wherever the limited slope should be negative.

Minimal repro on CPU with FP64, GCC 13.3:

Based on suite test 3EE4302C but replace the case.json with:

{
  "run_time_info": "T",
  "m": 63,
  "n": 0,
  "p": 0,
  "dt": 0.001,
  "t_step_start": 0,
  "t_step_stop": 100,
  "t_step_save": 100,
  "num_patches": 3,
  "model_eqns": 2,
  "alt_soundspeed": "F",
  "num_fluids": 1,
  "mpp_lim": "F",
  "mixture_err": "F",
  "time_stepper": 3,
  "recon_type": 2,
  "weno_order": 0,
  "mapped_weno": "F",
  "null_weights": "F",
  "mp_weno": "F",
  "riemann_solver": 2,
  "wave_speeds": 1,
  "avg_state": 2,
  "format": 1,
  "precision": 2,
  "patch_icpp(1)%pres": 1.0,
  "patch_icpp(1)%alpha_rho(1)": "1.0+0.2*sin(2.0*pi*x)",
  "patch_icpp(1)%alpha(1)": 1.0,
  "patch_icpp(2)%pres": 1.0,
  "patch_icpp(2)%alpha_rho(1)": "1.0+0.2*sin(2.0*pi*x)",
  "patch_icpp(2)%alpha(1)": 1.0,
  "patch_icpp(3)%pres": 1.0,
  "patch_icpp(3)%alpha_rho(1)": "1.0+0.2*sin(2.0*pi*x)",
  "patch_icpp(3)%alpha(1)": 1.0,
  "fluid_pp(1)%gamma": 2.5000000000000004,
  "fluid_pp(1)%eos": 2,
  "fluid_pp(1)%cv": 0.0,
  "fluid_pp(1)%qv": 0.0,
  "fluid_pp(1)%qvp": 0.0,
  "bubbles_euler": "F",
  "bubble_model": 3,
  "polytropic": "T",
  "polydisperse": "F",
  "thermal": 3,
  "patch_icpp(1)%r0": 1,
  "patch_icpp(1)%v0": 0,
  "patch_icpp(2)%r0": 1,
  "patch_icpp(2)%v0": 0,
  "patch_icpp(3)%r0": 1,
  "patch_icpp(3)%v0": 0,
  "qbmm": "F",
  "dist_type": 2,
  "poly_sigma": 0.3,
  "sigR": 0.1,
  "sigV": 0.1,
  "rhoRV": 0.0,
  "acoustic_source": "F",
  "num_source": 1,
  "acoustic(1)%loc(1)": 0.5,
  "acoustic(1)%mag": 0.2,
  "acoustic(1)%length": 0.25,
  "acoustic(1)%dir": 1.0,
  "acoustic(1)%npulse": 1,
  "acoustic(1)%pulse": 1,
  "rdma_mpi": "F",
  "bubbles_lagrange": "F",
  "lag_params%nBubs_glb": 1,
  "lag_params%solver_approach": 0,
  "lag_params%cluster_type": 2,
  "lag_params%pressure_corrector": "F",
  "lag_params%smooth_type": 1,
  "lag_params%epsilonb": 1.0,
  "lag_params%heatTransfer_model": "F",
  "lag_params%massTransfer_model": "F",
  "lag_params%valmaxvoid": 0.9,
  "x_domain%beg": 0.0,
  "x_domain%end": 1.0,
  "bc_x%beg": -1,
  "bc_x%end": -1,
  "patch_icpp(1)%geometry": 1,
  "patch_icpp(1)%x_centroid": 0.05,
  "patch_icpp(1)%length_x": 0.1,
  "patch_icpp(2)%x_centroid": 0.45,
  "patch_icpp(2)%length_x": 0.7,
  "patch_icpp(3)%x_centroid": 0.9,
  "patch_icpp(3)%length_x": 0.2,
  "patch_icpp(1)%vel(1)": 1.0,
  "patch_icpp(2)%geometry": 1,
  "patch_icpp(2)%vel(1)": 1.0,
  "patch_icpp(3)%geometry": 1,
  "patch_icpp(3)%vel(1)": 1.0,
  "muscl_order": 2,
  "muscl_lim": 5,
  "muscl_eps": 0,
  "parallel_io": "F",
  "prim_vars_wrt": "T",
  "cons_vars_wrt": "T"
}

./mfc.sh build -t pre_process simulation --no-mpi --no-gpu -j "$(nproc)"
./mfc.sh run superbee/case.json -t pre_process simulation --no-mpi --no-gpu -n 1 -j "$(nproc)"

See the density value in cons.1.00.000100.dat. At t=0.1, it should follow 1+0.2*sin(2*pi*(x-0.1)). The file labels cell-centered values with the right cell boundary, so subtract 1/128 from its coordinate. Right now, the density RMS error is 0.00568472.
The existing golden test result is wrong for this.

Fix is to restore the sign with if (slopeL < 0._wp) slope = -slope, as minmod and MC already do.