MFC does not work when multi-fluid and chemistry are set at the same time

[ChenHuangwei]

Hello！I am simulating the interaction between reactive shock waves and droplets. When I simulate the single fluid and turn on chemistry, it works well. However, when I set the droplet, it was wrong. The error reporting information and case file are given below.

import json, argparse
import cantera as ct

parser = argparse.ArgumentParser(prog="oran_detonation_6670", formatter_class=argparse.ArgumentDefaultsHelpFormatter)

parser.add_argument(
    "--mfc",
    type=json.loads,
    default="{}",
    metavar="DICT",
    help="MFC's toolchain's internal state.",
)
parser.add_argument(
    "--no-chem",
    dest="chemistry",
    default=True,
    action="store_false",
    help="Disable chemistry.",
)
parser.add_argument("--scale", type=float, default=1, help="Scale.")

args = parser.parse_args()

ctfile = "oran.yaml"
sol_ig = ct.Solution(ctfile)
sol_ig.TPX = 2000, 5000000, "H2:2,O2:1,AR:7"

sol_amb = ct.Solution(ctfile)
sol_amb.TPX = 300, 101325, "H2:2,O2:1,AR:7"

gamma_amb = sol_amb.cp / sol_amb.cv

# water droplet
rho_w = 1000
gam_w = 6.12
pi_w = 3.43e8
eps = 1e-9

L = 0.6
Nx = 6000
W = 0.06
Ny = 600

