cartesian_uniform_2d

Kind
grid
Family
cartesian finite_difference
Operators
D D2 Dxy grad_norm varcoeff_laplacian
Boundary conditions
dirichlet neumann periodic zero_flux zero_gradient
Source
grids/cartesian_uniform_2d/grid.esm
Tags
esd:grid family:cartesian ndim:2 spacing:uniform axis:x axis:y esd:stencil family:finite_difference grid:cartesian_uniform_2d op:D order:2 stencil_width:3 op:D2 order:4 stencil_width:6 stencil_width:5 order:6 stencil_width:8 stencil_width:7 op:Dxy op:grad_norm axes:x,y order:1 op:varcoeff_laplacian esd:rule bc:dirichlet spec:9.6.8 bc:periodic bc:zero_gradient bc:neumann bc:zero_flux

Uniform 2-D Cartesian grid with consumer-supplied extent: NX by NY cell-centered points, cell centers (x_i, y_j) = (x0 + (i - 1/2) dx, y0 + (j - 1/2) dy). The x and y index sets are the iteration axes; NX and NY (the cell counts) remain the esm-spec 9.7 metaparameters, bound at the import edge or the loader API. The real-valued geometry is NOT baked into this file: x0/dx (x-axis) and y0/dy (y-axis) are free names in the coordinate template bodies (and in this grid’s stencils/rules), resolving in the consuming component’s scope at evaluation — the same keyed-factor contract grids/cartesian_uniform_1d uses for x0/dx and grids/mpas uses for areaCell/dvEdge. A consuming model MUST define all four as ordinary real-valued model variables: for a box [ax, bx] x [ay, by], x0 = ax, dx = (bx - ax)/NX, y0 = ay, dy = (by - ay)/NY. Spell each spacing as an observed dividing by the METAPARAMETER NAME (dx = {op: /, args: [bx - ax, NX]}), which esm-spec 9.7.6 substitutes as an integer literal at load, so a loader-API rebinding of NX/NY — a convergence sweep — keeps the spacing consistent automatically; a model that binds NX/NY at the import edge closes them before they reach the model’s own scope and must spell the matching literals instead. BUILD-TIME SCOPE CAVEAT (AGENTS.md): the x_coord/y_coord free-name templates resolve in RUNTIME positions (rule bodies, equation RHS) but NOT in ic equations or §6.6.5 test references, which are evaluated build-time and cellwise with the consumer’s parameters out of scope — there, coordinates must fold from literals plus the metaparameter names (x_i = ax + (i - 1/2)((bx - ax)/NX)).

Metaparameters

NameTypeDefaultDescription
NXinteger64Number of cells along x.
NYinteger64Number of cells along y.

All sizes are load-time metaparameters (esm-spec §9.7.6): rebind them at the import edge (or through the loader API) and the same files serve every resolution.

Consumer-supplied free names

dx dy gB gL gR gT qXL qXR qYL qYR x0 y0

The geometry is consumer-supplied so the same files serve any domain extent, not just any resolution. For a domain [a, b] the consuming model defines x0 = a (a parameter) and dx = (b − a)/N (an observed whose expression divides by the metaparameter name N, so a loader-API rebinding of N — a convergence sweep — keeps dx consistent). A model that instead closes N at the import edge spells the matching literal ({op: /, args: [1, 8]} for N = 8 on the unit interval). See the extent proof problems/heat_1d_zero_grad_nonunit.esm (x ∈ [−1.5, 2.5], observed order 2.00).

Index sets

NameKindSize
xintervalNX
yintervalNY

Geometry templates

TemplateDefinitionDescription
x_coordx_coord[i] = x0 + (i − 0.5)·dx for i ∈ xCell-center x-coordinates x_i = x0 + (i - 1/2) dx, shaped over the x index set. x0 and dx are free names supplied by the consuming model (see metadata.description for the contract).
y_coordy_coord[j] = y0 + (j − 0.5)·dy for j ∈ yCell-center y-coordinates y_j = y0 + (j - 1/2) dy, shaped over the y index set. y0 and dy are free names supplied by the consuming model (see metadata.description for the contract).

Stencils

Interior-only, match-less templates importing grid.esm; the rules below wrap them with boundary-condition face regions into complete rewrite rules.

central_D1x_interior

op:D order:2 stencil width 3 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D1x_interior.esm

central_D1x_interior(f)[i,j] = (f[i+1, j] − f[i−1, j]) / (2·dx) for i ∈ [2, NX−1], j ∈ y

Interior-only 3-point centered first derivative along x on cartesian_uniform_2d: (f[i+1,j] - f[i-1,j]) / (2 dx) over i in [2, NX-1], all j. Match-less named template; the boundary-condition layer (rules/central_D1x_periodic.esm, rules/central_D1x_zero_grad_bc.esm) wraps it with the x-edge faces. dx is a free name resolving to the consuming model’s x-spacing variable (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). O(dx^2). The per-axis first derivative D(f,x) is lowered by its own rule; a consumer needing D(f,y) imports central_D1y_interior’s rule instead.

References

central_D1y_interior

op:D order:2 stencil width 3 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D1y_interior.esm

central_D1y_interior(f)[i,j] = (f[i, j+1] − f[i, j−1]) / (2·dy) for i ∈ x, j ∈ [2, NY−1]

Interior-only 3-point centered first derivative along y on cartesian_uniform_2d: (f[i,j+1] - f[i,j-1]) / (2 dy) over j in [2, NY-1], all i. Match-less named template; the boundary-condition layer (rules/central_D1y_periodic.esm, rules/central_D1y_zero_grad_bc.esm) wraps it with the y-edge faces. dy is a free name resolving to the consuming model’s y-spacing variable (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). O(dy^2). The per-axis first derivative D(f,y) is lowered by its own rule; a consumer needing D(f,x) imports central_D1x_interior’s rule instead.

References

central_D2x_4_dirichlet_face1

op:D2 order:4 stencil width 6 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_4_dirichlet_face1.esm

central_D2x_4_dirichlet_face1(g, f1, f2, f3, f4, f5) = (2560·g + −4095·f1 + 1848·f2 + −378·f3 + 72·f4 + −7·f5)/(756·dx·dx)

Fourth-order one-sided second-derivative closure along x for the cell nearest a Dirichlet (fixed-value) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_4_dirichlet_bc) applies it at the first interior x-column, tiled over all j. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_4_dirichlet_face1 (same Fornberg 1988 offset-node weights), rewritten for the x-axis so the free name is the x-spacing dx. On this cell-centered grid the x-wall sits half a cell outside the first center, so the wall value g and the five nearest cells along x sit at offsets {0, 1/2, 3/2, 5/2, 7/2, 9/2}*dx from the wall. The unique degree-5 polynomial through those six data, differentiated twice at the first center (offset 1/2 dx), gives D2 = (2560 g - 4095 f1 + 1848 f2 - 378 f3 + 72 f4 - 7 f5) / (756 dx^2), where f1..f5 are the five cells (at fixed j) in order of increasing distance from the x-wall and g is the wall value. Exact for polynomials up to degree 5, so the local truncation error is O(dx^4) — matching the interior 5-point stencil central_D2x_4_interior. The six weights sum to zero, so a steady linear-in-x profile carrying the wall value is discretely exact. Reused unchanged at the right x-wall by reflection (the second derivative is even under x -> -x), with f1..f5 mapped to the cells in increasing distance from the right wall and g bound to the right wall value. dx is a free name resolving to the consuming model’s x-spacing variable (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_4_dirichlet_face2

op:D2 order:4 stencil width 6 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_4_dirichlet_face2.esm

central_D2x_4_dirichlet_face2(g, f1, f2, f3, f4, f5) = (−256·g + 1323·f1 + −2100·f2 + 1134·f3 + −108·f4 + 7·f5)/(756·dx·dx)

Fourth-order one-sided second-derivative closure along x for the SECOND cell from a Dirichlet (fixed-value) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_4_dirichlet_bc) applies it at the second interior x-column, tiled over all j, whose centered 5-point stencil would reach one cell outside the domain. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_4_dirichlet_face2, rewritten for the x-axis so the free name is dx. The same six data as central_D2x_4_dirichlet_face1 (the wall value g at offset 0 and the five nearest cells f1..f5 at offsets {1/2,3/2,5/2,7/2,9/2}*dx) fit a degree-5 polynomial; differentiating it twice at the SECOND center (offset 3/2 dx) gives D2 = (-256 g + 1323 f1 - 2100 f2 + 1134 f3 - 108 f4 + 7 f5) / (756 dx^2). Exact for polynomials up to degree 5, so the local truncation error is O(dx^4); the six weights sum to zero, so a steady linear-in-x profile is discretely exact. Reused unchanged at the right x-wall’s second cell by reflection, with f1..f5 mapped to the cells in increasing distance from the right wall and g bound to the right wall value. dx is a free name resolving to the consuming model’s x-spacing variable (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_4_interior

op:D2 order:4 stencil width 5 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_4_interior.esm

central_D2x_4_interior(f)[i,j] = (−f[i−2, j] + 16·f[i−1, j] − 30·f[i, j] + 16·f[i+1, j] − f[i+2, j]) / (12·dx·dx) for i ∈ [3, NX−2], j ∈ y

Interior-only 5-point fourth-order centered second derivative along x on cartesian_uniform_2d: (-f[i-2,j] + 16 f[i-1,j] - 30 f[i,j] + 16 f[i+1,j] - f[i+2,j]) / (12 dx^2) over i in [3, NX-2], all j. The per-axis lift of the 1-D operator grids/cartesian_uniform_1d/stencils/central_D2_4_interior with the y axis a passive free index (the same coefficients -1,16,-30,16,-1 / 12, applied along x at fixed j). Match-less named template; the boundary-condition layer (rules/central_D2x_4_dirichlet_bc.esm and siblings) wraps it with the x-edge fourth-order face closures. dx is a free name resolving to the consuming model’s x-spacing variable (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). O(dx^4). The 2-D fourth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y), each half lowered by its own per-axis fourth-order rule.

References

central_D2x_4_neumann_face1

op:D2 order:4 stencil width 6 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_4_neumann_face1.esm

central_D2x_4_neumann_face1(q, f1, f2, f3, f4, f5) = (−1.92×10⁴·q·dx + −1.5285×10⁴·f1 + 7544·f2 + 1.2546×10⁴·f3 + −5784·f4 + 979·f5)/(2.0268×10⁴·dx·dx)

Fourth-order one-sided second-derivative closure along x for the cell nearest a Neumann (fixed-flux) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_4_neumann_bc) applies it at the first interior x-column, tiled over all j, against a prescribed wall gradient. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_4_neumann_face1, rewritten for the x-axis so the free name is dx. The datum is the wall derivative q = df/dx at the x-wall (offset 0) together with the five nearest cells f1..f5 (fixed j) at offsets {1/2,3/2,5/2,7/2,9/2}*dx. The unique degree-5 Hermite-type polynomial matching those six data differentiated twice at the first center (offset 1/2 dx) gives D2 = (-19200 q dx - 15285 f1 + 7544 f2 + 12546 f3 - 5784 f4 + 979 f5) / (20268 dx^2). Exact for polynomials up to degree 5, so the local truncation error is O(dx^4) — matching the interior 5-point stencil central_D2x_4_interior. The five value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-x with the matching wall flux is discretely exact. Written for a LEFT (low-index) x-wall where q = +df/dx points into the domain; reused unchanged at the right x-wall by reflection (the second derivative is even under x -> -x while df/dx is odd), so the rule binds q -> -qXR there with f1..f5 mapped to the cells in increasing distance from the right wall. dx is a free name resolving to the consuming model’s x-spacing variable (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_4_neumann_face2

op:D2 order:4 stencil width 6 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_4_neumann_face2.esm

central_D2x_4_neumann_face2(q, f1, f2, f3, f4, f5) = (1920·q·dx + 2.6019×10⁴·f1 + −5.21×10⁴·f2 + 2.8134×10⁴·f3 + −2124·f4 + 71·f5)/(2.0268×10⁴·dx·dx)

Fourth-order one-sided second-derivative closure along x for the SECOND cell from a Neumann (fixed-flux) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_4_neumann_bc) applies it at the second interior x-column, tiled over all j, whose centered 5-point stencil would reach one cell outside the domain. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_4_neumann_face2, rewritten for the x-axis so the free name is dx. The same six data as central_D2x_4_neumann_face1 (the wall derivative q = df/dx at offset 0 and the five nearest cells f1..f5 at offsets {1/2,3/2,5/2,7/2,9/2}*dx) determine the degree-5 Hermite-type polynomial; differentiating it twice at the SECOND center (offset 3/2 dx) gives D2 = (1920 q dx + 26019 f1 - 52100 f2 + 28134 f3 - 2124 f4 + 71 f5) / (20268 dx^2). Exact for polynomials up to degree 5, so the local truncation error is O(dx^4); the five value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-x with the matching wall flux is discretely exact. Written for a LEFT (low-index) x-wall where q = +df/dx points into the domain; reused unchanged at the right x-wall’s second cell by reflection (the second derivative even, df/dx odd), so the rule binds q -> -qXR there with f1..f5 mapped to the cells in increasing distance from the right wall. dx is a free name resolving to the consuming model’s x-spacing variable (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_6_dirichlet_face1

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_6_dirichlet_face1.esm

central_D2x_6_dirichlet_face1(g, f1, f2, f3, f4, f5, f6, f7) = (8.6016×10⁴·g + −1.36136×10⁵·f1 + 5.6628×10⁴·f2 + −3861·f3 + −5720·f4 + 4290·f5 + −1404·f6 + 187·f7)/(2.574×10⁴·dx·dx)

Sixth-order one-sided second-derivative closure along x for the FIRST cell from a Dirichlet (fixed-value) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_6_dirichlet_bc) applies it at the first interior column, tiled over all j. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_dirichlet_face1 (same Fornberg 1988 offset-node weights), rewritten for the x-axis so the free name is dx. The wall value g at offset 0 and the seven nearest cells f1..f7 (fixed j) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dx fit a degree-7 polynomial; differentiating it twice at the first center (offset 1/2 dx) gives D2 = (86016 g - 136136 f1 + 56628 f2 - 3861 f3 - 5720 f4 + 4290 f5 - 1404 f6 + 187 f7) / (25740 dx^2). Exact for polynomials up to degree 7, so the local truncation error is O(dx^6) — matching the interior 7-point stencil central_D2x_6_interior. The eight weights sum to zero, so a steady linear-in-x profile carrying the wall value is discretely exact. Reused unchanged at the far x-wall by reflection (the second derivative is even), with f1..f7 mapped to the cells in increasing distance from that wall. dx is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_6_dirichlet_face2

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_6_dirichlet_face2.esm

central_D2x_6_dirichlet_face2(g, f1, f2, f3, f4, f5, f6, f7) = (−2048·g + 1.1349×10⁴·f1 + −1.7745×10⁴·f2 + 8073·f3 + 1170·f4 + −1105·f5 + 351·f6 + −45·f7)/(7020·dx·dx)

Sixth-order one-sided second-derivative closure along x for the SECOND cell from a Dirichlet (fixed-value) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_6_dirichlet_bc) applies it at the second interior column, tiled over all j. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_dirichlet_face2 (same Fornberg 1988 offset-node weights), rewritten for the x-axis so the free name is dx. The wall value g at offset 0 and the seven nearest cells f1..f7 (fixed j) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dx fit a degree-7 polynomial; differentiating it twice at the second center (offset 3/2 dx) gives D2 = (-2048 g + 11349 f1 - 17745 f2 + 8073 f3 + 1170 f4 - 1105 f5 + 351 f6 - 45 f7) / (7020 dx^2). Exact for polynomials up to degree 7, so the local truncation error is O(dx^6) — matching the interior 7-point stencil central_D2x_6_interior. The eight weights sum to zero, so a steady linear-in-x profile carrying the wall value is discretely exact. Reused unchanged at the far x-wall by reflection (the second derivative is even), with f1..f7 mapped to the cells in increasing distance from that wall. dx is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_6_dirichlet_face3

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_6_dirichlet_face3.esm

central_D2x_6_dirichlet_face3(g, f1, f2, f3, f4, f5, f6, f7) = (4096·g + −1.7589×10⁴·f1 + 1.21836×10⁵·f2 + −2.16216×10⁵·f3 + 1.2012×10⁵·f4 + −1.3585×10⁴·f5 + 1404·f6 + −66·f7)/(7.722×10⁴·dx·dx)

Sixth-order one-sided second-derivative closure along x for the THIRD cell from a Dirichlet (fixed-value) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_6_dirichlet_bc) applies it at the third interior column, tiled over all j. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_dirichlet_face3 (same Fornberg 1988 offset-node weights), rewritten for the x-axis so the free name is dx. The wall value g at offset 0 and the seven nearest cells f1..f7 (fixed j) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dx fit a degree-7 polynomial; differentiating it twice at the third center (offset 5/2 dx) gives D2 = (4096 g - 17589 f1 + 121836 f2 - 216216 f3 + 120120 f4 - 13585 f5 + 1404 f6 - 66 f7) / (77220 dx^2). Exact for polynomials up to degree 7, so the local truncation error is O(dx^6) — matching the interior 7-point stencil central_D2x_6_interior. The eight weights sum to zero, so a steady linear-in-x profile carrying the wall value is discretely exact. Reused unchanged at the far x-wall by reflection (the second derivative is even), with f1..f7 mapped to the cells in increasing distance from that wall. dx is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_6_interior

op:D2 order:6 stencil width 7 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_6_interior.esm

central_D2x_6_interior(f)[i,j] = (2·f[i−3, j] − 27·f[i−2, j] + 270·f[i−1, j] − 490·f[i, j] + 270·f[i+1, j] − 27·f[i+2, j] + 2·f[i+3, j]) / (180·dx·dx) for i ∈ [4, NX−3], j ∈ y

Interior-only 7-point sixth-order centered second derivative along x on cartesian_uniform_2d: (2 f[i-3,j] - 27 f[..] + 270 f[..] - 490 f[i,j] + 270 f[..] - 27 f[..] + 2 f[..]) / (180 dx^2) over x in [4, NX-3], all y. The per-axis lift of the 1-D operator grids/cartesian_uniform_1d/stencils/central_D2_6_interior with the y axis a passive free index (the same coefficients 2,-27,270,-490,270,-27,2 / 180, applied along x at fixed y). Match-less named template; the boundary-condition layer (rules/central_D2x_6_dirichlet_bc.esm and siblings) wraps it with the x-edge sixth-order face closures. dx is a free name resolving to the consuming model’s x-spacing variable (grids/cartesian_uniform_2d/grid.esm). O(dx^6). The 2-D sixth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y), each half lowered by its own per-axis sixth-order rule.

References

central_D2x_6_neumann_face1

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_6_neumann_face1.esm

central_D2x_6_neumann_face1(q, f1, f2, f3, f4, f5, f6, f7) = (−1.354752×10⁷·q·dx + −4.382308×10⁶·f1 + −1.80976×10⁷·f2 + 4.529776×10⁷·f3 + −3.595516×10⁷·f4 + 1.735677×10⁷·f5 + −4.804812×10⁶·f6 + 5.85341×10⁵·f7)/(1.585242×10⁷·dx·dx)

Sixth-order one-sided second-derivative closure along x for the FIRST cell from a Neumann (fixed-flux) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_6_neumann_bc) applies it at the first interior column, tiled over all j, against a prescribed wall gradient. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_neumann_face1, rewritten for the x-axis so the free name is dx. The datum is the wall derivative q = df/dx at the wall (offset 0) together with the seven nearest cells f1..f7 (fixed j) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dx. The degree-7 Hermite-type polynomial matching those eight data differentiated twice at the first center (offset 1/2 dx) gives D2 = (-13547520 q dx - 4382308 f1 - 18097596 f2 + 45297765 f3 - 35955160 f4 + 17356770 f5 - 4804812 f6 + 585341 f7) / (15852420 dx^2). Exact for polynomials up to degree 7, so the local truncation error is O(dx^6) — matching the interior 7-point stencil central_D2x_6_interior. The seven value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-x with the matching wall flux is discretely exact. Written for a LEFT (low-index) x-wall where q = +df/dx points into the domain; reused unchanged at the far x-wall by reflection (the second derivative even, df/dx odd), so the rule binds q -> -qXR there with f1..f7 mapped to the cells in increasing distance from that wall. dx is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_6_neumann_face2

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_6_neumann_face2.esm

central_D2x_6_neumann_face2(q, f1, f2, f3, f4, f5, f6, f7) = (1.18272×10⁶·q·dx + 1.869115×10⁷·f1 + −3.544678×10⁷·f2 + 1.406812×10⁷·f3 + 5.47347×10⁶·f4 + −3.779905×10⁶·f5 + 1.136601×10⁶·f6 + −1.42665×10⁵·f7)/(1.585242×10⁷·dx·dx)

Sixth-order one-sided second-derivative closure along x for the SECOND cell from a Neumann (fixed-flux) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_6_neumann_bc) applies it at the second interior column, tiled over all j, against a prescribed wall gradient. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_neumann_face2, rewritten for the x-axis so the free name is dx. The datum is the wall derivative q = df/dx at the wall (offset 0) together with the seven nearest cells f1..f7 (fixed j) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dx. The degree-7 Hermite-type polynomial matching those eight data differentiated twice at the second center (offset 3/2 dx) gives D2 = (1182720 q dx + 18691149 f1 - 35446775 f2 + 14068125 f3 + 5473470 f4 - 3779905 f5 + 1136601 f6 - 142665 f7) / (15852420 dx^2). Exact for polynomials up to degree 7, so the local truncation error is O(dx^6) — matching the interior 7-point stencil central_D2x_6_interior. The seven value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-x with the matching wall flux is discretely exact. Written for a LEFT (low-index) x-wall where q = +df/dx points into the domain; reused unchanged at the far x-wall by reflection (the second derivative even, df/dx odd), so the rule binds q -> -qXR there with f1..f7 mapped to the cells in increasing distance from that wall. dx is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_6_neumann_face3

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_6_neumann_face3.esm

central_D2x_6_neumann_face3(q, f1, f2, f3, f4, f5, f6, f7) = (−2.1504×10⁵·q·dx + −2.349569×10⁶·f1 + 2.417076×10⁷·f2 + −4.363002×10⁷·f3 + 2.414452×10⁷·f4 + −2.555285×10⁶·f5 + 2.25684×10⁵·f6 + −6086·f7)/(1.585242×10⁷·dx·dx)

Sixth-order one-sided second-derivative closure along x for the THIRD cell from a Neumann (fixed-flux) x-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2x_6_neumann_bc) applies it at the third interior column, tiled over all j, against a prescribed wall gradient. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_neumann_face3, rewritten for the x-axis so the free name is dx. The datum is the wall derivative q = df/dx at the wall (offset 0) together with the seven nearest cells f1..f7 (fixed j) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dx. The degree-7 Hermite-type polynomial matching those eight data differentiated twice at the third center (offset 5/2 dx) gives D2 = (-215040 q dx - 2349569 f1 + 24170756 f2 - 43630020 f3 + 24144520 f4 - 2555285 f5 + 225684 f6 - 6086 f7) / (15852420 dx^2). Exact for polynomials up to degree 7, so the local truncation error is O(dx^6) — matching the interior 7-point stencil central_D2x_6_interior. The seven value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-x with the matching wall flux is discretely exact. Written for a LEFT (low-index) x-wall where q = +df/dx points into the domain; reused unchanged at the far x-wall by reflection (the second derivative even, df/dx odd), so the rule binds q -> -qXR there with f1..f7 mapped to the cells in increasing distance from that wall. dx is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2x_interior

op:D2 order:2 stencil width 3 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2x_interior.esm

central_D2x_interior(f)[i,j] = (f[i+1, j] − 2·f[i, j] + f[i−1, j]) / dx² for i ∈ [2, NX−1], j ∈ y

Interior-only 3-point centered second derivative along x on cartesian_uniform_2d: (f[i+1,j] - 2 f[i,j] + f[i-1,j]) / dx^2 over i in [2, NX-1], all j. Match-less named template; the boundary-condition layer (rules/central_D2x_neumann_bc.esm) wraps it with the x-edge faces. dx is a free name resolving to the consuming model’s x-spacing variable (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). O(dx^2). The 2-D Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y), each half lowered by its own per-axis rule.

References

central_D2y_4_dirichlet_face1

op:D2 order:4 stencil width 6 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_4_dirichlet_face1.esm

central_D2y_4_dirichlet_face1(g, f1, f2, f3, f4, f5) = (2560·g + −4095·f1 + 1848·f2 + −378·f3 + 72·f4 + −7·f5)/(756·dy·dy)

Fourth-order one-sided second-derivative closure along y for the cell nearest a Dirichlet (fixed-value) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_4_dirichlet_bc) applies it at the first interior y-row, tiled over all i. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_4_dirichlet_face1 (same Fornberg 1988 offset-node weights), rewritten for the y-axis so the free name is the y-spacing dy. On this cell-centered grid the y-wall sits half a cell outside the first center, so the wall value g and the five nearest cells along y sit at offsets {0, 1/2, 3/2, 5/2, 7/2, 9/2}*dy from the wall. The unique degree-5 polynomial through those six data, differentiated twice at the first center (offset 1/2 dy), gives D2 = (2560 g - 4095 f1 + 1848 f2 - 378 f3 + 72 f4 - 7 f5) / (756 dy^2), where f1..f5 are the five cells (at fixed i) in order of increasing distance from the y-wall and g is the wall value. Exact for polynomials up to degree 5, so the local truncation error is O(dy^4) — matching the interior 5-point stencil central_D2y_4_interior. The six weights sum to zero, so a steady linear-in-y profile carrying the wall value is discretely exact. Reused unchanged at the top y-wall by reflection (the second derivative is even under y -> -y), with f1..f5 mapped to the cells in increasing distance from the top wall and g bound to the top wall value. dy is a free name resolving to the consuming model’s y-spacing variable (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_4_dirichlet_face2

op:D2 order:4 stencil width 6 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_4_dirichlet_face2.esm

central_D2y_4_dirichlet_face2(g, f1, f2, f3, f4, f5) = (−256·g + 1323·f1 + −2100·f2 + 1134·f3 + −108·f4 + 7·f5)/(756·dy·dy)

Fourth-order one-sided second-derivative closure along y for the SECOND cell from a Dirichlet (fixed-value) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_4_dirichlet_bc) applies it at the second interior y-row, tiled over all i, whose centered 5-point stencil would reach one cell outside the domain. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_4_dirichlet_face2, rewritten for the y-axis so the free name is dy. The same six data as central_D2y_4_dirichlet_face1 (the wall value g at offset 0 and the five nearest cells f1..f5 at offsets {1/2,3/2,5/2,7/2,9/2}*dy) fit a degree-5 polynomial; differentiating it twice at the SECOND center (offset 3/2 dy) gives D2 = (-256 g + 1323 f1 - 2100 f2 + 1134 f3 - 108 f4 + 7 f5) / (756 dy^2). Exact for polynomials up to degree 5, so the local truncation error is O(dy^4); the six weights sum to zero, so a steady linear-in-y profile is discretely exact. Reused unchanged at the top y-wall’s second cell by reflection, with f1..f5 mapped to the cells in increasing distance from the top wall and g bound to the top wall value. dy is a free name resolving to the consuming model’s y-spacing variable (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_4_interior

op:D2 order:4 stencil width 5 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_4_interior.esm

central_D2y_4_interior(f)[i,j] = (−f[i, j−2] + 16·f[i, j−1] − 30·f[i, j] + 16·f[i, j+1] − f[i, j+2]) / (12·dy·dy) for i ∈ x, j ∈ [3, NY−2]

Interior-only 5-point fourth-order centered second derivative along y on cartesian_uniform_2d: (-f[i,j-2] + 16 f[i,j-1] - 30 f[i,j] + 16 f[i,j+1] - f[i,j+2]) / (12 dy^2) over j in [3, NY-2], all i. The per-axis lift of the 1-D operator grids/cartesian_uniform_1d/stencils/central_D2_4_interior with the x axis a passive free index (the same coefficients -1,16,-30,16,-1 / 12, applied along y at fixed i). Match-less named template; the boundary-condition layer (rules/central_D2y_4_dirichlet_bc.esm and siblings) wraps it with the y-edge fourth-order face closures. dy is a free name resolving to the consuming model’s y-spacing variable (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). O(dy^4). The 2-D fourth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y), each half lowered by its own per-axis fourth-order rule.

