[drcole17]

Adds a temperature-dependent surface tension coefficient so the CSF capillary model can drive thermocapillary (Marangoni) flow:

$$\sigma(T) = \sigma + \frac{\mathrm{d}\sigma}{\mathrm{d}T},(T - T_\mathrm{ref})$$

Selected by sigma_model = 1 with sigma_T_ref and sigma_dTdT. The default sigma_model = 0 keeps the existing constant-sigma path; no existing golden moved.

Approach

s_get_capillary evaluates σ(T) per cell over the ghost-cell bounds (c_sigma, after the color-gradient reconstruction), and s_compute_capillary_source_flux averages it to each face before building the capillary stress tensor and the energy term. Nothing else in the CSF model changes: the tangential Marangoni force comes out of the face-to-face variation of the existing stress tensor, so there is no separate source term to keep consistent with it.

Temperature source: the q_T_sf that Fourier conduction (#1906) already maintains. sigma_model = 1 therefore requires some fluid_pp(i)%k_therm > 0, enforced in case_validator.py. Reasons:

• Under heat_conduction, q_T_sf is already allocated, filled in the ghost cells at every BC type, and MPI-exchanged in s_populate_variables_buffers, which runs before s_get_capillary. Reading it adds one argument to s_get_capillary and nothing to the boundary or MPI code. A separate temperature for σ would duplicate that ~36-site ghost fill or leave σ wrong next to boundaries and rank faces.
• It is the same mixture temperature conduction diffuses, so σ and the heat flux can never disagree about T.
• Thermocapillary migration is a conduction-coupled problem: the imposed gradient has to be sustained against the drop's own advection, which needs bulk conduction and isothermal walls (bc = -15/-16) anyway.
• k_therm is already rejected with chemistry, so σ(T) never has to choose between conduction's temperature and chemistry's.

The cost is that a purely advected (zero-diffusivity) temperature field cannot be expressed; a small k_therm is required. The Ma → 0 benchmark below is reached through large diffusivity, not frozen T.

Parameters go through the registry (definitions.py, SIM_GPU_DECL_VARS), so the declarations, namelist binding, broadcast, and GPU_DECLARE are generated. sigma_model carries choices/value_labels and a when_value dependency requiring both other parameters; the validator also rejects sigma_model without surface_tension and sigma_T_ref/sigma_dTdT without sigma_model = 1.

Verification

Regression. New golden E39EFF77 (2D -> 2 Fluid(s) -> capillary=T -> model_eqns=3 -> sigma_model=1): the existing 2D capillary case with conduction on and patches at different temperatures, so σ differs across the interface. Rerunning the trace with σ held constant (sigma_dTdT = 0) misses the golden by up to 14.6% relative (y-momentum at step 50), while an unmodified rerun matches it exactly, so the golden pins the σ(T) path and not just the constant one. No existing golden changed. The golden was generated with Intel ifx; built with gfortran 15 + Intel MPI, all three capillary=T tests (3859AF63, AA9C637F, E39EFF77) also pass against their goldens.

Full suite (./mfc.sh test --no-debug, Intel ifx 2026.1 + Intel MPI, CPU): 723 passed, 39 skipped. (13 tests, none involving surface tension, failed on the first pass for local reasons: the default srun launcher had no PMI for Intel MPI on the development host, and a concurrent validate rewrote a Cantera IC cache. All 13 pass on rerun with -- -b mpirun.) ./mfc.sh precheck: all checks pass.

Validation: examples/2D_thermocapillary_migration, against the σ(T)-only benchmarks of Samareh, Mostaghimi & Moreau, Int. J. Heat Mass Transfer 73 (2014) 616–626. Intel ifx + Intel MPI, CPU, 64 cells across the box, --case-optimization:

Case	Samareh	Reference	MFC	Wall time
TC1, Ma → 0, identical fluids (case_Ma_0p001.py)	§4.1.1, Fig 5	terminal v/v_YGB ≈ 0.80	plateau 0.857 at t/t_r = 0.6; 0.876 at t/t_r = 2	3.7 h, 6 ranks
TC2, Re = 5, Ma = 20, distinct fluids (case_Ma_20.py)	§4.1.2, Fig 7 (Nas & Tryggvason)	peak U* ≈ 0.13	peak 0.134 at t* = 3.7	18.5 min, 8 ranks

The drop migrates toward the hot wall in both. TC2's U* decays steadily after the peak and is 0.071 at t* = 15 when the run ends, so the late-time value is not converged and is not claimed. TC1 sits 7% above Samareh at 12.8 cells per D (TC2 has 32) and then creeps up at +0.005 per t_r. An earlier, pre-#1906 implementation of this feature with its own conduction measured 0.81, 0.89 and 0.87 at 12.8, 25.6 and 51.2 cells per D, so TC1 is not grid-converged at this resolution in either version. The gap between the two at the same grid (0.857 vs 0.81) is not isolated: both the conduction model and the y walls (reflective -2 → slip -15, which isothermal walls require) changed. Finer TC1 grids were not rerun here; the time step is conduction-limited (dt ∝ Δx²), so the next grid costs ~16× (~60 h on 6 ranks).

Notes for reviewers

• c_sigma is real(wp), not stp, like the gL_x/gR_x arrays beside it in the same module. Under --mixed (stp = half) σ would be quantized at ~1e-3 relative, which is the size of the cell-to-cell σ difference that is the Marangoni force in TC1. It is one scalar field per rank.
• c_sigma is allocated (1,1,1) when sigma_model = 0, the existing pattern for optional module arrays, so the GPU declare and deallocate stay unconditional.
• s_get_capillary gained a q_T_sf argument. It has one call site (m_rhs.fpp).
• The TC2 initial condition floors both volume fractions at 1e-9. With each phase exactly empty on the far side of the interface, the run reached a cell holding α₁ = 1.04e-16 with α₁ρ₁ = 0 at step 14846, which pressure relaxation's non-physical-density check (State-dependent equations of state: Mie-Gruneisen (linear and cubic Hugoniot), JWL and Vinet, with temperature, per-phase evaluation on the 5- and 6-equation models #1811) rejects. The peak U* was the same (0.1340 at t* = 3.74) in the aborted run.

Known limitations

• The temperature must come from conduction (above).
• σ(T) is only physical while it stays positive, |T − T_ref| < σ/|dσ/dT|. Documented in case.md, not checked at runtime.
• MFC's compressibility enters the benchmarks: the closed box rings acoustically (the velocity is read as the oscillation mean), and in TC2 each fluid's absolute density stratifies as ~1/T, which Samareh's incompressible solver does not have. TC1's identical fluids avoid the latter.
• GPU: compiles under NVHPC 25.5 (the nvcr.io/nvidia/nvhpc:25.5-devel-cuda12.9 container, MFC_CUDA_CC=70,75,80,89,90) for --gpu acc and --gpu mp, and with --case-optimization --gpu acc on benchmarks/5eq_rk3_weno3_hllc. Not run on a GPU; c_sigma and the three parameters are device-declared and updated in s_initialize_gpu_vars.
• Only a 2D golden. The 3D id == 3 face average is the same one-line pattern as x and y but has no trace.

Size

Two commits. The feature (b7971324) is about 100 changed lines of source, toolchain and docs: 50 in src (mostly m_surface_tension.fpp), the rest the registry entries, validator rule and case.md. The golden adds 216 lines. The validation example (6cda6754, 895 lines including the measurement script and a schematic) is separate and can be dropped from this PR if you would rather not carry it.

Pre-existing issue found (not addressed here)

s_initialize_surface_tension_module allocates c_divs(j)%sf for j = 1, num_dims + 1, but s_finalize_surface_tension_module deallocates only j = 1, num_dims, so the |∇c| field c_divs(num_dims + 1)%sf (host and device) is never freed. Harmless at program exit; left alone to keep this diff to σ(T).

---

Contribution Policy

Written with Claude Code. The physics is checked against a published benchmark, as above. It is a port of a feature developed and validated on a fork, rebuilt onto upstream conduction (#1906).

Acknowledgement

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

Left unticked deliberately: that attestation is the author's to make, not the tool's.

🤖 Generated with Claude Code

https://claude.ai/code/session_01VtsGVvCfsDJz5dmtSk9dKt

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_surface_tension.fpp	292	+29
src/simulation/m_global_parameters.fpp	798	+3

Directory	Lines	Diff
simulation	28563	+32
total	47622	+32

[github-actions]

Claude Code Review

Head SHA: 6cda675

Files changed:

• 17
• docs/documentation/case.md
• examples/2D_thermocapillary_migration/case_Ma_20.py
• examples/2D_thermocapillary_migration/run.py
• src/simulation/m_global_parameters.fpp
• src/simulation/m_rhs.fpp
• src/simulation/m_start_up.fpp
• src/simulation/m_surface_tension.fpp
• toolchain/mfc/case_validator.py
• toolchain/mfc/params/definitions.py
• toolchain/mfc/test/cases.py

Findings:
No issues found. Checked in step 4: src/simulation/m_global_parameters.fpp and src/simulation/m_surface_tension.fpp, plus the m_rhs.fpp, m_start_up.fpp and toolchain hunks in the diff. m_rhs.fpp has no context file, so the claim that q_T_sf is allocated and filled (including ghost cells) before s_get_capillary runs when sigma_model = 1 was not checked. The reviewer should confirm that the conduction-temperature field is populated at that point in the RHS.

[codecov]

Codecov Report

❌ Patch coverage is 69.23077% with 8 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.88%. Comparing base (3ef5e1e) to head (6cda675).

Files with missing lines	Patch %	Lines
src/simulation/m_surface_tension.fpp	63.63%	4 Missing and 4 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2102      +/-   ##
==========================================
- Coverage   61.89%   61.88%   -0.01%     
==========================================
  Files          86       86              
  Lines       22888    22906      +18     
  Branches     3372     3377       +5     
==========================================
+ Hits        14166    14176      +10     
- Misses       6220     6224       +4     
- Partials     2502     2506       +4     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.