Consistent NEMO viscous Jacobian: distinct left block and subtractive edge assembly by bellonarts · Pull Request #2885 · su2code/SU2

bellonarts · 2026-09-05T20:58:07Z

Proposed Changes

The gradient-corrected NEMO viscous routine returns the right Jacobian block twice, and the solver assembles viscous Jacobians with additive signs despite subtracting the flux at node i and adding it at node j. This patch returns the distinct left and right blocks and uses UpdateBlocksSub for the matching assembly signs. The residual equations are unchanged.

For the controlled constant density and viscosity momentum test, Newtonian stress gives raw flux derivatives of -8/15 at i and +8/15 at j. The opposite signs inside numerics come from the velocity difference; they do not include the residual assembly sign. The assembled blocks must therefore be (-Ji, -Jj, +Ji, +Jj). The general viscous Jacobian remains a TSL approximation.

Tests

The numerics unit test checks the independently derived values and central differences of the returned momentum flux at both endpoints, using two perturbation sizes. The unused gradient reconstruction helper is removed.
The solver assembly test checks the signs and placement of the supplied blocks, with explicit finite checks. It is an integration test, not an independent proof of the full Jacobian.
The existing ion_gy regression now runs through iteration 99 using test_iter, the existing mesh and the same stored restart. Its residual threshold is lowered to prevent the previous stop at iteration 10. The duplicate case and config are removed, and the output restart has a separate filename so the input seed stays unchanged.
The installed unit-test container now gets the Mutation++ library path required by the new NEMO coverage.

In the earlier captured two-rank comparison, the final log10 energy residual was 8.78 with the original code and 4.54 with the patch. The patched trajectory stayed bounded but still reported nonphysical states. This is regression stability evidence, not convergence or physical validation. The existing 0.01 tolerance is retained.

In the focused Linux run, both MPI and NoMPI builds passed the full normal unit suite: 36 cases and 73,848 assertions each, including MLPCpp. super_cat passed, and both fresh two-rank ion_gy repeats passed at iteration 99 with identical printed rows and unchanged input seeds. Native checks also confirmed that restoring either original bug makes its relevant test fail.

The cone references were refreshed from identical results in the fork and upstream runs, keeping the existing 1e-5 tolerance. The original regression comparisons failed only those stored cone references. The focused solver runs used cd8a39d; the final commit 50e3343 changes only the two reference arrays and the unit container's library path. Compiled source, unit tests and solver configs are identical.

Final upstream CI now passes the installed normal unit suite (36 cases, 73,848 assertions), AD/DD unit suites, and all 157 serial regression tests. The parallel cone, super_cat, and ion_gy at iteration 99 also pass; the rest of the parallel suite is still running. Both cone comparisons have zero printed difference.

The overlapping visc_cone, super_cat, and ion_gy expectations must be refreshed if another NEMO PR merges first.

Related Work

Possibly relevant to #2717. The inviscid case in #2026 is outside this patch's scope. Earlier matched Mach 7.95 cone runs showed the original code aborting at iteration 1,406 and this patch completing 3,000 bounded iterations; those runs did not establish convergence.

PR Checklist

I am submitting my contribution to the develop branch.
My contribution is commented and consistent with SU2 style.
I have added a test case that demonstrates my contribution.
Documentation is not applicable; no user-facing configuration option or API is added.
Pre-commit formatting checks pass.
The focused Linux checks and final upstream results listed above are verified.

…an and assemble it subtractively CAvgGradCorrected_NEMO::ComputeResidual returned the j-side Jacobian block in both slots, while its sibling CAvgGrad_NEMO returns the distinct blocks. CNEMONSSolver::Viscous_Residual subtracts the viscous residual at iPoint and adds it at jPoint but assembled the Jacobian with UpdateBlocks, the additive convention. With NEMO's viscous Jacobians being flux derivatives (dF/dU_i = -dF/dU_j), UpdateBlocksSub(J_i, J_j) assembles exactly dR/dU in all four blocks; the two hunks are one fix and stay in one commit. The residual is unchanged, so converged solutions do not change; the implicit operator does. Add the distinct-block unit test on a 2-D AIR-5 edge (NEMOViscousFixture) and an assembly-level test that builds a NEMO Navier-Stokes solver on a box mesh, calls Viscous_Residual and checks that the assembled residual and Jacobian blocks equal the derivative of the subtracted and added edge fluxes; both fail on develop and pass with this change. Rebaseline the regression vectors that move as a consequence: visc_cone (serial and parallel), super_cat and ion_gy (parallel). Add a regression case, ion_gy_march, that marches the ionized cylinder 100 iterations from its stored restart with the residual gate disabled: develop diverges (energy residual above 8 and a wall heat flux of order -1e10 by iteration 99) while the corrected operator holds the restart state. Give the two cylinder configurations distinct output restart names so that neither overwrites the restart the other reads. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

Keep git diff --check clean before upstream PR submission.

Update the serial and two-rank viscous-cone aarch64 expectations from GitHub Actions run 33998938091 at exact pre-update head 077f73e. Artifact digest: sha256:e7eedd1985e4e87830bdb740aa1b2a02e2f5047f0afe1a18bb88d8dce9d0b243. No x86 expectation or solver source changes.

Reseal PRs su2code#2883-su2code#2885 at heads b8504d8, 2925b4a, and 8bfba2f, switch their NEMO jobs from harvest to strict vector verification, and apply the public regression drivers' 1e-5 default when tol remains unset. S2/S3 stay harvest-only because they preserve the pre-vector-update cumulative stages.

pcarruscag · 2026-09-07T16:06:56Z

    ion_gy.cfg_file = "cyl_ion_gy.cfg"
    ion_gy.test_iter = 10
