Isothermal immersed boundaries without chemistry (thermal_bc = 1 on the k_therm path)

[BCKim55]

Contribution Policy

We do not accept pull requests generated primarily by AI without genuine understanding or real-world usage context.

All contributions are expected to demonstrate:

• A clear understanding of the codebase
• Alignment with product direction
• Thoughtful reasoning behind changes
• Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.

---

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

---

PR template credit: junegunn

---

Summary

patch_ib%thermal_bc = 1 (prescribed Twall) now also works with chemistry = F, on the Fourier conduction path of #1906 (fluid_pp%k_therm). Until now an IB without chemistry was always adiabatic. Non-reacting single-gas cases (heated/cooled particle beds, shock–particle) need a fixed surface temperature without turning on the chemistry machinery.

Discussed on the CM3C Slack and at the CM3C weekly meeting with @sbryngelson and @danieljvickers.

Implementation

In s_ibm_correct_state, for chemistry = F and thermal_bc = 1 only:

• The image-point temperature comes from s_phase_temperature, and the ghost density from the new s_phase_density in m_eos (the inverse of s_phase_temperature at fixed pressure). The ghost-point kernel holds no EOS coefficients.
• Ghost temperature T_g = T_w + θ (T_w − T_IP), with θ = 1 (the mirror) unless the temperature limiter cuts it. The pressure and velocity are unchanged.
• q_T_sf at the ghost is therefore T_g, and the existing m_conduction face flux gives the wall heat flux. There is no new kernel.
• s_blend_ghost_state is split into s_blend_ghost_temperature (which takes its bounds as arguments) and s_blend_ghost_species. With chemistry the bounds are the NASA window, unchanged. Without chemistry the only bound is positivity, since there is no fit range. Each bound is compared before dividing, so T_hi = huge() cannot overflow.
• thermal_bc = 0 without chemistry is untouched (the existing mirror is already adiabatic).
• Validator: without chemistry, thermal_bc = 1 requires k_therm > 0, num_fluids = 1 and Twall > 0, and any thermal_bc other than 0 or 1 is rejected. k_therm already limits the EOS to ideal or stiffened gas, so the state-dependent families never reach s_phase_density. The 200–5000 K Twall window now applies only with chemistry.

No new parameter: heat_conduction (any k_therm > 0) acts as the switch. A temperature flag separate from chemistry, as Dan suggested, can come later.

Verification

A cylinder (2D) or sphere (3D) held at Twall in quiescent air (ideal gas, constant k) in a periodic box, run to τ = αt/R² = 1. This is pure conduction: no gravity, and no flow except the expansion or contraction the heating drives.

Reference. A 1-D radial solution of the low-Mach equations in the same closed volume (Majda & Sethian 1985): variable density ρ = p₀/(R T), the induced radial velocity from ∇·u, and the uniform thermodynamic pressure p₀(t) of the closed box. These are the equations MFC solves as Ma → 0. As ΔT → 0 the solver reduces to the constant-density solution and, with an open far field, to Carslaw & Jaeger. Errors below are for the absorbed heat (MFC: change of the total gas energy in the closed box) and the rms gas temperature.

Reference solution: equations, numerics and sources

Radial domain R < r < R_out, where R_out gives the same gas volume as the periodic box (d = 2 cylinder, d = 3 sphere), with an adiabatic, impermeable outer wall:

$$\rho c_p\left(\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial r}\right)=\frac{1}{r^{d-1}}\frac{\partial}{\partial r}\left(r^{d-1}k\frac{\partial T}{\partial r}\right)+\frac{dp_0}{dt},\qquad \rho=\frac{p_0}{R_g T}$$

$$\frac{1}{r^{d-1}}\frac{\partial\left(r^{d-1}u\right)}{\partial r}=\frac{(\gamma-1),\nabla\cdot(k\nabla T)-dp_0/dt}{\gamma,p_0},\qquad \frac{dp_0}{dt}=\frac{\gamma-1}{V_{gas}},\dot Q_{gas}$$

with T(R, t) = T_w, u(R) = u(R_out) = 0, T(r, 0) = T_∞ and p_0(0) = p_∞. Here \dot Q_gas is the heat flow from the wall into the gas. The absorbed heat follows from the thermodynamic pressure, H(t) = V_gas (p_0(0) − p_0(t))/(γ − 1), which is the change of the total gas energy that MFC measures.