References

central_D2y_4_neumann_face1

op:D2 order:4 stencil width 6 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_4_neumann_face1.esm

central_D2y_4_neumann_face1(q, f1, f2, f3, f4, f5) = (−1.92×10⁴·q·dy + −1.5285×10⁴·f1 + 7544·f2 + 1.2546×10⁴·f3 + −5784·f4 + 979·f5)/(2.0268×10⁴·dy·dy)

Fourth-order one-sided second-derivative closure along y for the cell nearest a Neumann (fixed-flux) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_4_neumann_bc) applies it at the first interior y-row, tiled over all i, against a prescribed wall gradient. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_4_neumann_face1, rewritten for the y-axis so the free name is dy. The datum is the wall derivative q = df/dy at the y-wall (offset 0) together with the five nearest cells f1..f5 (fixed i) at offsets {1/2,3/2,5/2,7/2,9/2}*dy. The unique degree-5 Hermite-type polynomial matching those six data differentiated twice at the first center (offset 1/2 dy) gives D2 = (-19200 q dy - 15285 f1 + 7544 f2 + 12546 f3 - 5784 f4 + 979 f5) / (20268 dy^2). Exact for polynomials up to degree 5, so the local truncation error is O(dy^4) — matching the interior 5-point stencil central_D2y_4_interior. The five value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-y with the matching wall flux is discretely exact. Written for a BOTTOM (low-index) y-wall where q = +df/dy points into the domain; reused unchanged at the top y-wall by reflection (the second derivative is even under y -> -y while df/dy is odd), so the rule binds q -> -qYR there with f1..f5 mapped to the cells in increasing distance from the top wall. dy is a free name resolving to the consuming model’s y-spacing variable (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_4_neumann_face2

op:D2 order:4 stencil width 6 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_4_neumann_face2.esm

central_D2y_4_neumann_face2(q, f1, f2, f3, f4, f5) = (1920·q·dy + 2.6019×10⁴·f1 + −5.21×10⁴·f2 + 2.8134×10⁴·f3 + −2124·f4 + 71·f5)/(2.0268×10⁴·dy·dy)

Fourth-order one-sided second-derivative closure along y for the SECOND cell from a Neumann (fixed-flux) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_4_neumann_bc) applies it at the second interior y-row, tiled over all i, whose centered 5-point stencil would reach one cell outside the domain. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_4_neumann_face2, rewritten for the y-axis so the free name is dy. The same six data as central_D2y_4_neumann_face1 (the wall derivative q = df/dy at offset 0 and the five nearest cells f1..f5 at offsets {1/2,3/2,5/2,7/2,9/2}*dy) determine the degree-5 Hermite-type polynomial; differentiating it twice at the SECOND center (offset 3/2 dy) gives D2 = (1920 q dy + 26019 f1 - 52100 f2 + 28134 f3 - 2124 f4 + 71 f5) / (20268 dy^2). Exact for polynomials up to degree 5, so the local truncation error is O(dy^4); the five value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-y with the matching wall flux is discretely exact. Written for a BOTTOM (low-index) y-wall where q = +df/dy points into the domain; reused unchanged at the top y-wall’s second cell by reflection (the second derivative even, df/dy odd), so the rule binds q -> -qYR there with f1..f5 mapped to the cells in increasing distance from the top wall. dy is a free name resolving to the consuming model’s y-spacing variable (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_6_dirichlet_face1

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_6_dirichlet_face1.esm

central_D2y_6_dirichlet_face1(g, f1, f2, f3, f4, f5, f6, f7) = (8.6016×10⁴·g + −1.36136×10⁵·f1 + 5.6628×10⁴·f2 + −3861·f3 + −5720·f4 + 4290·f5 + −1404·f6 + 187·f7)/(2.574×10⁴·dy·dy)

Sixth-order one-sided second-derivative closure along y for the FIRST cell from a Dirichlet (fixed-value) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_6_dirichlet_bc) applies it at the first interior row, tiled over all i. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_dirichlet_face1 (same Fornberg 1988 offset-node weights), rewritten for the y-axis so the free name is dy. The wall value g at offset 0 and the seven nearest cells f1..f7 (fixed i) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dy fit a degree-7 polynomial; differentiating it twice at the first center (offset 1/2 dy) gives D2 = (86016 g - 136136 f1 + 56628 f2 - 3861 f3 - 5720 f4 + 4290 f5 - 1404 f6 + 187 f7) / (25740 dy^2). Exact for polynomials up to degree 7, so the local truncation error is O(dy^6) — matching the interior 7-point stencil central_D2y_6_interior. The eight weights sum to zero, so a steady linear-in-y profile carrying the wall value is discretely exact. Reused unchanged at the far y-wall by reflection (the second derivative is even), with f1..f7 mapped to the cells in increasing distance from that wall. dy is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_6_dirichlet_face2

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_6_dirichlet_face2.esm

central_D2y_6_dirichlet_face2(g, f1, f2, f3, f4, f5, f6, f7) = (−2048·g + 1.1349×10⁴·f1 + −1.7745×10⁴·f2 + 8073·f3 + 1170·f4 + −1105·f5 + 351·f6 + −45·f7)/(7020·dy·dy)

Sixth-order one-sided second-derivative closure along y for the SECOND cell from a Dirichlet (fixed-value) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_6_dirichlet_bc) applies it at the second interior row, tiled over all i. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_dirichlet_face2 (same Fornberg 1988 offset-node weights), rewritten for the y-axis so the free name is dy. The wall value g at offset 0 and the seven nearest cells f1..f7 (fixed i) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dy fit a degree-7 polynomial; differentiating it twice at the second center (offset 3/2 dy) gives D2 = (-2048 g + 11349 f1 - 17745 f2 + 8073 f3 + 1170 f4 - 1105 f5 + 351 f6 - 45 f7) / (7020 dy^2). Exact for polynomials up to degree 7, so the local truncation error is O(dy^6) — matching the interior 7-point stencil central_D2y_6_interior. The eight weights sum to zero, so a steady linear-in-y profile carrying the wall value is discretely exact. Reused unchanged at the far y-wall by reflection (the second derivative is even), with f1..f7 mapped to the cells in increasing distance from that wall. dy is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_6_dirichlet_face3

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_6_dirichlet_face3.esm

central_D2y_6_dirichlet_face3(g, f1, f2, f3, f4, f5, f6, f7) = (4096·g + −1.7589×10⁴·f1 + 1.21836×10⁵·f2 + −2.16216×10⁵·f3 + 1.2012×10⁵·f4 + −1.3585×10⁴·f5 + 1404·f6 + −66·f7)/(7.722×10⁴·dy·dy)

Sixth-order one-sided second-derivative closure along y for the THIRD cell from a Dirichlet (fixed-value) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_6_dirichlet_bc) applies it at the third interior row, tiled over all i. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_dirichlet_face3 (same Fornberg 1988 offset-node weights), rewritten for the y-axis so the free name is dy. The wall value g at offset 0 and the seven nearest cells f1..f7 (fixed i) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dy fit a degree-7 polynomial; differentiating it twice at the third center (offset 5/2 dy) gives D2 = (4096 g - 17589 f1 + 121836 f2 - 216216 f3 + 120120 f4 - 13585 f5 + 1404 f6 - 66 f7) / (77220 dy^2). Exact for polynomials up to degree 7, so the local truncation error is O(dy^6) — matching the interior 7-point stencil central_D2y_6_interior. The eight weights sum to zero, so a steady linear-in-y profile carrying the wall value is discretely exact. Reused unchanged at the far y-wall by reflection (the second derivative is even), with f1..f7 mapped to the cells in increasing distance from that wall. dy is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_6_interior

op:D2 order:6 stencil width 7 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_6_interior.esm

central_D2y_6_interior(f)[i,j] = (2·f[i, j−3] − 27·f[i, j−2] + 270·f[i, j−1] − 490·f[i, j] + 270·f[i, j+1] − 27·f[i, j+2] + 2·f[i, j+3]) / (180·dy·dy) for i ∈ x, j ∈ [4, NY−3]

Interior-only 7-point sixth-order centered second derivative along y on cartesian_uniform_2d: (2 f[i,j-3] - 27 f[..] + 270 f[..] - 490 f[i,j] + 270 f[..] - 27 f[..] + 2 f[..]) / (180 dy^2) over y in [4, NY-3], all x. The per-axis lift of the 1-D operator grids/cartesian_uniform_1d/stencils/central_D2_6_interior with the x axis a passive free index (the same coefficients 2,-27,270,-490,270,-27,2 / 180, applied along y at fixed x). Match-less named template; the boundary-condition layer (rules/central_D2y_6_dirichlet_bc.esm and siblings) wraps it with the y-edge sixth-order face closures. dy is a free name resolving to the consuming model’s y-spacing variable (grids/cartesian_uniform_2d/grid.esm). O(dy^6). The 2-D sixth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y), each half lowered by its own per-axis sixth-order rule.

References

central_D2y_6_neumann_face1

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_6_neumann_face1.esm

central_D2y_6_neumann_face1(q, f1, f2, f3, f4, f5, f6, f7) = (−1.354752×10⁷·q·dy + −4.382308×10⁶·f1 + −1.80976×10⁷·f2 + 4.529776×10⁷·f3 + −3.595516×10⁷·f4 + 1.735677×10⁷·f5 + −4.804812×10⁶·f6 + 5.85341×10⁵·f7)/(1.585242×10⁷·dy·dy)

Sixth-order one-sided second-derivative closure along y for the FIRST cell from a Neumann (fixed-flux) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_6_neumann_bc) applies it at the first interior row, tiled over all i, against a prescribed wall gradient. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_neumann_face1, rewritten for the y-axis so the free name is dy. The datum is the wall derivative q = df/dy at the wall (offset 0) together with the seven nearest cells f1..f7 (fixed i) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dy. The degree-7 Hermite-type polynomial matching those eight data differentiated twice at the first center (offset 1/2 dy) gives D2 = (-13547520 q dy - 4382308 f1 - 18097596 f2 + 45297765 f3 - 35955160 f4 + 17356770 f5 - 4804812 f6 + 585341 f7) / (15852420 dy^2). Exact for polynomials up to degree 7, so the local truncation error is O(dy^6) — matching the interior 7-point stencil central_D2y_6_interior. The seven value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-y with the matching wall flux is discretely exact. Written for a BOTTOM (low-index) y-wall where q = +df/dy points into the domain; reused unchanged at the far y-wall by reflection (the second derivative even, df/dy odd), so the rule binds q -> -qYR there with f1..f7 mapped to the cells in increasing distance from that wall. dy is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_6_neumann_face2

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_6_neumann_face2.esm

central_D2y_6_neumann_face2(q, f1, f2, f3, f4, f5, f6, f7) = (1.18272×10⁶·q·dy + 1.869115×10⁷·f1 + −3.544678×10⁷·f2 + 1.406812×10⁷·f3 + 5.47347×10⁶·f4 + −3.779905×10⁶·f5 + 1.136601×10⁶·f6 + −1.42665×10⁵·f7)/(1.585242×10⁷·dy·dy)

Sixth-order one-sided second-derivative closure along y for the SECOND cell from a Neumann (fixed-flux) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_6_neumann_bc) applies it at the second interior row, tiled over all i, against a prescribed wall gradient. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_neumann_face2, rewritten for the y-axis so the free name is dy. The datum is the wall derivative q = df/dy at the wall (offset 0) together with the seven nearest cells f1..f7 (fixed i) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dy. The degree-7 Hermite-type polynomial matching those eight data differentiated twice at the second center (offset 3/2 dy) gives D2 = (1182720 q dy + 18691149 f1 - 35446775 f2 + 14068125 f3 + 5473470 f4 - 3779905 f5 + 1136601 f6 - 142665 f7) / (15852420 dy^2). Exact for polynomials up to degree 7, so the local truncation error is O(dy^6) — matching the interior 7-point stencil central_D2y_6_interior. The seven value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-y with the matching wall flux is discretely exact. Written for a BOTTOM (low-index) y-wall where q = +df/dy points into the domain; reused unchanged at the far y-wall by reflection (the second derivative even, df/dy odd), so the rule binds q -> -qYR there with f1..f7 mapped to the cells in increasing distance from that wall. dy is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_6_neumann_face3

op:D2 order:6 stencil width 8 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_6_neumann_face3.esm

central_D2y_6_neumann_face3(q, f1, f2, f3, f4, f5, f6, f7) = (−2.1504×10⁵·q·dy + −2.349569×10⁶·f1 + 2.417076×10⁷·f2 + −4.363002×10⁷·f3 + 2.414452×10⁷·f4 + −2.555285×10⁶·f5 + 2.25684×10⁵·f6 + −6086·f7)/(1.585242×10⁷·dy·dy)

Sixth-order one-sided second-derivative closure along y for the THIRD cell from a Neumann (fixed-flux) y-wall on cartesian_uniform_2d, cell-centered. Match-less named scalar template; the boundary-condition rule (rules/central_D2y_6_neumann_bc) applies it at the third interior row, tiled over all i, against a prescribed wall gradient. Identical algebra to the 1-D closure grids/cartesian_uniform_1d/stencils/central_D2_6_neumann_face3, rewritten for the y-axis so the free name is dy. The datum is the wall derivative q = df/dy at the wall (offset 0) together with the seven nearest cells f1..f7 (fixed i) at offsets {1/2,3/2,5/2,7/2,9/2,11/2,13/2}*dy. The degree-7 Hermite-type polynomial matching those eight data differentiated twice at the third center (offset 5/2 dy) gives D2 = (-215040 q dy - 2349569 f1 + 24170756 f2 - 43630020 f3 + 24144520 f4 - 2555285 f5 + 225684 f6 - 6086 f7) / (15852420 dy^2). Exact for polynomials up to degree 7, so the local truncation error is O(dy^6) — matching the interior 7-point stencil central_D2y_6_interior. The seven value-weights sum to zero and the derivative term supplies the exact wall slope, so a steady ramp-in-y with the matching wall flux is discretely exact. Written for a BOTTOM (low-index) y-wall where q = +df/dy points into the domain; reused unchanged at the far y-wall by reflection (the second derivative even, df/dy odd), so the rule binds q -> -qYR there with f1..f7 mapped to the cells in increasing distance from that wall. dy is a free name (grids/cartesian_uniform_2d/grid.esm).

References

central_D2y_interior

op:D2 order:2 stencil width 3 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_D2y_interior.esm

central_D2y_interior(f)[i,j] = (f[i, j+1] − 2·f[i, j] + f[i, j−1]) / dy² for i ∈ x, j ∈ [2, NY−1]

Interior-only 3-point centered second derivative along y on cartesian_uniform_2d: (f[i,j+1] - 2 f[i,j] + f[i,j-1]) / dy^2 over j in [2, NY-1], all i. Match-less named template; the boundary-condition layer (rules/central_D2y_neumann_bc.esm) wraps it with the y-edge faces. dy is a free name resolving to the consuming model’s y-spacing variable (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). O(dy^2). The 2-D Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y), each half lowered by its own per-axis rule.

References

central_Dxy_interior

op:Dxy order:2 stencil width 3 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/central_Dxy_interior.esm

central_Dxy_interior(f)[i,j] = (f[i+1, j+1] − f[i+1, j−1] − (f[i−1, j+1] − f[i−1, j−1])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [2, NY−1]

Interior-only 4-point centered MIXED second derivative d^2f/dx dy along the x and y axes on cartesian_uniform_2d: (f[i+1,j+1] - f[i+1,j-1] - f[i-1,j+1] + f[i-1,j-1]) / (4 dx dy) over i in [2, NX-1] and j in [2, NY-1] (the doubly-interior block, both axes one cell in from either wall so all four corner samples f[i+-1,j+-1] exist). Match-less named template; the boundary-condition layer (rules/central_Dxy_periodic.esm) wraps it with the x- and y-edge and corner faces. This scheme is the tensor product of the two centered first differences (central_D1x_interior applied to central_D1y_interior), so it inherits their O(dx^2 + dy^2) accuracy; it is the off-diagonal term of a full anisotropic diffusion tensor K = [[Kxx, Kxy], [Kxy, Kyy]], the 2 Kxy d^2f/dx dy contribution that the diagonal central_D2x/central_D2y stencils do not supply. dx and dy are free names resolving to the consuming model’s x- and y-spacing variables (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). O(dx^2 + dy^2).

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (mixed partial derivatives as products of one-dimensional centered differences).

godunov_norm_D2_interior

op:grad_norm order:1 stencil width 3 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/godunov_norm_D2_interior.esm

godunov_norm_D2_interior(u)[i,j] = sqrt(max((u[i, j] − u[i−1, j]) / dx, 0)² + min((u[i+1, j] − u[i, j]) / dx, 0)² + max((u[i, j] − u[i, j−1]) / dy, 0)² + min((u[i, j+1] − u[i, j]) / dy, 0)²) for i ∈ [2, NX−1], j ∈ [2, NY−1]

Interior-only first-order Godunov (Rouy-Tourin / Osher-Sethian) upwind gradient magnitude |grad u| = sqrt(u_x^2 + u_y^2) on cartesian_uniform_2d over the doubly-interior block [2, NX-1] x [2, NY-1] where the 5-point stencil u[i,j], u[i-1,j], u[i+1,j], u[i,j-1], u[i,j+1] fits without reaching a boundary. Per-axis backward D-_x=(u[i,j]-u[i-1,j])/dx, D-_y=(u[i,j]-u[i,j-1])/dy and forward D+_x=(u[i+1,j]-u[i,j])/dx, D+_y=(u[i,j+1]-u[i,j])/dy; the norm is sqrt( max(D-_x,0)^2 + min(D+_x,0)^2 + max(D-_y,0)^2 + min(D+_y,0)^2 ), the per-axis entropy (upwind) fix of Osher-Sethian (1988) / Rouy-Tourin (1992). Match-less named template extracting the doubly-interior aggregate of grids/cartesian_uniform_2d/rules/godunov_norm_D2_periodic verbatim so a boundary-condition layer (rules/godunov_norm_D2_dirichlet_bc.esm) wraps it with prescribed-value edge and corner faces via apply_expression_template, exactly as central_D2x_dirichlet_bc reuses central_D2x_interior. EXACT on a linear field u = a x + b y + c (D-_x=D+_x=a, D-_y=D+_y=b, norm = sqrt(a^2+b^2) = |grad u|). First order O(dx + dy) on smooth data, uniform through the |grad| kinks. dx, dy are free names resolving to the consuming model’s spacing variables (grids/cartesian_uniform_2d/grid.esm consumer-supplied geometry contract).

References

  • Osher, S., & Sethian, J. A. (1988). Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79(1), 12-49. The max/min upwind (entropy) construction of |grad phi| for the level-set Hamilton-Jacobi term. doi:10.1016/0021-9991(88)90002-2
  • Rouy, E., & Tourin, A. (1992). A viscosity solutions approach to shape-from-shading. SIAM Journal on Numerical Analysis, 29(3), 867-884. The Godunov Hamiltonian sqrt(max(D-,0)^2 + min(D+,0)^2) as the monotone, first-order-convergent discretization of the eikonal |grad u|. doi:10.1137/0729053

varcoeff_laplacian_x_interior

op:varcoeff_laplacian order:2 stencil width 3 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/varcoeff_laplacian_x_interior.esm

varcoeff_laplacian_x_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i+1, j]))·(u[i+1, j] − u[i, j]) − (0.5·(k[i−1, j] + k[i, j]))·(u[i, j] − u[i−1, j])) / dx² for i ∈ [2, NX−1], j ∈ y

Interior-only conservative finite-volume x-half of the variable-coefficient / nonlinear Laplacian div(k grad u) on cartesian_uniform_2d: the x-flux-divergence Lx u[i,j] = ( kf[i+1/2,j] (u[i+1,j]-u[i,j]) - kf[i-1/2,j] (u[i,j]-u[i-1,j]) ) / dx^2 over i in [2, NX-1], all j, the difference of the two x-face-normal fluxes divided by the uniform cell measure dx^2. This is the exact 2-D per-axis lift of grids/cartesian_uniform_1d/stencils/varcoeff_laplacian_interior with the y index a free (full-rank) spectator: the same arithmetic-mean face algebra runs independently on every y-row j. The x-face conductivity is the ARITHMETIC MEAN of the two adjacent cell values along x, kf[i+1/2,j] = (k[i,j]+k[i+1,j])/2; this is the standard second-order harmonic-vs-arithmetic choice for a smooth coefficient (both are O(dx^2)-accurate at the face for a differentiable k, and the arithmetic mean keeps the operator affine in k so it serves the nonlinear case k=k(u) with no reformulation). Both the state u and the coefficient k are consumer-supplied fields of shape [x, y] passed as template parameters (bound by the matching rule from the compound D(kD(u,x),x) the consumer writes); the SAME rule serves the variable-coefficient case (k=k(x,y) a given field) and the nonlinear case (k declared as an observed k=f(u)) because the body only reads k by index at the x-faces. This interior aggregate loops over i and j; a coefficient k supplied as an aggregate observed must therefore loop over DIFFERENT indices (the library convention is gi/gj, as the grid geometry observeds use) so its inlining into this aggregate is not captured — see the rule (rules/varcoeff_laplacian_x_zero_flux_bc) authoring contract. Match-less named template; the boundary-condition layer (rules/varcoeff_laplacian_x_zero_flux_bc et al.) wraps it into the complete rewrite rule and adds the two x-edge faces. dx is a free name resolving to the consuming model’s x-spacing variable (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(k*D(u,y),y), each half lowered by its own per-axis rule; this stencil is the interior of the first half. Reduces to the constant-coefficient 3-point (k/dx^2)(u[i+1,j]-2u[i,j]+u[i-1,j]) when k is spatially constant.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient / conservative form) and Chapter 3 (multidimensional operators by dimension-wise assembly).

varcoeff_laplacian_y_interior

op:varcoeff_laplacian order:2 stencil width 3 interior-only (match-less)

Source: grids/cartesian_uniform_2d/stencils/varcoeff_laplacian_y_interior.esm

varcoeff_laplacian_y_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i, j+1]))·(u[i, j+1] − u[i, j]) − (0.5·(k[i, j−1] + k[i, j]))·(u[i, j] − u[i, j−1])) / dy² for i ∈ x, j ∈ [2, NY−1]

Interior-only conservative finite-volume y-half of the variable-coefficient / nonlinear Laplacian div(k grad u) on cartesian_uniform_2d: the y-flux-divergence Ly u[i,j] = ( kf[i,j+1/2] (u[i,j+1]-u[i,j]) - kf[i,j-1/2] (u[i,j]-u[i,j-1]) ) / dy^2 over j in [2, NY-1], all i, the difference of the two y-face-normal fluxes divided by the uniform cell measure dy^2. This is the exact 2-D per-axis lift of grids/cartesian_uniform_1d/stencils/varcoeff_laplacian_interior with the x index a free (full-rank) spectator: the same arithmetic-mean face algebra runs independently on every x-column i. The y-face conductivity is the ARITHMETIC MEAN of the two adjacent cell values along y, kf[i,j+1/2] = (k[i,j]+k[i,j+1])/2; this is the standard second-order harmonic-vs-arithmetic choice for a smooth coefficient (both are O(dy^2)-accurate at the face for a differentiable k, and the arithmetic mean keeps the operator affine in k so it serves the nonlinear case k=k(u) with no reformulation). Both the state u and the coefficient k are consumer-supplied fields of shape [x, y] passed as template parameters (bound by the matching rule from the compound D(kD(u,y),y) the consumer writes); the SAME rule serves the variable-coefficient case (k=k(x,y) a given field) and the nonlinear case (k declared as an observed k=f(u)) because the body only reads k by index at the y-faces. This interior aggregate loops over i and j; a coefficient k supplied as an aggregate observed must therefore loop over DIFFERENT indices (the library convention is gi/gj, as the grid geometry observeds use) so its inlining into this aggregate is not captured — see the rule (rules/varcoeff_laplacian_y_zero_flux_bc) authoring contract. Match-less named template; the boundary-condition layer (rules/varcoeff_laplacian_y_zero_flux_bc et al.) wraps it into the complete rewrite rule and adds the two y-edge faces. dy is a free name resolving to the consuming model’s y-spacing variable (the grid’s consumer-supplied geometry contract; see grids/cartesian_uniform_2d/grid.esm). The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(k*D(u,y),y), each half lowered by its own per-axis rule; this stencil is the interior of the second half. Reduces to the constant-coefficient 3-point (k/dy^2)(u[i,j+1]-2u[i,j]+u[i,j-1]) when k is spatially constant.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient / conservative form) and Chapter 3 (multidimensional operators by dimension-wise assembly).

Rules

Complete auto-applied rewrite rules on spatial D: imported stencil + boundary conditions in one makearray (esm-spec §9.6.8). Import a rule and every matching derivative in your model lowers through it.

central_D1x_dirichlet_bc

op:D order:2 bc:dirichlet axes:x,y stencil width 3

Source: grids/cartesian_uniform_2d/rules/central_D1x_dirichlet_bc.esm

Rewrites ∂f/∂x — match pattern:

