[drcole17]

Fixes two missing terms in 2D axisymmetric (cyl_coord, grid_geometry == 2) surface tension: the azimuthal (hoop) capillary stress, and the capillary flux divergence on the axis row. Without them a static axisymmetric drop reaches at most half its Laplace pressure (0.052 of 0.1), develops 6.6× capillary-velocity currents where the interface meets the axis, and aborts within a third of a capillary time.

The capillary stress $\boldsymbol\Omega = -\sigma\left(|\nabla c|,\mathbf I - \frac{\nabla c\otimes\nabla c}{|\nabla c|}\right)$ enters the momentum flux. In axisymmetric coordinates its radial divergence is

$$(\nabla\cdot\boldsymbol\Omega)_r = \partial_x\Omega_{xr} + \frac{1}{r}\partial_r(r,\Omega_{rr}) - \frac{\Omega_{\theta\theta}}{r},\qquad \Omega_{\theta\theta} = -\sigma|\nabla c|$$

since $\nabla c$ has no azimuthal component. Upstream supplies the first two terms through the face fluxes and the generic $-F/r$ geometric source, but never $-\Omega_{\theta\theta}/r$ — the second principal curvature of an axisymmetric interface. Separately, the axis branch of m_rhs skips the face-flux divergence of the axis cell (k = 0) and reads tau_Re_vf there instead, which only the viscous and conduction axis sources fill; the capillary contribution in that row was dropped.

Approach

Two kernels in m_surface_tension.fpp, both gated on surface_tension .and. grid_geometry == 2:

• s_compute_capillary_hoop_source adds $-\sigma\overline{|\nabla c|}/r$ to the radial momentum RHS for k >= 1 (or k >= 0 off-axis). $\overline{|\nabla c|}$ is averaged over the same two faces, with the same capillary_cutoff, as the $-F/r$ geometric source it pairs with, so $-(\Omega_{rr}-\Omega_{\theta\theta})/r$ is a per-face normal force.
• s_compute_capillary_axis_source fills tau_Re_vf(mom%beg:E) at k = -1, 0, 1 with the cell-centered radial row of $\boldsymbol\Omega$ and its energy flux $\boldsymbol\Omega\cdot\mathbf u + \sigma|\nabla c|u_r$, following s_compute_conduction_axis_source. The axis cell's hoop term ($+\sigma|\nabla c|$) is added to tau_Re_vf(mom%beg+1) at k = 0 only, because m_rhs uses that single value in its $-\tau/r$ axis geometric term. It accumulates on top of the viscous and conduction axis sources, and zeroes the fields itself when neither ran.

m_rhs.fpp allocates tau_Re_vf(mom%beg:E) for axisymmetric surface tension (previously viscous or conduction only), calls the axis kernel after the conduction one, and widens the two existing axis-row k = 0 loops to include it. The energy equation needs no hoop term: $\nabla\cdot(\boldsymbol\Omega\cdot\mathbf u)$ is a vector divergence, which the existing geometric source already completes.

Cartesian and 3D cylindrical code paths are unchanged.

Verification

Static-sphere Laplace gate. $R=1$ drop centered on the axis in the quarter domain $[0,2.5]^2$ (symmetry at $x=0$), $\sigma=0.05$, $\rho_l/\rho_g = 10$, 6-equation model, HLLC, WENO5, $\Delta t = 0.125/N$, 4 ranks (2×2, so both axis and off-axis ranks run the new code). Target $\Delta p = 2\sigma/R = 0.1$; $t_c=\sqrt{\rho_l R^3/\sigma}=4.47$, $U_c=\sqrt{\sigma/(\rho_l R)}=0.224$.

Run	$\Delta p$	Peak $\lvert\mathbf u\rvert$	Aborts at
master, N=50	0.016–0.052 (saved frames $t$ = 0.5, 1.0)	1.47 (6.6 $U_c$), at the axis crossing	$t=1.37$ (0.31 $t_c$)
master, N=50, 0.4 Δt	—	—	$t=1.36$
master, N=100	—	—	$t=0.25$ (0.06 $t_c$)
this PR, N=50	0.0996–0.1009 for $t>6$ (mean 0.1002)	0.23 (1.0 $U_c$) mid-latitude, $t$=6–18; 0.06–0.07 after $t=18$	$t=35.7$ (8.0 $t_c$)
this PR, N=100	0.030–0.164 ringing ($t\le3$), 0.108–0.109 ($t$=3–4)	0.22	$t=4.11$ (0.92 $t_c$)
this PR, N=100, 0.5 Δt	—	—	$t=3.74$
Cartesian cylinder N=50, master and this PR	0.0498–0.0502 (target $\sigma/R = 0.05$)	0.026–0.061	—

Flat interface crossing the axis. A planar $x=0.5$ interface has zero curvature and must stay at rest. Max $|v_r|$ near the axis after 50 steps: Cartesian 1.6e-12, master axisymmetric 0.293, this PR 7.2e-12. This is the dropped axis-row flux in isolation.

Cartesian inertness. The Cartesian N=50 gate run is bit-identical to master across all 4212 output files.

New golden AC1AA22C (2D -> Axisymmetric -> model_eqns=3 -> capillary=T): the flat axis-crossing interface above plus an $r=0.8$ interface that carries the hoop stress, $\sigma=0.01$. It fails on master (max error 2.27e-1) and passes here bitwise. The existing capillary goldens 3859AF63 and AA9C637F are unchanged; no existing golden file changed.

Full suite (gfortran 15, MPI, CPU): 723 passed, 0 failed, 39 skipped. Run without -a, because h5dump is not installed on the development host, so the post-process NaN/Inf check was not exercised locally.

