Miura Weno scheme - #1396
Draft
jcanton wants to merge 22 commits into
Draft
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>
|
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
1 participant
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