{"op": "D", "args": ["f"], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorcentral_D1x_interior(f)[i,j] = (f[i+1, j] − f[i−1, j]) / (2·dx) for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face(f[2, j] + f[1, j] − 2·gL) / (2·dx) for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face(2·gR − f[NX, j] − f[NX−1, j]) / (2·dx) for i ∈ [NX, NX], j ∈ y

Second-order centered first derivative along x with inhomogeneous DIRICHLET (fixed-value) conditions on the two x-walls, over the full 2-D field f[x,y] on cartesian_uniform_2d, cell-centered — the first-derivative companion of central_D2x_dirichlet_bc (same mirror ghost, different stencil). Interior i in [2, NX-1] (all j): the imported 3-point centered x-stencil (f[i+1,j]-f[i-1,j])/(2 dx). Faces i=1 and i=NX: mirror-ghost elimination against a prescribed wall value. gL and gR are consumer-supplied free names (the same free-name contract as the grid geometry dx; see grids/cartesian_uniform_2d/grid.esm) giving the field value on the left x-wall x=x0 and the right x-wall x=x0+NXdx, constant along that wall (all j). On this cell-centered grid the wall sits half a cell outside the first center, so linear reconstruction to the wall gives the ghost f[0,j]=2gL-f[1,j] (and f[NX+1,j]=2gR-f[NX,j]); substituting into the centered stencil (f[2,j]-f[0,j])/(2 dx) and (f[NX+1,j]-f[NX-1,j])/(2 dx) yields D1x[1,j] = (f[2,j]+f[1,j]-2gL)/(2 dx) and D1x[NX,j] = (2*gR-f[NX,j]-f[NX-1,j])/(2 dx). A consuming model declares gL/gR as ordinary real parameters; a parameter not overridden falls back to its default (esm-spec §6: default, or 0), so gL=gR=0 is the homogeneous-Dirichlet x-wall and both wall terms drop out. The gL/gR terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). dx is the grid’s consumer-supplied x-spacing free name. The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring central_D2x_dirichlet_bc and central_D1x_zero_grad_bc, so every binding’s simulator fills the (1 x NY) regions without rank promotion. FACE TRUNCATION: D1x[1,j] - u_x(x_1,y_j) = (dx/8) u_xx(x0,y_j) + O(dx^2), so the mirror-ghost face is O(dx) locally for a generic field (SUPERCONVERGENT to O(dx^2) when u_xx=0 at the wall). WELL-POSED USE: a fixed value pins the field at the wall, the DIFFUSION operator’s natural boundary condition, so this gradient closure belongs in a well-posed advection-DIFFUSION setting (or flow into the wall), NOT a pure first-order advection with a fixed value at an OUTFLOW wall, where the pinned value reflects (a Kreiss/GKS-unstable mode); the MMS driver problems/advection_2d_x_dirichlet_central.esm exercises it as the advective x-term of an advection-diffusion decay problem whose physical x-diffusion supplies the boundary dissipation, confirming 2nd order. Matches plain D(f, wrt: x) at default priority, so compound rules (e.g. a second-derivative rule at priority 10) fire first, match-scoped with where: {f: {shape: [x, y]}} (esm-spec §9.6.1) so it fires only on the first x-derivative of a bare field declared over this 2-D grid’s [x, y] index sets (distinct from the 1-D shape-[x] rules); under import-edge renaming (esm-spec §9.7.7) the wrt literal and the where shape follow the renamed axes together. A model imports exactly one first-derivative rule for x (this Dirichlet rule, the zero-gradient rule, or the periodic variant), never more than one.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (cell-centered Dirichlet via mirror ghost) and Chapter 10 (well-posed boundary conditions for advection-diffusion).
  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0

Conformance fixtures

Convergence — case advection_2d_x_dirichlet_central

Error norms read from the committed golden tests/conformance/convergence/advection_2d_x_dirichlet_central/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
165.690e-046.559e-03
321.442e-041.981.829e-031.84
643.617e-052.004.809e-041.93
1289.050e-062.001.232e-041.96

Expected order: 2 (± 0.2).

Convergence of advection_2d_x_dirichlet_central (log-log)

central_D1x_periodic

op:D order:2 bc:periodic axes:x,y stencil width 3

Source: grids/cartesian_uniform_2d/rules/central_D1x_periodic.esm

Rewrites ∂f/∂x — match pattern:

{"op": "D", "args": ["f"], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorcentral_D1x_interior(f)[i,j] = (f[i+1, j] − f[i−1, j]) / (2·dx) for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face(f[2, j] − f[NX, j]) / (2·dx) for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face(f[1, j] − f[NX−1, j]) / (2·dx) for i ∈ [NX, NX], j ∈ y

Second-order centered first derivative along x with periodic (wrap) boundary on cartesian_uniform_2d, over the full 2-D field f[x,y]. Interior i in [2, NX-1] (all j): the imported x-stencil (f[i+1,j]-f[i-1,j])/(2 dx). The two x-face columns wrap periodically: i=1 uses the ghost f[0,j]=f[NX,j], giving (f[2,j]-f[NX,j])/(2 dx); i=NX uses f[NX+1,j]=f[1,j], giving (f[1,j]-f[NX-1,j])/(2 dx). The scheme is the same O(dx^2) centered difference at every column and conserves sum_x(f) exactly (the wrap telescopes). The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) rather than the esm-spec 9.6.8 worked example’s rank-reduced form, mirroring grids/latlon/rules/central_D_lon_periodic so every binding’s simulator fills the (1 x NY) regions without rank promotion. dx is the grid’s consumer-supplied x-spacing free name; the periodic wrap is dx-scale-invariant, so the rule serves any x-extent unchanged. Matches plain D(f, wrt: x) at default priority, so compound rules (e.g. a second-derivative rule at priority 10) fire first, match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the first x-derivative of a bare field declared over this grid’s [x, y] index sets (distinct from the 1-D shape-[x] rules); under import-edge renaming (esm-spec 9.7.7) the wrt literal and the where shape follow the renamed axes together, so two instances coexist scoped to their own meshes. O(dx^2). A model imports exactly one first-derivative rule for x (this centered periodic rule or the zero-gradient variant), never both.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains).

Conformance fixtures

Convergence — case advection_2d_x_periodic_central

Error norms read from the committed golden tests/conformance/convergence/advection_2d_x_periodic_central/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
161.578e-021.568e-02
323.964e-031.993.964e-031.98
649.924e-042.009.919e-042.00
1282.482e-042.002.482e-042.00

Expected order: 2 (± 0.2).

Convergence of advection_2d_x_periodic_central (log-log)

central_D1x_zero_grad_bc

op:D order:2 bc:zero_gradient axes:x,y stencil width 3

Source: grids/cartesian_uniform_2d/rules/central_D1x_zero_grad_bc.esm

Rewrites ∂f/∂x — match pattern:

{"op": "D", "args": ["f"], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorcentral_D1x_interior(f)[i,j] = (f[i+1, j] − f[i−1, j]) / (2·dx) for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face(−3·f[1, j] + 4·f[2, j] − f[3, j]) / (2·dx) for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face((f[NX−2, j] − 4·f[NX−1, j]) + 3·f[NX, j]) / (2·dx) for i ∈ [NX, NX], j ∈ y

Second-order centered first derivative along x with one-sided (zero-gradient frame) boundary columns on cartesian_uniform_2d, over the full 2-D field f[x,y]. Interior i in [2, NX-1] (all j): the imported x-stencil (f[i+1,j]-f[i-1,j])/(2 dx). The two x-face columns use 2nd-order one-sided differences that reach only into the interior (no ghost, no periodic wrap): i=1 gives (-3 f[1,j]+4 f[2,j]-f[3,j])/(2 dx) and i=NX gives (f[NX-2,j]-4 f[NX-1,j]+3 f[NX,j])/(2 dx), so the discretized derivative is fully defined with its BC and stays O(dx^2) at every column (esm-spec 9.6.8). The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring grids/cartesian_uniform_2d/rules/central_D2y_neumann_bc so every binding’s simulator fills the (1 x NY) regions without rank promotion. dx is the grid’s consumer-supplied x-spacing free name. The one-sided faces are O(dx^2) locally and preserve global 2nd-order convergence for boundary-compatible solutions — those whose odd x-derivatives vanish at the walls, so the discrete derivative matches the zero-gradient closure superconvergently (the MMS driver problems/advection_2d_x_zero_grad_central.esm uses cos(pi x), whose x-derivative is zero on both walls). Matches plain D(f, wrt: x) at default priority, so compound rules (e.g. a second-derivative rule at priority 10) fire first, match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the first x-derivative of a bare field declared over this grid’s [x, y] index sets (distinct from the 1-D shape-[x] rules); under import-edge renaming (esm-spec 9.7.7) the wrt literal and the where shape follow the renamed axes together. A model imports exactly one first-derivative rule for x (this zero-gradient rule or the periodic variant), never both.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (one-sided boundary treatment).
  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0

Conformance fixtures

Convergence — case advection_2d_x_zero_grad_central

Error norms read from the committed golden tests/conformance/convergence/advection_2d_x_zero_grad_central/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
166.276e-046.034e-04
321.527e-042.041.513e-042.00
643.799e-052.013.782e-052.00
1289.494e-062.009.458e-062.00

Expected order: 2 (± 0.2).

Convergence of advection_2d_x_zero_grad_central (log-log)

central_D1y_dirichlet_bc

op:D order:2 bc:dirichlet axes:x,y stencil width 3

Source: grids/cartesian_uniform_2d/rules/central_D1y_dirichlet_bc.esm

Rewrites ∂f/∂y — match pattern:

{"op": "D", "args": ["f"], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorcentral_D1y_interior(f)[i,j] = (f[i, j+1] − f[i, j−1]) / (2·dy) for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face(f[i, 2] + f[i, 1] − 2·gB) / (2·dy) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face(2·gT − f[i, NY] − f[i, NY−1]) / (2·dy) for i ∈ x, j ∈ [NY, NY]

Second-order centered first derivative along y with inhomogeneous DIRICHLET (fixed-value) conditions on the two y-walls, over the full 2-D field f[x,y] on cartesian_uniform_2d, cell-centered — the first-derivative companion of central_D2y_dirichlet_bc (same mirror ghost, different stencil). Interior j in [2, NY-1] (all i): the imported 3-point centered y-stencil (f[i,j+1]-f[i,j-1])/(2 dy). Faces j=1 and j=NY: mirror-ghost elimination against a prescribed wall value. gB and gT are consumer-supplied free names (the same free-name contract as the grid geometry dy; see grids/cartesian_uniform_2d/grid.esm) giving the field value on the BOTTOM y-wall y=y0 (the j=1 side) and the TOP y-wall y=y0+NYdy (the j=NY side), constant along that wall (all i). On this cell-centered grid the wall sits half a cell outside the first center, so linear reconstruction to the wall gives the ghost f[i,0]=2gB-f[i,1] (and f[i,NY+1]=2gT-f[i,NY]); substituting into the centered stencil (f[i,2]-f[i,0])/(2 dy) and (f[i,NY+1]-f[i,NY-1])/(2 dy) yields D1y[i,1] = (f[i,2]+f[i,1]-2gB)/(2 dy) and D1y[i,NY] = (2*gT-f[i,NY]-f[i,NY-1])/(2 dy). A consuming model declares gB/gT as ordinary real parameters; a parameter not overridden falls back to its default (esm-spec §6: default, or 0), so gB=gT=0 is the homogeneous-Dirichlet y-wall and both wall terms drop out. The gB/gT terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). dy is the grid’s consumer-supplied y-spacing free name. The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring central_D2y_dirichlet_bc and central_D1y_zero_grad_bc, so every binding’s simulator fills the (NX x 1) regions without rank promotion. FACE TRUNCATION: D1y[i,1] - u_y(x_i,y_1) = (dy/8) u_yy(x_i,y0) + O(dy^2), so the mirror-ghost face is O(dy) locally for a generic field (SUPERCONVERGENT to O(dy^2) when u_yy=0 at the wall). WELL-POSED USE: a fixed value pins the field at the wall, the DIFFUSION operator’s natural boundary condition, so this gradient closure belongs in a well-posed advection-DIFFUSION setting (or flow into the wall), NOT a pure first-order advection with a fixed value at an OUTFLOW wall, where the pinned value reflects (a Kreiss/GKS-unstable mode); the MMS driver problems/advection_2d_y_dirichlet_central.esm exercises it as the advective y-term of an advection-diffusion decay problem whose physical y-diffusion supplies the boundary dissipation, confirming 2nd order. Matches plain D(f, wrt: y) at default priority, so compound rules (e.g. a second-derivative rule at priority 10) fire first, match-scoped with where: {f: {shape: [x, y]}} (esm-spec §9.6.1) so it fires only on the first y-derivative of a bare field declared over this 2-D grid’s [x, y] index sets; under import-edge renaming (esm-spec §9.7.7) the wrt literal and the where shape follow the renamed axes together. A model imports exactly one first-derivative rule for y (this Dirichlet rule, the zero-gradient rule, or the periodic variant), never more than one.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (cell-centered Dirichlet via mirror ghost) and Chapter 10 (well-posed boundary conditions for advection-diffusion).
  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0

Conformance fixtures

Convergence — case advection_2d_y_dirichlet_central

Error norms read from the committed golden tests/conformance/convergence/advection_2d_y_dirichlet_central/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
165.690e-046.559e-03
321.442e-041.981.829e-031.84
643.617e-052.004.809e-041.93
1289.050e-062.001.232e-041.96

Expected order: 2 (± 0.2).

Convergence of advection_2d_y_dirichlet_central (log-log)

central_D1y_periodic

op:D order:2 bc:periodic axes:x,y stencil width 3

Source: grids/cartesian_uniform_2d/rules/central_D1y_periodic.esm

Rewrites ∂f/∂y — match pattern:

{"op": "D", "args": ["f"], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorcentral_D1y_interior(f)[i,j] = (f[i, j+1] − f[i, j−1]) / (2·dy) for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face(f[i, 2] − f[i, NY]) / (2·dy) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face(f[i, 1] − f[i, NY−1]) / (2·dy) for i ∈ x, j ∈ [NY, NY]

Second-order centered first derivative along y with periodic (wrap) boundary on cartesian_uniform_2d, over the full 2-D field f[x,y]. Interior j in [2, NY-1] (all i): the imported y-stencil (f[i,j+1]-f[i,j-1])/(2 dy). The two y-face rows wrap periodically: j=1 uses the ghost f[i,0]=f[i,NY], giving (f[i,2]-f[i,NY])/(2 dy); j=NY uses f[i,NY+1]=f[i,1], giving (f[i,1]-f[i,NY-1])/(2 dy). The scheme is the same O(dy^2) centered difference at every row and conserves sum_y(f) exactly (the wrap telescopes). The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring grids/cartesian_uniform_2d/rules/central_D2y_neumann_bc so every binding’s simulator fills the (NX x 1) regions without rank promotion — the rank-reduced face form of the esm-spec 9.6.8 worked example only fills leading-axis faces portably. dy is the grid’s consumer-supplied y-spacing free name; the periodic wrap is dy-scale-invariant, so the rule serves any y-extent unchanged. Matches plain D(f, wrt: y) at default priority, so compound rules (e.g. a second-derivative rule at priority 10) fire first, match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the first y-derivative of a bare field declared over this grid’s [x, y] index sets; under import-edge renaming (esm-spec 9.7.7) the wrt literal and the where shape follow the renamed axes together, so two instances coexist scoped to their own meshes. O(dy^2). A model imports exactly one first-derivative rule for y (this centered periodic rule or the zero-gradient variant), never both.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains).

Conformance fixtures

Convergence — case advection_2d_y_periodic_central

Error norms read from the committed golden tests/conformance/convergence/advection_2d_y_periodic_central/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
161.578e-021.568e-02
323.964e-031.993.964e-031.98
649.924e-042.009.919e-042.00
1282.482e-042.002.482e-042.00

Expected order: 2 (± 0.2).

Convergence of advection_2d_y_periodic_central (log-log)

central_D1y_zero_grad_bc

op:D order:2 bc:zero_gradient axes:x,y stencil width 3

Source: grids/cartesian_uniform_2d/rules/central_D1y_zero_grad_bc.esm

Rewrites ∂f/∂y — match pattern:

{"op": "D", "args": ["f"], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorcentral_D1y_interior(f)[i,j] = (f[i, j+1] − f[i, j−1]) / (2·dy) for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face(−3·f[i, 1] + 4·f[i, 2] − f[i, 3]) / (2·dy) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face((f[i, NY−2] − 4·f[i, NY−1]) + 3·f[i, NY]) / (2·dy) for i ∈ x, j ∈ [NY, NY]

Second-order centered first derivative along y with one-sided (zero-gradient frame) boundary rows on cartesian_uniform_2d, over the full 2-D field f[x,y]. Interior j in [2, NY-1] (all i): the imported y-stencil (f[i,j+1]-f[i,j-1])/(2 dy). The two y-face rows use 2nd-order one-sided differences that reach only into the interior (no ghost, no periodic wrap): j=1 gives (-3 f[i,1]+4 f[i,2]-f[i,3])/(2 dy) and j=NY gives (f[i,NY-2]-4 f[i,NY-1]+3 f[i,NY])/(2 dy), so the discretized derivative is fully defined with its BC and stays O(dy^2) at every row (esm-spec 9.6.8). The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring grids/cartesian_uniform_2d/rules/central_D2y_neumann_bc so every binding’s simulator fills the (NX x 1) regions without rank promotion — the rank-reduced face form of the esm-spec 9.6.8 worked example only fills leading-axis faces portably. dy is the grid’s consumer-supplied y-spacing free name. The one-sided faces are O(dy^2) locally and preserve global 2nd-order convergence for boundary-compatible solutions — those whose odd y-derivatives vanish at the walls, so the discrete derivative matches the zero-gradient closure superconvergently (the MMS driver problems/advection_2d_y_zero_grad_central.esm uses cos(pi y), whose y-derivative is zero on both walls). Matches plain D(f, wrt: y) at default priority, so compound rules (e.g. a second-derivative rule at priority 10) fire first, match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the first y-derivative of a bare field declared over this grid’s [x, y] index sets; under import-edge renaming (esm-spec 9.7.7) the wrt literal and the where shape follow the renamed axes together. A model imports exactly one first-derivative rule for y (this zero-gradient rule or the periodic variant), never both.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (one-sided boundary treatment).
  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0

Conformance fixtures

Convergence — case advection_2d_y_zero_grad_central

Error norms read from the committed golden tests/conformance/convergence/advection_2d_y_zero_grad_central/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
166.276e-046.034e-04
321.527e-042.041.513e-042.00
643.799e-052.013.782e-052.00
1289.494e-062.009.458e-062.00

Expected order: 2 (± 0.2).

Convergence of advection_2d_y_zero_grad_central (log-log)

central_D2x_4_dirichlet_bc

op:D2 order:4 bc:dirichlet axes:x,y stencil width 5 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_4_dirichlet_bc.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [3, NX−2], y ∈ [1, NY]interiorcentral_D2x_4_interior(f)[i,j] = (−f[i−2, j] + 16·f[i−1, j] − 30·f[i, j] + 16·f[i+1, j] − f[i+2, j]) / (12·dx·dx) for i ∈ [3, NX−2], j ∈ y
x = 1, y ∈ [1, NY]boundary facecentral_D2x_4_dirichlet_face1(gL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j]) for i ∈ [1, 1], j ∈ y
x = 2, y ∈ [1, NY]boundary facecentral_D2x_4_dirichlet_face2(gL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j]) for i ∈ [2, 2], j ∈ y
x = NX−1, y ∈ [1, NY]boundary facecentral_D2x_4_dirichlet_face2(gR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j]) for i ∈ [NX−1, NX−1], j ∈ y
x = NX, y ∈ [1, NY]boundary facecentral_D2x_4_dirichlet_face1(gR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j]) for i ∈ [NX, NX], j ∈ y

The x-half of a fourth-order 2-D Laplacian on cartesian_uniform_2d with inhomogeneous DIRICHLET (fixed-value) conditions on the two x-walls, over the full 2-D field f[x,y], cell-centered. Interior i in [3, NX-2] (all j): the imported 5-point fourth-order x-stencil central_D2x_4_interior, (-f[i-2,j] + 16 f[i-1,j] - 30 f[i,j] + 16 f[i+1,j] - f[i+2,j]) / (12 dx^2), O(dx^4). The four near-x-wall columns i in {1, 2, NX-1, NX}, whose centered 5-point stencils would reach outside the domain, are closed with the imported one-sided fourth-order face templates central_D2x_4_dirichlet_face1 (nearest-wall column) and central_D2x_4_dirichlet_face2 (second column), each applied as a FULL-RANK aggregate over all j (singleton output range on x): the unique degree-5 polynomial through the prescribed x-wall value and the five nearest cells (at fixed j), differentiated twice at that cell. The wall values gL (left wall x=x0) and gR (right wall x=x0+NX*dx) are consumer-supplied free names bound into the face templates (the g parameter); same free-name contract as the grid geometry dx (grids/cartesian_uniform_2d/grid.esm). A consuming model declares gL/gR as ordinary real parameters, constant along the wall (all j); an un-overridden parameter falls back to its default (esm-spec 6: default, or 0), so gL=gR=0 is the homogeneous-Dirichlet case. gL/gR appear only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). The right x-wall reuses the same face weights by reflection (the second derivative is even under x -> -x), so f1..f5 there map to cells NX, NX-1, NX-2, NX-3, NX-4 (increasing distance from the right wall) and g binds to gR. Because BOTH the interior stencil and all four boundary closures are O(dx^4), the x-operator is uniformly fourth order with no boundary reduction; a steady linear-in-x profile carrying the wall values is discretely exact at every column (all closure weight-sums are zero). The full fourth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the first term. Compound match on D(D(f,x),x) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes (distinct from the 1-D shape-[x] rules); under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. Requires NX >= 5. A model imports exactly one second-derivative rule for x.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (one-sided high-order closures) and Chapter 3 (multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case heat_2d_x_dirichlet_o4

Error norms read from the committed golden tests/conformance/convergence/heat_2d_x_dirichlet_o4/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
89.756e-072.500e-06
166.667e-083.871.601e-073.96
324.225e-093.981.007e-083.99
642.647e-104.006.301e-104.00

Expected order: 4 (± 0.2).

Convergence of heat_2d_x_dirichlet_o4 (log-log)

central_D2x_4_neumann_bc

op:D2 order:4 bc:neumann axes:x,y stencil width 5 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_4_neumann_bc.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [3, NX−2], y ∈ [1, NY]interiorcentral_D2x_4_interior(f)[i,j] = (−f[i−2, j] + 16·f[i−1, j] − 30·f[i, j] + 16·f[i+1, j] − f[i+2, j]) / (12·dx·dx) for i ∈ [3, NX−2], j ∈ y
x = 1, y ∈ [1, NY]boundary facecentral_D2x_4_neumann_face1(qXL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j]) for i ∈ [1, 1], j ∈ y
x = 2, y ∈ [1, NY]boundary facecentral_D2x_4_neumann_face2(qXL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j]) for i ∈ [2, 2], j ∈ y
x = NX−1, y ∈ [1, NY]boundary facecentral_D2x_4_neumann_face2(−qXR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j]) for i ∈ [NX−1, NX−1], j ∈ y
x = NX, y ∈ [1, NY]boundary facecentral_D2x_4_neumann_face1(−qXR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j]) for i ∈ [NX, NX], j ∈ y

The x-half of a fourth-order 2-D Laplacian on cartesian_uniform_2d with inhomogeneous NEUMANN (fixed-flux) conditions on the two x-walls, over the full 2-D field f[x,y], cell-centered. Interior i in [3, NX-2] (all j): the imported 5-point fourth-order x-stencil central_D2x_4_interior, (-f[i-2,j] + 16 f[i-1,j] - 30 f[i,j] + 16 f[i+1,j] - f[i+2,j]) / (12 dx^2), O(dx^4). The four near-x-wall columns i in {1, 2, NX-1, NX}, whose centered 5-point stencils would reach outside the domain, are closed with the imported one-sided fourth-order face templates central_D2x_4_neumann_face1 (nearest-wall column) and central_D2x_4_neumann_face2 (second column), each applied as a FULL-RANK aggregate over all j (singleton output range on x): the unique degree-5 Hermite-type polynomial matching the prescribed x-wall derivative and the five nearest cells (at fixed j), differentiated twice at that cell. The wall gradients qXL (df/dx at the left wall x=x0, +x direction) and qXR (df/dx at the right wall x=x0+NX*dx, +x direction) are consumer-supplied free names; same free-name contract as the grid geometry dx (grids/cartesian_uniform_2d/grid.esm). A consuming model declares qXL/qXR as ordinary real parameters, constant along the wall (all j); an un-overridden parameter falls back to its default (esm-spec 6: default, or 0), so qXL=qXR=0 is the homogeneous (zero-flux) case and the rule reduces exactly to a zero-gradient fourth-order closure. qXL/qXR appear only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). The left face templates take the wall derivative datum q as their +x-into-domain slope, so the left columns bind q -> qXL. The right x-wall reuses the SAME face weights by reflection (the second derivative is even under x -> -x while df/dx is odd), so the right columns bind q -> -qXR and map f1..f5 to cells NX, NX-1, NX-2, NX-3, NX-4 (increasing distance from the right wall). Because BOTH the interior stencil and all four boundary closures are O(dx^4), the x-operator is uniformly fourth order with no boundary reduction; a steady ramp-in-x carrying the wall fluxes is discretely exact at every column (the value-weight sums are zero and the derivative term supplies the exact wall slope). The full fourth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the first term and central_D2y_4_neumann_bc the second, and the two per-axis makearrays superpose to the correct 9-point interior and one-sided edges/corners. Compound match on D(D(f,x),x) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. Requires NX >= 5. A model imports exactly one second-derivative rule for x.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (Neumann one-sided high-order closures) and Chapter 3 (multidimensional heat by dimension-wise second differences).

Conformance fixtures

Convergence — case heat_2d_neumann_flux_o4

Error norms read from the committed golden tests/conformance/convergence/heat_2d_neumann_flux_o4/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
84.366e-061.413e-05
162.077e-074.391.006e-063.81
329.348e-094.475.036e-084.32
645.271e-104.151.688e-094.90

Expected order: 4 (± 0.5).

Convergence of heat_2d_neumann_flux_o4 (log-log)

central_D2x_4_periodic

op:D2 order:4 bc:periodic axes:x,y stencil width 5 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_4_periodic.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [3, NX−2], y ∈ [1, NY]interiorcentral_D2x_4_interior(f)[i,j] = (−f[i−2, j] + 16·f[i−1, j] − 30·f[i, j] + 16·f[i+1, j] − f[i+2, j]) / (12·dx·dx) for i ∈ [3, NX−2], j ∈ y
x = 1, y ∈ [1, NY]boundary face(−f[NX−1, j] + 16·f[NX, j] − 30·f[1, j] + 16·f[2, j] − f[3, j]) / (12·dx·dx) for i ∈ [1, 1], j ∈ y
x = 2, y ∈ [1, NY]boundary face(−f[NX, j] + 16·f[1, j] − 30·f[2, j] + 16·f[3, j] − f[4, j]) / (12·dx·dx) for i ∈ [2, 2], j ∈ y
x = NX−1, y ∈ [1, NY]boundary face(−f[NX−3, j] + 16·f[NX−2, j] − 30·f[NX−1, j] + 16·f[NX, j] − f[1, j]) / (12·dx·dx) for i ∈ [NX−1, NX−1], j ∈ y
x = NX, y ∈ [1, NY]boundary face(−f[NX−2, j] + 16·f[NX−1, j] − 30·f[NX, j] + 16·f[1, j] − f[2, j]) / (12·dx·dx) for i ∈ [NX, NX], j ∈ y

The x-half of a fourth-order 2-D Laplacian on cartesian_uniform_2d with a PERIODIC x-boundary: fourth-order centered second derivative along x over the full 2-D field f[x,y], the x direction wrapped as a circle. Interior i in [3, NX-2] (all j): the imported 5-point x-stencil central_D2x_4_interior, (-f[i-2,j] + 16 f[i-1,j] - 30 f[i,j] + 16 f[i+1,j] - f[i+2,j]) / (12 dx^2). The four x-face columns i in {1, 2, NX-1, NX} carry the SAME 5-point weights with the out-of-range x-neighbors wrapped periodically (i-2, i-1 wrap past 1 to NX-1, NX; i+1, i+2 wrap past NX to 1, 2), so the whole x-operator is the exact circulant fourth-order Laplacian along x — O(dx^4) with no boundary order reduction, the 2-D per-axis analogue of grids/cartesian_uniform_1d/rules/central_D2_4_periodic. The four face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring central_D2x_periodic so every binding’s simulator fills the (1 x NY) columns without rank promotion. dx is the grid’s consumer-supplied x-spacing free name; the periodic wrap is dx-scale-invariant. Requires NX >= 5 for the 5-point stencil to have distinct points. The full fourth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the first term. Compound match on D(D(f,x),x) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes (distinct from the 1-D shape-[x] rules); under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. A model imports exactly one second-derivative rule for x.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains).

Conformance fixtures

Convergence — case heat_2d_x_periodic_o4

