Bug: pTg phase depletion violates volume and energy closure

[aliencaocao]

After pTg projects a reacting phase mass to zero,

MFC/src/common/m_phase_change.fpp

Lines 124 to 170 in 511cda5

	! Check if pTg-equilibrium needed; only partial densities require updating
	if ((relax_model == 6) .and. ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
	& l) > mixM*rM) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, &
	& l) > mixM*rM)) .and. (pS < pCr) .and. (TS < TCr)) then
	! Solve pTg-equilibrium directly on the actual reacting masses. The Newton solver projects
	! the liquid mass onto [0, mT], so it recovers the single-phase limits itself (ml -> 0 for
	! all-vapor, ml -> mT for all-liquid). The former overheated-vapor / subcooled-liquid pT
	! shortcuts were removed: their pT states differ O(1) from the pTg equilibrium, so the
	! sub-ULP shortcut/pTg branch decision flipped across backends (CPU vs GPU) near a phase
	! boundary and destroyed cross-backend reproducibility.
	q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m1
	q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m2
	call s_infinite_ptg_relaxation_k(j, k, l, pS, rhoe, q_cons_vf, TS)
	end if
	! Calculations AFTER equilibrium
	$:GPU_LOOP(parallelism='[seq]')
	do i = 1, num_fluids
	! entropy
	sk(i) = cvs(i)*log((TS**isentrope_n(i))/((pS + isentrope_B(i))**(isentrope_n(i) - 1.0_wp))) + qvps(i)
	! enthalpy
	hk(i) = isentrope_n(i)*cvs(i)*TS + qvs(i)
	! Gibbs-free energy
	gk(i) = hk(i) - TS*sk(i)
	! densities
	rhok(i) = f_sg_thermal(pS, TS, isentrope_n(i), isentrope_B(i), cvs(i))
	! internal energy
	ek(i) = (pS + isentrope_n(i)*isentrope_B(i))/(pS + isentrope_B(i))*cvs(i)*TS + qvs(i)
	end do
	! calculating volume fractions, internal energies, and total entropy
	rhos = 0.0_wp
	$:GPU_LOOP(parallelism='[seq]')
	do i = 1, num_fluids
	! volume fractions
	q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l) = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)/rhok(i)
	! alpha*rho*e
	if (model_eqns == model_eqns_6eq) then
	q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
	& l)*ek(i)

reconstructs fractions and phasic energies without restoring pT equilibrium for the final masses. Gibbs equality no longer constrains an absent phase. This affects relax_model=6 depletion states, producing invalid fractions and inconsistent six-equation phasic energies.

Minimal repro:
Based on D21F4F38 but replace case.json:

{
  "run_time_info": "T",
  "m": 299,
  "n": 0,
  "p": 0,
  "dt": 1e-07,
  "t_step_start": 0,
  "t_step_stop": 50,
  "t_step_save": 50,
  "num_patches": 3,
  "model_eqns": 3,
  "alt_soundspeed": "F",
  "num_fluids": 2,
  "mpp_lim": "F",
  "mixture_err": "F",
  "time_stepper": 3,
  "recon_type": 1,
  "weno_order": 5,
  "weno_eps": 1e-16,
  "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": 437550.0,
  "patch_icpp(1)%alpha_rho(1)": 0.008406131093,
  "patch_icpp(1)%alpha(1)": 8.7149e-06,
  "patch_icpp(2)%pres": 96602.0,
  "patch_icpp(2)%alpha_rho(1)": 0.040265879299999995,
  "patch_icpp(2)%alpha(1)": 3.6749e-05,
  "patch_icpp(3)%pres": 96602.0,
  "patch_icpp(3)%alpha_rho(1)": 0.040265879299999995,
  "patch_icpp(3)%alpha(1)": 3.6749e-05,
  "fluid_pp(1)%gamma": 0.7409,
  "fluid_pp(1)%eos": 1,
  "fluid_pp(1)%cv": 1816,
  "fluid_pp(1)%qv": -1167000,
  "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)": 606.15,
  "patch_icpp(2)%geometry": 1,
  "patch_icpp(2)%vel(1)": 10.0,
  "patch_icpp(3)%geometry": 1,
  "patch_icpp(3)%vel(1)": 10.0,
  "relax": "T",
  "relax_model": 6,
  "palpha_eps": 0.01,
  "ptgalpha_eps": 0.01,
  "fluid_pp(1)%pi_inf": 1740900000.0,
  "fluid_pp(2)%gamma": 2.3266,
  "fluid_pp(2)%eos": 2,
  "fluid_pp(2)%cv": 1040,
  "fluid_pp(2)%qv": 2030000,
  "fluid_pp(2)%qvp": -23400,
  "patch_icpp(1)%alpha(2)": 0.9999912851,
  "patch_icpp(1)%alpha_rho(2)": 2.31317984069332,
  "patch_icpp(2)%alpha(2)": 0.999963251,
  "patch_icpp(2)%alpha_rho(2)": 0.5802786745553,
  "patch_icpp(3)%alpha(2)": 0.999963251,
  "patch_icpp(3)%alpha_rho(2)": 0.5802786745553,
  "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 phase/case.json -t pre_process simulation --no-mpi --no-gpu -n 1 -j "$(nproc)"

See the last columns of cons.5.00.000050.dat and cons.6.00.000050.dat. Currently it produces alpha_liquid=0, alpha_vapor=1.051528, their sum must be one. The existin golden file contain wrong data.

liquid mass = 0
liquid volume fraction = 0
vapor volume fraction = 1.05152801179084

The vapor volume fraction should not be more than 1.

Fix: if either reacting mass is exactly zero after pTg relaxation, call the existing pT solver before recovering fractions and phasic energies.