Skip to content

Discretisation

The solver works on a uniform rectangular grid. The domain is \([0,\ \text{aspect}] \times [0,\ 1]\), where aspect \(= h_x / h_z\) is the cell aspect ratio (Settings::aspectRatio). FDBase::calculate_geometry places \(n_x \times n_z\) points on it, giving spacings

\[ h_z = \frac{1}{n_z - 1}, \qquad h_x = \text{aspect}\;h_z. \]

Stack builds \(n_x = n_z\) at every level of the hierarchy (coarsest level minimumResolution on both axes, each finer level \(n \to 2(n-1)+1\)), so the point counts are independent of aspect -- it only scales \(h_x\). With aspect \(= 1\) the cells are square on the unit square. FDArray caches the derivative pre-factors that go with these spacings (\(1/2h_x\), \(1/h_x^2\), \(1/4h_xh_z\), and the \(z\) equivalents) so the stencil methods are a handful of multiply--adds.

Interior stencils

At an interior point FDArray uses standard second-order central differences.

Method Approximates Stencil
dx, dz \(\partial_x f\), \(\partial_z f\) \(\dfrac{f_{i+1,j} - f_{i-1,j}}{2h_x}\)
dxx, dzz \(\partial_{xx} f\), \(\partial_{zz} f\) \(\dfrac{f_{i-1,j} - 2f_{i,j} + f_{i+1,j}}{h_x^2}\)
dxz \(\partial_{xz} f\) \(\dfrac{f_{i+1,j+1} - f_{i+1,j-1} - f_{i-1,j+1} + f_{i-1,j-1}}{4h_xh_z}\)

gradient and gradient_magnitude build on dx / dz; the discrete Laplacian used by the examples is dxx + dzz, the 5-point operator

\[ (\Delta_h u)_{i,j} = \frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{h_x^2} + \frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{h_z^2}, \]

second-order accurate in the (square) cell size \(h\). tests/test_poisson_mms.cpp confirms the observed order on a manufactured solution: 2.00.

One-sided stencils at the boundary

A central stencil reads points that do not exist on the first and last rows and columns. Each stencil method switches to a one-sided form there:

  • First derivatives (dx, dz): three-point one-sided, \((\pm 3f_{0} \mp 4f_{1} \pm f_{2})/2h\) -- second order.
  • Second derivatives (dxx, dzz): four-point one-sided, \((2f_0 - 5f_1 + 4f_2 - f_3)/h^2\) -- second order, but it reaches three points into the grid.
  • Mixed derivative (dxz): a 9-region switch covering the four edges and four corners, each region using the one-sided form in whichever direction runs off the grid.

The *u variants (dxu, dxxu, dzz, ...) are the same stencils evaluated on the constant field \(u \equiv 1\); the solver uses them to assemble the diagonal of the operator for the relaxation update without a second pass over the data.

tests/test_stencils_mixed.cpp pins every dxz region against a constant field (which must map to zero) and a bilinear field (which must be exact).

Quadrature

FDArray::calculate_flux integrates the field over the domain with an extended Simpson rule (end weights \(17/48,\ 59/48,\ 43/48,\ 49/48\)), applied along rows and then down the resulting column of row-integrals. For the viscoplastic example this is the dimensionless down-channel flux \(q\).