Error norms read from the committed golden tests/conformance/convergence/heat_2d_x_periodic_o4/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
81.580e-041.403e-04
161.029e-053.949.701e-063.85
326.497e-073.996.216e-073.96
644.071e-084.003.909e-083.99

Expected order: 4 (± 0.2).

Convergence of heat_2d_x_periodic_o4 (log-log)

central_D2x_6_dirichlet_bc

op:D2 order:6 bc:dirichlet axes:x,y stencil width 7 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_6_dirichlet_bc.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [4, NX−3], y ∈ [1, NY]interiorcentral_D2x_6_interior(f)[i,j] = (2·f[i−3, j] − 27·f[i−2, j] + 270·f[i−1, j] − 490·f[i, j] + 270·f[i+1, j] − 27·f[i+2, j] + 2·f[i+3, j]) / (180·dx·dx) for i ∈ [4, NX−3], j ∈ y
x = 1, y ∈ [1, NY]boundary facecentral_D2x_6_dirichlet_face1(gL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j], f[6, j], f[7, j]) for i ∈ [1, 1], j ∈ y
x = 2, y ∈ [1, NY]boundary facecentral_D2x_6_dirichlet_face2(gL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j], f[6, j], f[7, j]) for i ∈ [2, 2], j ∈ y
x = 3, y ∈ [1, NY]boundary facecentral_D2x_6_dirichlet_face3(gL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j], f[6, j], f[7, j]) for i ∈ [3, 3], j ∈ y
x = NX−2, y ∈ [1, NY]boundary facecentral_D2x_6_dirichlet_face3(gR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j], f[NX−5, j], f[NX−6, j]) for i ∈ [NX−2, NX−2], j ∈ y
x = NX−1, y ∈ [1, NY]boundary facecentral_D2x_6_dirichlet_face2(gR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j], f[NX−5, j], f[NX−6, j]) for i ∈ [NX−1, NX−1], j ∈ y
x = NX, y ∈ [1, NY]boundary facecentral_D2x_6_dirichlet_face1(gR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j], f[NX−5, j], f[NX−6, j]) for i ∈ [NX, NX], j ∈ y

The x-half of a sixth-order 2-D Laplacian on cartesian_uniform_2d with inhomogeneous DIRICHLET conditions on the two x-walls, over the full 2-D field f[x,y], cell-centered. Interior x in [4, NX-3] (all y): the imported 7-point sixth-order x-stencil central_D2x_6_interior, O(dx^6). The six near-x-wall lines x in {1,2,3,NX-2,NX-1,NX}, whose centered 7-point stencils would reach outside the domain, are closed with the imported one-sided sixth-order face templates central_D2x_6_dirichlet_face1/face2/face3 (nearest, second, third column), each applied as a FULL-RANK aggregate over all y (singleton output range on x): the unique degree-7 polynomial matching the prescribed wall datum and the seven nearest cells (at fixed y), differentiated twice at that cell. The wall values gL and gR are consumer-supplied free names bound into the face templates (the g parameter); an un-overridden parameter defaults to 0, so gL=gR=0 is the homogeneous-Dirichlet case. The far wall reuses the same face weights by reflection (the second derivative is even), so f1..f7 there map to cells NX, NX-1, …, NX-6. The wall free names appear only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). Because BOTH the interior stencil and all six boundary closures are O(dx^6), the x-operator is uniformly sixth order with no boundary reduction; a steady linear-in-x profile carrying the wall values is discretely exact at every line (all closure weight-sums are zero). The full sixth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the x-term. Compound match on D(D(f,x),x) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. Requires NX >= 7. A model imports exactly one second-derivative rule for x.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (one-sided high-order closures) and Chapter 3 (multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case heat_2d_x_dirichlet_o6

Error norms read from the committed golden tests/conformance/convergence/heat_2d_x_dirichlet_o6/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
87.259e-081.873e-07
164.874e-107.229.892e-107.56
326.520e-126.221.558e-115.99
641.422e-135.521.044e-123.90

Expected order: 6 (± 0.6).

Convergence of heat_2d_x_dirichlet_o6 (log-log)

central_D2x_6_neumann_bc

op:D2 order:6 bc:neumann axes:x,y stencil width 7 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_6_neumann_bc.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [4, NX−3], y ∈ [1, NY]interiorcentral_D2x_6_interior(f)[i,j] = (2·f[i−3, j] − 27·f[i−2, j] + 270·f[i−1, j] − 490·f[i, j] + 270·f[i+1, j] − 27·f[i+2, j] + 2·f[i+3, j]) / (180·dx·dx) for i ∈ [4, NX−3], j ∈ y
x = 1, y ∈ [1, NY]boundary facecentral_D2x_6_neumann_face1(qXL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j], f[6, j], f[7, j]) for i ∈ [1, 1], j ∈ y
x = 2, y ∈ [1, NY]boundary facecentral_D2x_6_neumann_face2(qXL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j], f[6, j], f[7, j]) for i ∈ [2, 2], j ∈ y
x = 3, y ∈ [1, NY]boundary facecentral_D2x_6_neumann_face3(qXL, f[1, j], f[2, j], f[3, j], f[4, j], f[5, j], f[6, j], f[7, j]) for i ∈ [3, 3], j ∈ y
x = NX−2, y ∈ [1, NY]boundary facecentral_D2x_6_neumann_face3(−qXR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j], f[NX−5, j], f[NX−6, j]) for i ∈ [NX−2, NX−2], j ∈ y
x = NX−1, y ∈ [1, NY]boundary facecentral_D2x_6_neumann_face2(−qXR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j], f[NX−5, j], f[NX−6, j]) for i ∈ [NX−1, NX−1], j ∈ y
x = NX, y ∈ [1, NY]boundary facecentral_D2x_6_neumann_face1(−qXR, f[NX, j], f[NX−1, j], f[NX−2, j], f[NX−3, j], f[NX−4, j], f[NX−5, j], f[NX−6, j]) for i ∈ [NX, NX], j ∈ y

The x-half of a sixth-order 2-D Laplacian on cartesian_uniform_2d with inhomogeneous NEUMANN conditions on the two x-walls, over the full 2-D field f[x,y], cell-centered. Interior x in [4, NX-3] (all y): the imported 7-point sixth-order x-stencil central_D2x_6_interior, O(dx^6). The six near-x-wall lines x in {1,2,3,NX-2,NX-1,NX}, whose centered 7-point stencils would reach outside the domain, are closed with the imported one-sided sixth-order face templates central_D2x_6_neumann_face1/face2/face3 (nearest, second, third column), each applied as a FULL-RANK aggregate over all y (singleton output range on x): the unique degree-7 polynomial matching the prescribed wall datum and the seven nearest cells (at fixed y), differentiated twice at that cell. The wall gradients qXL (df/dx at the low wall) and qXR (df/dx at the high wall) are consumer-supplied free names; an un-overridden parameter defaults to 0, so qXL=qXR=0 is the homogeneous (zero-flux) case and the rule reduces exactly to a zero-gradient sixth-order closure. The low face templates take q as the +x-into-domain slope, so the low rows bind q -> qXL. The far wall reuses the SAME face weights by reflection (the second derivative even, df/dx odd), so the far rows bind q -> -qXR and map f1..f7 to cells NX, NX-1, …, NX-6. The wall free names appear only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). Because BOTH the interior stencil and all six boundary closures are O(dx^6), the x-operator is uniformly sixth order with no boundary reduction; a steady ramp-in-x carrying the wall fluxes is discretely exact at every line (all closure value-weight sums are zero). The full sixth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the x-term. Compound match on D(D(f,x),x) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. Requires NX >= 7. A model imports exactly one second-derivative rule for x.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (one-sided high-order closures) and Chapter 3 (multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case heat_2d_neumann_flux_o6

Error norms read from the committed golden tests/conformance/convergence/heat_2d_neumann_flux_o6/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
83.229e-071.019e-06
165.963e-095.762.827e-085.17
326.581e-116.504.193e-106.08
646.494e-136.664.322e-126.60

Expected order: 6 (± 0.6).

Convergence of heat_2d_neumann_flux_o6 (log-log)

central_D2x_6_periodic

op:D2 order:6 bc:periodic axes:x,y stencil width 7 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_6_periodic.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [4, NX−3], y ∈ [1, NY]interiorcentral_D2x_6_interior(f)[i,j] = (2·f[i−3, j] − 27·f[i−2, j] + 270·f[i−1, j] − 490·f[i, j] + 270·f[i+1, j] − 27·f[i+2, j] + 2·f[i+3, j]) / (180·dx·dx) for i ∈ [4, NX−3], j ∈ y
x = 1, y ∈ [1, NY]boundary face(2·f[NX−2, j] − 27·f[NX−1, j] + 270·f[NX, j] − 490·f[1, j] + 270·f[2, j] − 27·f[3, j] + 2·f[4, j]) / (180·dx·dx) for i ∈ [1, 1], j ∈ y
x = 2, y ∈ [1, NY]boundary face(2·f[NX−1, j] − 27·f[NX, j] + 270·f[1, j] − 490·f[2, j] + 270·f[3, j] − 27·f[4, j] + 2·f[5, j]) / (180·dx·dx) for i ∈ [2, 2], j ∈ y
x = 3, y ∈ [1, NY]boundary face(2·f[NX, j] − 27·f[1, j] + 270·f[2, j] − 490·f[3, j] + 270·f[4, j] − 27·f[5, j] + 2·f[6, j]) / (180·dx·dx) for i ∈ [3, 3], j ∈ y
x = NX−2, y ∈ [1, NY]boundary face(2·f[NX−5, j] − 27·f[NX−4, j] + 270·f[NX−3, j] − 490·f[NX−2, j] + 270·f[NX−1, j] − 27·f[NX, j] + 2·f[1, j]) / (180·dx·dx) for i ∈ [NX−2, NX−2], j ∈ y
x = NX−1, y ∈ [1, NY]boundary face(2·f[NX−4, j] − 27·f[NX−3, j] + 270·f[NX−2, j] − 490·f[NX−1, j] + 270·f[NX, j] − 27·f[1, j] + 2·f[2, j]) / (180·dx·dx) for i ∈ [NX−1, NX−1], j ∈ y
x = NX, y ∈ [1, NY]boundary face(2·f[NX−3, j] − 27·f[NX−2, j] + 270·f[NX−1, j] − 490·f[NX, j] + 270·f[1, j] − 27·f[2, j] + 2·f[3, j]) / (180·dx·dx) for i ∈ [NX, NX], j ∈ y

The x-half of a sixth-order 2-D Laplacian on cartesian_uniform_2d with a PERIODIC x-boundary: sixth-order centered second derivative along x over the full 2-D field f[x,y], the x direction wrapped as a circle. Interior x in [4, NX-3] (all y): the imported 7-point x-stencil central_D2x_6_interior, (2 f[..-3] - 27 f[..-2] + 270 f[..-1] - 490 f[..] + 270 f[..+1] - 27 f[..+2] + 2 f[..+3]) / (180 dx^2). The six x-face lines x in {1,2,3,NX-2,NX-1,NX} carry the SAME 7-point weights with out-of-range x-neighbors wrapped periodically, so the whole x-operator is the exact circulant sixth-order Laplacian along x — O(dx^6) with no boundary order reduction, the 2-D per-axis analogue of grids/cartesian_uniform_1d/rules/central_D2_6_periodic. The six face values are FULL-RANK aggregates (singleton output range on the x axis, all y) so every binding’s simulator fills the lines without rank promotion. Requires NX >= 7 for the 7-point stencil to have distinct points. dx is the grid’s consumer-supplied x-spacing free name; the periodic wrap is dx-scale-invariant. The full sixth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the x-term. Compound match on D(D(f,x),x) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. A model imports exactly one second-derivative rule for x.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains).

Conformance fixtures

Convergence — case heat_2d_x_periodic_o6

Error norms read from the committed golden tests/conformance/convergence/heat_2d_x_periodic_o6/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
81.503e-051.335e-05
162.524e-075.902.380e-075.81
324.016e-095.973.842e-095.95
646.303e-115.996.052e-115.99

Expected order: 6 (± 0.3).

Convergence of heat_2d_x_periodic_o6 (log-log)

central_D2x_dirichlet_bc

op:D2 order:2 bc:dirichlet axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_dirichlet_bc.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorcentral_D2x_interior(f)[i,j] = (f[i+1, j] − 2·f[i, j] + f[i−1, j]) / dx² for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face((f[2, j] − 3·f[1, j]) + 2·gL) / dx² for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face((f[NX−1, j] − 3·f[NX, j]) + 2·gR) / dx² for i ∈ [NX, NX], j ∈ y

The x-half of a 2-D Laplacian on cartesian_uniform_2d with inhomogeneous DIRICHLET (fixed-value) conditions on the two x-walls, over the full 2-D field f[x,y], cell-centered. Interior i in [2, NX-1] (all j): the imported x-stencil (f[i+1,j] - 2 f[i,j] + f[i-1,j]) / dx^2. Faces i=1 and i=NX: mirror-ghost elimination against a prescribed wall value. gL and gR are consumer-supplied free names (the same free-name contract as the grid geometry dx; see grids/cartesian_uniform_2d/grid.esm) giving the field value on the left x-wall x=x0 and the right x-wall x=x0+NXdx, constant along that wall (all j). On this cell-centered grid the wall sits half a cell outside the first center, so linear reconstruction to the wall gives the ghost f[0,j]=2gL-f[1,j] (and f[NX+1,j]=2gR-f[NX,j]); substituting into the 3-point stencil yields D2x[1,j] = (f[2,j]-3f[1,j]+2gL)/dx^2 and D2x[NX,j] = (f[NX-1,j]-3f[NX,j]+2*gR)/dx^2. A consuming model declares gL/gR as ordinary real parameters; a parameter not overridden falls back to its default (esm-spec §6: default, or 0), so gL=gR=0 is the homogeneous-Dirichlet x-wall. The gL/gR terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). dx is the grid’s consumer-supplied x-spacing free name. The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring central_D2x_neumann_bc. The full Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the first term. Compound match on D(D(f,x),x) at priority 10 (esm-spec §9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec §9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes (distinct from the 1-D shape-[x] rules); under import-edge renaming (esm-spec §9.7.7) the wrt literals and the where shape follow the renamed axes together. Interior O(dx^2); the one-sided mirror-ghost faces preserve global 2nd-order convergence for Dirichlet-compatible solutions.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (boundary conditions; cell-centered Dirichlet via mirror ghost) and Chapter 3 (multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case advection_2d_x_dirichlet_central

Error norms read from the committed golden tests/conformance/convergence/advection_2d_x_dirichlet_central/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
165.690e-046.559e-03
321.442e-041.981.829e-031.84
643.617e-052.004.809e-041.93
1289.050e-062.001.232e-041.96

Expected order: 2 (± 0.2).

Convergence of advection_2d_x_dirichlet_central (log-log)

Convergence — case anisotropic_diffusion_2d_dirichlet

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_dirichlet/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.657e-032.377e-03
326.589e-042.015.932e-042.00
641.644e-042.001.481e-042.00
1284.107e-052.003.703e-052.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_dirichlet (log-log)

Convergence — case heat_2d_x_dirichlet

Error norms read from the committed golden tests/conformance/convergence/heat_2d_x_dirichlet/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
163.167e-053.121e-05
327.925e-062.007.837e-061.99
641.982e-062.001.962e-062.00
1284.954e-072.004.905e-072.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_x_dirichlet (log-log)

central_D2x_neumann_bc

op:D2 order:2 bc:neumann axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_neumann_bc.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorcentral_D2x_interior(f)[i,j] = (f[i+1, j] − 2·f[i, j] + f[i−1, j]) / dx² for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face(f[i+1, j] − f[i, j] − qXL·dx) / dx² for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face((f[i−1, j] − f[i, j]) + qXR·dx) / dx² for i ∈ [NX, NX], j ∈ y

The x-half of a 2-D Laplacian on cartesian_uniform_2d: second-order centered second derivative along x with inhomogeneous Neumann (fixed-flux) conditions on the two x-walls, over the full 2-D field f[x,y]. Interior i in [2, NX-1] (all j): the imported x-stencil. x-left face i=1: ghost elimination (f[i,j]-f[0,j])/dx = qXL gives D2x[1,j] = (f[2,j]-f[1,j]-qXLdx)/dx^2. x-right face i=NX: (f[NX+1,j]-f[NX,j])/dx = qXR gives D2x[NX,j] = (f[NX-1,j]-f[NX,j]+qXRdx)/dx^2. qXL/qXR are consumer-supplied free names (df/dx at the left/right walls, +x direction) declared as ordinary real parameters; an un-overridden parameter falls back to its default (esm-spec §6: default, or 0), so qXL=qXR=0 is the homogeneous (zero-gradient) x-wall. dx is the grid’s consumer-supplied spacing free name. The full Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the first term and central_D2y_neumann_bc the second, and the two per-axis makearrays superpose to the correct 5-point interior and one-sided edges/corners. Compound match on D(D(f,x),x) at priority 10 (esm-spec §9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec §9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes (distinct from the 1-D shape-[x] rules); under import-edge renaming (esm-spec §9.7.7) the wrt literals and the where shape follow the renamed axes together. Interior O(dx^2); the one-sided faces preserve global 2nd-order convergence for Neumann-compatible solutions.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 and Chapter 3 (Neumann conditions; multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case anisotropic_diffusion_2d_neumann

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_neumann/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
169.382e-036.067e-03
322.347e-032.001.520e-032.00
645.868e-042.003.802e-042.00
1281.467e-042.009.507e-052.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_neumann (log-log)

Convergence — case heat_2d_neumann_flux

Error norms read from the committed golden tests/conformance/convergence/heat_2d_neumann_flux/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
81.300e-042.380e-04
163.261e-052.006.150e-051.95
328.159e-062.001.550e-051.99
642.040e-062.003.883e-062.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_neumann_flux (log-log)

central_D2x_periodic

op:D2 order:2 bc:periodic axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2x_periodic.esm

Rewrites ∂²f/∂x² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorcentral_D2x_interior(f)[i,j] = (f[i+1, j] − 2·f[i, j] + f[i−1, j]) / dx² for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face(f[2, j] − 2·f[1, j] + f[NX, j]) / dx² for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face(f[1, j] − 2·f[NX, j] + f[NX−1, j]) / dx² for i ∈ [NX, NX], j ∈ y

The x-half of a 2-D Laplacian on cartesian_uniform_2d with a PERIODIC x-boundary: second-order centered second derivative along x over the full 2-D field f[x,y], the x direction wrapped as a circle. Interior i in [2, NX-1] (all j): the imported x-stencil (f[i+1,j] - 2 f[i,j] + f[i-1,j]) / dx^2. The two x-face columns wrap periodically: i=1 uses the ghost f[0,j]=f[NX,j], giving (f[2,j] - 2 f[1,j] + f[NX,j]) / dx^2; i=NX uses f[NX+1,j]=f[1,j], giving (f[1,j] - 2 f[NX,j] + f[NX-1,j]) / dx^2. The scheme is the same O(dx^2) centered second difference at every column, so the whole x-operator is the exact circulant symmetric Laplacian along x — no boundary order reduction, the natural second-order analogue of central_D1x_periodic. The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring central_D2x_neumann_bc so every binding’s simulator fills the (1 x NY) columns without rank promotion. dx is the grid’s consumer-supplied x-spacing free name; the periodic wrap is dx-scale-invariant, so the rule serves any x-extent unchanged. The full Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the first term. Compound match on D(D(f,x),x) at priority 10 (esm-spec §9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec §9.6.1) so it fires only on the xx-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes (distinct from the 1-D shape-[x] rules); under import-edge renaming (esm-spec §9.7.7) the wrt literals and the where shape follow the renamed axes together, so two instances coexist scoped to their own meshes. O(dx^2). A model imports exactly one second-derivative rule for x (this periodic rule or a Neumann/Dirichlet variant), never more than one.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains) and Section 2.12 (second differences).

Conformance fixtures

Convergence — case anisotropic_diffusion_2d_periodic

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_periodic/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.489e-022.396e-02
326.326e-031.986.266e-031.94
641.589e-031.991.585e-031.98
1283.976e-042.003.973e-042.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_periodic (log-log)

Convergence — case heat_2d_x_periodic

Error norms read from the committed golden tests/conformance/convergence/heat_2d_x_periodic/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
165.049e-044.760e-04
321.267e-041.991.212e-041.97
643.170e-052.003.044e-051.99
1287.927e-062.007.617e-062.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_x_periodic (log-log)

central_D2y_4_dirichlet_bc

op:D2 order:4 bc:dirichlet axes:x,y stencil width 5 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_4_dirichlet_bc.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [3, NY−2]interiorcentral_D2y_4_interior(f)[i,j] = (−f[i, j−2] + 16·f[i, j−1] − 30·f[i, j] + 16·f[i, j+1] − f[i, j+2]) / (12·dy·dy) for i ∈ x, j ∈ [3, NY−2]
x ∈ [1, NX], y = 1boundary facecentral_D2y_4_dirichlet_face1(gB, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5]) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = 2boundary facecentral_D2y_4_dirichlet_face2(gB, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5]) for i ∈ x, j ∈ [2, 2]
x ∈ [1, NX], y = NY−1boundary facecentral_D2y_4_dirichlet_face2(gT, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4]) for i ∈ x, j ∈ [NY−1, NY−1]
x ∈ [1, NX], y = NYboundary facecentral_D2y_4_dirichlet_face1(gT, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4]) for i ∈ x, j ∈ [NY, NY]

The y-half of a fourth-order 2-D Laplacian on cartesian_uniform_2d with inhomogeneous DIRICHLET (fixed-value) conditions on the two y-walls, over the full 2-D field f[x,y], cell-centered. Interior j in [3, NY-2] (all i): the imported 5-point fourth-order y-stencil central_D2y_4_interior, (-f[i,j-2] + 16 f[i,j-1] - 30 f[i,j] + 16 f[i,j+1] - f[i,j+2]) / (12 dy^2), O(dy^4). The four near-y-wall rows j in {1, 2, NY-1, NY}, whose centered 5-point stencils would reach outside the domain, are closed with the imported one-sided fourth-order face templates central_D2y_4_dirichlet_face1 (nearest-wall row) and central_D2y_4_dirichlet_face2 (second row), each applied as a FULL-RANK aggregate over all i (singleton output range on y): the unique degree-5 polynomial through the prescribed y-wall value and the five nearest cells (at fixed i), differentiated twice at that cell. The wall values gB (bottom wall y=y0) and gT (top wall y=y0+NY*dy) are consumer-supplied free names bound into the face templates (the g parameter); same free-name contract as the grid geometry dy (grids/cartesian_uniform_2d/grid.esm). A consuming model declares gB/gT as ordinary real parameters, constant along the wall (all i); an un-overridden parameter falls back to its default (esm-spec 6: default, or 0), so gB=gT=0 is the homogeneous-Dirichlet case. gB/gT appear only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). The top y-wall reuses the same face weights by reflection (the second derivative is even under y -> -y), so f1..f5 there map to cells NY, NY-1, NY-2, NY-3, NY-4 (increasing distance from the top wall) and g binds to gT. Because BOTH the interior stencil and all four boundary closures are O(dy^4), the y-operator is uniformly fourth order with no boundary reduction; a steady linear-in-y profile carrying the wall values is discretely exact at every row (all closure weight-sums are zero). The full fourth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the second term. Compound match on D(D(f,y),y) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. Requires NY >= 5. A model imports exactly one second-derivative rule for y.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (one-sided high-order closures) and Chapter 3 (multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case heat_2d_y_dirichlet_o4

Error norms read from the committed golden tests/conformance/convergence/heat_2d_y_dirichlet_o4/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
89.756e-072.500e-06
166.667e-083.871.601e-073.96
324.225e-093.981.007e-083.99
642.647e-104.006.301e-104.00

Expected order: 4 (± 0.2).

Convergence of heat_2d_y_dirichlet_o4 (log-log)

central_D2y_4_neumann_bc

op:D2 order:4 bc:neumann axes:x,y stencil width 5 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_4_neumann_bc.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [3, NY−2]interiorcentral_D2y_4_interior(f)[i,j] = (−f[i, j−2] + 16·f[i, j−1] − 30·f[i, j] + 16·f[i, j+1] − f[i, j+2]) / (12·dy·dy) for i ∈ x, j ∈ [3, NY−2]
x ∈ [1, NX], y = 1boundary facecentral_D2y_4_neumann_face1(qYL, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5]) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = 2boundary facecentral_D2y_4_neumann_face2(qYL, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5]) for i ∈ x, j ∈ [2, 2]
x ∈ [1, NX], y = NY−1boundary facecentral_D2y_4_neumann_face2(−qYR, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4]) for i ∈ x, j ∈ [NY−1, NY−1]
x ∈ [1, NX], y = NYboundary facecentral_D2y_4_neumann_face1(−qYR, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4]) for i ∈ x, j ∈ [NY, NY]

The y-half of a fourth-order 2-D Laplacian on cartesian_uniform_2d with inhomogeneous NEUMANN (fixed-flux) conditions on the two y-walls, over the full 2-D field f[x,y], cell-centered. Interior j in [3, NY-2] (all i): the imported 5-point fourth-order y-stencil central_D2y_4_interior, (-f[i,j-2] + 16 f[i,j-1] - 30 f[i,j] + 16 f[i,j+1] - f[i,j+2]) / (12 dy^2), O(dy^4). The four near-y-wall rows j in {1, 2, NY-1, NY}, whose centered 5-point stencils would reach outside the domain, are closed with the imported one-sided fourth-order face templates central_D2y_4_neumann_face1 (nearest-wall row) and central_D2y_4_neumann_face2 (second row), each applied as a FULL-RANK aggregate over all i (singleton output range on y): the unique degree-5 Hermite-type polynomial matching the prescribed y-wall derivative and the five nearest cells (at fixed i), differentiated twice at that cell. The wall gradients qYL (df/dy at the bottom wall y=y0, +y direction) and qYR (df/dy at the top wall y=y0+NY*dy, +y direction) are consumer-supplied free names; same free-name contract as the grid geometry dy (grids/cartesian_uniform_2d/grid.esm). A consuming model declares qYL/qYR as ordinary real parameters, constant along the wall (all i); an un-overridden parameter falls back to its default (esm-spec 6: default, or 0), so qYL=qYR=0 is the homogeneous (zero-flux) case and the rule reduces exactly to a zero-gradient fourth-order closure. qYL/qYR appear only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). The bottom face templates take the wall derivative datum q as their +y-into-domain slope, so the bottom rows bind q -> qYL. The top y-wall reuses the SAME face weights by reflection (the second derivative is even under y -> -y while df/dy is odd), so the top rows bind q -> -qYR and map f1..f5 to cells NY, NY-1, NY-2, NY-3, NY-4 (increasing distance from the top wall). Because BOTH the interior stencil and all four boundary closures are O(dy^4), the y-operator is uniformly fourth order with no boundary reduction; a steady ramp-in-y carrying the wall fluxes is discretely exact at every row (the value-weight sums are zero and the derivative term supplies the exact wall slope). The full fourth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): central_D2x_4_neumann_bc lowers the first term and this rule the second, and the two per-axis makearrays superpose to the correct 9-point interior and one-sided edges/corners. Compound match on D(D(f,y),y) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. Requires NY >= 5. A model imports exactly one second-derivative rule for y.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (Neumann one-sided high-order closures) and Chapter 3 (multidimensional heat by dimension-wise second differences).

Conformance fixtures

Convergence — case heat_2d_neumann_flux_o4