case = {
    # Logistics
    "run_time_info": "T",
    # Computational Domain Parameters
    "x_domain%beg": 0.0,
    "x_domain%end": L,
    "y_domain%beg": 0.0,
    "y_domain%end": W,
    "m": Nx,
    "n": Ny,
    "p": 0,
    # time set
             "cfl_adap_dt": "T",
             "cfl_target": 0.1,
             "n_start": 0,
             "t_save": 1e-7,
             "t_stop": 5e-6,

    #"parallel_io": "F" if args.mfc.get("mpi", True) else "F",
    # Simulation Algorithm Parameters
    "model_eqns": 2,  # Use the 5 equation model
    "num_fluids": 2,  # One fluids: air
    "alt_soundspeed": "F",
    "num_patches": 6,
    "mpp_lim": "F",  #Mixture physical parameters limits
    "mixture_err": "T",  # Mixture properties correction
    "time_stepper": 1,  # Runge–Kutta order [1-3]

            "relax": "F",
            "relax_model": 5,
            "palpha_eps": 1.0e-8,
            "ptgalpha_eps": 1.0e-2,
    
    "recon_type": 2,  # Reconstruction Type: [1] WENO; [2] MUSCL
    
    "weno_order": 3,   # WENO order [1,3,5] 
    "weno_eps": 1e-16,  #WENO perturbation (avoid division by zero)
    "weno_avg": "F",
    "mapped_weno": "T",  #WENO-M (WENO with mapping of nonlinear weights)
    "mp_weno": "F",   #Monotonicity preserving WENO
    
    "muscl_order": 2,
    "muscl_lim": 4,
    #"int_comp": "T",

    "riemann_solver": 2,   #iemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact*; [4] HLLD (only for MHD) [5] Laxf
    "wave_speeds": 1, #[1] Direct (Batten et al. [5]); [2] Pressure-velocity* (Toro [51])
    "avg_state": 2,  #[1] Roe average*; [2] Arithmetic mean
    "viscous": "F",
    
    # set BC
    "bc_x%beg": -2,
    "bc_x%end": -2,
    "bc_y%beg": -2,
    "bc_y%end": -2,
    
    # Chemistry
    #"chemistry": "F" if not args.chemistry else "T",
    "chemistry": "T",
    "chem_params%diffusion": "F",
    "chem_params%reactions": "T",
    "cantera_file": ctfile,
    # Formatted Database Files Structure Parameters
            "format": 1,
            "precision": 2,
            "prim_vars_wrt": "T",
            #"E_wrt": "T",
            "parallel_io": "T",
            "c_wrt": "T",
            "chem_wrt_T": "T",
            #"schlieren_wrt": "T",

     # Patch 1: amb
    "patch_icpp(1)%geometry": 3,
    "patch_icpp(1)%x_centroid": L/2,
    "patch_icpp(1)%y_centroid": W/2,
    "patch_icpp(1)%length_x": L,
    "patch_icpp(1)%length_y": W,
    "patch_icpp(1)%vel(1)": 0.0e00,
    "patch_icpp(1)%vel(2)": 0.0e00,
    "patch_icpp(1)%pres": sol_amb.P,
    "patch_icpp(1)%alpha_rho(1)": (1 - eps) *sol_amb.density,
    "patch_icpp(1)%alpha_rho(2)": eps*rho_w,
    "patch_icpp(1)%alpha(1)": 1 - eps,
    "patch_icpp(1)%alpha(2)": eps,

    # Patch 2: ig
    "patch_icpp(2)%geometry": 2,
    "patch_icpp(2)%alter_patch(1)": "T",
    "patch_icpp(2)%x_centroid": 0.0,
    "patch_icpp(2)%y_centroid": 0.0075,
    "patch_icpp(2)%radius": 0.0075,
    "patch_icpp(2)%length_x": 0.012,
    "patch_icpp(2)%length_y": 0.012,
    "patch_icpp(2)%vel(1)": 10.0e00,
    "patch_icpp(2)%vel(2)": 0.0e00,
    "patch_icpp(2)%pres": sol_ig.P,
    "patch_icpp(2)%alpha_rho(1)": (1 - eps) *sol_ig.density,
    "patch_icpp(2)%alpha_rho(2)": eps*rho_w,
    "patch_icpp(2)%alpha(1)": 1 - eps,
    "patch_icpp(2)%alpha(2)": eps,

    # Patch 3: ig
    "patch_icpp(3)%geometry": 2,
    "patch_icpp(3)%alter_patch(1)": "T",
    "patch_icpp(3)%x_centroid": 0.0,
    "patch_icpp(3)%y_centroid": 0.0225,
    "patch_icpp(3)%radius": 0.0075,
    "patch_icpp(3)%length_x": 0.012,
    "patch_icpp(3)%length_y": 0.012,
    "patch_icpp(3)%vel(1)": 10.0e00,
    "patch_icpp(3)%vel(2)": 0.0e00,
    "patch_icpp(3)%pres": sol_ig.P,
    "patch_icpp(3)%alpha_rho(1)": (1 - eps) *sol_ig.density,
    "patch_icpp(3)%alpha_rho(2)": eps*rho_w,
    "patch_icpp(3)%alpha(1)": 1 - eps,
    "patch_icpp(3)%alpha(2)": eps,

    # Patch 4: ig
    "patch_icpp(4)%geometry": 2,
    "patch_icpp(4)%alter_patch(1)": "T",
    "patch_icpp(4)%x_centroid": 0.0,
    "patch_icpp(4)%y_centroid": 0.0375,
    "patch_icpp(4)%radius": 0.0075,
    "patch_icpp(4)%length_x": 0.012,
    "patch_icpp(4)%length_y": 0.012,
    "patch_icpp(4)%vel(1)": 10.0e00,
    "patch_icpp(4)%vel(2)": 0.0e00,
    "patch_icpp(4)%pres": sol_ig.P,
    "patch_icpp(4)%alpha_rho(1)": sol_ig.density,
    "patch_icpp(4)%alpha_rho(2)": eps*rho_w,
    "patch_icpp(4)%alpha(1)": 1 - eps,
    "patch_icpp(4)%alpha(2)": eps,

    # Patch 5: ig
    "patch_icpp(5)%geometry": 2,
    "patch_icpp(5)%alter_patch(1)": "T",
    "patch_icpp(5)%x_centroid": 0.0,
    "patch_icpp(5)%y_centroid": 0.0525,
    "patch_icpp(5)%radius": 0.0075,
    "patch_icpp(5)%length_x": 0.012,
    "patch_icpp(5)%length_y": 0.012,
    "patch_icpp(5)%vel(1)": 10.0e00,
    "patch_icpp(5)%vel(2)": 0.0e00,
    "patch_icpp(5)%pres": sol_ig.P,
    "patch_icpp(5)%alpha_rho(1)": (1 - eps) *sol_ig.density,
    "patch_icpp(5)%alpha_rho(2)": eps*rho_w,
    "patch_icpp(5)%alpha(1)": 1 - eps,
    "patch_icpp(5)%alpha(2)": eps,

    # Patch 6: droplet
    "patch_icpp(6)%geometry": 2,
    "patch_icpp(6)%x_centroid": 0.2,
    "patch_icpp(6)%y_centroid": 0.03,
    "patch_icpp(6)%radius": 0.015,
    "patch_icpp(6)%length_x": 0.06,
    "patch_icpp(6)%length_y": 0.06,
    "patch_icpp(6)%alter_patch(1)": "T",
    "patch_icpp(6)%vel(1)": 0.0e00,
    "patch_icpp(6)%vel(2)": 0.0e00,
    "patch_icpp(6)%pres": 101325.0,
    "patch_icpp(6)%alpha_rho(1)": eps*sol_amb.density,
    "patch_icpp(6)%alpha_rho(2)": (1 - eps) * rho_w,
    "patch_icpp(6)%alpha(1)": eps,
    "patch_icpp(6)%alpha(2)": 1 - eps,

    
    # Fluids Physical Parameters
    "fluid_pp(1)%gamma": 1.0 / (1.4 - 1.0),
    "fluid_pp(1)%pi_inf": 0,
            
    "fluid_pp(2)%gamma": 1.0e00 / (gam_w - 1.0e00),
    "fluid_pp(2)%pi_inf": pi_w * gam_w / (gam_w - 1.0e00),
}

