Miura Weno scheme - #1396
Draft
jcanton wants to merge 77 commits into
Draft
Miura Weno scheme#1396jcanton wants to merge 77 commits into
jcanton wants to merge 77 commits into
Conversation
Fix from Andreas Jocksch: the distance-vector allocation and the dimension passed on to compute_z_lsq_mat_c/compute_lsq_pseudoinv must be the actual stencil size, not the allocation size. Latent on main (both are 3 for the linear reconstruction) but required for 9-point stencils. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Row 8 listed cell 3 as an outer neighbor, but the outer neighbors of
direct neighbor 9 are {6, 12} per c2e2c_table (cell 3 is not adjacent to
any neighbor of cell 8): the entry must be 6.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Port the torus branch of the candidate least-squares setup from ICON (mo_intp_coeffs_lsq_bln.f90, icon-exclaim branch transport_ajocksch): 9-point stencil construction with orientation swap, torus moments and translated (hat) moments, the 27 quadratic and 3 linear candidate pseudoinverses with the live zero patterns and l_weights_s correction, plus scattering of the Fortran-ordered coefficients onto the C2E2C and C2E2C2E2C connectivity slots for the runtime stencils. Pure init-time numpy; unit tests run on the SimpleGrid tables and a synthetic periodic equilateral-triangle torus patch, including exact quadratic/linear reconstruction and an SVD conditioning gate (worst singular value ratio 2.3e-2 across all cells and candidates). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Two runtime gt4py stencils consuming the Task 1 WENO coefficients, one candidate slice at a time (a later task loops them 27x): - reconstruct_quadratic_coefficients_weno_candidate: per-candidate quadratic reconstruction on cells, reassembling the 9-point pseudoinverse dot products from the direct (C2E2C) and butterfly (C2E2C2E2C) row sets. - accumulate_weno_candidate_flux_weights: per-candidate WENO smoothness weighting on edges, gathering the upwind cell's coefficients via E2C and accumulating the weighted sums and weight sum. Faithful to mo_advection_hflux.f90 2447-2511. StencilTests with numpy references pass on embedded and gtfn_cpu on SimpleGrid. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Cell-domain gt4py stencil for the ihadv_tracer=102 reconstruction (mo_advection_hflux.f90 upwind_hflux_miura_weno, 983-1019): blends the 3 linear least-squares candidates -- candidate i drops direct neighbour i -- by the inverse-square smoothness weight of each candidate gradient. The zonal/meridional candidate pseudoinverses are slices of the Task-1 (n_cells, 3, 2, 3) array over the C2E2C rows. The conservative branch (llsq_lin_consv) is not ported (Fortran default off), so the constant coefficient stays p_cc. StencilTest with a literal-Fortran numpy reference passes on embedded and gtfn_cpu. Adds a numpy consistency unit test: on a linear field the 3 candidates coincide, so the WENO blend reduces to the plain LSQ gradient. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add the AdvectionWenoLinearState (6 zonal/meridional candidate pseudoinverses over C2E2C) and the SecondOrderMiuraWeno tracer flux, which mirrors SecondOrderMiura -- same cell/edge domains and reused flux stencil and limiter hook -- but swaps the SVD linear reconstruction for the WENO blend. Adds the HorizontalAdvectionType enum values LINEAR_2ND_ORDER_WENO (102) and QUADRATIC_3RD_ORDER_WENO (103); the factory builds SecondOrderMiuraWeno for 102 (raising a ValueError when the state is missing) and raises NotImplementedError for the not-yet-wired 103. Unit tests cover the from_fortran_dict int->enum mapping and the factory error paths. The existing LINEAR_2ND_ORDER path is unchanged. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
When the configured horizontal advection type is WENO (102/103), compute the init-time linear WENO pseudoinverses from the grid C2E2C table, cell centers and torus domain size, wrap them into an AdvectionWenoLinearState and pass it to the advection factory. Guarded to a torus grid on a single rank (the sphere coefficients are not ported and the coefficient halo exchange is not implemented); the non-WENO path is byte-identical. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
New analytical initial condition for a tracer-advection-only torus experiment: constant zonal wind u0 in a hydrostatically balanced base state (gauss3d-style) plus a constant circular qv disc. Since the dycore is off, the IC also prescribes the advection driving fields the dycore would normally provide: vn_traj = vn, mass_flx_me from edge-interpolated rho and ddqz_z_full_e, mass_flx_ic = 0 and airmass_now/new = rho * ddqz_z_full. The driver now allocates the tracer advection states before the IC and threads the target fields through initial_condition.create as a plain dataclass of fields (keeping icon4py.model.common free of advection imports); the other ICs ignore it. While at it, allocate mass_flx_ic with the KDim + 1 extension: the advection granule reads it at KDim + 1, and the previous full-level allocation made gtfn read out of bounds (garbage in the bottom tracer level). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Add the active tracer fields (qv, ...) to the driver's output pipeline: a tracer_state_to_dataarrays assembler using the common tracer CF metadata, an output_variables helper that extends the default variable list with the active tracer names, and the corresponding store call in the driver. Runs without tracers are unchanged (empty tracer state adds nothing and the variable list stays the default). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Integration test for the tracer-advection-only experiment on the torus_50000mx5000m_res500m grid (no serialized data): builds the ExperimentConfig programmatically with dycore and diffusion off, verifies the grid is fully periodic, runs the driver for 24 steps at CFL ~ 0.3 and checks that the initial output frame equals the sampled disc, that total tracer mass is conserved to 1e-12 and that the final qv matches the analytically translated disc. Observed relative L2 errors: 0.328 for linear WENO (ihadv=102) and 0.481 for plain miura (ihadv=2); the tolerances are frozen at ~1.5x these values, while a non-moving tracer produces sqrt(2) ~ 1.41. xfails on the embedded backend like the neighboring driver tests (concat_where in the metrics factory). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Reads the driver's icon4py_output_*.nc and the *_ugrid.nc companion and writes two figures into the output directory: tripcolor panels of the initial, final, analytically translated and difference qv fields (periodic seam triangles masked), and time series of the relative L2 error versus the translated disc and of the total tracer mass drift. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Port of prep_gauss_quadrature_q_miura3 (mo_advection_quadrature.f90
561-703): 4-point Gauss-Legendre quadrature of the monomials
{1, x, y, x^2, y^2, xy} over the parallelogram departure region. The
constant Jacobian (a quarter of the parallelogram area) is folded into
the weights z_wgt = 0.25 * wgt, so the sums are the area AVERAGES of
the monomials and, unlike the plain/cubic variants, there is no
Jacobian determinant and no dreg_area output. Argument conventions
mirror prepare_numerical_quadrature_for_cubic_reconstruction; the test
uses the live gaussq_2d_o2 shape functions and weights against a
literal-port numpy reference.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Port of the final flux write of upwind_hflux_miura3_weno (mo_advection_hflux.f90 2514-2521): the accumulated candidate-weighted coefficients are divided by the accumulated smoothness weights (before the dot product, as in the Fortran), dotted with the miura3 quadrature vector and multiplied by the mass flux. No division by the departure-region area: the miura3 quadrature already averages. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Port of calculate_planar_distance_at_edge (mo_intp_coeffs.f90 2278-2406): torus-periodic separation vectors from the edge midpoint to the E2C cell circumcenters and E2V vertices, rotated into the edge-local (primal normal, dual normal) frame - ICON's pos_on_tplane_e components 1:4, including the edge-vertex components 3:4 that compute_ffsl_backtrajectory needs and that no icon4py code produced so far. Unlike interpolation_fields.compute_pos_on_tplane_e_x_y_torus this is the full projection without the equilateral shortcut. The remaining backtrajectory inputs need no new code: primal/dual_normal_cell equal the per-edge normals on both E2C slots on the plane torus (cvec2gvec identity) and lvn_sys_pos is computed at runtime by the existing counterclockwise-indicator stencil. The synthetic torus patch moves from the weno unit tests to the shared tests utils and grows edge tables (E2V, E2C, edge centers, primal/dual normals, tangent orientation) with ICON's conventions: the primal normal points from cell 1 to cell 2, the stored tangent is tangent_orientation * normalize(v2 - v1). Unit tests check the derived closed forms on the equilateral patch: cell offsets at -/+ the apothem a/(2*sqrt(3)) with zero tangential part (matching the existing equilateral shortcut), vertex offsets purely tangential at +-a/2 and ordered by the tangent orientation. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Closed-form test of the never-wired backtrajectory stencil on the synthetic equilateral torus patch with constant vn = c > 0, vt = 0 and the new init-time torus geometry: the departure-region vertices must equal the edge vertices and their images translated upstream by c*dt along the edge normal, expressed relative to the upwind (cell 1) circumcenter in the global frame, ordered counterclockwise per lvn_sys_pos (btraj_dreg, mo_advection_traj.f90 363-773). The expected values are derived from the unwrapped patch triangles only - independent of the edge-local-frame arrays under test - and the output loop is additionally checked to be counterclockwise with shoelace area edge_length * c * dt. Sign mutations (flipped primal normal, inverted lvn_sys_pos) make the test fail, so it pins down the frame conventions of both the geometry and the stencil. Runs on embedded and gtfn_cpu. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Literal numpy port of the upwind_hflux_miura3_weno cell loop (f90 2443-2546) compared against the composed gt4py stencils (reconstruct x27, accumulate x27, flux) on the torus patch with real candidate pseudoinverses. Gated three ways: hand-derived constant-reconstruction closed form on the reference, 1e-12 reference-vs-pipeline agreement with live and synthetic candidate weights, and a polynomial-only two-group closed form on an exactly quadratic field. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Runtime assembly of upwind_hflux_miura3_weno (mo_advection_hflux.f90 2033-2620): per-step lvn_sys_pos, ffsl backtrajectory over the full time step, quadratic Gauss quadrature, a Python loop over the 27 candidate reconstruction/accumulate launches with L_WEIGHTS_S, and the WENO flux. AdvectionWenoQuadraticState carries the scattered candidate pseudoinverses, moments, cell area and the torus ffsl geometry. p_vt is threaded from SemiLagrangian through the SemiLagrangianTracerFlux interface; the second order schemes ignore it. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
_construct_weno_quadratic_state assembles the 9-point stencil, torus moments, 27 candidate pseudoinverses and the ffsl backtrajectory geometry from the grid geometry sources; the linear state is now built only for 102 and the quadratic one only for 103 (previously 103 built the linear state before the factory raised). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Observed relative L2 error 0.297 (vs 0.328 for 102 and 0.481 for 2); tolerance frozen at ~1.5x observed. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
- config: honor run_nml["ldynamics"]; with dynamics off, leave nonhydrostatic/diffusion None so tracer advection runs from the IC-prescribed mass fluxes instead of silently running the dycore over the frozen prescription. - tracer blob test: use the squared disc formula (dx**2+dy**2 <= r**2) to match the tracer_blob IC bit-for-bit, so the frame-0 array-equal check can't flake on a boundary cell. - interpolation_factory: TODO guarding the lsq_high_ord -> lsq_dim_stencil mapping (WENO namelists set lsq_high_ord=2 while this assumes 3). - advection tests: add level markers (integration on the miura3 pipeline, unit on the ffsl uniform-flow stencil test). - driver_io: rename create_io_monitor's local variable so it no longer shadows the module-level output_variables function. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Upstream fixed the torus periodicity itself, more broadly than the fork did: rectangular independently wrapped periods are now the default rather than an opt-in mode, the same correction reaches StretchedTorusGridSpec, the coupled lattice stays reachable as periodic_layout='skew' including for odd row counts, and the layout is recorded in the grid identity. Drop the git pin and take 0.8.0 from the index. The option is renamed, so 'icon_periodicity=True' becomes periodic_layout='rectangular'. It is now the default but is passed explicitly: 'skew' is exactly the convention that does not survive the per-axis minimum image icon4py applies on a torus. Grids are unchanged: the convergence test reports the same per level errors as the fork did, bit for bit, so the release reproduces the same mesh with the same numbering. Regenerated grids also keep every coordinate inside one period, no edge misreconstructed under a per-axis minimum image at any of the four refinement levels, and a discrete divergence of 5e-18 to 5e-17. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Ar8TP7p4QVPryDErtCigzg
Main renamed two packages under the WENO work: `model/atmosphere/advection` -> `tracer_advection` (#1392) and `model/standalone_driver` -> `model/driver` (#1283). Git moved most of the branch's files itself; the new test files and `standalone_driver/config.py` were relocated by hand. Beyond the renames: - `#1379` fixed the `lsq_dim_stencil` / `lsq_dim_c` mix-up upstream, which supersedes the branch's own fix and the TODO in `interpolation_factory`. `compute_lsq_weights_c` also lost its now-unused stencil-size argument. - `#1391`/`#1361` turned the driver config declarative: `tracer_blob` is now registered in `IC_CONFIG` and in the config-union map, and the WENO driver test builds its `DriverConfig` through `make_initial`. - `#1404` removed `PrognosticState.tracer`, so `tracer_blob` takes the tracer state explicitly and `_store_output` is passed `tracers.current`. - `#1404` also made the tracer prep state share the dycore's accumulated buffers. `assemble_driver_states` keeps that path and only accepts pre-allocated states, which the driver now builds solely for the prescribing `tracer_blob` initial condition. 🤖 Written by an agent on behalf of @jcanton
Brings the idealized linear-advection convergence study (the "cylinder", i.e. `circle_2d`, and `gaussian_2d`) onto the WENO branch, as Andreas asked. The two branches had grown separate machinery for the same thing: an initial condition prescribing the advection driving fields the dycore would otherwise accumulate. The convergence branch's is kept, since it is the more general of the two, and the WENO one is dropped: - `AdvectionPrepAdvState` lives in `common/states/adv_states.py`, so an initial condition in `common` can fill it without importing the advection package. - `link_prep_adv_to_dycore` builds the dycore's `PrepAdvection` from those buffers, replacing `initialize_prep_tracer_advection`, which shared them the other way round. Same identity invariant, so its unit tests are rewritten rather than dropped. - `tracer_blob` takes `adv_prep_adv_state` like the linear-advection initial conditions, and `TracerAdvectionPrescription` is gone. Its airmass fields were already redundant: the driver recomputes the airmass from rho every step, including with the dycore off. Also here: relative times serialize as floats now (the convergence study needs a non-integer dtime), so the two config fixtures that pinned integer seconds are updated; `linear_advection_tests.py` moves to `model/driver/tests/driver/` after the `#1283` rename; and the two `*2d.yaml` experiment configs that no test references are dropped in favour of the `*_2d.yaml` ones. 🤖 Written by an agent on behalf of @jcanton
🤖 Written by an agent on behalf of @jcanton
Picks up the work that landed on that branch after the local checkout: the plotting helpers are dropped (`plot_utils`, `enable_plot`, the matplotlib dependency), the convergence grid family and CFL are retuned, the `circle_2d` bands are honest about a discontinuity (degraded first order in L1, zeroth in L-infinity), and `link_prep_adv_to_dycore` gets its own unit tests, which replace the ones written here for the same invariant. Two fixes on top: the new test annotated a fixture with a `TYPE_CHECKING`-only import, and called `data_allocation.get_allocator`, which does not exist — allocators come from `model_backends.get_allocator`. 🤖 Written by an agent on behalf of @jcanton
Zalesak flux-corrected transport, from `hflx_limiter_mo` in `mo_advection_hlimit.f90`.
All six stencils it needs were already in the package, imported by nothing; this wires
them into a `Monotonic` limiter beside `PositiveDefinite`.
`beta_fct` is per-scheme, not per-config: Fortran passes `opt_beta_fct` only from the
quadratic-reconstruction schemes (`mo_advection_hflux.f90:3083` and `:4810`), and the
linear ones fall through to the routine's own default of 1, which is strictly monotonic.
The namelist value arrives as `AdvectionConfig.monotonic_limiter_boost_factor`, range
checked to `[1, 2)` as `mo_advection_nml.f90:336` does.
The boundary-interpolation repair stage (`f90:361-378`) is built only on a limited-area
grid, mirroring the Fortran guard; the zero-initialisation of `r_p`/`r_m` stands in for
`f90:390-397`, which zeroes exactly the rows the factor stencil never writes.
The limiter needs the mass flux and the new density, so `apply_flux_limiter` and
`compute_tracer_flux` grow those arguments and `FiniteVolume.run` forwards them.
Measured on the tracer disc, relative to the initial tracer range, over 24 steps:
scheme limiter overshoot undershoot
miura none 2.7e-01 3.6e-01
miura WENO none 7.0e-02 1.0e-02
miura3 WENO none 1.3e-01 8.2e-02
miura monotonic 2.1e-14 3.0e-17
miura3 WENO monotonic 4.2e-02 7.5e-18
The unlimited rows assert that they *do* leave the range, so the monotonic rows cannot
quietly stop proving anything.
🤖 Written by an agent on behalf of @jcanton
`upwind_hflux_miura3` for `lsq_high_ord=2`, live path only. It is
`ThirdOrderMiuraWeno` with the 27-candidate loop collapsed to the single
full-stencil reconstruction, so it shares the backtrajectory, the quadrature and
the reconstruction stencil and differs only in the last step.
Almost everything already existed. What is new:
- `compute_lsq_pseudoinverse_quadratic`, the full-stencil quadratic pseudoinverse.
It is the block the WENO path already computed for its candidates 0-2, lifted
out and shared, so the two cannot drift apart. `compute_lsq_coeffs` in
`interpolation_fields` is deliberately untouched: it describes ICON's
`lsq_lin_set`, and the torus high-order path lives in `weno_least_squares`.
- `compute_horizontal_tracer_flux_from_quadratic_coefficients`, the one genuinely
missing stencil (`mo_advection_hflux.f90:4764-4770`): the upwind cell's six
coefficients dotted with the quadrature vector, times the mass flux.
- `AdvectionQuadraticState`, and a shared `_quadratic_reconstruction_inputs` in the
driver so the two quadratic states are built from one code path.
`reconstruct_quadratic_coefficients_weno_candidate` is renamed to
`reconstruct_quadratic_coefficients_svd`: both schemes use it, and only the WENO
one passes a candidate.
The cubic variant (`lsq_high_ord=3`) is not implemented — it needs 10 unknowns,
the third-order moments and its own quadrature.
Measured on the tracer disc, relative L2 error against the translated disc after
24 steps, which orders the schemes the way the theory says it should:
miura 0.370
miura3 0.326
miura3 WENO 0.299
The decisive check is a unit test rather than that number: applied to the exact
cell averages of a random quadratic, the pseudoinverse recovers that quadratic's
derivative coefficients to 1e-10, which is the property the whole scheme rests on.
🤖 Written by an agent on behalf of @jcanton
`upwind_hflux_miura_cycl`, mo_advection_hflux.f90 2219-2458. No new stencils: the
five it needs were already in the package, imported by nothing.
Unlike the other Miura variants this is a `FiniteVolume` rather than a
`SemiLagrangianTracerFlux`, because each substep advances its own tracer and
density before the next reconstruction, so the flux and the update interleave.
The mass flux is constant over the step, so the tangential velocity, the
backtrajectory (built with `0.5 * dtime / n`) and the mass-flux divergence are
computed once; the divergence reuse is what the stencil's `nsub` argument gates.
Two fidelity points worth stating, both from the Fortran:
- the substep limiter is always the positive definite one, whatever the
configured limiter is: `f90:2338` calls `hflx_limiter_pd` for both `ifluxl_sm`
and `ifluxl_m`. The configured limiter only decides whether to limit at all.
- the step flux is the plain mean of the substep fluxes, each evaluated with the
full mass flux, and the caller then applies it over the full step, which is
what makes it equal to the sum of substep flux times substep length.
`n_substeps` comes from the namelist's `nadv_substeps` (default 3) and is
restricted to 2 or 3, the two averaging stencils ICON provides.
Measured on the tracer disc, relative L2 against the translated disc, unlimited:
miura 0.483
miura, 2 substeps 0.438
miura, 3 substeps 0.433
Subcycling improves on plain miura and improves further with more substeps, which
is the expected ordering. Mass is conserved to 1e-12 over the 24 steps.
Not yet validated against reference data — the evidence is the Fortran read
line by line, mass conservation, and that ordering.
🤖 Written by an agent on behalf of @jcanton
… two-stage form
This is the half of Andreas Jocksch's `upwind_hflux_miura_cell` reformulation that
carries over to GT4Py: fusing the least-squares reconstruction into the flux kernel
so the three coefficient fields are never written to memory.
The other half, scattering from the cell to its three edges, does not carry over and
would not buy anything if it did. GT4Py has no scatter, and the E2C gather we already
use performs one reconstruction per cell and one read per edge — the same arithmetic
and the same traffic as his cell loop, without the scatter, and it stays correct under
MPI, whereas his cell range is interior-only and drops the halo exchange.
Measured on `icon_benchmark_regional` (MCH_OPR_R19B08_DOMAIN01, 80 levels, gtfn_cpu):
Name (time in ms) Min Max Mean StdDev
fused 398.6145 426.3306 413.4167 12.2363
two_stage 774.9087 1215.0837 979.8793 170.0719
Roughly 2x, and the distributions do not overlap — the slowest fused round is faster
than the fastest two-stage one. The fused form is also bit-identical to the two-stage
one, which the accompanying test asserts with assert_array_equal rather than a
tolerance, since this is meant to be a pure schedule change.
Nothing uses the fused kernel yet; this commit is the experiment the decision rests on.
🤖 Written by an agent on behalf of @jcanton
The benchmark said to: about 2x faster on `icon_benchmark_regional` with gtfn_cpu, and bit-identical, so this is a pure schedule change. Evidence that it really is: the seven tracer-disc cases produce L2 errors identical to the last digit before and after the swap, and the stencil-level test asserts `assert_array_equal` rather than a tolerance. The subcycled and WENO schemes still use their two-stage forms. The WENO one is where the same fusion should be worth most — 54 launches per tracer per step, each writing six cell fields and reading them straight back — but that is a bigger change than this one and wants its own measurement. 🤖 Written by an agent on behalf of @jcanton
The building block FFSL needs. Not FFSL itself. Correcting an earlier assumption while I was at it: FFSL is not *required* to be cubic — `upwind_hflux_ffsl` and `hflux_ffsl_hybrid` both branch on `lsq_high_ord` 1, 2 or 3, and the `== 2` branch calls the same `recon_lsq_cell_q_svd` that MIURA3 already gave us. What forces cubic here is the icon4py side: the FFSL quadrature and flux stencils that already exist were written for ICON's default of 3, so they take ten coefficients and ten quadrature components. Choosing 2 instead would mean writing six-coefficient variants of three stencils. What this adds: - `compute_lsq_moments_torus(cubic=True)`, extending the moment set from `[x, y, x^2, y^2, xy]` to `+ [x^3, y^3, x^2 y, x y^2]`. ICON's ordering, so the cubic set appends to the quadratic one rather than reordering it — a test pins exactly that, since a silent reordering would be very hard to find later. - the matching cubic shifts in `compute_lsq_moments_hat`. - `compute_lsq_pseudoinverse_cubic`. ICON gives this set `wgt_exp = 0`, so unlike the quadratic one the rows are unweighted, and with 9 unknowns over the same 9-point stencil the system is square rather than overdetermined. The `fxxx`/`fyyy` integrands follow the `delx`/`dely` form the Fortran actually uses, not the algebraically equal MAPLE form quoted beside it, which it says is "marginally more accurate" but does not use. Three tests, the last being the one that matters: the analytic polygon integrals agree with a quadrature exact for cubics (the existing 3-point midpoint rule is only exact to degree 2, so this needed a 4-point rule), and the pseudoinverse applied to the exact cell averages of a random cubic recovers that cubic's nine derivative coefficients. Still missing for FFSL: the patch1/patch2 cell gather. `compute_ffsl_flux_area_list` emits those as absolute cell indices, and GT4Py cannot index a field by a runtime index, so they need to become a relative slot into an E2C2E2C offset. 🤖 Written by an agent on behalf of @jcanton
All three affordable schemes land inside the bands miura already used, so they share
them rather than each getting a hand-tuned range. Those bands are principled rather
than fitted -- a discontinuity converges at a degraded first order in L1 and not at
all in the maximum norm -- which is why one set covers every scheme:
circle_2d, coarse family, spacings 200/100/50
scheme L1 stderr Linf stderr
miura_subcycled 0.7026 0.0084 -0.0527 0.0575
miura3 0.7234 0.0084 -0.0595 0.0536
miura WENO 0.7305 0.0091 -0.0745 0.0296
asserted L1 in [0.0, 1.0] Linf in [-0.1, 0.1]
The miura3 WENO rows are kept but skipped, with the cost recorded in the skip reason.
Its 27-candidate loop is two stencil launches per candidate per cell, so a row of this
study costs about two orders of magnitude more than the single-launch schemes: the
other three finish in minutes each, while that one was killed after 20 hours without
producing a rate. It stays covered end to end by the tracer-disc test.
gaussian_2d rows still assert nothing pending their measurement.
🤖 Written by an agent on behalf of @jcanton
FFSL splits each edge's departure region into three patches. Patch 0 sits in one of the edge's two E2C cells, but patches 1 and 2 sit in cells ICON identifies by ABSOLUTE index (`butterfly_idx`), and GT4Py cannot index a field by a runtime index. They have to become a relative slot into a connectivity instead. This adds that connectivity; the FFSL scheme itself is still to come. The slot order is not ours to invent, and it is not an artefact of loop order either. `mo_model_domimp_setup.f90:380-443` assigns each of the four cells by geometry: which of the edge's two cells it flanks, and which of the edge's two vertices it shares. So `slot = 2 * side + vertex`, and the construction is reproducible from `e2c`, `e2v`, `c2e2c` and `c2v` alone. That the order is geometric is also what makes it testable, which matters here more than usual: a permuted slot order still gathers real cells and still produces plausible fluxes, so it would not announce itself. The test asserts the defining property directly -- the cell in each slot must both neighbour `e2c[side]` and contain `e2v[vertex]` -- and a further test permutes the slots to confirm that check actually fails when it should. The SimpleGrid table is generated from the construction rather than hand-written, and pinned by a test that the literal still equals the algorithm. Its hand-written `c2e2c2e2c` neighbour does NOT match `grid_manager`'s ordering, which has cost debugging time before; this one cannot drift the same way. Proven usable from GT4Py, not merely declared: a runtime-slot gather written as a chain of `where()` over the four compile-time slots reproduces the numpy reference exactly, on both embedded and gtfn_cpu. 🤖 Written by an agent on behalf of @jcanton
…lute index
With the E2C2E2C offset in place, `compute_ffsl_flux_area_list` can hand the two
outer patches a relative slot, which GT4Py can gather with, rather than the absolute
cell index it inherited from the Fortran, which GT4Py cannot.
The substitution collapses much further than expected. ICON picks the cell with
MERGE(butterfly_idx(je,jb,1,p), butterfly_idx(je,jb,2,p), lvn_pos). Since the slot is
2 * side + vertex, and patch 1 is the vertex-0 wing while patch 2 is the vertex-1
wing, that MERGE is only a choice of side -- so the lookup becomes two constants, the
eight butterfly index/block input fields disappear entirely, and the program drops
from 38 arguments to 28.
The test pins the reduction rather than the arithmetic, because that is where a port
like this goes wrong silently: a flipped sign convention or a swapped patch-to-vertex
mapping would gather a real cell and produce plausible fluxes.
The empty-patch sentinel is inlined as a literal in the field operator, because gtfn
does not fold a module-level constant referenced inside one ("Symbols not found") --
the same limitation `_WENO_EPS` already carries in
accumulate_weno_candidate_flux_weights. The name stays for callers and tests, so the
value is still defined once.
🤖 Written by an agent on behalf of @jcanton
… 0 side Three follow-ups an adversarial review of the previous two commits turned up. `test_local_connectivity` is parametrized over a HARDCODED list of field offsets rather than over the registered ones, so a new connectivity is silently not localization-tested. E2C2E2C joins the list; all four ranks pass, which also confirms the derived table is built from the LOCAL tables and its slots hold local indices. `test_icon.py` needed nothing, but that was worth checking rather than assuming: its skip-value tests walk `dims.local_dims()`, so they picked the new dimension up by themselves, and they pass on the real MCH and global grids -- which is what justifies listing E2C2E2CDim in CONNECTIVITIES_ON_BOUNDARIES. The slot test gains the `vn == 0` case. Fortran's `lvn_pos` is `p_vn >= 0` (mo_advection_geometry.f90:768), not `> 0`, so zero normal velocity takes side 0; nothing was pinning which way that went. 🤖 Written by an agent on behalf of @jcanton
Completes the convergence study for every scheme the cost allows.
gaussian_2d, coarse family, spacings 200/100/50
scheme L1 stderr Linf stderr
miura_subcycled 2.2149 0.0050 2.5507 0.0727
miura3 3.0391 0.0234 2.9975 0.0383
miura WENO 2.3809 0.0758 1.3071 0.0463
miura3 earns its formal order and gets a principled band: third order in both
norms, to within the same tolerance miura's second order uses. That is the
strongest evidence so far that the quadratic reconstruction is right, since a
scheme can produce plausible error magnitudes while converging at the wrong rate.
The other two get bands centred on the measurement, at the width the discontinuous
rows already use, and `_measured()` says so in its name and docstring. Neither has
a clean formal rate on a smooth profile -- WENO blending and subcycling both change
the effective order -- and dressing a fitted band up as a derived one would be the
worse outcome.
The whole sweep is now 3 rows in 37 minutes per profile, against the 20 hours the
single skipped row was still not finished in.
🤖 Written by an agent on behalf of @jcanton
No new reference data was needed: the MCH serialization already contains it. The ntracer legend at the top of this file records that tracer 2 runs miura with itype_hlimit = 3, and the savepoints for it have been on disk all along -- only tracers 1 and 4 had parametrize rows. So the monotonic limiter is no longer resting on a careful Fortran read plus property tests. It reproduces ICON's own output, at the tolerances this harness already used: rtol 1e-10 / atol 1e-11 on hfl_tracer, which is exactly the field the limiter rewrites, and atol 1e-16 on the new tracer. Tracer 0 is the other untested one, ihadv_tracer = 52 (FFSL_HYB_MCYCL). Its reference data is equally available, and it is what a future FFSL port should be held to. 🤖 Written by an agent on behalf of @jcanton
* origin/main: Use multiple output domains for compute_rho_theta_pgrad_and_update_vn (#1405)
… plots The branch had moved on with four fixes: the standalone -> driver move, a wrong import, and two float-time fixes in test_config_io and test_driver_states. The test_driver_states conflict was the same bug fixed two ways. Both branches hit the fixture annotation being evaluated at definition time; here it was made a runtime import, there `from __future__ import annotations` was added so annotations stay lazy. Theirs is kept, since it is the file's own idiom. The 44 convergence plots at the repo root are untracked and deleted, and *.pdf is ignored. They were generated output from a plotting helper that origin/advection_convergence has since removed, so nothing produces them any more; the ignore rule is there so a local run of an older revision cannot commit them again. Nothing outside those files referenced a pdf. 🤖 Written by an agent on behalf of @jcanton
|
When developing, you can test your changes on CSCS CI before merge with the You can pass options to override pipeline variables, for example:
Available options are:
For each option, See The Merging Once your PR is approved and ready for merging, add it to the merge queue. The Optional Tests To run benchmarks you can use:
For more detailed information please look at CI in the EXCLAIM universe. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
(Claude) port of Andreas' Weno advection scheme