Error norms read from the committed golden tests/conformance/convergence/heat_2d_neumann_flux_o4/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
84.366e-061.413e-05
162.077e-074.391.006e-063.81
329.348e-094.475.036e-084.32
645.271e-104.151.688e-094.90

Expected order: 4 (± 0.5).

Convergence of heat_2d_neumann_flux_o4 (log-log)

central_D2y_4_periodic

op:D2 order:4 bc:periodic axes:x,y stencil width 5 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_4_periodic.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [3, NY−2]interiorcentral_D2y_4_interior(f)[i,j] = (−f[i, j−2] + 16·f[i, j−1] − 30·f[i, j] + 16·f[i, j+1] − f[i, j+2]) / (12·dy·dy) for i ∈ x, j ∈ [3, NY−2]
x ∈ [1, NX], y = 1boundary face(−f[i, NY−1] + 16·f[i, NY] − 30·f[i, 1] + 16·f[i, 2] − f[i, 3]) / (12·dy·dy) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = 2boundary face(−f[i, NY] + 16·f[i, 1] − 30·f[i, 2] + 16·f[i, 3] − f[i, 4]) / (12·dy·dy) for i ∈ x, j ∈ [2, 2]
x ∈ [1, NX], y = NY−1boundary face(−f[i, NY−3] + 16·f[i, NY−2] − 30·f[i, NY−1] + 16·f[i, NY] − f[i, 1]) / (12·dy·dy) for i ∈ x, j ∈ [NY−1, NY−1]
x ∈ [1, NX], y = NYboundary face(−f[i, NY−2] + 16·f[i, NY−1] − 30·f[i, NY] + 16·f[i, 1] − f[i, 2]) / (12·dy·dy) for i ∈ x, j ∈ [NY, NY]

The y-half of a fourth-order 2-D Laplacian on cartesian_uniform_2d with a PERIODIC y-boundary: fourth-order centered second derivative along y over the full 2-D field f[x,y], the y direction wrapped as a circle. Interior j in [3, NY-2] (all i): the imported 5-point y-stencil central_D2y_4_interior, (-f[i,j-2] + 16 f[i,j-1] - 30 f[i,j] + 16 f[i,j+1] - f[i,j+2]) / (12 dy^2). The four y-face rows j in {1, 2, NY-1, NY} carry the SAME 5-point weights with the out-of-range y-neighbors wrapped periodically (j-2, j-1 wrap past 1 to NY-1, NY; j+1, j+2 wrap past NY to 1, 2), so the whole y-operator is the exact circulant fourth-order Laplacian along y — O(dy^4) with no boundary order reduction, the 2-D per-axis analogue of grids/cartesian_uniform_1d/rules/central_D2_4_periodic. The four face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring central_D2y_periodic so every binding’s simulator fills the (NX x 1) rows without rank promotion. dy is the grid’s consumer-supplied y-spacing free name; the periodic wrap is dy-scale-invariant. Requires NY >= 5 for the 5-point stencil to have distinct points. The full fourth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the second term. Compound match on D(D(f,y),y) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. A model imports exactly one second-derivative rule for y.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains).

Conformance fixtures

Convergence — case heat_2d_y_periodic_o4

Error norms read from the committed golden tests/conformance/convergence/heat_2d_y_periodic_o4/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
81.580e-041.403e-04
161.029e-053.949.701e-063.85
326.497e-073.996.216e-073.96
644.071e-084.003.909e-083.99

Expected order: 4 (± 0.2).

Convergence of heat_2d_y_periodic_o4 (log-log)

central_D2y_6_dirichlet_bc

op:D2 order:6 bc:dirichlet axes:x,y stencil width 7 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_6_dirichlet_bc.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [4, NY−3]interiorcentral_D2y_6_interior(f)[i,j] = (2·f[i, j−3] − 27·f[i, j−2] + 270·f[i, j−1] − 490·f[i, j] + 270·f[i, j+1] − 27·f[i, j+2] + 2·f[i, j+3]) / (180·dy·dy) for i ∈ x, j ∈ [4, NY−3]
x ∈ [1, NX], y = 1boundary facecentral_D2y_6_dirichlet_face1(gB, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5], f[i, 6], f[i, 7]) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = 2boundary facecentral_D2y_6_dirichlet_face2(gB, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5], f[i, 6], f[i, 7]) for i ∈ x, j ∈ [2, 2]
x ∈ [1, NX], y = 3boundary facecentral_D2y_6_dirichlet_face3(gB, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5], f[i, 6], f[i, 7]) for i ∈ x, j ∈ [3, 3]
x ∈ [1, NX], y = NY−2boundary facecentral_D2y_6_dirichlet_face3(gT, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4], f[i, NY−5], f[i, NY−6]) for i ∈ x, j ∈ [NY−2, NY−2]
x ∈ [1, NX], y = NY−1boundary facecentral_D2y_6_dirichlet_face2(gT, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4], f[i, NY−5], f[i, NY−6]) for i ∈ x, j ∈ [NY−1, NY−1]
x ∈ [1, NX], y = NYboundary facecentral_D2y_6_dirichlet_face1(gT, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4], f[i, NY−5], f[i, NY−6]) for i ∈ x, j ∈ [NY, NY]

The y-half of a sixth-order 2-D Laplacian on cartesian_uniform_2d with inhomogeneous DIRICHLET conditions on the two y-walls, over the full 2-D field f[x,y], cell-centered. Interior y in [4, NY-3] (all x): the imported 7-point sixth-order y-stencil central_D2y_6_interior, O(dy^6). The six near-y-wall lines y in {1,2,3,NY-2,NY-1,NY}, whose centered 7-point stencils would reach outside the domain, are closed with the imported one-sided sixth-order face templates central_D2y_6_dirichlet_face1/face2/face3 (nearest, second, third row), each applied as a FULL-RANK aggregate over all x (singleton output range on y): the unique degree-7 polynomial matching the prescribed wall datum and the seven nearest cells (at fixed x), differentiated twice at that cell. The wall values gB and gT are consumer-supplied free names bound into the face templates (the g parameter); an un-overridden parameter defaults to 0, so gB=gT=0 is the homogeneous-Dirichlet case. The far wall reuses the same face weights by reflection (the second derivative is even), so f1..f7 there map to cells NY, NY-1, …, NY-6. The wall free names appear only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). Because BOTH the interior stencil and all six boundary closures are O(dy^6), the y-operator is uniformly sixth order with no boundary reduction; a steady linear-in-y profile carrying the wall values is discretely exact at every line (all closure weight-sums are zero). The full sixth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the y-term. Compound match on D(D(f,y),y) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. Requires NY >= 7. A model imports exactly one second-derivative rule for y.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (one-sided high-order closures) and Chapter 3 (multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case heat_2d_y_dirichlet_o6

Error norms read from the committed golden tests/conformance/convergence/heat_2d_y_dirichlet_o6/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
87.259e-081.873e-07
164.874e-107.229.892e-107.56
326.520e-126.221.558e-115.99
641.422e-135.521.044e-123.90

Expected order: 6 (± 0.6).

Convergence of heat_2d_y_dirichlet_o6 (log-log)

central_D2y_6_neumann_bc

op:D2 order:6 bc:neumann axes:x,y stencil width 7 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_6_neumann_bc.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [4, NY−3]interiorcentral_D2y_6_interior(f)[i,j] = (2·f[i, j−3] − 27·f[i, j−2] + 270·f[i, j−1] − 490·f[i, j] + 270·f[i, j+1] − 27·f[i, j+2] + 2·f[i, j+3]) / (180·dy·dy) for i ∈ x, j ∈ [4, NY−3]
x ∈ [1, NX], y = 1boundary facecentral_D2y_6_neumann_face1(qYL, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5], f[i, 6], f[i, 7]) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = 2boundary facecentral_D2y_6_neumann_face2(qYL, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5], f[i, 6], f[i, 7]) for i ∈ x, j ∈ [2, 2]
x ∈ [1, NX], y = 3boundary facecentral_D2y_6_neumann_face3(qYL, f[i, 1], f[i, 2], f[i, 3], f[i, 4], f[i, 5], f[i, 6], f[i, 7]) for i ∈ x, j ∈ [3, 3]
x ∈ [1, NX], y = NY−2boundary facecentral_D2y_6_neumann_face3(−qYR, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4], f[i, NY−5], f[i, NY−6]) for i ∈ x, j ∈ [NY−2, NY−2]
x ∈ [1, NX], y = NY−1boundary facecentral_D2y_6_neumann_face2(−qYR, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4], f[i, NY−5], f[i, NY−6]) for i ∈ x, j ∈ [NY−1, NY−1]
x ∈ [1, NX], y = NYboundary facecentral_D2y_6_neumann_face1(−qYR, f[i, NY], f[i, NY−1], f[i, NY−2], f[i, NY−3], f[i, NY−4], f[i, NY−5], f[i, NY−6]) for i ∈ x, j ∈ [NY, NY]

The y-half of a sixth-order 2-D Laplacian on cartesian_uniform_2d with inhomogeneous NEUMANN conditions on the two y-walls, over the full 2-D field f[x,y], cell-centered. Interior y in [4, NY-3] (all x): the imported 7-point sixth-order y-stencil central_D2y_6_interior, O(dy^6). The six near-y-wall lines y in {1,2,3,NY-2,NY-1,NY}, whose centered 7-point stencils would reach outside the domain, are closed with the imported one-sided sixth-order face templates central_D2y_6_neumann_face1/face2/face3 (nearest, second, third row), each applied as a FULL-RANK aggregate over all x (singleton output range on y): the unique degree-7 polynomial matching the prescribed wall datum and the seven nearest cells (at fixed x), differentiated twice at that cell. The wall gradients qYL (df/dy at the low wall) and qYR (df/dy at the high wall) are consumer-supplied free names; an un-overridden parameter defaults to 0, so qYL=qYR=0 is the homogeneous (zero-flux) case and the rule reduces exactly to a zero-gradient sixth-order closure. The low face templates take q as the +y-into-domain slope, so the low rows bind q -> qYL. The far wall reuses the SAME face weights by reflection (the second derivative even, df/dy odd), so the far rows bind q -> -qYR and map f1..f7 to cells NY, NY-1, …, NY-6. The wall free names appear only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). Because BOTH the interior stencil and all six boundary closures are O(dy^6), the y-operator is uniformly sixth order with no boundary reduction; a steady ramp-in-y carrying the wall fluxes is discretely exact at every line (all closure value-weight sums are zero). The full sixth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the y-term. Compound match on D(D(f,y),y) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. Requires NY >= 7. A model imports exactly one second-derivative rule for y.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (one-sided high-order closures) and Chapter 3 (multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case heat_2d_neumann_flux_o6

Error norms read from the committed golden tests/conformance/convergence/heat_2d_neumann_flux_o6/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
83.229e-071.019e-06
165.963e-095.762.827e-085.17
326.581e-116.504.193e-106.08
646.494e-136.664.322e-126.60

Expected order: 6 (± 0.6).

Convergence of heat_2d_neumann_flux_o6 (log-log)

central_D2y_6_periodic

op:D2 order:6 bc:periodic axes:x,y stencil width 7 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_6_periodic.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [4, NY−3]interiorcentral_D2y_6_interior(f)[i,j] = (2·f[i, j−3] − 27·f[i, j−2] + 270·f[i, j−1] − 490·f[i, j] + 270·f[i, j+1] − 27·f[i, j+2] + 2·f[i, j+3]) / (180·dy·dy) for i ∈ x, j ∈ [4, NY−3]
x ∈ [1, NX], y = 1boundary face(2·f[i, NY−2] − 27·f[i, NY−1] + 270·f[i, NY] − 490·f[i, 1] + 270·f[i, 2] − 27·f[i, 3] + 2·f[i, 4]) / (180·dy·dy) for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = 2boundary face(2·f[i, NY−1] − 27·f[i, NY] + 270·f[i, 1] − 490·f[i, 2] + 270·f[i, 3] − 27·f[i, 4] + 2·f[i, 5]) / (180·dy·dy) for i ∈ x, j ∈ [2, 2]
x ∈ [1, NX], y = 3boundary face(2·f[i, NY] − 27·f[i, 1] + 270·f[i, 2] − 490·f[i, 3] + 270·f[i, 4] − 27·f[i, 5] + 2·f[i, 6]) / (180·dy·dy) for i ∈ x, j ∈ [3, 3]
x ∈ [1, NX], y = NY−2boundary face(2·f[i, NY−5] − 27·f[i, NY−4] + 270·f[i, NY−3] − 490·f[i, NY−2] + 270·f[i, NY−1] − 27·f[i, NY] + 2·f[i, 1]) / (180·dy·dy) for i ∈ x, j ∈ [NY−2, NY−2]
x ∈ [1, NX], y = NY−1boundary face(2·f[i, NY−4] − 27·f[i, NY−3] + 270·f[i, NY−2] − 490·f[i, NY−1] + 270·f[i, NY] − 27·f[i, 1] + 2·f[i, 2]) / (180·dy·dy) for i ∈ x, j ∈ [NY−1, NY−1]
x ∈ [1, NX], y = NYboundary face(2·f[i, NY−3] − 27·f[i, NY−2] + 270·f[i, NY−1] − 490·f[i, NY] + 270·f[i, 1] − 27·f[i, 2] + 2·f[i, 3]) / (180·dy·dy) for i ∈ x, j ∈ [NY, NY]

The y-half of a sixth-order 2-D Laplacian on cartesian_uniform_2d with a PERIODIC y-boundary: sixth-order centered second derivative along y over the full 2-D field f[x,y], the y direction wrapped as a circle. Interior y in [4, NY-3] (all x): the imported 7-point y-stencil central_D2y_6_interior, (2 f[..-3] - 27 f[..-2] + 270 f[..-1] - 490 f[..] + 270 f[..+1] - 27 f[..+2] + 2 f[..+3]) / (180 dy^2). The six y-face lines y in {1,2,3,NY-2,NY-1,NY} carry the SAME 7-point weights with out-of-range y-neighbors wrapped periodically, so the whole y-operator is the exact circulant sixth-order Laplacian along y — O(dy^6) with no boundary order reduction, the 2-D per-axis analogue of grids/cartesian_uniform_1d/rules/central_D2_6_periodic. The six face values are FULL-RANK aggregates (singleton output range on the y axis, all x) so every binding’s simulator fills the lines without rank promotion. Requires NY >= 7 for the 7-point stencil to have distinct points. dy is the grid’s consumer-supplied y-spacing free name; the periodic wrap is dy-scale-invariant. The full sixth-order Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the y-term. Compound match on D(D(f,y),y) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. A model imports exactly one second-derivative rule for y.

References

  • Fornberg, B. (1988). Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of Computation 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains).

Conformance fixtures

Convergence — case heat_2d_y_periodic_o6

Error norms read from the committed golden tests/conformance/convergence/heat_2d_y_periodic_o6/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
81.503e-051.335e-05
162.524e-075.902.380e-075.81
324.016e-095.973.842e-095.95
646.303e-115.996.052e-115.99

Expected order: 6 (± 0.3).

Convergence of heat_2d_y_periodic_o6 (log-log)

central_D2y_dirichlet_bc

op:D2 order:2 bc:dirichlet axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_dirichlet_bc.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorcentral_D2y_interior(f)[i,j] = (f[i, j+1] − 2·f[i, j] + f[i, j−1]) / dy² for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face((f[i, 2] − 3·f[i, 1]) + 2·gB) / dy² for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face((f[i, NY−1] − 3·f[i, NY]) + 2·gT) / dy² for i ∈ x, j ∈ [NY, NY]

The y-half of a 2-D Laplacian on cartesian_uniform_2d with inhomogeneous DIRICHLET (fixed-value) conditions on the two y-walls, over the full 2-D field f[x,y], cell-centered. Interior j in [2, NY-1] (all i): the imported y-stencil (f[i,j+1] - 2 f[i,j] + f[i,j-1]) / dy^2. Faces j=1 and j=NY: mirror-ghost elimination against a prescribed wall value. gB and gT are consumer-supplied free names (the same free-name contract as the grid geometry dy; see grids/cartesian_uniform_2d/grid.esm) giving the field value on the BOTTOM y-wall y=y0 (the j=1 side) and the TOP y-wall y=y0+NYdy (the j=NY side), constant along that wall (all i). The names deliberately read gB/gT (bottom/top) so the y-wall assignment is unambiguous, the y-analogue of gL/gR on the x-walls. On this cell-centered grid the wall sits half a cell outside the first center, so linear reconstruction to the wall gives the ghost f[i,0]=2gB-f[i,1] (and f[i,NY+1]=2gT-f[i,NY]); substituting into the 3-point stencil yields D2y[i,1] = (f[i,2]-3f[i,1]+2gB)/dy^2 and D2y[i,NY] = (f[i,NY-1]-3f[i,NY]+2*gT)/dy^2. A consuming model declares gB/gT as ordinary real parameters; a parameter not overridden falls back to its default (esm-spec §6: default, or 0), so gB=gT=0 is the homogeneous-Dirichlet y-wall. The gB/gT terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). dy is the grid’s consumer-supplied y-spacing free name. The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring central_D2y_neumann_bc. The full Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the second term. Compound match on D(D(f,y),y) at priority 10 (esm-spec §9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec §9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec §9.7.7) the wrt literals and the where shape follow the renamed axes together. Interior O(dy^2); the one-sided mirror-ghost faces preserve global 2nd-order convergence for Dirichlet-compatible solutions.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (boundary conditions; cell-centered Dirichlet via mirror ghost) and Chapter 3 (multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case advection_2d_y_dirichlet_central

Error norms read from the committed golden tests/conformance/convergence/advection_2d_y_dirichlet_central/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
165.690e-046.559e-03
321.442e-041.981.829e-031.84
643.617e-052.004.809e-041.93
1289.050e-062.001.232e-041.96

Expected order: 2 (± 0.2).

Convergence of advection_2d_y_dirichlet_central (log-log)

Convergence — case anisotropic_diffusion_2d_dirichlet

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_dirichlet/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.657e-032.377e-03
326.589e-042.015.932e-042.00
641.644e-042.001.481e-042.00
1284.107e-052.003.703e-052.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_dirichlet (log-log)

Convergence — case heat_2d_y_dirichlet

Error norms read from the committed golden tests/conformance/convergence/heat_2d_y_dirichlet/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
163.167e-053.121e-05
327.925e-062.007.837e-061.99
641.982e-062.001.962e-062.00
1284.954e-072.004.905e-072.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_y_dirichlet (log-log)

central_D2y_neumann_bc

op:D2 order:2 bc:neumann axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_neumann_bc.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorcentral_D2y_interior(f)[i,j] = (f[i, j+1] − 2·f[i, j] + f[i, j−1]) / dy² for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face(f[i, j+1] − f[i, j] − qYL·dy) / dy² for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face((f[i, j−1] − f[i, j]) + qYR·dy) / dy² for i ∈ x, j ∈ [NY, NY]

The y-half of a 2-D Laplacian on cartesian_uniform_2d: second-order centered second derivative along y with inhomogeneous Neumann (fixed-flux) conditions on the two y-walls, over the full 2-D field f[x,y]. Interior j in [2, NY-1] (all i): the imported y-stencil. y-bottom face j=1: ghost elimination (f[i,j]-f[i,0])/dy = qYL gives D2y[i,1] = (f[i,2]-f[i,1]-qYLdy)/dy^2. y-top face j=NY: (f[i,NY+1]-f[i,NY])/dy = qYR gives D2y[i,NY] = (f[i,NY-1]-f[i,NY]+qYRdy)/dy^2. qYL/qYR are consumer-supplied free names (df/dy at the bottom/top walls, +y direction) declared as ordinary real parameters; an un-overridden parameter falls back to its default (esm-spec §6: default, or 0), so qYL=qYR=0 is the homogeneous (zero-gradient) y-wall. dy is the grid’s consumer-supplied spacing free name. The full Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): central_D2x_neumann_bc lowers the first term and this rule the second, and the two per-axis makearrays superpose to the correct 5-point interior and one-sided edges/corners. Compound match on D(D(f,y),y) at priority 10 (esm-spec §9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec §9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec §9.7.7) the wrt literals and the where shape follow the renamed axes together. Interior O(dy^2); the one-sided faces preserve global 2nd-order convergence for Neumann-compatible solutions.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 and Chapter 3 (Neumann conditions; multidimensional Poisson/heat by dimension-wise second differences).

Conformance fixtures

Convergence — case anisotropic_diffusion_2d_neumann

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_neumann/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
169.382e-036.067e-03
322.347e-032.001.520e-032.00
645.868e-042.003.802e-042.00
1281.467e-042.009.507e-052.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_neumann (log-log)

Convergence — case heat_2d_neumann_flux

Error norms read from the committed golden tests/conformance/convergence/heat_2d_neumann_flux/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
81.300e-042.380e-04
163.261e-052.006.150e-051.95
328.159e-062.001.550e-051.99
642.040e-062.003.883e-062.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_neumann_flux (log-log)

central_D2y_periodic

op:D2 order:2 bc:periodic axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_D2y_periodic.esm

Rewrites ∂²f/∂y² at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "y"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorcentral_D2y_interior(f)[i,j] = (f[i, j+1] − 2·f[i, j] + f[i, j−1]) / dy² for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face(f[i, 2] − 2·f[i, 1] + f[i, NY]) / dy² for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face(f[i, 1] − 2·f[i, NY] + f[i, NY−1]) / dy² for i ∈ x, j ∈ [NY, NY]

The y-half of a 2-D Laplacian on cartesian_uniform_2d with a PERIODIC y-boundary: second-order centered second derivative along y over the full 2-D field f[x,y], the y direction wrapped as a circle. Interior j in [2, NY-1] (all i): the imported y-stencil (f[i,j+1] - 2 f[i,j] + f[i,j-1]) / dy^2. The two y-face rows wrap periodically: j=1 uses the ghost f[i,0]=f[i,NY], giving (f[i,2] - 2 f[i,1] + f[i,NY]) / dy^2; j=NY uses f[i,NY+1]=f[i,1], giving (f[i,1] - 2 f[i,NY] + f[i,NY-1]) / dy^2. The scheme is the same O(dy^2) centered second difference at every row, so the whole y-operator is the exact circulant symmetric Laplacian along y — no boundary order reduction, the natural second-order analogue of central_D1y_periodic. The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring central_D2y_neumann_bc so every binding’s simulator fills the (NX x 1) rows without rank promotion. dy is the grid’s consumer-supplied y-spacing free name; the periodic wrap is dy-scale-invariant, so the rule serves any y-extent unchanged. The full Laplacian is assembled in the consuming equation as D(D(f,x),x) + D(D(f,y),y): this rule lowers the second term. Compound match on D(D(f,y),y) at priority 10 (esm-spec §9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec §9.6.1) so it fires only on the yy-second-derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec §9.7.7) the wrt literals and the where shape follow the renamed axes together, so two instances coexist scoped to their own meshes. O(dy^2). A model imports exactly one second-derivative rule for y (this periodic rule or a Neumann/Dirichlet variant), never more than one.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains) and Section 2.12 (second differences).

Conformance fixtures

Convergence — case anisotropic_diffusion_2d_periodic

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_periodic/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.489e-022.396e-02
326.326e-031.986.266e-031.94
641.589e-031.991.585e-031.98
1283.976e-042.003.973e-042.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_periodic (log-log)

Convergence — case heat_2d_y_periodic

Error norms read from the committed golden tests/conformance/convergence/heat_2d_y_periodic/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
165.049e-044.760e-04
321.267e-041.991.212e-041.97
643.170e-052.003.044e-051.99
1287.927e-062.007.617e-062.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_y_periodic (log-log)

central_Dxy_dirichlet_bc

op:Dxy order:2 bc:dirichlet axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_Dxy_dirichlet_bc.esm