for i in range(len(sol_amb.Y)):
    case[f"chem_wrt_Y({i + 1})"] = "T"
    case[f"patch_icpp(1)%Y({i+1})"] = sol_amb.Y[i]
    case[f"patch_icpp(2)%Y({i+1})"] = sol_ig.Y[i]
    case[f"patch_icpp(3)%Y({i+1})"] = sol_ig.Y[i]
    case[f"patch_icpp(4)%Y({i+1})"] = sol_ig.Y[i]
    case[f"patch_icpp(5)%Y({i+1})"] = sol_ig.Y[i]
    case[f"patch_icpp(6)%Y({i+1})"] = sol_amb.Y[i]
        
if __name__ == "__main__":
    print(json.dumps(case))

[sbryngelson]

Root cause: multi-fluid + chemistry coupling is not implemented

The NaNs are expected with the current code — MFC's chemistry path assumes a single reacting ideal-gas fluid and is not coupled to the multi-fluid (5-equation, model_eqns = 2) model. Single-fluid + chemistry works; adding the water droplet (num_fluids = 2, a stiffened-gas fluid) breaks two things in src/common/m_variables_conversion.fpp:

1. Partial densities are overwritten with the total density (chemistry branch of s_convert_conservative_to_primitive_variables):

if (chemistry) then
    rho_K = 0._wp
    do i = eqn_idx%species%beg, eqn_idx%species%end
        rho_K = rho_K + max(0._wp, qK_cons_vf(i)%sf(j, k, l))   ! total density = Σ(ρYₖ)
    end do
    do i = 1, eqn_idx%cont%end          ! i = 1 .. num_fluids
        qK_prim_vf(i)%sf(j, k, l) = rho_K   ! sets ALL partial densities to the total
    end do
    ...

For num_fluids = 1 this is correct. For num_fluids = 2 it overwrites both partial densities (gas and water) with the species-summed total density, destroying the multi-fluid partial-density / volume-fraction structure.

2. The Cantera temperature inversion is fed physically inconsistent state (s_compute_pressure, chemistry branch):

Y_rs(:)  = rhoYks(:)/rho
e_Per_Kg = energy/rho            ! total energy / total density
call get_temperature(e_Per_Kg - Pdyn_Per_Kg, T_guess, Y_rs, .true., T)
call get_pressure(rho, T, Y_rs, pres)

Temperature is recovered from the total internal energy using Cantera's ideal-gas EOS, with pi_inf ignored entirely. The water droplet is a stiffened gas with pi_inf ≈ 4×10⁸, so its internal energy is enormous. Cantera then tries to find an ideal-gas temperature that matches that energy, which diverges / returns garbage → NaN, and it propagates into the "NaN(s) in timestep output" you see.

In other words, there is no consistent formulation in the code today for "fluid 1 reacts via Cantera, fluid 2 is an inert stiffened gas." This is an unimplemented feature, not a one-line bug, and there is currently no input-validation guard to reject the combination (so it silently NaNs rather than erroring out).

Suggested path forward

• Short term: add a validation check (toolchain/mfc/case_validator.py, mirroring the existing mhd and num_fluids != 1 guard, plus a Fortran-side @:PROHIBIT in src/common/m_checker_common.fpp) that rejects chemistry with num_fluids > 1 with a clear message, so users aren't left debugging NaNs.
• Long term: proper coupling of a reacting gas phase with inert (stiffened-gas) fluids — separating the gas-phase density/energy used for the Cantera temperature/pressure inversion from the total mixture, and preserving partial densities and volume fractions through the chemistry conversion path.

@DimAdam-01 — flagging for you as the chemistry owner / assignee.