Numerics: finite volume with 1500 cells clustered at the wall, BDF time integration. Grid-converged to 1e-6 in H, and gas mass conserved to 1e-10 (ΔT = 6 K) and 1e-5 (ΔT = 300 K). As ΔT → 0 it matches the constant-density solution to 1e-4. That solution, in turn, reproduces the Carslaw & Jaeger sphere closed form and cylinder integral to 1e-5 when given an open far field and dp_0/dt = 0.

• A. Majda & J. Sethian, "The derivation and numerical solution of the equations for zero Mach number combustion," Combust. Sci. Technol. 42, 185–205 (1985).
• H. S. Carslaw & J. C. Jaeger, Conduction of Heat in Solids, 2nd ed., Oxford University Press (1959): §9.10 (sphere) and §13.5 (cylinder).

Small ΔT (wall 294 K, gas 300 K):

case	absorbed heat at τ = 1	temperature rms / ΔT
2D, D/Δx = 10 / 20 / 40 / 80	+4.2 / +2.6 / +1.5 / +0.8 %	0.67 / 0.37 / 0.24 / 0.12 %
3D, D/Δx = 10 / 20 / 40 (MI300A)	+7.8 / +3.7 / +2.0 %	1.25 / 0.53 / 0.29 %

Large ΔT, compressibility (wall 600 K, gas 300 K; the near-wall density halves):

D/Δx	heat vs constant-density reference	heat vs low-Mach reference	temperature rms / ΔT (low-Mach)
10	−9.6 %	+3.4 %	0.70 %
20	−10.8 %	+2.0 %	0.35 %
40	−11.4 %	+1.3 %	0.22 %

MFC converges to the variable-density solution and not to the constant-density one, so the density change and the flow it drives are captured. Both errors are first order, as expected for a ghost-cell IB.

Other checks

• This branch reproduces an earlier implementation (the same scheme with explicit ideal-gas algebra) to round-off: ρ and E agree to 1e-13.
• GPU, Cray ftn OpenMP offload on MI300A, matches CPU: 2D D/Δx = 20 absorbed heat 1.7632e-6 J/m on both, and 3D D/Δx = 10 identical to the earlier GPU run.
• Limiter: gas at 700 K, wall at 293 K. The mirror would be negative; the run stays stable with no NaN.
• Regression: all IBM tests and the reacting-surface example (F52F0D4C) pass unchanged. One new test, 2D -> 1 Fluid(s) -> IBM -> Isothermal (5BBDFB29), whose golden holds the isothermal ghost densities.

Relation to open PRs

• IBM: extrapolate hot/reacting-wall ghost state from where the image point is actually sampled #1969, Add ib_surface_wrt: wall temperature and gasified mass flux at each IB surface point #1956, Add thermal_bc = 3: a lumped solid whose wall temperature follows its heat balance (after #1956) #1957 change the same s_blend_ghost_state call (and Add thermal_bc = 3: a lumped solid whose wall temperature follows its heat balance (after #1956) #1957 the Twall window check and the thermal_bc docs), so whichever lands second needs a small rebase; I'll rebase this one.
• IBM: extrapolate hot/reacting-wall ghost state from where the image point is actually sampled #1969's theta_max = |levelset|/d applies to this path too: once it lands, the non-chemistry branch should pass it to s_blend_ghost_temperature, which should also improve the heat convergence above.
• With Add thermal_bc = 3: a lumped solid whose wall temperature follows its heat balance (after #1956) #1957, thermal_bc = 3 without chemistry is rejected by the validator here rather than silently ignored.
• With Add patch_ib%thermal_layer: start a hot IB in its conduction solution #1955, thermal_layer becomes reachable without chemistry; its initial density ρ ∝ 1/T at fixed pressure holds for ideal gas but not for stiffened gas.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_ibm.fpp	1857	+16
src/common/m_eos.fpp	525	+7

Directory	Lines	Diff
common	10468	+7
simulation	28461	+16
total	47462	+23

[github-actions]

Claude Code Review

Head SHA: 4442ef8

The review produced no output (Claude step: failure). Run and transcript: https://github.com/MFlowCode/MFC/actions/runs/38328154585