Rewrites ∂∂f/∂x/∂y at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [2, NY−1]interiorcentral_Dxy_interior(f)[i,j] = (f[i+1, j+1] − f[i+1, j−1] − (f[i−1, j+1] − f[i−1, j−1])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [2, NY−1]
x = 1, y ∈ [2, NY−1]boundary face(f[2, j+1] + f[1, j+1] − (f[2, j−1] + f[1, j−1])) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [2, NY−1]
x = NX, y ∈ [2, NY−1]boundary face(f[NX, j−1] + f[NX−1, j−1] − (f[NX, j+1] + f[NX−1, j+1])) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [2, NY−1]
x ∈ [2, NX−1], y = 1boundary face(f[i+1, 2] + f[i+1, 1] − (f[i−1, 2] + f[i−1, 1])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [1, 1]
x ∈ [2, NX−1], y = NYboundary face(f[i−1, NY] + f[i−1, NY−1] − (f[i+1, NY] + f[i+1, NY−1])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [NY, NY]
x = 1, y = 1boundary face(f[2, 2] + f[2, 1] + f[1, 2] + f[1, 1]) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [1, 1]
x = 1, y = NYboundary face(0 − (f[2, NY] + f[2, NY−1] + f[1, NY] + f[1, NY−1])) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [NY, NY]
x = NX, y = 1boundary face(0 − (f[NX, 2] + f[NX, 1] + f[NX−1, 2] + f[NX−1, 1])) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [1, 1]
x = NX, y = NYboundary face(f[NX, NY] + f[NX, NY−1] + f[NX−1, NY] + f[NX−1, NY−1]) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [NY, NY]

Second-order centered MIXED second derivative d^2f/dx dy on cartesian_uniform_2d with DIRICHLET (fixed-value) walls on all four sides, over the full 2-D field f[x,y], cell-centered. This is the off-diagonal 2 Kxy d^2f/dx dy term of a full anisotropic diffusion tensor K = [[Kxx, Kxy], [Kxy, Kyy]] that the diagonal central_D2x/central_D2y rules do not supply; the non-periodic wall analogue of central_Dxy_periodic. INTERIOR block [2,NX-1]x[2,NY-1]: the imported central_Dxy_interior 4-point stencil (f[i+1,j+1]-f[i+1,j-1]-f[i-1,j+1]+f[i-1,j-1])/(4 dx dy), which reaches diagonally to all four corner neighbours, so the wall closure needs ghosts across whichever axis is out of the domain (one axis on the four edge strips, two axes plus the diagonal corner ghost on the four corner cells). GHOST RECONSTRUCTION: on this cell-centered grid the wall sits half a cell outside the first center, so the Dirichlet value g_wall gives the mirror (odd-reflection) ghost f_ghost = 2 g_wall - f_interior (left f[0,j]=2 gL-f[1,j], right f[NX+1,j]=2 gR-f[NX,j], bottom f[i,0]=2 gB-f[i,1], top f[i,NY+1]=2 gT-f[i,NY]); the diagonal corner ghost is that reconstruction applied in both directions (for the homogeneous-compatible case f[0,0]=f[1,1], etc.). CANCELLATION OF THE WALL VALUE (why no free names appear): the mixed derivative differences ALONG the wall (the tangential direction), and the prescribed value is constant along each wall (scalar gL/gR/gB/gT, matching central_D2x_dirichlet_bc / central_D2y_dirichlet_bc). At the left x-wall edge (i=1) the stencil’s antisymmetric y-difference of the ghost column is -f[0,j+1]+f[0,j-1] = -(2 gL-f[1,j+1])+(2 gL-f[1,j-1]) = f[1,j+1]-f[1,j-1]: the 2 gL drops IDENTICALLY. The same tangential-difference cancellation removes gR, gB, gT on the other three edges and every corner (the diagonal-ghost constant terms cancel against the two edge-ghost constant terms). So this Dirichlet cross closure carries NO wall free names - it encodes only the structural fact that df/dy=0 along the x-walls and df/dx=0 along the y-walls (both implied by f being constant along the wall), which is exactly the boundary information the mixed derivative needs; the specific wall values are irrelevant to d^2f/dx dy. (A wall datum that VARIES along the wall would add a tangential-derivative term this scalar-wall closure omits, consistent with the scalar-wall D2 rules; out of scope.) RESULTING FACE FORMULAS (all over 4 dx dy): left edge i=1 -> (f[2,j+1]-f[2,j-1])+(f[1,j+1]-f[1,j-1]); right edge i=NX -> (f[NX,j-1]-f[NX,j+1])+(f[NX-1,j-1]-f[NX-1,j+1]); bottom edge j=1 -> (f[i+1,2]-f[i-1,2])+(f[i+1,1]-f[i-1,1]); top edge j=NY -> (f[i-1,NY]-f[i+1,NY])+(f[i-1,NY-1]-f[i+1,NY-1]); corner (1,1) -> f[2,2]+f[2,1]+f[1,2]+f[1,1]; corner (1,NY) -> -(f[2,NY]+f[2,NY-1]+f[1,NY]+f[1,NY-1]); corner (NX,1) -> -(f[NX,2]+f[NX,1]+f[NX-1,2]+f[NX-1,1]); corner (NX,NY) -> f[NX,NY]+f[NX,NY-1]+f[NX-1,NY]+f[NX-1,NY-1]. Each edge/corner is the mirror-ghost centered cross difference, O(dx^2+dy^2) like the interior, and is exactly the odd-reflected form of the interior stencil (the odd extension is exact for a solution that is odd about the wall, e.g. sin(pi x) about x=0). The makearray tiles [x,y] as the 3x3 block partition {1}u[2,NX-1]u{NX} x {1}u[2,NY-1]u{NY} (interior + four edge strips + four corner cells), mirroring central_Dxy_periodic and central_D2x/D2y_dirichlet_bc so every binding’s simulator fills the strips and corners without rank promotion. dx/dy are the grid’s consumer-supplied spacing free names (grids/cartesian_uniform_2d/grid.esm). WELL-POSEDNESS: the isolated cross operator is symmetric indefinite (eigenvalues span -N^2..+N^2), not a stable time-evolution operator on its own; the accompanying MMS (problems/anisotropic_diffusion_2d_dirichlet.esm) embeds it in the positive-definite tensor K=[[1,1/2],[1/2,1]] so L = D2x+Dxy+D2y is dissipative. With homogeneous Dirichlet walls the discrete energy estimate’s boundary flux vanishes (the test function is zero on the wall), so the combined operator is provably negative definite and the parabolic evolution is well-posed - unlike the Neumann cross closure, whose walls leave a residual conormal boundary flux. Compound match on D(D(f,x),y) at priority 10 (esm-spec 9.6.3); the consumer writes d^2u/dx dy as D(D(u,x),y) (inner wrt x, outer wrt y), and the priority-10 compound fires before any first-derivative rule could rewrite the inner D. Match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the mixed second derivative of a bare field declared over this 2-D grid’s [x,y] axes; under import-edge renaming (esm-spec 9.7.7) the two wrt literals and the where shape follow the renamed axes together, so two instances coexist scoped to their own meshes. O(dx^2+dy^2).

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (mixed partial derivatives; cell-centered ghost reconstruction for Dirichlet/Neumann walls) and Chapter 3 (multidimensional operators by tensor products of one-dimensional differences).
  • Morton, K. W. and Mayers, D. F. (2005). Numerical Solution of Partial Differential Equations, 2nd ed. Cambridge University Press. Section 3 (boundary closures by ghost cells; tangential vs normal derivatives at walls).

Conformance fixtures

Convergence — case anisotropic_diffusion_2d_dirichlet

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_dirichlet/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.657e-032.377e-03
326.589e-042.015.932e-042.00
641.644e-042.001.481e-042.00
1284.107e-052.003.703e-052.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_dirichlet (log-log)

central_Dxy_neumann_bc

op:Dxy order:2 bc:neumann axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_Dxy_neumann_bc.esm

Rewrites ∂∂f/∂x/∂y at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [2, NY−1]interiorcentral_Dxy_interior(f)[i,j] = (f[i+1, j+1] − f[i+1, j−1] − (f[i−1, j+1] − f[i−1, j−1])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [2, NY−1]
x = 1, y ∈ [2, NY−1]boundary face(f[2, j+1] + f[1, j−1] − (f[2, j−1] + f[1, j+1])) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [2, NY−1]
x = NX, y ∈ [2, NY−1]boundary face(f[NX, j+1] + f[NX−1, j−1] − (f[NX, j−1] + f[NX−1, j+1])) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [2, NY−1]
x ∈ [2, NX−1], y = 1boundary face(f[i+1, 2] + f[i−1, 1] − (f[i+1, 1] + f[i−1, 2])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [1, 1]
x ∈ [2, NX−1], y = NYboundary face(f[i+1, NY] + f[i−1, NY−1] − (f[i+1, NY−1] + f[i−1, NY])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [NY, NY]
x = 1, y = 1boundary face(f[2, 2] + f[1, 1] − (f[2, 1] + f[1, 2])) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [1, 1]
x = 1, y = NYboundary face(f[2, NY] + f[1, NY−1] − (f[2, NY−1] + f[1, NY])) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [NY, NY]
x = NX, y = 1boundary face(f[NX, 2] + f[NX−1, 1] − (f[NX, 1] + f[NX−1, 2])) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [1, 1]
x = NX, y = NYboundary face(f[NX, NY] + f[NX−1, NY−1] − (f[NX, NY−1] + f[NX−1, NY])) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [NY, NY]

Second-order centered MIXED second derivative d^2f/dx dy on cartesian_uniform_2d with NEUMANN (fixed-flux) walls on all four sides, over the full 2-D field f[x,y], cell-centered. This is the off-diagonal 2 Kxy d^2f/dx dy term of a full anisotropic diffusion tensor K = [[Kxx, Kxy], [Kxy, Kyy]] that the diagonal central_D2x/central_D2y rules do not supply; the fixed-flux wall analogue of central_Dxy_periodic and the Neumann counterpart of central_Dxy_dirichlet_bc. INTERIOR block [2,NX-1]x[2,NY-1]: the imported central_Dxy_interior 4-point stencil (f[i+1,j+1]-f[i+1,j-1]-f[i-1,j+1]+f[i-1,j-1])/(4 dx dy), which reaches diagonally, so the wall closure needs ghosts across whichever axis is out of the domain (one axis on the edges, two axes plus the diagonal corner ghost on the corners). GHOST RECONSTRUCTION: the prescribed normal derivative q_wall gives the even-reflection ghost from the one-sided flux relation (left (f[1,j]-f[0,j])/dx=qXL -> f[0,j]=f[1,j]-qXL dx; right f[NX+1,j]=f[NX,j]+qXR dx; bottom f[i,0]=f[i,1]-qYL dy; top f[i,NY+1]=f[i,NY]+qYR dy); the diagonal corner ghost applies the flux relation in both directions and is UNAMBIGUOUS (both orders give f[0,0]=f[1,1]-qXL dx-qYL dy). CANCELLATION OF THE WALL FLUX (why no free names appear): the mixed derivative differences ALONG the wall, and q_wall is constant along each wall (scalar qXL/qXR/qYL/qYR, matching central_D2x_neumann_bc / central_D2y_neumann_bc). At the left x-wall edge (i=1) the antisymmetric y-difference of the ghost column is -f[0,j+1]+f[0,j-1] = -(f[1,j+1]-qXL dx)+(f[1,j-1]-qXL dx) = -f[1,j+1]+f[1,j-1]: the qXL dx drops IDENTICALLY. The same tangential-difference cancellation removes qXR, qYL, qYR on the other edges and every corner (the diagonal-ghost -qXL dx-qYL dy terms cancel against the two edge-ghost flux terms). So this Neumann cross closure carries NO wall free names either; it encodes only the EVEN reflection of the field across each wall (the fixed-flux continuation), the sign-flipped counterpart of the odd-reflection Dirichlet closure. RESULTING FACE FORMULAS (all over 4 dx dy), the interior stencil with each out-of-range sample replaced by its EVEN mirror (reflected index, same sign): left edge i=1 -> (f[2,j+1]-f[2,j-1])-(f[1,j+1]-f[1,j-1]); right edge i=NX -> (f[NX,j+1]-f[NX,j-1])-(f[NX-1,j+1]-f[NX-1,j-1]); bottom edge j=1 -> (f[i+1,2]-f[i-1,2])-(f[i+1,1]-f[i-1,1]); top edge j=NY -> (f[i+1,NY]-f[i-1,NY])-(f[i+1,NY-1]-f[i-1,NY-1]); corner (1,1) -> (f[2,2]+f[1,1])-(f[2,1]+f[1,2]); corner (1,NY) -> (f[2,NY]+f[1,NY-1])-(f[2,NY-1]+f[1,NY]); corner (NX,1) -> (f[NX,2]+f[NX-1,1])-(f[NX,1]+f[NX-1,2]); corner (NX,NY) -> (f[NX,NY]+f[NX-1,NY-1])-(f[NX,NY-1]+f[NX-1,NY]). Each face is the even-reflected centered cross difference, O(dx^2+dy^2) like the interior, exact for a solution even about the wall (e.g. cos(pi x) about x=0). The makearray tiles [x,y] as the 3x3 block partition {1}u[2,NX-1]u{NX} x {1}u[2,NY-1]u{NY}, mirroring central_Dxy_periodic and central_D2x/D2y_neumann_bc. dx/dy are the grid’s consumer-supplied spacing free names (grids/cartesian_uniform_2d/grid.esm). WELL-POSEDNESS CAVEAT: the isolated cross operator is symmetric indefinite; embedding it in the positive-definite tensor K=[[1,1/2],[1/2,1]] (L = D2x+Dxy+D2y) makes the DIRICHLET-walled operator dissipative, but with per-axis Neumann walls (df/dx=0 on x-walls, df/dy=0 on y-walls) the discrete energy estimate retains a residual conormal boundary flux n.K grad u = Kxy u_tangential on each wall that the imposed normal condition does not kill, so the combined operator is NOT guaranteed negative definite and the pure time-evolution MMS is only conditionally well-posed. See problems/anisotropic_diffusion_2d_neumann.esm for the manufactured-solution setup and the observed stability at the tested resolutions. Compound match on D(D(f,x),y) at priority 10 (esm-spec 9.6.3), match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the mixed second derivative of a bare field declared over this 2-D grid’s [x,y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and the where shape follow the renamed axes together. O(dx^2+dy^2).

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (mixed partial derivatives; cell-centered ghost reconstruction for Dirichlet/Neumann walls) and Chapter 3 (multidimensional operators by tensor products of one-dimensional differences).
  • Morton, K. W. and Mayers, D. F. (2005). Numerical Solution of Partial Differential Equations, 2nd ed. Cambridge University Press. Section 3 (boundary closures by ghost cells; tangential vs normal derivatives at walls).

Conformance fixtures

Convergence — case anisotropic_diffusion_2d_neumann

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_neumann/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
169.382e-036.067e-03
322.347e-032.001.520e-032.00
645.868e-042.003.802e-042.00
1281.467e-042.009.507e-052.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_neumann (log-log)

central_Dxy_periodic

op:Dxy order:2 bc:periodic axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/central_Dxy_periodic.esm

Rewrites ∂∂f/∂x/∂y at priority 10 — match pattern:

{"op": "D", "args": [{"op": "D", "args": ["f"], "wrt": "x"}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when f is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [2, NY−1]interiorcentral_Dxy_interior(f)[i,j] = (f[i+1, j+1] − f[i+1, j−1] − (f[i−1, j+1] − f[i−1, j−1])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [2, NY−1]
x = 1, y ∈ [2, NY−1]boundary face(f[2, j+1] − f[2, j−1] − (f[NX, j+1] − f[NX, j−1])) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [2, NY−1]
x = NX, y ∈ [2, NY−1]boundary face(f[1, j+1] − f[1, j−1] − (f[NX−1, j+1] − f[NX−1, j−1])) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [2, NY−1]
x ∈ [2, NX−1], y = 1boundary face(f[i+1, 2] − f[i+1, NY] − (f[i−1, 2] − f[i−1, NY])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [1, 1]
x ∈ [2, NX−1], y = NYboundary face(f[i+1, 1] − f[i+1, NY−1] − (f[i−1, 1] − f[i−1, NY−1])) / (4·(dx·dy)) for i ∈ [2, NX−1], j ∈ [NY, NY]
x = 1, y = 1boundary face(f[2, 2] − f[2, NY] − (f[NX, 2] − f[NX, NY])) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [1, 1]
x = 1, y = NYboundary face(f[2, 1] − f[2, NY−1] − (f[NX, 1] − f[NX, NY−1])) / (4·(dx·dy)) for i ∈ [1, 1], j ∈ [NY, NY]
x = NX, y = 1boundary face(f[1, 2] − f[1, NY] − (f[NX−1, 2] − f[NX−1, NY])) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [1, 1]
x = NX, y = NYboundary face(f[1, 1] − f[1, NY−1] − (f[NX−1, 1] − f[NX−1, NY−1])) / (4·(dx·dy)) for i ∈ [NX, NX], j ∈ [NY, NY]

Second-order centered MIXED second derivative d^2f/dx dy on cartesian_uniform_2d with a DOUBLY-PERIODIC boundary (both x and y wrapped as circles), over the full 2-D field f[x,y]. This is the off-diagonal term of a full anisotropic diffusion tensor K = [[Kxx, Kxy], [Kxy, Kyy]] — the 2 Kxy d^2f/dx dy contribution the diagonal central_D2x/central_D2y rules do not supply. The 4-point stencil (f[i+1,j+1] - f[i+1,j-1] - f[i-1,j+1] + f[i-1,j-1]) / (4 dx dy) reaches diagonally to all four corner neighbours, so the periodic closure needs BOTH x-walls AND both y-walls to wrap, plus the four corners where the two wraps meet. The makearray therefore tiles [x,y] as a 3x3 block partition {1}u[2,NX-1]u{NX} x {1}u[2,NY-1]u{NY}: the doubly-interior block [2,NX-1]x[2,NY-1] holds the imported central_Dxy_interior stencil; the four edge strips wrap the single out-of-range axis (i=1 uses the ghost column NX for i-1, i=NX uses column 1 for i+1; j=1 uses row NY for j-1, j=NY uses row 1 for j+1); the four corner cells wrap both axes. Every cell uses the identical O(dx^2 + dy^2) centered cross difference (the wrap is exact on the periodic torus, no boundary order reduction), so the whole operator is the circulant tensor product D1x (x) D1y — the periodic cross-derivative analogue of central_D1x_periodic/central_D1y_periodic. The edge and corner faces are FULL-RANK aggregates (singleton output range on the wrapped axis) mirroring central_D2x_periodic/central_D2y_periodic so every binding’s simulator fills the strips and corner cells without rank promotion. dx/dy are the grid’s consumer-supplied spacing free names; the periodic wrap is scale-invariant, so the rule serves any x/y extent unchanged. Compound match on D(D(f,x),y) at priority 10 (esm-spec 9.6.3) — the consumer writes d^2u/dx dy as D(D(u,x),y) (inner wrt x, outer wrt y); the priority-10 compound fires before any first-derivative rule could rewrite the inner D. Match-scoped with where: {f: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the mixed second derivative of a bare field declared over this 2-D grid’s [x, y] axes; under import-edge renaming (esm-spec 9.7.7) the two wrt literals and the where shape follow the renamed axes together, so two instances coexist scoped to their own meshes. A one-sided (non-periodic) closure of the cross derivative is a follow-up; periodic is the cleanest self-contained boundary for a diagonal-reaching cross stencil (no ghost-gradient bookkeeping on four walls plus four corners), which is why the accompanying MMS (problems/anisotropic_diffusion_2d_periodic.esm) is doubly periodic. NOTE the isolated cross operator is symmetric indefinite (its discrete eigenvalues span -N^2..+N^2), so it is not a well-posed time-evolution operator on its own; the MMS embeds it in a positive-definite diffusion tensor, whose combined operator is dissipative — see the problem file for the derivation.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 10.4 (centered differences on periodic domains) and Section 2.12 (mixed partial derivatives).

Conformance fixtures

Convergence — case anisotropic_diffusion_2d_periodic

Error norms read from the committed golden tests/conformance/convergence/anisotropic_diffusion_2d_periodic/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.489e-022.396e-02
326.326e-031.986.266e-031.94
641.589e-031.991.585e-031.98
1283.976e-042.003.973e-042.00

Expected order: 2 (± 0.2).

Convergence of anisotropic_diffusion_2d_periodic (log-log)

godunov_norm_D2_dirichlet_bc

op:grad_norm order:1 bc:dirichlet axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/godunov_norm_D2_dirichlet_bc.esm

Rewrites √(∂u/∂x² + ∂u/∂y²) at priority 10 — match pattern:

{"op": "sqrt", "args": [{"op": "+", "args": [{"op": "^", "args": [{"op": "D", "args": ["u"], "wrt": "x"}, 2]}, {"op": "^", "args": [{"op": "D", "args": ["u"], "wrt": "y"}, 2]}]}]}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [2, NY−1]interiorgodunov_norm_D2_interior(u)[i,j] = sqrt(max((u[i, j] − u[i−1, j]) / dx, 0)² + min((u[i+1, j] − u[i, j]) / dx, 0)² + max((u[i, j] − u[i, j−1]) / dy, 0)² + min((u[i, j+1] − u[i, j]) / dy, 0)²) for i ∈ [2, NX−1], j ∈ [2, NY−1]
x = 1, y ∈ [2, NY−1]boundary facesqrt(max((u[1, j] − (2·gL − u[1, j])) / dx, 0)² + min((u[2, j] − u[1, j]) / dx, 0)² + max((u[1, j] − u[1, j−1]) / dy, 0)² + min((u[1, j+1] − u[1, j]) / dy, 0)²) for i ∈ [1, 1], j ∈ [2, NY−1]
x = NX, y ∈ [2, NY−1]boundary facesqrt(max((u[NX, j] − u[NX−1, j]) / dx, 0)² + min((2·gR − u[NX, j] − u[NX, j]) / dx, 0)² + max((u[NX, j] − u[NX, j−1]) / dy, 0)² + min((u[NX, j+1] − u[NX, j]) / dy, 0)²) for i ∈ [NX, NX], j ∈ [2, NY−1]
x ∈ [2, NX−1], y = 1boundary facesqrt(max((u[i, 1] − u[i−1, 1]) / dx, 0)² + min((u[i+1, 1] − u[i, 1]) / dx, 0)² + max((u[i, 1] − (2·gB − u[i, 1])) / dy, 0)² + min((u[i, 2] − u[i, 1]) / dy, 0)²) for i ∈ [2, NX−1], j ∈ [1, 1]
x ∈ [2, NX−1], y = NYboundary facesqrt(max((u[i, NY] − u[i−1, NY]) / dx, 0)² + min((u[i+1, NY] − u[i, NY]) / dx, 0)² + max((u[i, NY] − u[i, NY−1]) / dy, 0)² + min((2·gT − u[i, NY] − u[i, NY]) / dy, 0)²) for i ∈ [2, NX−1], j ∈ [NY, NY]
x = 1, y = 1boundary facesqrt(max((u[1, 1] − (2·gL − u[1, 1])) / dx, 0)² + min((u[2, 1] − u[1, 1]) / dx, 0)² + max((u[1, 1] − (2·gB − u[1, 1])) / dy, 0)² + min((u[1, 2] − u[1, 1]) / dy, 0)²) for i ∈ [1, 1], j ∈ [1, 1]
x = 1, y = NYboundary facesqrt(max((u[1, NY] − (2·gL − u[1, NY])) / dx, 0)² + min((u[2, NY] − u[1, NY]) / dx, 0)² + max((u[1, NY] − u[1, NY−1]) / dy, 0)² + min((2·gT − u[1, NY] − u[1, NY]) / dy, 0)²) for i ∈ [1, 1], j ∈ [NY, NY]
x = NX, y = 1boundary facesqrt(max((u[NX, 1] − u[NX−1, 1]) / dx, 0)² + min((2·gR − u[NX, 1] − u[NX, 1]) / dx, 0)² + max((u[NX, 1] − (2·gB − u[NX, 1])) / dy, 0)² + min((u[NX, 2] − u[NX, 1]) / dy, 0)²) for i ∈ [NX, NX], j ∈ [1, 1]
x = NX, y = NYboundary facesqrt(max((u[NX, NY] − u[NX−1, NY]) / dx, 0)² + min((2·gR − u[NX, NY] − u[NX, NY]) / dx, 0)² + max((u[NX, NY] − u[NX, NY−1]) / dy, 0)² + min((2·gT − u[NX, NY] − u[NX, NY]) / dy, 0)²) for i ∈ [NX, NX], j ∈ [NY, NY]

First-order Godunov (Rouy-Tourin / Osher-Sethian) upwind discretization of the gradient magnitude |grad u| = sqrt(u_x^2 + u_y^2) for the non-negative-speed eikonal term psi_t = -S|grad psi| (S >= 0) on cartesian_uniform_2d, with an inhomogeneous DIRICHLET (state-constraint) boundary u = g prescribed on all four rectangle walls – the natural boundary condition for eikonal / level-set problems |grad u| = f with u given on the boundary. It is the Dirichlet-closed sibling of grids/cartesian_uniform_2d/rules/godunov_norm_D2_periodic and the 2-D companion of godunov_norm_D1_dirichlet_bc: the interior Godunov Hamiltonian and its per-axis entropy (max/min) upwind selection are UNCHANGED – only the out-of-range neighbor on each boundary strip changes from a periodic wrap to a prescribed-value (mirror-ghost) source. COMPOUND MATCH (esm-spec 9.6.3, the canonical 2-D gradient-magnitude compound-precedence example): matches the whole compound sqrt((D(u, wrt: x))^2 + (D(u, wrt: y))^2) rooted at the outer sqrt at priority 10, IDENTICAL to godunov_norm_D2_periodic’s match, firing before any plain first-derivative rule could lower an inner D. Match-scoped with where {u: {shape: [x, y]}} (esm-spec 9.6.1). THE SCHEME: the makearray tiles [x, y] as the same 3x3 block partition {1} u [2,NX-1] u {NX} x {1} u [2,NY-1] u {NY} as central_D2x/y and the periodic sibling. The doubly-interior block [2,NX-1] x [2,NY-1] is a single apply_expression_template of the imported stencils/godunov_norm_D2_interior aggregate. The four edge strips wrap their single out-of-range axis with a mirror ghost against the prescribed wall value: on cell-centered walls half a cell outside the first/last center, the ghost is u[0,j]=2gL-u[1,j] (left), u[NX+1,j]=2gR-u[NX,j] (right), u[i,0]=2gB-u[i,1] (bottom), u[i,NY+1]=2gT-u[i,NY] (top); the four corner cells wrap BOTH axes with the corresponding pair of ghosts. Every cell evaluates the identical per-axis max/min entropy norm sqrt(max(D-_x,0)^2 + min(D+_x,0)^2 + max(D-_y,0)^2 + min(D+_y,0)^2); the edge and corner faces are FULL-RANK aggregates (singleton output range on the wrapped axis) mirroring central_D2x_dirichlet_bc so every simulator fills the strips and corners without rank promotion. gL, gR are consumer-supplied free names giving the value on the left x-wall x=x0 and right x-wall x=x0+NXdx (constant along that wall, all j); gB, gT give the value on the bottom y-wall y=y0 and top y-wall y=y0+NYdy (constant along that wall, all i); a parameter not overridden defaults to 0 (esm-spec 6), so gL=gR=gB=gT=0 is the homogeneous case. The g* terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). EXACTNESS: the scheme ANNIHILATES a constant field exactly (every difference is 0 to the bit, the null-space property |grad(const)| = 0) and reproduces |grad u| exactly for any field whose wall traces are captured by the scalar wall values (e.g. u = c constant, or the manufactured g = c + sin(pi x) sin(pi y) whose four walls all equal c). A general non-constant affine field u = a x + b y + c has wall traces that VARY along each wall (the left wall trace c + b*y depends on y), which a scalar wall value cannot represent exactly; such array-valued wall data is out of scope for the constant-along-wall Dirichlet free-name convention shared with central_D2x/y_dirichlet_bc, so 2-D exactness is stated for constant / constant-wall-trace fields. First order O(dx + dy) in space on smooth data; the mirror ghost is 2nd-order-consistent so the boundary retains the interior order for Dirichlet-compatible solutions. dx, dy are free names (consumer-supplied geometry contract). A model imports exactly one gradient-magnitude rule for the [x, y] pair. IMPLEMENTATION NOTE: the interior Godunov Hamiltonian’s sqrt(sum of squares) / min-max is non-smooth and the 2-D lowered RHS is large; the eikonal MMS driver problems/godunov_norm_2d_eikonal_dirichlet_mms.esm verifies this rule in Julia (Tsit5) and Python (LSODA); the Rust diffsol integrator lacks a step-size fail-fast guard and wedges on this RHS (see the case manifests’ blocked_upstream_bindings, by analogy to the periodic sibling).

References

  • Osher, S., & Sethian, J. A. (1988). Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79(1), 12-49. The max/min upwind (entropy) construction of |grad phi| for the level-set Hamilton-Jacobi term. doi:10.1016/0021-9991(88)90002-2
  • Rouy, E., & Tourin, A. (1992). A viscosity solutions approach to shape-from-shading. SIAM Journal on Numerical Analysis, 29(3), 867-884. The Godunov Hamiltonian sqrt(max(D-,0)^2 + min(D+,0)^2) as the monotone, first-order-convergent discretization of the eikonal |grad u|. doi:10.1137/0729053
  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.12 (cell-centered Dirichlet via mirror ghost).

Conformance fixtures

Convergence — case godunov_norm_2d_eikonal_dirichlet_mms

Error norms read from the committed golden tests/conformance/convergence/godunov_norm_2d_eikonal_dirichlet_mms/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.604e-026.653e-02
321.344e-020.953.502e-020.93
646.838e-030.981.819e-020.95
1283.449e-030.999.366e-030.96

Expected order: 1 (± 0.2).

Convergence of godunov_norm_2d_eikonal_dirichlet_mms (log-log)

godunov_norm_D2_periodic

op:grad_norm order:1 bc:periodic axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/godunov_norm_D2_periodic.esm

Rewrites √(∂u/∂x² + ∂u/∂y²) at priority 10 — match pattern:

{"op": "sqrt", "args": [{"op": "+", "args": [{"op": "^", "args": [{"op": "D", "args": ["u"], "wrt": "x"}, 2]}, {"op": "^", "args": [{"op": "D", "args": ["u"], "wrt": "y"}, 2]}]}]}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [2, NY−1]interiorsqrt(max((u[i, j] − u[i−1, j]) / dx, 0)² + min((u[i+1, j] − u[i, j]) / dx, 0)² + max((u[i, j] − u[i, j−1]) / dy, 0)² + min((u[i, j+1] − u[i, j]) / dy, 0)²) for i ∈ [2, NX−1], j ∈ [2, NY−1]
x = 1, y ∈ [2, NY−1]boundary facesqrt(max((u[1, j] − u[NX, j]) / dx, 0)² + min((u[2, j] − u[1, j]) / dx, 0)² + max((u[1, j] − u[1, j−1]) / dy, 0)² + min((u[1, j+1] − u[1, j]) / dy, 0)²) for i ∈ [1, 1], j ∈ [2, NY−1]
x = NX, y ∈ [2, NY−1]boundary facesqrt(max((u[NX, j] − u[NX−1, j]) / dx, 0)² + min((u[1, j] − u[NX, j]) / dx, 0)² + max((u[NX, j] − u[NX, j−1]) / dy, 0)² + min((u[NX, j+1] − u[NX, j]) / dy, 0)²) for i ∈ [NX, NX], j ∈ [2, NY−1]
x ∈ [2, NX−1], y = 1boundary facesqrt(max((u[i, 1] − u[i−1, 1]) / dx, 0)² + min((u[i+1, 1] − u[i, 1]) / dx, 0)² + max((u[i, 1] − u[i, NY]) / dy, 0)² + min((u[i, 2] − u[i, 1]) / dy, 0)²) for i ∈ [2, NX−1], j ∈ [1, 1]
x ∈ [2, NX−1], y = NYboundary facesqrt(max((u[i, NY] − u[i−1, NY]) / dx, 0)² + min((u[i+1, NY] − u[i, NY]) / dx, 0)² + max((u[i, NY] − u[i, NY−1]) / dy, 0)² + min((u[i, 1] − u[i, NY]) / dy, 0)²) for i ∈ [2, NX−1], j ∈ [NY, NY]
x = 1, y = 1boundary facesqrt(max((u[1, 1] − u[NX, 1]) / dx, 0)² + min((u[2, 1] − u[1, 1]) / dx, 0)² + max((u[1, 1] − u[1, NY]) / dy, 0)² + min((u[1, 2] − u[1, 1]) / dy, 0)²) for i ∈ [1, 1], j ∈ [1, 1]
x = 1, y = NYboundary facesqrt(max((u[1, NY] − u[NX, NY]) / dx, 0)² + min((u[2, NY] − u[1, NY]) / dx, 0)² + max((u[1, NY] − u[1, NY−1]) / dy, 0)² + min((u[1, 1] − u[1, NY]) / dy, 0)²) for i ∈ [1, 1], j ∈ [NY, NY]
x = NX, y = 1boundary facesqrt(max((u[NX, 1] − u[NX−1, 1]) / dx, 0)² + min((u[1, 1] − u[NX, 1]) / dx, 0)² + max((u[NX, 1] − u[NX, NY]) / dy, 0)² + min((u[NX, 2] − u[NX, 1]) / dy, 0)²) for i ∈ [NX, NX], j ∈ [1, 1]
x = NX, y = NYboundary facesqrt(max((u[NX, NY] − u[NX−1, NY]) / dx, 0)² + min((u[1, NY] − u[NX, NY]) / dx, 0)² + max((u[NX, NY] − u[NX, NY−1]) / dy, 0)² + min((u[NX, 1] − u[NX, NY]) / dy, 0)²) for i ∈ [NX, NX], j ∈ [NY, NY]

First-order Godunov (Rouy-Tourin / Osher-Sethian) upwind discretization of the gradient magnitude |grad u| = sqrt(u_x^2 + u_y^2) for a non-negative-speed Hamilton-Jacobi term (psi_t = -S|grad psi|, S >= 0; the level-set / eikonal norm) on cartesian_uniform_2d, with a DOUBLY-PERIODIC boundary (both x and y wrapped as circles). This is the full 2-D scheme of archive/discretizations/finite_difference/godunov_norm_1st_uniform_cartesian_2d.json expressed in the D convention that grids/cartesian_uniform_2d already uses (central_D1x_, central_D2x_), and the 2-D companion of grids/cartesian_uniform_1d/rules/godunov_norm_D1_periodic (the 1-D reduction keeps only the x term; this rule sums the x AND y axes). COMPOUND MATCH (esm-spec 9.6.3 compound precedence): the rule matches the whole compound sqrt((D(u, wrt: x))^2 + (D(u, wrt: y))^2) rooted at the outer sqrt — cartesian’s plain-D first-derivative idiom squared and summed under a square root, the exact structural analogue of the archive’s applies_to sqrt(grad(u,x)^2 + grad(u,y)^2). This 2-D gradient magnitude is precisely the canonical compound-precedence example esm-spec 9.6.3 (and 9.6.8) names. It carries priority 10 (>= the D1/D2 rules’ 10) so that, under the section-9.6.3 outermost-first fixpoint, it fires on the entire sqrt(D_x^2 + D_y^2) term BEFORE any plain first-derivative rule (central_D1x_periodic / central_D1y_periodic) could lower the inner D(u, wrt: x) or D(u, wrt: y); the whole compound is replaced by this one makearray in a single step and the engine does not descend into the fresh body. Match-scoped with where: {u: {shape: [x, y]}} (esm-spec 9.6.1) so it fires only on the gradient magnitude of a bare field declared over this grid’s [x, y] index sets; under import-edge renaming (esm-spec 9.7.7) the two wrt literals and the where shape follow the renamed axes together, so two cartesian instances coexist scoped to their own grids. Because the constraint requires a bare shaped variable, a consumer forming the gradient magnitude of a COMPOUND inner expression must first bind it to a declared [x,y]-shaped observed. THE SCHEME: with per-axis backward differences D-_x = (u[i,j] - u[i-1,j]) / dx, D-_y = (u[i,j] - u[i,j-1]) / dy and forward differences D+_x = (u[i+1,j] - u[i,j]) / dx, D+_y = (u[i,j+1] - u[i,j]) / dy, the discrete norm is sqrt( max(D-_x,0)^2 + min(D+_x,0)^2 + max(D-_y,0)^2 + min(D+_y,0)^2 ). The per-axis max/min is the entropy (upwind) fix of Osher & Sethian (1988) and Rouy & Tourin (1992): for the non-negative-speed term it selects the one-sided difference along the incoming characteristic per axis per point (a centered difference is unstable for this hyperbolic operator). It is EXACT on a linear field: for u = a x + b y + c the interior gives D-_x = D+_x = a and D-_y = D+_y = b, so sqrt(max(a,0)^2 + min(a,0)^2 + max(b,0)^2 + min(b,0)^2) = sqrt(a^2 + b^2) = |grad u| to rounding, reproducing the continuous gradient magnitude with no truncation error (the low-order analogue of the null-space property that annihilates constants, a = b = 0 -> 0). First order in space, O(dx + dy), on smooth data; it retains first-order accuracy uniformly through the |grad| kinks (points where a gradient component changes sign, or where grad u vanishes) that the entropy fix is designed to resolve. The makearray tiles [x, y] as a 3x3 block partition {1}u[2,NX-1]u{NX} x {1}u[2,NY-1]u{NY} — a 9-region closure identical in shape to grids/cartesian_uniform_2d/rules/central_Dxy_periodic, because the Godunov 5-point stencil reaches i+-1 AND j+-1: the doubly-interior block [2,NX-1]x[2,NY-1] holds the Godunov-norm aggregate over the full 5-point stencil (u[i,j], u[i-1,j], u[i+1,j], u[i,j-1], u[i,j+1]); the four edge strips wrap their single out-of-range axis (i=1 uses the ghost column NX for D-_x, i=NX uses column 1 for D+_x; j=1 uses row NY for D-_y, j=NY uses row 1 for D+_y) with the in-range axis interior; the four corner cells wrap BOTH axes. Every cell evaluates the identical max/min entropy norm; the periodic wrap is exact on the torus, so the operator is self-contained and testable on a closed 2-torus. The edge and corner faces are FULL-RANK aggregates (singleton output range on the wrapped axis) mirroring central_Dxy_periodic so every binding’s simulator fills the strips and corner cells without rank promotion. A doubly-periodic linear ramp is a contradiction (a non-zero ramp cannot wrap), so the wrap faces carry the seam of a ramp rather than sqrt(a^2+b^2); the exactness property is the interior statement, and the periodic gate that ships with this rule uses the constant/interior form (problems/godunov_norm_2d_periodic_exact.esm). dx/dy are free names resolving to the consuming model’s spacing variables (grids/cartesian_uniform_2d/grid.esm consumer-supplied geometry contract); the periodic wrap is dx/dy-scale-invariant, so the rule serves any x/y extent unchanged. A model imports exactly one gradient-magnitude rule for the [x, y] pair.

References

  • Osher, S., & Sethian, J. A. (1988). Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79(1), 12-49. The max/min upwind (entropy) construction of |grad phi| for the level-set Hamilton-Jacobi term. doi:10.1016/0021-9991(88)90002-2
  • Rouy, E., & Tourin, A. (1992). A viscosity solutions approach to shape-from-shading. SIAM Journal on Numerical Analysis, 29(3), 867-884. The Godunov Hamiltonian sqrt(sum_d max(D-_d,0)^2 + min(D+_d,0)^2) as the monotone, first-order-convergent discretization of the eikonal |grad u|. doi:10.1137/0729053

Conformance fixtures

Convergence — case godunov_norm_2d_eikonal_mms

Error norms read from the committed golden tests/conformance/convergence/godunov_norm_2d_eikonal_mms/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
161.558e-013.467e-01
327.807e-021.001.742e-010.99
643.940e-020.998.720e-021.00
1281.986e-020.994.361e-021.00

Expected order: 1 (± 0.15).

Convergence of godunov_norm_2d_eikonal_mms (log-log)

varcoeff_laplacian_x_dirichlet_bc

op:varcoeff_laplacian order:2 bc:dirichlet axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/varcoeff_laplacian_x_dirichlet_bc.esm

Rewrites ∂(k·∂u/∂x)/∂x at priority 10 — match pattern:

{"op": "D", "args": [{"op": "*", "args": ["k", {"op": "D", "args": ["u"], "wrt": "x"}]}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]; k is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorvarcoeff_laplacian_x_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i+1, j]))·(u[i+1, j] − u[i, j]) − (0.5·(k[i−1, j] + k[i, j]))·(u[i, j] − u[i−1, j])) / dx² for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face((0.5·(k[1, j] + k[2, j]))·(u[2, j] − u[1, j]) − (2·k[1, j])·(u[1, j] − gL)) / dx² for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face((0.5·(k[NX−1, j] + k[NX, j]))·(u[NX−1, j] − u[NX, j]) + (2·k[NX, j])·(gR − u[NX, j])) / dx² for i ∈ [NX, NX], j ∈ y

The x-half of a 2-D conservative finite-volume VARIABLE-COEFFICIENT / NONLINEAR Laplacian div(k grad u) on cartesian_uniform_2d with inhomogeneous DIRICHLET (fixed-value) conditions on the two x-walls, over the full 2-D field u[x,y], cell-centered. The operator carries a coefficient, so a bare laplacian match cannot express it; the rule matches the CONSERVATIVE COMPOUND the consumer writes, D( {op:, args:[k, D(u,wrt:x)]}, wrt:x ) at priority 10, so under esm-spec 9.6.3 outermost-first compound precedence it fires on the whole compound before any plain-D rule could lower the inner derivative. This is the 2-D per-axis lift of grids/cartesian_uniform_1d/rules/varcoeff_laplacian_dirichlet_bc: the operated axis x gets an interior band plus two mirror-ghost x-faces, while the y axis is full-rank (all rows j treated identically). BOTH operands are match wildcards bound to consumer fields: u (the state) and k (the coefficient, the LEFT factor of the product); matching is structural and non-commutative, so the consumer MUST write k * D(u, wrt:x). Match-scoped with where {u:{shape:[x,y]}, k:{shape:[x,y]}} (esm-spec 9.6.1) so it fires only on this 2-D grid’s declared [x,y]-shaped fields; under import-edge renaming (esm-spec 9.7.7) the wrt literals and both where shapes follow the renamed axes together. ONE rule serves the VARIABLE-COEFFICIENT case (k=k(x,y)) and the NONLINEAR case (k=f(u)); the body only reads k by index at the faces. AUTHORING CONTRACT for k (identical to varcoeff_laplacian_x_zero_flux_bc): k is a declared [x,y]-shaped field; an aggregate-observed coefficient MUST loop over indices DISTINCT from the stencil’s interior indices i,j (the library convention is gi/gj), because it is inlined into the stencil’s [2,NX-1]x[1,NY] aggregate and a shared loop-index name is captured. Face conductivity is the ARITHMETIC MEAN kf[i+1/2,j]=(k[i,j]+k[i+1,j])/2 (O(dx^2) at the face for a differentiable coefficient, affine in k so it composes with the nonlinear dual). Interior i in [2,NX-1] (all j): the imported conservative x-flux-difference stencil. Faces i=1 and i=NX: mirror-ghost elimination against a prescribed wall value combined with the boundary-face conductivity. On this cell-centered grid the wall sits half a cell outside the first/last x-center, so linear reconstruction to the wall gives u[0,j]=2gL-u[1,j] (and u[NX+1,j]=2gR-u[NX,j]); the left x-face flux F[1/2,j]=kf[1/2,j](u[1,j]-u[0,j])/dx uses the wall conductivity k[1,j] (the arithmetic-mean face value degenerates to the first cell value under the even ghost k[0,j]=k[1,j]), so F[1/2,j]=k[1,j]2(u[1,j]-gL)/dx and Lx u[1,j] = ( (k[1,j]+k[2,j])/2 (u[2,j]-u[1,j]) - 2 k[1,j] (u[1,j]-gL) ) / dx^2; symmetrically at the right x-wall F[NX+1/2,j]=k[NX,j]2(gR-u[NX,j])/dx and Lx u[NX,j] = ( (k[NX-1,j]+k[NX,j])/2 (u[NX-1,j]-u[NX,j]) + 2 k[NX,j] (gR-u[NX,j]) ) / dx^2. gL and gR are consumer-supplied free names (the same free-name contract as the grid geometry dx; see grids/cartesian_uniform_2d/grid.esm) giving the field value on the left x-wall x=x0 and the right x-wall x=x0+NXdx, constant along that wall (all j); a consuming model declares them as ordinary real parameters, and a parameter not overridden falls back to its default (esm-spec 6: default, or 0), so gL=gR=0 is the homogeneous-Dirichlet x-wall. The gL/gR terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring central_D2x_dirichlet_bc. dx is the grid’s consumer-supplied x-spacing free name. The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(k*D(u,y),y): this rule lowers the first term. Reduces to central_D2x_dirichlet_bc when k is spatially constant (k=1 gives the mirror-ghost residuals (u[2,j]-3u[1,j]+2gL)/dx^2 and (u[NX-1,j]-3u[NX,j]+2gR)/dx^2). Interior O(dx^2); the one-sided mirror-ghost faces are O(dx) locally, preserving global 2nd-order L2 convergence for Dirichlet-compatible solutions. A model imports exactly one second-derivative/Laplacian rule for x.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient conservative form) and Section 2.12 (cell-centered Dirichlet via mirror ghost).

Conformance fixtures

Convergence — case heat_2d_varcoeff_x_dirichlet

Error norms read from the committed golden tests/conformance/convergence/heat_2d_varcoeff_x_dirichlet/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.366e-042.108e-04
325.687e-052.065.193e-052.02
641.408e-052.011.290e-052.01
1283.513e-062.003.218e-062.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_varcoeff_x_dirichlet (log-log)

varcoeff_laplacian_x_neumann_bc

op:varcoeff_laplacian order:2 bc:neumann axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/varcoeff_laplacian_x_neumann_bc.esm

Rewrites ∂(k·∂u/∂x)/∂x at priority 10 — match pattern:

{"op": "D", "args": [{"op": "*", "args": ["k", {"op": "D", "args": ["u"], "wrt": "x"}]}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]; k is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorvarcoeff_laplacian_x_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i+1, j]))·(u[i+1, j] − u[i, j]) − (0.5·(k[i−1, j] + k[i, j]))·(u[i, j] − u[i−1, j])) / dx² for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face((0.5·(k[1, j] + k[2, j]))·(u[2, j] − u[1, j])) / dx² − qXL/dx for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face((0.5·(k[NX−1, j] + k[NX, j]))·(u[NX−1, j] − u[NX, j])) / dx² + qXR/dx for i ∈ [NX, NX], j ∈ y

The x-half of a 2-D conservative finite-volume VARIABLE-COEFFICIENT / NONLINEAR Laplacian div(k grad u) on cartesian_uniform_2d with inhomogeneous NEUMANN (fixed-flux) conditions on the two x-walls, over the full 2-D field u[x,y], cell-centered. The operator carries a coefficient, so a bare laplacian match cannot express it; the rule matches the CONSERVATIVE COMPOUND the consumer writes, D( {op:, args:[k, D(u,wrt:x)]}, wrt:x ) at priority 10, so under esm-spec 9.6.3 outermost-first compound precedence it fires on the whole compound before any plain-D rule could lower the inner derivative. This is the 2-D per-axis lift of grids/cartesian_uniform_1d/rules/varcoeff_laplacian_neumann_bc: the operated axis x gets an interior band plus two prescribed-flux x-faces, while the y axis is full-rank (all rows j treated identically). BOTH operands are match wildcards bound to consumer fields: u (the state) and k (the coefficient, the LEFT factor of the product); matching is structural and non-commutative, so the consumer MUST write k * D(u, wrt:x). Match-scoped with where {u:{shape:[x,y]}, k:{shape:[x,y]}} (esm-spec 9.6.1) so it fires only on this 2-D grid’s declared [x,y]-shaped fields; under import-edge renaming (esm-spec 9.7.7) the wrt literals and both where shapes follow the renamed axes together. ONE rule serves the VARIABLE-COEFFICIENT case (k=k(x,y)) and the NONLINEAR case (k=f(u)); the body only reads k by index at the faces. AUTHORING CONTRACT for k (identical to varcoeff_laplacian_x_zero_flux_bc): k is a declared [x,y]-shaped field; an aggregate-observed coefficient MUST loop over indices DISTINCT from the stencil’s interior indices i,j (the library convention is gi/gj), because it is inlined into the stencil’s [2,NX-1]x[1,NY] aggregate and a shared loop-index name is captured. Face conductivity is the ARITHMETIC MEAN kf[i+1/2,j]=(k[i,j]+k[i+1,j])/2. Interior i in [2,NX-1] (all j): the imported conservative x-flux-difference stencil. Faces i=1 and i=NX: the domain-boundary x-face flux is PRESCRIBED DIRECTLY – the natural Neumann treatment for the finite-volume form, requiring no ghost geometry and no boundary conductivity. qXL and qXR are consumer-supplied free names (the same free-name contract as the grid geometry dx; see grids/cartesian_uniform_2d/grid.esm) giving the PHYSICAL FLUX (k du/dx, in the +x direction) at the left x-wall x=x0 and the right x-wall x=x0+NXdx respectively, constant along that wall (all j). The left cell balance Lx u[1,j] = (F[3/2,j]-F[1/2,j])/dx with F[1/2,j]=qXL gives Lx u[1,j] = (k[1,j]+k[2,j])/2 (u[2,j]-u[1,j]) / dx^2 - qXL/dx; the right cell balance Lx u[NX,j] = (F[NX+1/2,j]-F[NX-1/2,j])/dx with F[NX+1/2,j]=qXR gives Lx u[NX,j] = (k[NX-1,j]+k[NX,j])/2 (u[NX-1,j]-u[NX,j]) / dx^2 + qXR/dx. A consuming model declares qXL/qXR as ordinary real parameters; a parameter not overridden falls back to its default (esm-spec 6: default, or 0), so the homogeneous case is the zero default and this rule REDUCES EXACTLY to varcoeff_laplacian_x_zero_flux_bc when qXL=qXR=0 – the prescribed-flux terms qXL/dx and qXR/dx vanish and the two face values become byte-identical to that rule’s dropped-flux faces (kf[3/2,j](u[2,j]-u[1,j])/dx^2 at i=1, kf[NX-1/2,j](u[NX-1,j]-u[NX,j])/dx^2 at i=NX), the flux at each x-wall dropped to zero. The qXL/qXR terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring central_D2x_neumann_bc. dx is the grid’s consumer-supplied x-spacing free name. The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(kD(u,y),y): this rule lowers the first term. Interior O(dx^2); the prescribed-flux faces impose the exact boundary flux, preserving global 2nd-order L2 convergence for Neumann-compatible solutions. A model imports exactly one second-derivative/Laplacian rule for x.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient conservative form) and Section 2.12 (Neumann boundary conditions by prescribing the domain-face flux).

Conformance fixtures

Convergence — case heat_2d_varcoeff_x_neumann

Error norms read from the committed golden tests/conformance/convergence/heat_2d_varcoeff_x_neumann/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
161.293e-038.851e-04
323.147e-042.042.165e-042.03
647.814e-052.015.354e-052.02
1281.950e-052.001.338e-052.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_varcoeff_x_neumann (log-log)

varcoeff_laplacian_x_periodic

op:varcoeff_laplacian order:2 bc:periodic axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/varcoeff_laplacian_x_periodic.esm

Rewrites ∂(k·∂u/∂x)/∂x at priority 10 — match pattern:

{"op": "D", "args": [{"op": "*", "args": ["k", {"op": "D", "args": ["u"], "wrt": "x"}]}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]; k is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorvarcoeff_laplacian_x_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i+1, j]))·(u[i+1, j] − u[i, j]) − (0.5·(k[i−1, j] + k[i, j]))·(u[i, j] − u[i−1, j])) / dx² for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face((0.5·(k[1, j] + k[2, j]))·(u[2, j] − u[1, j]) − (0.5·(k[NX, j] + k[1, j]))·(u[1, j] − u[NX, j])) / dx² for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face((0.5·(k[NX, j] + k[1, j]))·(u[1, j] − u[NX, j]) − (0.5·(k[NX−1, j] + k[NX, j]))·(u[NX, j] − u[NX−1, j])) / dx² for i ∈ [NX, NX], j ∈ y

