Bug: LF exports volume fraction instead of interface velocity

[aliencaocao]

In

MFC/src/simulation/m_riemann_solver_lf.fpp

Lines 282 to 289 in 511cda5

	! Advection flux and source: interface velocity for volume fraction transport
	$:GPU_LOOP(parallelism='[seq]')
	do i = eqn_idx%adv%beg, eqn_idx%adv%end
	flux_rsx_vf(${SF('')}$, i) = (qL_prim_rsx_vf(${SF('')}$, i) - qR_prim_rsx_vf(${SF(' + 1')}$, &
	& i))*s_M*s_P/(s_M - s_P)
	flux_src_rsx_vf(${SF('')}$, i) = (s_M*qR_prim_rsx_vf(${SF(' + 1')}$, &
	& i) - s_P*qL_prim_rsx_vf(${SF('')}$, i))/(s_M - s_P)
	end do

, LF exports an averaged fraction where the RHS expects face-normal velocity, and omits the advective fraction flux. Identical states alpha=0.25, u=2 return (flux, velocity)=(0, 0.25) instead of (0.5, 2). This corrupts multifluid fraction transport and its source/boundary terms.

Minimal repro on GCC 13.3 on CPU with FP64:

Based on 96AF00B3 but replace case.json with:

{
  "run_time_info": "T",
  "m": 299,
  "n": 0,
  "p": 0,
  "dt": 0.0005,
  "t_step_start": 0,
  "t_step_stop": 50,
  "t_step_save": 50,
  "num_patches": 3,
  "model_eqns": 2,
  "alt_soundspeed": "F",
  "num_fluids": 2,
  "mpp_lim": "F",
  "mixture_err": "T",
  "time_stepper": 3,
  "recon_type": 1,
  "weno_order": 5,
  "weno_eps": 1e-16,
  "mapped_weno": "F",
  "null_weights": "F",
  "mp_weno": "F",
  "riemann_solver": 5,
  "wave_speeds": 1,
  "avg_state": 2,
  "format": 1,
  "precision": 2,
  "patch_icpp(1)%pres": 1.0,
  "patch_icpp(1)%alpha_rho(1)": 0.81,
  "patch_icpp(1)%alpha(1)": 0.9,
  "patch_icpp(2)%pres": 0.5,
  "patch_icpp(2)%alpha_rho(1)": 0.25,
  "patch_icpp(2)%alpha(1)": 0.5,
  "patch_icpp(3)%pres": 0.1,
  "patch_icpp(3)%alpha_rho(1)": 0.08,
  "patch_icpp(3)%alpha(1)": 0.2,
  "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": -3,
  "bc_x%end": -3,
  "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)": 0.0,
  "patch_icpp(2)%geometry": 1,
  "patch_icpp(2)%vel(1)": 0.0,
  "patch_icpp(3)%geometry": 1,
  "patch_icpp(3)%vel(1)": 0.0,
  "fluid_pp(2)%gamma": 2.5,
  "fluid_pp(2)%eos": 2,
  "patch_icpp(1)%alpha_rho(2)": 0.19,
  "patch_icpp(1)%alpha(2)": 0.1,
  "patch_icpp(2)%alpha_rho(2)": 0.25,
  "patch_icpp(2)%alpha(2)": 0.5,
  "patch_icpp(3)%alpha_rho(2)": 0.0225,
  "patch_icpp(3)%alpha(2)": 0.8,
  "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 lf/case.json -t pre_process simulation --no-mpi --no-gpu -n 1 -j "$(nproc)"

Check the last columns of lf/D/cons.5.00.000050.dat and cons.6.00.000050.dat (the two volume fractions). Their sum must be one (as per https://comp-physics.group/papers/wilfong-mfc-26.pdf page 3 section 2.1.1 eqn 4). The existing golden test passes despite this physical error because the reference golden.txt is wrong. It has

alpha1 = 0.52599988041605
alpha2 = 0.06994055693542
sum    = 0.59594043735147 != 1

The fix is to include the advective LF fraction flux, export the shared face-normal velocity, and zero the radial fraction geometric source. I have one ready, waiting for confirmation on the issue before I PR.