./mfc.sh precheck: all checks pass.

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.

Notes for reviewers

• tau_Re_vf is now allocated for axisymmetric surface tension without viscous or conduction, under the same condition as its deallocation. Only tau_Re_vf(mom%beg:E) is used, as for conduction.
• Call order matters. s_compute_capillary_axis_source accumulates into tau_Re_vf, so it must follow the viscous and conduction axis sources, which assign. When neither is active it zeroes the fields itself.
• k_beg in the hoop kernel tests bc_y%beg == BC_REFLECTIVE to mirror m_rhs's bc_y%beg == -2 .or. -14 axis branch; for 2D axisymmetric runs the validator already forces -2 at the axis, and off-axis MPI ranks see a processor boundary and start at k = 0.
• Why face-averaged rather than cell-centered. A cell-centered hoop term, $-\sigma|\nabla c|{j,k}/r_k$, gives the same $\Delta p$ but does not pair face-for-face with the $-\Omega{rr}/r$ geometric source. It aborted at 2.55 $t_c$ (N=50) and 0.34 $t_c$ (N=100), against 8.0 and 0.92 $t_c$ for the face-averaged form here.
• Net +98 lines of Fortran (+92 in m_surface_tension.fpp, +6 in m_rhs.fpp), plus a 26-line test case.

Known limitations

• Parasitic currents remain. At N=50 they peak near 1.0 $U_c$ at mid-latitude and reach 0.08–0.25 $U_c$ (0.018–0.055) where the interface crosses the axis. The Cartesian cylinder at the same resolution peaks at 0.27 $U_c$.
• Static drops still eventually abort: N=50 at 8.0 $t_c$; N=100 at 0.92 $t_c$, and also with half the time step, so it is not a time-step limit. No precursor shows in the saved fields.
• Translating drops are not fixed. In the fork this change came from, a translating axisymmetric drop still dies near $t/t_c \approx 0.21$. Not re-run on this branch.
• The real fix for all three is a balanced-force curvature discretization, which is out of scope here.
• 3D cylindrical + surface tension is untouched and unvalidated. By inspection, the shared m_rhs axis branch drops its axis-row capillary flux the same way; not run.
• The existing step-interface capillary layout (sigma = 1) blows up in Cartesian coordinates too, which is why the new test sets $\sigma = 0.01$.

Pre-existing issues found (not addressed here)

• s_finalize_surface_tension_module deallocates c_divs(1:num_dims)%sf but allocates c_divs(1:num_dims+1)%sf, so the $|\nabla c|$ field leaks.

---

Contribution Policy

Written with Claude Code. Reimplemented from fixes developed and exercised in my fork's axisymmetric droplet runs, with the signs re-derived and every claim above re-measured on this branch except where marked.

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	327	+64
src/simulation/m_rhs.fpp	1972	+4

Directory	Lines	Diff
simulation	28599	+68
total	47658	+68

[github-actions]

Claude Code Review

Head SHA: aad40cc

Files changed:

• 3
• src/simulation/m_rhs.fpp
• src/simulation/m_surface_tension.fpp
• toolchain/mfc/test/cases.py

Findings:

• medium src/simulation/m_surface_tension.fpp: In s_compute_capillary_axis_source and s_compute_capillary_hoop_source, w3 is hard-coded to 0 (s_compute_capillary_stress_tensor(sigma, w1, w2, 0._wp, ...)), but normW is c_divs(num_dims+1), which includes c_divs(3)**2 when p > 0. For 3D cylindrical grids (grid_geometry == 2, p > 0) the tensor then mixes a norm containing the z/theta gradient with components that omit it. The hoop term also uses that full normW. The existing face-flux routine passes gL_x(...,3) when p > 0. I could not confirm that the new paths are valid for 3D cylindrical. Either pass c_divs(3) when p > 0, or add a validator or test that restricts the new terms to the 2D axisymmetric case. Confidence: moderate.
• medium src/simulation/m_surface_tension.fpp: The zeroing of tau_Re_vf(mom%beg:E) in the axis source depends on .not. (viscous .or. heat_conduction). The docstring requires the new routine to run after the viscous and conduction axis sources. When only heat_conduction is on, or only viscous, I could not check that the other source initialises the full mom%beg:E range at k = -1..1. m_rhs.fpp is not in the context directory, so this is unverified. If either routine skips momentum or energy rows, the new code would accumulate onto uninitialised memory. Please confirm both routines zero all of mom%beg:E for k = -1:1, or zero unconditionally and let the other sources add.
• medium toolchain/mfc/test/cases.py: The new case redefines model_eqns in both the define_case_d modifier list and the parameter dict, and stack is unchanged by it. Check that the case's trace and UUID are stable and that grid_geometry = 2 is actually active here, since the other alter_2d cases are Cartesian. The comment says the case is meant to hit the axis and hoop paths, but nothing in the hunk sets cylindrical geometry. If it is not set elsewhere in alter_2d, the new routines are never exercised. Confidence: low to moderate.

[codecov]

Codecov Report

❌ Patch coverage is 76.31579% with 9 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.92%. Comparing base (3ef5e1e) to head (aad40cc).

Files with missing lines	Patch %	Lines
src/simulation/m_surface_tension.fpp	80.64%	0 Missing and 6 partials ⚠️
src/simulation/m_rhs.fpp	57.14%	0 Missing and 3 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2103      +/-   ##
==========================================
+ Coverage   61.89%   61.92%   +0.03%     
==========================================
  Files          86       86              
  Lines       22888    22922      +34     
  Branches     3372     3380       +8     
==========================================
+ Hits        14166    14195      +29     
  Misses       6220     6220              
- Partials     2502     2507       +5     

☔ 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.