The x-half of a 2-D conservative finite-volume VARIABLE-COEFFICIENT / NONLINEAR Laplacian div(k grad u) on cartesian_uniform_2d with a PERIODIC x-boundary, over the full 2-D field u[x,y]. The operator carries a coefficient, so a bare laplacian match cannot express it; the rule matches the CONSERVATIVE COMPOUND the consumer writes, D( {op:, args:[k, D(u,wrt:x)]}, wrt:x ) – outer-D of (k times the inner x-gradient) – at priority 10, so under esm-spec 9.6.3 outermost-first compound precedence it fires on the whole compound before any plain-D rule could lower the inner derivative (exactly as central_D2x_periodic fires on D(D(u,x),x)). This is the 2-D per-axis lift of grids/cartesian_uniform_1d/rules/varcoeff_laplacian_periodic: the operated axis x gets an interior band plus two wrapped x-faces, while the y axis is full-rank (all rows j treated identically). BOTH operands are match wildcards bound to consumer fields: u (the state, inner D’s argument) and k (the coefficient, the LEFT factor of the product); matching is structural and non-commutative, so the consumer MUST write the coefficient on the LEFT: k * D(u, wrt:x). Match-scoped with where {u:{shape:[x,y]}, k:{shape:[x,y]}} (esm-spec 9.6.1) so it fires only on the xx-conservative-compound of bare fields declared over this 2-D grid’s [x,y] axes (distinct from the 1-D shape-[x] varcoeff rules); under import-edge renaming (esm-spec 9.7.7) the wrt literals and both where shapes follow the renamed axes together, so two instances coexist scoped to their own meshes. ONE rule serves two duals: the VARIABLE-COEFFICIENT case (k=k(x,y), a given periodic field) and the NONLINEAR case (k declared as an observed k=f(u)); the body only reads k by index at the faces, so a state-dependent k just makes the expanded stencil nonlinear in u with no reformulation (the keyed-factor field contract, like the grid’s dx). AUTHORING CONTRACT for k: it is a declared [x,y]-shaped field; when supplied as an aggregate observed it MUST loop over indices DISTINCT from the stencil’s interior indices i,j (the library convention is gi/gj), because the coefficient observed is inlined into the stencil’s [2,NX-1]x[1,NY] aggregate at evaluation and a shared loop-index name is captured (silently wrong interior values). Face conductivity is the ARITHMETIC MEAN kf[i+1/2,j]=(k[i,j]+k[i+1,j])/2 (O(dx^2) at the face for a differentiable coefficient, and affine in k so it composes with the nonlinear dual). Interior i in [2,NX-1] (all j): the imported conservative x-flux-difference stencil. The two x-face columns i=1 and i=NX carry the same conservative x-flux difference with the out-of-range x-neighbor AND its coefficient wrapped periodically: at i=1 the left x-face flux uses the ghost u[0,j]=u[NX,j] and k[0,j]=k[NX,j], giving kf[1/2,j]=(k[NX,j]+k[1,j])/2 and Lx u[1,j] = ( (k[1,j]+k[2,j])/2 (u[2,j]-u[1,j]) - (k[NX,j]+k[1,j])/2 (u[1,j]-u[NX,j]) ) / dx^2; at i=NX the right x-face flux uses u[NX+1,j]=u[1,j] and k[NX+1,j]=k[1,j], giving kf[NX+1/2,j]=(k[NX,j]+k[1,j])/2 and Lx u[NX,j] = ( (k[NX,j]+k[1,j])/2 (u[1,j]-u[NX,j]) - (k[NX-1,j]+k[NX,j])/2 (u[NX,j]-u[NX-1,j]) ) / dx^2. The SAME wrapped x-face flux (k[NX,j]+k[1,j])/2 (u[1,j]-u[NX,j])/dx is the left face of column 1 and the right face of column NX, so the discrete x-operator is exactly conservative on each row (the circulant variable-coefficient x-Laplacian) and its row sums vanish. The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring central_D2x_neumann_bc so every binding’s simulator fills the (1 x NY) columns without rank promotion. dx is the grid’s consumer-supplied x-spacing free name; the periodic wrap is dx-scale-invariant, so the rule serves any x-extent unchanged. The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(k*D(u,y),y): this rule lowers the first term. Reduces to central_D2x_periodic when k is spatially constant. Interior and periodic faces are all O(dx^2), so the whole x-operator is 2nd order with no boundary reduction. A model imports exactly one second-derivative/Laplacian rule for x.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient conservative form) and Section 10.4 (centered differences on periodic domains).

Conformance fixtures

Convergence — case heat_2d_varcoeff_x_periodic

Error norms read from the committed golden tests/conformance/convergence/heat_2d_varcoeff_x_periodic/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
164.046e-028.114e-04
329.975e-032.022.028e-042.00
642.485e-032.015.069e-052.00
1286.207e-042.001.267e-052.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_varcoeff_x_periodic (log-log)

varcoeff_laplacian_x_zero_flux_bc

op:varcoeff_laplacian order:2 bc:zero_flux axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/varcoeff_laplacian_x_zero_flux_bc.esm

Rewrites ∂(k·∂u/∂x)/∂x at priority 10 — match pattern:

{"op": "D", "args": [{"op": "*", "args": ["k", {"op": "D", "args": ["u"], "wrt": "x"}]}], "wrt": "x"}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]; k is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [2, NX−1], y ∈ [1, NY]interiorvarcoeff_laplacian_x_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i+1, j]))·(u[i+1, j] − u[i, j]) − (0.5·(k[i−1, j] + k[i, j]))·(u[i, j] − u[i−1, j])) / dx² for i ∈ [2, NX−1], j ∈ y
x = 1, y ∈ [1, NY]boundary face((0.5·(k[1, j] + k[2, j]))·(u[2, j] − u[1, j])) / dx² for i ∈ [1, 1], j ∈ y
x = NX, y ∈ [1, NY]boundary face((0.5·(k[NX−1, j] + k[NX, j]))·(u[NX−1, j] − u[NX, j])) / dx² for i ∈ [NX, NX], j ∈ y