-    ion_gy.test_vals = [-11.629873, -4.165562, -4.702662, -4.950351, -5.146155, -4.993878, -6.893332, 5.990109, 5.990004, -0.014849, 0.000000, 90090.000000]
+    ion_gy.test_vals = [-12.682344, -4.169796, -4.714319, -5.658742, -5.712605, -5.269125, -7.961372, 3.359972, 1.705292, -0.014849, 0.000000, 90286.000000]
    test_list.append(ion_gy)

+    # Ionized cylinder marched from its restart with the residual gate disabled: the
+    # corrected NEMO viscous Jacobian holds the restart state, the previous one diverged.
+    ion_gy_march = TestCase('ion_gy_march')
+    ion_gy_march.cfg_dir = "nonequilibrium/visc_cylinder"
+    ion_gy_march.cfg_file = "cyl_ion_gy_march.cfg"
+    ion_gy_march.test_iter = 99
+    ion_gy_march.test_vals = [-11.662039, -4.203178, -4.868257, -5.462497, -5.232052, -4.960881, -6.951391, 4.541901, 4.552855, -0.014861, 0.000001, 90357.000000]
+    ion_gy_march.tol = 0.01
+    test_list.append(ion_gy_march)


This is the same case but with more iterations, right?
If so just modify the case and remove the new config you added, since it's mostly duplicated.
Configs do not have to be copied to change the number of iterations, the test framework modified ITER based on test_iter.

Done. I extended the existing ion_gy test to iteration 99 and removed the duplicate case and config. I also lowered the residual threshold so it reaches that iteration and kept the restart input protected.

pcarruscag · 2026-09-07T16:15:55Z

+  void set_edge_linear_velocity_gradient() {
+    gradient_i = su2double(0.0);
+    gradient_j = su2double(0.0);
+    constexpr su2double distance_squared = 25.0;
+    for (unsigned short component = 0; component < nDim; ++component) {
+      const su2double jump = primitive_j[VEL_INDEX + component] - primitive_i[VEL_INDEX + component];
+      for (unsigned short dimension = 0; dimension < nDim; ++dimension) {
+        const su2double edge_component = coord_j[dimension] - coord_i[dimension];
+        gradient_i(VEL_INDEX + component, dimension) = jump * edge_component / distance_squared;
+        gradient_j(VEL_INDEX + component, dimension) = jump * edge_component / distance_squared;
+      }
+    }
+  }


Same comment from the previous PR about not duplicating code in unit tests so that the comparison is not code against code, the proof that the code is correct should be independent.

Removed the unused gradient helper. The test now checks values derived from Newtonian stress and finite differences of the actual momentum flux at both endpoints. The separate assembly test is described as checking signs and block placement.

pcarruscag · 2026-09-07T16:16:34Z

      LinSysRes.SubtractBlock(iPoint, residual);
      LinSysRes.AddBlock(jPoint, residual);
      if (implicit) {
-        Jacobian.UpdateBlocks(iEdge, iPoint, jPoint, residual.jacobian_i, residual.jacobian_j);


I always assumed the NEMO folks had it like this because the signs were flipped inside the numerics...

I checked this inside numerics. Those signs come from differentiating the velocity difference; they do not include the residual assembly sign. In the simple momentum test the raw derivatives are -8/15 at i and +8/15 at j. Since the solver applies -F at i and +F at j, the blocks need to be (-Ji, -Jj, +Ji, +Jj).

Reuse ion_gy for 100 iterations and remove the duplicate configuration. Check the momentum flux derivatives independently and clarify the assembly test coverage. Validated the focused tests and both restored-bug controls.

bellonarts and others added 3 commits September 3, 2026 15:42

Remove trailing whitespace from ion cylinder march config

077f73e

Keep git diff --check clean before upstream PR submission.

bellonarts mentioned this pull request Sep 7, 2026

Make the NEMO total-energy under-relaxation check reachable #2884

Open

6 tasks

pcarruscag reviewed Sep 7, 2026

View reviewed changes

bellonarts added 2 commits September 7, 2026 15:41

Refine NEMO viscous Jacobian regression tests

cd8a39d

Reuse ion_gy for 100 iterations and remove the duplicate configuration. Check the momentum flux derivatives independently and clarify the assembly test coverage. Validated the focused tests and both restored-bug controls.

Refresh cone references and load Mutation++ in unit CI

50e3343

Navigation Menu

Uh oh!

Consistent NEMO viscous Jacobian: distinct left block and subtractive edge assembly - #2885

Consistent NEMO viscous Jacobian: distinct left block and subtractive edge assembly#2885
bellonarts wants to merge 5 commits into
su2code:developfrom
babybluechips:fix_nemo_viscous_jacobian

bellonarts commented Sep 5, 2026 •

edited

Loading

Uh oh!

pcarruscag Sep 7, 2026

Uh oh!

bellonarts Sep 7, 2026

Uh oh!

pcarruscag Sep 7, 2026

Uh oh!

bellonarts Sep 7, 2026

Uh oh!

pcarruscag Sep 7, 2026

Uh oh!

bellonarts Sep 7, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

bellonarts commented Sep 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed Changes

Tests

Related Work

PR Checklist

Uh oh!

pcarruscag Sep 7, 2026

Choose a reason for hiding this comment

Uh oh!

bellonarts Sep 7, 2026

Choose a reason for hiding this comment

Uh oh!

pcarruscag Sep 7, 2026

Choose a reason for hiding this comment

Uh oh!

bellonarts Sep 7, 2026

Choose a reason for hiding this comment

Uh oh!

pcarruscag Sep 7, 2026

Choose a reason for hiding this comment

Uh oh!

bellonarts Sep 7, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

bellonarts commented Sep 5, 2026 •

edited

Loading