FluidNumerics / FluidNumerics/SELF

Test acceptance criteria: entropy bounds are too weak; add manufactured solutions and convergence tests

Open
#177 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
Fortran
Stars
92
Forks
13
Avg merge
20h 38m
Merged PRs (30d)
7

Description

## Problem

Nearly every model-level integration test in the suite accepts or rejects on an entropy bound and nothing else. Counting the assertion conditions across `advection_diffusion_*`, `lineareuler*`, `linear_shallow_water_2d_*`, `ec_advection_*`, `esatmo*` and `burgers1d_*`:

| assertion | occurrences |
|---|---|
| `if(ef > e0)` — entropy did not grow | 33 |
| `if(ef /= ef)` — entropy is not NaN | 12 |
| a quantitative error against a known state (`maxerr > tolerance`) | 2 (`esatmo3d_motionless`, `esatmo3d_hydrostaticbalance`) |
| an order-of-accuracy check | 0 |

A monotonicity bound on a dissipative scheme is a *stability* statement. It says almost nothing about correctness, and it is satisfied by a large space of wrong answers:

- **The trivial solution.** Entropy of zeros is zero, and `0 > 0` is false.
- **A frozen solution.** If `dSdt` is identically zero the bound holds exactly.
- **Wrong wave speed or a mis-scaled flux coefficient.** Still dissipative, still bounded, arrives in the wrong place at the wrong time.
- **Order reduction.** A mortar projection or boundary treatment that silently drops to first order still dissipates.
- **An over-absorbing boundary or sponge.** More dissipation passes the bound more comfortably.

## This is not hypothetical

While fixing review feedback on #176 we found that `advection_diffusion_2d_rk3_pickup` and `advection_diffusion_2d_rk3_pickup_mpi` had been passing on **every GPU build while integrating an all-zero solution**.

`Read_DGModel{2,3}D_t` never published the restored solution to the device (`Read_DGModel1D_t` always has), and nothing in `ForwardStep` uploads it, so a GPU restart silently discarded the pickup file. The test then read the pickup, computed `e0` = entropy of zeros = 0, stepped zeros, got `ef` = 0, and `ef > e0` was false. Green, every time, asserting nothing. The bug was found only because a new test in #176 happened to compare restored values directly, and it failed on the MI210 runner with a mismatch of exactly `c0`.

That is a GPU restart path being completely broken for all 2D and 3D models, sitting behind a passing test suite. The acceptance criterion could not distinguish "restarts correctly" from "restarts with nothing at all".

## Proposal

Three tiers, cheapest first. The first is a small change that closes the specific hole above; the other two are the real fix.

### 1. Non-vacuity guards (cheap, immediate)

Any test that asserts a bound on a norm should also assert the solution is not trivial — e.g. `maxval(abs(solution%interior))` is within a stated factor of the expected amplitude, and `e0` itself is above a floor. This alone would have failed the GPU pickup tests on the first build after the regression. It is a few lines per test and needs no new infrastructure.

### 2. Method of manufactured solutions

The infrastructure already exists: `source2d`/`source3d` are overridable pure functions and `SourceMethod` is a type-bound procedure, so a test can subclass a model, install the MMS forcing, and assert a pointwise error against the manufactured field. `esatmo3d_motionless` is the closest existing precedent (steady exact state, absolute tolerance) and shows the shape of it.

Good first candidates, because their exact solutions are already written down in the docs and examples:

- **Linear Euler 2D/3D** — the plane-wave solution used in `examples/linear_euler2d_planewave_propagation.f90` is an exact solution of the discretised system's PDE. It is currently only an example, not an assertion.
- **Advection–diffusion** — the standard Gaussian/sinusoid with analytic decay.
- **Linear shallow water** — geostrophic or gravity-wave modes.

### 3. Convergence tests

For a spectral element method the strongest cheap statement is **p-convergence**: hold the mesh fixed, sweep `controlDegree`, and assert the error against an exact or manufactured solution falls at the expected rate — exponentially for smooth data, or at the design algebraic order where the scheme is limited. `controlDegree` is already a parameter in every test, so a sweep is a loop, not new machinery.

h-convergence is also available through the structured mesh generator (`nxPerTile`/`nTileX`), and matters specifically for the paths where an entropy bound is weakest:

- mortar/nonconforming interfaces (`lineareuler{2,3}d_mortar_soundwave`) — the place order reduction is most likely to hide;
- AMR prolongation/restriction (`lineareuler{2,3}d_amr_soundwave`);
- boundary conditions, where a wrong exterior state is dissipative and therefore invisible today.

Note there is precedent at the operator level: `twopointvectordivergence_*` and `mappedtwopointvectordivergence_*` assert exactness on polynomial fields. The gap is that nothing equivalent exists one level up, at the solver.

### Keep the entropy bounds

They are cheap, they run everywhere, and they are the right tool for what they actually test — stability, and entropy-conservation claims like `ec_advection_2d_entropy_conservation`. The proposal is to stop treating them as the *primary* acceptance criterion, not to remove them.

## Suggested sequencing

1. Add non-vacuity guards across the existing model tests (small, mechanical, immediately protective).
2. Add one MMS or exact-solution convergence test per model family, starting with linear Euler 2D since the plane-wave solution is already in the tree.
3. Extend to the mortar/AMR/boundary paths, where the current criterion is weakest.
4. Write the expectation down in `CONTRIBUTING.md` so new model tests start here rather than copying `if(ef > e0)`.

Worth deciding early: whether convergence tests run in the default `ctest` suite or behind a label, given the coverage job's runtime ceiling (documented in `docs/Learning/AdaptiveMeshRefinement.md` as having been hit before). A p-sweep at low degree on a small mesh is cheap; an h-sweep is not.

Context: #176.

Contributor guide

Open the contributing guide

Research direction

Start by reviewing the existing model-level integration tests named in the issue and the plane-wave solution in examples/linear_euler2d_planewave_propagation.f90, with esatmo3d_motionless as the closest assertion precedent. Confirm the initial scope and runtime policy, then define completion as non-vacuity protection plus a selected MMS or convergence test, while documenting the expectation in CONTRIBUTING.md.

Written by the indexing model from the issue text.

Assessment

Tech stack
fortran
Domain
testing
Issue type
Feature
Difficulty
5/5
Estimated time
Over a week
Activity status
Active
Clarity
Needs clarification
Newbie friendliness
35/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.