The x-half of a 2-D conservative finite-volume VARIABLE-COEFFICIENT / NONLINEAR Laplacian div(k grad u) on cartesian_uniform_2d with zero-flux (homogeneous Neumann) conditions on the two x-walls, over the full 2-D field u[x,y], cell-centered. The operator carries a coefficient, so a bare laplacian match cannot express it; the rule matches the CONSERVATIVE COMPOUND the consumer writes, D( {op:, args:[k, D(u,wrt:x)]}, wrt:x ) at priority 10, so under esm-spec 9.6.3 outermost-first compound precedence it fires on the whole compound before any plain-D rule could lower the inner derivative (exactly as central_D2x fires on D(D(u,x),x)). This is the 2-D per-axis lift of grids/cartesian_uniform_1d/rules/varcoeff_laplacian_zero_flux_bc: the operated axis x gets an interior band plus two dropped-flux x-faces, while the y axis is full-rank (all rows j treated identically). BOTH operands are match wildcards bound to consumer fields: u (the state) and k (the coefficient, the LEFT factor of the product); matching is structural and non-commutative, so the consumer MUST write k * D(u, wrt:x). Match-scoped with where {u:{shape:[x,y]}, k:{shape:[x,y]}} (esm-spec 9.6.1) so it fires only on this 2-D grid’s declared [x,y]-shaped fields; under import-edge renaming (esm-spec 9.7.7) the wrt literals and both where shapes follow the renamed axes together. ONE rule serves the VARIABLE-COEFFICIENT case (k=k(x,y)) and the NONLINEAR case (k=f(u)); the body only reads k by index at the faces. AUTHORING CONTRACT for k: k is a declared [x,y]-shaped field; an aggregate-observed coefficient MUST loop over indices DISTINCT from the stencil’s interior indices i,j (the library convention is gi/gj), because it is inlined into the stencil’s [2,NX-1]x[1,NY] aggregate and a shared loop-index name is captured (silently wrong interior values). Face conductivity is the ARITHMETIC MEAN kf[i+1/2,j]=(k[i,j]+k[i+1,j])/2 (O(dx^2) at the face for a differentiable coefficient, affine in k so it composes with the nonlinear dual). Interior i in [2,NX-1] (all j): the imported conservative x-flux-difference stencil. The x-boundary faces impose zero flux exactly (natural for the finite-volume form – no ghost geometry): at i=1 the left x-domain-face flux is dropped, giving Lx u[1,j] = kf[3/2,j](u[2,j]-u[1,j])/dx^2 = (k[1,j]+k[2,j])/2 (u[2,j]-u[1,j])/dx^2; at i=NX the right x-domain-face flux is dropped, giving Lx u[NX,j] = -kf[NX-1/2,j](u[NX,j]-u[NX-1,j])/dx^2 = (k[NX-1,j]+k[NX,j])/2 (u[NX-1,j]-u[NX,j])/dx^2. This is exactly varcoeff_laplacian_x_neumann_bc at qXL=qXR=0. The face values are FULL-RANK aggregates (singleton output range on the x axis, all j) mirroring central_D2x_neumann_bc. dx is the grid’s consumer-supplied x-spacing free name; the zero-flux closure is dx-scale-invariant, so the rule serves any x-extent unchanged. The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(k*D(u,y),y): this rule lowers the first term. Reduces to the constant-coefficient zero-gradient x-half when k is spatially constant. Interior O(dx^2); the one-sided zero-flux faces are O(dx) locally, which preserves global 2nd-order L2 convergence for zero-flux-compatible solutions (flux k u_x vanishing at both x-walls). A model imports exactly one second-derivative/Laplacian rule for x.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient conservative form) and Section 2.12 (Neumann boundary conditions by dropping the domain-face flux).

Conformance fixtures

Convergence — case heat_2d_varcoeff_x_zero_flux

Error norms read from the committed golden tests/conformance/convergence/heat_2d_varcoeff_x_zero_flux/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
166.726e-032.584e-03
321.692e-031.996.522e-041.99
644.236e-042.001.630e-042.00
1281.060e-042.004.078e-052.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_varcoeff_x_zero_flux (log-log)

varcoeff_laplacian_y_dirichlet_bc

op:varcoeff_laplacian order:2 bc:dirichlet axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/varcoeff_laplacian_y_dirichlet_bc.esm

Rewrites ∂(k·∂u/∂y)/∂y at priority 10 — match pattern:

{"op": "D", "args": [{"op": "*", "args": ["k", {"op": "D", "args": ["u"], "wrt": "y"}]}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]; k is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorvarcoeff_laplacian_y_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i, j+1]))·(u[i, j+1] − u[i, j]) − (0.5·(k[i, j−1] + k[i, j]))·(u[i, j] − u[i, j−1])) / dy² for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face((0.5·(k[i, 1] + k[i, 2]))·(u[i, 2] − u[i, 1]) − (2·k[i, 1])·(u[i, 1] − gB)) / dy² for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face((0.5·(k[i, NY−1] + k[i, NY]))·(u[i, NY−1] − u[i, NY]) + (2·k[i, NY])·(gT − u[i, NY])) / dy² for i ∈ x, j ∈ [NY, NY]

The y-half of a 2-D conservative finite-volume VARIABLE-COEFFICIENT / NONLINEAR Laplacian div(k grad u) on cartesian_uniform_2d with inhomogeneous DIRICHLET (fixed-value) conditions on the two y-walls, over the full 2-D field u[x,y], cell-centered. The operator carries a coefficient, so a bare laplacian match cannot express it; the rule matches the CONSERVATIVE COMPOUND the consumer writes, D( {op:, args:[k, D(u,wrt:y)]}, wrt:y ) at priority 10, so under esm-spec 9.6.3 outermost-first compound precedence it fires on the whole compound before any plain-D rule could lower the inner derivative. This is the 2-D per-axis lift of grids/cartesian_uniform_1d/rules/varcoeff_laplacian_dirichlet_bc: the operated axis y gets an interior band plus two mirror-ghost y-faces, while the x axis is full-rank (all columns i treated identically). BOTH operands are match wildcards bound to consumer fields: u (the state) and k (the coefficient, the LEFT factor of the product); matching is structural and non-commutative, so the consumer MUST write k * D(u, wrt:y). Match-scoped with where {u:{shape:[x,y]}, k:{shape:[x,y]}} (esm-spec 9.6.1) so it fires only on this 2-D grid’s declared [x,y]-shaped fields; under import-edge renaming (esm-spec 9.7.7) the wrt literals and both where shapes follow the renamed axes together. ONE rule serves the VARIABLE-COEFFICIENT case (k=k(x,y)) and the NONLINEAR case (k=f(u)); the body only reads k by index at the faces. AUTHORING CONTRACT for k (identical to varcoeff_laplacian_y_zero_flux_bc): k is a declared [x,y]-shaped field; an aggregate-observed coefficient MUST loop over indices DISTINCT from the stencil’s interior indices i,j (the library convention is gi/gj), because it is inlined into the stencil’s [1,NX]x[2,NY-1] aggregate and a shared loop-index name is captured. Face conductivity is the ARITHMETIC MEAN kf[i,j+1/2]=(k[i,j]+k[i,j+1])/2 (O(dy^2) at the face for a differentiable coefficient, affine in k so it composes with the nonlinear dual). Interior j in [2,NY-1] (all i): the imported conservative y-flux-difference stencil. Faces j=1 and j=NY: mirror-ghost elimination against a prescribed wall value combined with the boundary-face conductivity. On this cell-centered grid the wall sits half a cell outside the first/last y-center, so linear reconstruction to the wall gives u[i,0]=2gB-u[i,1] (and u[i,NY+1]=2gT-u[i,NY]); the bottom y-face flux F[i,1/2]=kf[i,1/2](u[i,1]-u[i,0])/dy uses the wall conductivity k[i,1] (the arithmetic-mean face value degenerates to the first cell value under the even ghost k[i,0]=k[i,1]), so F[i,1/2]=k[i,1]2(u[i,1]-gB)/dy and Ly u[i,1] = ( (k[i,1]+k[i,2])/2 (u[i,2]-u[i,1]) - 2 k[i,1] (u[i,1]-gB) ) / dy^2; symmetrically at the top y-wall F[i,NY+1/2]=k[i,NY]2(gT-u[i,NY])/dy and Ly u[i,NY] = ( (k[i,NY-1]+k[i,NY])/2 (u[i,NY-1]-u[i,NY]) + 2 k[i,NY] (gT-u[i,NY]) ) / dy^2. gB and gT are consumer-supplied free names (the same free-name contract as the grid geometry dy; see grids/cartesian_uniform_2d/grid.esm) giving the field value on the BOTTOM y-wall y=y0 (the j=1 side) and the TOP y-wall y=y0+NYdy (the j=NY side), constant along that wall (all i). The names deliberately read gB/gT (bottom/top) so the y-wall assignment is unambiguous, the y-analogue of gL/gR on the x-walls. A consuming model declares them as ordinary real parameters, and a parameter not overridden falls back to its default (esm-spec 6: default, or 0), so gB=gT=0 is the homogeneous-Dirichlet y-wall. The gB/gT terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring central_D2y_dirichlet_bc. dy is the grid’s consumer-supplied y-spacing free name. The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(k*D(u,y),y): this rule lowers the second term. Reduces to central_D2y_dirichlet_bc when k is spatially constant (k=1 gives the mirror-ghost residuals (u[i,2]-3u[i,1]+2gB)/dy^2 and (u[i,NY-1]-3u[i,NY]+2gT)/dy^2). Interior O(dy^2); the one-sided mirror-ghost faces are O(dy) locally, preserving global 2nd-order L2 convergence for Dirichlet-compatible solutions. A model imports exactly one second-derivative/Laplacian rule for y.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient conservative form) and Section 2.12 (cell-centered Dirichlet via mirror ghost).

Conformance fixtures

Convergence — case heat_2d_varcoeff_y_dirichlet

Error norms read from the committed golden tests/conformance/convergence/heat_2d_varcoeff_y_dirichlet/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
162.366e-042.108e-04
325.687e-052.065.193e-052.02
641.408e-052.011.290e-052.01
1283.513e-062.003.218e-062.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_varcoeff_y_dirichlet (log-log)

varcoeff_laplacian_y_neumann_bc

op:varcoeff_laplacian order:2 bc:neumann axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/varcoeff_laplacian_y_neumann_bc.esm

Rewrites ∂(k·∂u/∂y)/∂y at priority 10 — match pattern:

{"op": "D", "args": [{"op": "*", "args": ["k", {"op": "D", "args": ["u"], "wrt": "y"}]}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]; k is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorvarcoeff_laplacian_y_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i, j+1]))·(u[i, j+1] − u[i, j]) − (0.5·(k[i, j−1] + k[i, j]))·(u[i, j] − u[i, j−1])) / dy² for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face((0.5·(k[i, 1] + k[i, 2]))·(u[i, 2] − u[i, 1])) / dy² − qYL/dy for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face((0.5·(k[i, NY−1] + k[i, NY]))·(u[i, NY−1] − u[i, NY])) / dy² + qYR/dy for i ∈ x, j ∈ [NY, NY]

The y-half of a 2-D conservative finite-volume VARIABLE-COEFFICIENT / NONLINEAR Laplacian div(k grad u) on cartesian_uniform_2d with inhomogeneous NEUMANN (fixed-flux) conditions on the two y-walls, over the full 2-D field u[x,y], cell-centered. The operator carries a coefficient, so a bare laplacian match cannot express it; the rule matches the CONSERVATIVE COMPOUND the consumer writes, D( {op:, args:[k, D(u,wrt:y)]}, wrt:y ) at priority 10, so under esm-spec 9.6.3 outermost-first compound precedence it fires on the whole compound before any plain-D rule could lower the inner derivative. This is the 2-D per-axis lift of grids/cartesian_uniform_1d/rules/varcoeff_laplacian_neumann_bc: the operated axis y gets an interior band plus two prescribed-flux y-faces, while the x axis is full-rank (all columns i treated identically). BOTH operands are match wildcards bound to consumer fields: u (the state) and k (the coefficient, the LEFT factor of the product); matching is structural and non-commutative, so the consumer MUST write k * D(u, wrt:y). Match-scoped with where {u:{shape:[x,y]}, k:{shape:[x,y]}} (esm-spec 9.6.1) so it fires only on this 2-D grid’s declared [x,y]-shaped fields; under import-edge renaming (esm-spec 9.7.7) the wrt literals and both where shapes follow the renamed axes together. ONE rule serves the VARIABLE-COEFFICIENT case (k=k(x,y)) and the NONLINEAR case (k=f(u)); the body only reads k by index at the faces. AUTHORING CONTRACT for k (identical to varcoeff_laplacian_y_zero_flux_bc): k is a declared [x,y]-shaped field; an aggregate-observed coefficient MUST loop over indices DISTINCT from the stencil’s interior indices i,j (the library convention is gi/gj), because it is inlined into the stencil’s [1,NX]x[2,NY-1] aggregate and a shared loop-index name is captured. Face conductivity is the ARITHMETIC MEAN kf[i,j+1/2]=(k[i,j]+k[i,j+1])/2. Interior j in [2,NY-1] (all i): the imported conservative y-flux-difference stencil. Faces j=1 and j=NY: the domain-boundary y-face flux is PRESCRIBED DIRECTLY – the natural Neumann treatment for the finite-volume form, requiring no ghost geometry and no boundary conductivity. qYL and qYR are consumer-supplied free names (the same free-name contract as the grid geometry dy; see grids/cartesian_uniform_2d/grid.esm) giving the PHYSICAL FLUX (k du/dy, in the +y direction) at the bottom y-wall y=y0 and the top y-wall y=y0+NYdy respectively, constant along that wall (all i). The bottom cell balance Ly u[i,1] = (F[i,3/2]-F[i,1/2])/dy with F[i,1/2]=qYL gives Ly u[i,1] = (k[i,1]+k[i,2])/2 (u[i,2]-u[i,1]) / dy^2 - qYL/dy; the top cell balance Ly u[i,NY] = (F[i,NY+1/2]-F[i,NY-1/2])/dy with F[i,NY+1/2]=qYR gives Ly u[i,NY] = (k[i,NY-1]+k[i,NY])/2 (u[i,NY-1]-u[i,NY]) / dy^2 + qYR/dy. A consuming model declares qYL/qYR as ordinary real parameters; a parameter not overridden falls back to its default (esm-spec 6: default, or 0), so the homogeneous case is the zero default and this rule REDUCES EXACTLY to varcoeff_laplacian_y_zero_flux_bc when qYL=qYR=0 – the prescribed-flux terms qYL/dy and qYR/dy vanish and the two face values become byte-identical to that rule’s dropped-flux faces (kf[i,3/2](u[i,2]-u[i,1])/dy^2 at j=1, kf[i,NY-1/2](u[i,NY-1]-u[i,NY])/dy^2 at j=NY), the flux at each y-wall dropped to zero. The qYL/qYR terms live only in the runtime rule body, never in build-time ic/reference positions (AGENTS.md build-time scope caveat). The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring central_D2y_neumann_bc. dy is the grid’s consumer-supplied y-spacing free name. The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(kD(u,y),y): this rule lowers the second term. Interior O(dy^2); the prescribed-flux faces impose the exact boundary flux, preserving global 2nd-order L2 convergence for Neumann-compatible solutions. A model imports exactly one second-derivative/Laplacian rule for y.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient conservative form) and Section 2.12 (Neumann boundary conditions by prescribing the domain-face flux).

Conformance fixtures

Convergence — case heat_2d_varcoeff_y_neumann

Error norms read from the committed golden tests/conformance/convergence/heat_2d_varcoeff_y_neumann/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
161.293e-038.851e-04
323.147e-042.042.165e-042.03
647.814e-052.015.354e-052.02
1281.950e-052.001.338e-052.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_varcoeff_y_neumann (log-log)

varcoeff_laplacian_y_periodic

op:varcoeff_laplacian order:2 bc:periodic axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/varcoeff_laplacian_y_periodic.esm

Rewrites ∂(k·∂u/∂y)/∂y at priority 10 — match pattern:

{"op": "D", "args": [{"op": "*", "args": ["k", {"op": "D", "args": ["u"], "wrt": "y"}]}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]; k is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorvarcoeff_laplacian_y_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i, j+1]))·(u[i, j+1] − u[i, j]) − (0.5·(k[i, j−1] + k[i, j]))·(u[i, j] − u[i, j−1])) / dy² for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face((0.5·(k[i, 1] + k[i, 2]))·(u[i, 2] − u[i, 1]) − (0.5·(k[i, NY] + k[i, 1]))·(u[i, 1] − u[i, NY])) / dy² for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face((0.5·(k[i, NY] + k[i, 1]))·(u[i, 1] − u[i, NY]) − (0.5·(k[i, NY−1] + k[i, NY]))·(u[i, NY] − u[i, NY−1])) / dy² for i ∈ x, j ∈ [NY, NY]

The y-half of a 2-D conservative finite-volume VARIABLE-COEFFICIENT / NONLINEAR Laplacian div(k grad u) on cartesian_uniform_2d with a PERIODIC y-boundary, over the full 2-D field u[x,y]. The operator carries a coefficient, so a bare laplacian match cannot express it; the rule matches the CONSERVATIVE COMPOUND the consumer writes, D( {op:, args:[k, D(u,wrt:y)]}, wrt:y ) – outer-D of (k times the inner y-gradient) – at priority 10, so under esm-spec 9.6.3 outermost-first compound precedence it fires on the whole compound before any plain-D rule could lower the inner derivative (exactly as central_D2y_periodic fires on D(D(u,y),y)). This is the 2-D per-axis lift of grids/cartesian_uniform_1d/rules/varcoeff_laplacian_periodic: the operated axis y gets an interior band plus two wrapped y-faces, while the x axis is full-rank (all columns i treated identically). BOTH operands are match wildcards bound to consumer fields: u (the state, inner D’s argument) and k (the coefficient, the LEFT factor of the product); matching is structural and non-commutative, so the consumer MUST write the coefficient on the LEFT: k * D(u, wrt:y). Match-scoped with where {u:{shape:[x,y]}, k:{shape:[x,y]}} (esm-spec 9.6.1) so it fires only on the yy-conservative-compound of bare fields declared over this 2-D grid’s [x,y] axes; under import-edge renaming (esm-spec 9.7.7) the wrt literals and both where shapes follow the renamed axes together, so two instances coexist scoped to their own meshes. ONE rule serves two duals: the VARIABLE-COEFFICIENT case (k=k(x,y), a given periodic field) and the NONLINEAR case (k declared as an observed k=f(u)); the body only reads k by index at the faces, so a state-dependent k just makes the expanded stencil nonlinear in u with no reformulation (the keyed-factor field contract, like the grid’s dy). AUTHORING CONTRACT for k: it is a declared [x,y]-shaped field; when supplied as an aggregate observed it MUST loop over indices DISTINCT from the stencil’s interior indices i,j (the library convention is gi/gj), because the coefficient observed is inlined into the stencil’s [1,NX]x[2,NY-1] aggregate at evaluation and a shared loop-index name is captured (silently wrong interior values). Face conductivity is the ARITHMETIC MEAN kf[i,j+1/2]=(k[i,j]+k[i,j+1])/2 (O(dy^2) at the face for a differentiable coefficient, and affine in k so it composes with the nonlinear dual). Interior j in [2,NY-1] (all i): the imported conservative y-flux-difference stencil. The two y-face rows j=1 and j=NY carry the same conservative y-flux difference with the out-of-range y-neighbor AND its coefficient wrapped periodically: at j=1 the bottom y-face flux uses the ghost u[i,0]=u[i,NY] and k[i,0]=k[i,NY], giving kf[i,1/2]=(k[i,NY]+k[i,1])/2 and Ly u[i,1] = ( (k[i,1]+k[i,2])/2 (u[i,2]-u[i,1]) - (k[i,NY]+k[i,1])/2 (u[i,1]-u[i,NY]) ) / dy^2; at j=NY the top y-face flux uses u[i,NY+1]=u[i,1] and k[i,NY+1]=k[i,1], giving kf[i,NY+1/2]=(k[i,NY]+k[i,1])/2 and Ly u[i,NY] = ( (k[i,NY]+k[i,1])/2 (u[i,1]-u[i,NY]) - (k[i,NY-1]+k[i,NY])/2 (u[i,NY]-u[i,NY-1]) ) / dy^2. The SAME wrapped y-face flux (k[i,NY]+k[i,1])/2 (u[i,1]-u[i,NY])/dy is the bottom face of row 1 and the top face of row NY, so the discrete y-operator is exactly conservative on each column (the circulant variable-coefficient y-Laplacian) and its column sums vanish. The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring central_D2y_neumann_bc so every binding’s simulator fills the (NX x 1) rows without rank promotion. dy is the grid’s consumer-supplied y-spacing free name; the periodic wrap is dy-scale-invariant, so the rule serves any y-extent unchanged. The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(k*D(u,y),y): this rule lowers the second term. Reduces to central_D2y_periodic when k is spatially constant. Interior and periodic faces are all O(dy^2), so the whole y-operator is 2nd order with no boundary reduction. A model imports exactly one second-derivative/Laplacian rule for y.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient conservative form) and Section 10.4 (centered differences on periodic domains).

Conformance fixtures

Convergence — case heat_2d_varcoeff_y_periodic

Error norms read from the committed golden tests/conformance/convergence/heat_2d_varcoeff_y_periodic/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
164.046e-028.114e-04
329.975e-032.022.028e-042.00
642.485e-032.015.069e-052.00
1286.207e-042.001.267e-052.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_varcoeff_y_periodic (log-log)

varcoeff_laplacian_y_zero_flux_bc

op:varcoeff_laplacian order:2 bc:zero_flux axes:x,y stencil width 3 priority 10

Source: grids/cartesian_uniform_2d/rules/varcoeff_laplacian_y_zero_flux_bc.esm

Rewrites ∂(k·∂u/∂y)/∂y at priority 10 — match pattern:

{"op": "D", "args": [{"op": "*", "args": ["k", {"op": "D", "args": ["u"], "wrt": "y"}]}], "wrt": "y"}
Match scope (esm-spec §9.6.1 where): fires only when u is a bare field shaped [x, y]; k is a bare field shaped [x, y]. Under import-edge renaming (§9.7.7) the wrt literal and this shape follow the renamed axis together, so the rule can be imported more than once (each instance scoped to its own grid) without a first-declared-wins collision. A consumer differentiating a compound inline expression must bind it to a declared shaped observed first.

Discretization — one makearray tiling the output axes; boundary conditions live in the face regions of this same rule (esm-spec §9.6.8; later regions overwrite earlier, §4.3.2):

RegionTreatmentValue
x ∈ [1, NX], y ∈ [2, NY−1]interiorvarcoeff_laplacian_y_interior(u, k)[i,j] = ((0.5·(k[i, j] + k[i, j+1]))·(u[i, j+1] − u[i, j]) − (0.5·(k[i, j−1] + k[i, j]))·(u[i, j] − u[i, j−1])) / dy² for i ∈ x, j ∈ [2, NY−1]
x ∈ [1, NX], y = 1boundary face((0.5·(k[i, 1] + k[i, 2]))·(u[i, 2] − u[i, 1])) / dy² for i ∈ x, j ∈ [1, 1]
x ∈ [1, NX], y = NYboundary face((0.5·(k[i, NY−1] + k[i, NY]))·(u[i, NY−1] − u[i, NY])) / dy² for i ∈ x, j ∈ [NY, NY]

The y-half of a 2-D conservative finite-volume VARIABLE-COEFFICIENT / NONLINEAR Laplacian div(k grad u) on cartesian_uniform_2d with zero-flux (homogeneous Neumann) conditions on the two y-walls, over the full 2-D field u[x,y], cell-centered. The operator carries a coefficient, so a bare laplacian match cannot express it; the rule matches the CONSERVATIVE COMPOUND the consumer writes, D( {op:, args:[k, D(u,wrt:y)]}, wrt:y ) at priority 10, so under esm-spec 9.6.3 outermost-first compound precedence it fires on the whole compound before any plain-D rule could lower the inner derivative (exactly as central_D2y fires on D(D(u,y),y)). This is the 2-D per-axis lift of grids/cartesian_uniform_1d/rules/varcoeff_laplacian_zero_flux_bc: the operated axis y gets an interior band plus two dropped-flux y-faces, while the x axis is full-rank (all columns i treated identically). BOTH operands are match wildcards bound to consumer fields: u (the state) and k (the coefficient, the LEFT factor of the product); matching is structural and non-commutative, so the consumer MUST write k * D(u, wrt:y). Match-scoped with where {u:{shape:[x,y]}, k:{shape:[x,y]}} (esm-spec 9.6.1) so it fires only on this 2-D grid’s declared [x,y]-shaped fields; under import-edge renaming (esm-spec 9.7.7) the wrt literals and both where shapes follow the renamed axes together. ONE rule serves the VARIABLE-COEFFICIENT case (k=k(x,y)) and the NONLINEAR case (k=f(u)); the body only reads k by index at the faces. AUTHORING CONTRACT for k: k is a declared [x,y]-shaped field; an aggregate-observed coefficient MUST loop over indices DISTINCT from the stencil’s interior indices i,j (the library convention is gi/gj), because it is inlined into the stencil’s [1,NX]x[2,NY-1] aggregate and a shared loop-index name is captured (silently wrong interior values). Face conductivity is the ARITHMETIC MEAN kf[i,j+1/2]=(k[i,j]+k[i,j+1])/2 (O(dy^2) at the face for a differentiable coefficient, affine in k so it composes with the nonlinear dual). Interior j in [2,NY-1] (all i): the imported conservative y-flux-difference stencil. The y-boundary faces impose zero flux exactly (natural for the finite-volume form – no ghost geometry): at j=1 the bottom y-domain-face flux is dropped, giving Ly u[i,1] = kf[i,3/2](u[i,2]-u[i,1])/dy^2 = (k[i,1]+k[i,2])/2 (u[i,2]-u[i,1])/dy^2; at j=NY the top y-domain-face flux is dropped, giving Ly u[i,NY] = -kf[i,NY-1/2](u[i,NY]-u[i,NY-1])/dy^2 = (k[i,NY-1]+k[i,NY])/2 (u[i,NY-1]-u[i,NY])/dy^2. This is exactly varcoeff_laplacian_y_neumann_bc at qYL=qYR=0. The face values are FULL-RANK aggregates (singleton output range on the y axis, all i) mirroring central_D2y_neumann_bc. dy is the grid’s consumer-supplied y-spacing free name; the zero-flux closure is dy-scale-invariant, so the rule serves any y-extent unchanged. The full variable-coefficient 2-D Laplacian is assembled in the consuming equation as D(kD(u,x),x) + D(k*D(u,y),y): this rule lowers the second term. Reduces to the constant-coefficient zero-gradient y-half when k is spatially constant. Interior O(dy^2); the one-sided zero-flux faces are O(dy) locally, which preserves global 2nd-order L2 convergence for zero-flux-compatible solutions (flux k u_y vanishing at both y-walls). A model imports exactly one second-derivative/Laplacian rule for y.

References

  • LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM. Section 2.15 (variable-coefficient conservative form) and Section 2.12 (Neumann boundary conditions by dropping the domain-face flux).

Conformance fixtures

Convergence — case heat_2d_varcoeff_y_zero_flux

Error norms read from the committed golden tests/conformance/convergence/heat_2d_varcoeff_y_zero_flux/golden/errors.json (binding: julia; nothing recomputed for display).

nL2_errorobserved orderLinf_errorobserved order
166.726e-032.584e-03
321.692e-031.996.522e-041.99
644.236e-042.001.630e-042.00
1281.060e-042.004.078e-052.00

Expected order: 2 (± 0.2).

Convergence of heat_2d_varcoeff_y_zero_flux (log-log)