Channels between two surfaces

A channel whose walls are not flat: a wavy bottom, a bump, a corrugated lid. The equations are written in the physical coordinates, exactly as for a flat box; Conatus maps the region between the two surfaces onto a rectangle and carries the change of variables into every operator.

A wavy wall

cn.MappedBox takes periodic axes and one vertical axis cn.between(bottom=..., top=...). Each surface is a number (a flat wall) or a function of the periodic coordinates, named after them (lambda x: ... in 2D, lambda x, y: ... in 3D). The boundaries are bottom, top and walls.

import numpy as np

import conatus as cn

eps = 0.01                                     # wall amplitude (wavelength 2 pi, gap 1)
box = cn.MappedBox(x=cn.periodic(2 * np.pi),
                   z=cn.between(bottom=lambda x: eps * np.cos(x), top=1.0))
m = cn.Model(box, name="Couette over a wavy wall")
m.velocity("u")
m.pressure("p")
m.equations("grad(p) - lap(u) = 0", "div(u) = 0")      # Stokes flow
m.at("bottom", u=0)                            # no slip on the wavy wall
m.at("top", u=(1, 0))                          # the upper plate slides at speed 1
sol = m.solve((16, 24))                        # (N_x, N_z)

U, X, Z = sol.grid("u"), sol.coordinate("x"), sol.coordinate("z")
print(sol.explain())
print("grid:", U.shape, "| z of the bottom row runs from", f"{Z[:, 0].min():.4f}",
      "to", f"{Z[:, 0].max():.4f}")

LBVP solution on a mapped domain: one coupled system, relative residual 2.24e-16 (sparse LU, scipy SuperLU)
grid: (2, 16, 24) | z of the bottom row runs from -0.0089 to 0.0111
max |u_x - (z + eps u_1)| = 4.30e-04 = 4.30 eps^2

sol.coordinate("z") is a full array: the height of a grid point depends on its horizontal position. The grid is \(x_j\) (uniform) times the Gauss–Chebyshev points \(\xi_l\) of the terrain-following coordinate.

For a small amplitude \(\varepsilon\) the flow is \(u = z + \varepsilon u_1 + O(\varepsilon^2)\), with \(u_1\) from the stream function \(\phi(z)\cos x\), \((\partial_z^2 - 1)^2\phi = 0\), and the no-slip condition moved from \(z = \varepsilon\cos x\) to \(z = 0\). The difference with the computed field is of second order in \(\varepsilon\) (it falls by 4.0 each time \(\varepsilon\) is halved; this is a test of the library):

# first order in eps: psi_1 = phi(z) cos x with (D^2 - 1)^2 phi = 0,
# phi(0) = phi(1) = 0, phi'(0) = -1 (no slip moved to z = 0), phi'(1) = 0
def rows(z):
    ch, sh = np.cosh(z), np.sinh(z)
    return np.array([ch, z * ch, sh, z * sh]), np.array([sh, ch + z * sh, ch, sh + z * ch])

co = np.linalg.solve([rows(0)[0], rows(1)[0], rows(0)[1], rows(1)[1]], [0, 0, -1, 0])
u1 = np.tensordot(co, rows(Z)[1], 1) * np.cos(X)
rem = np.abs(U[0] - Z - eps * u1).max()
print(f"max |u_x - (z + eps u_1)| = {rem:.2e} = {rem / eps**2:.2f} eps^2")

What the map adds

The vertical coordinate is replaced by \(\xi \in [-1, 1]\),

\[ z = b(x) + \bigl(t(x) - b(x)\bigr)\,\frac{1 + \xi}{2} \]

(Gal-Chen and Somerville, 1975), and every physical derivative is written by the chain rule,

\[ \partial_z = J\,\partial_\xi, \qquad \partial_x\big|_z = \partial_x\big|_\xi + (A + B\,\xi)\,\partial_\xi, \qquad J = \frac{2}{t - b},\; A = -J\,\frac{b' + t'}{2},\; B = -J\,\frac{t' - b'}{2}, \]

with the volume element \(\tfrac{t - b}{2}\,d\xi\,dx\) for integrals (the pressure gauge \(\int p\,dV = 0\) is over the physical domain). The metric is a polynomial of degree one in \(\xi\): along the vertical the operators stay banded, exactly. Along \(x\) the functions \(J\), \(A\), \(B\) multiply, so the Fourier modes are coupled: Conatus builds one sparse system for the whole domain (all modes, the metric inside the operator) and solves it directly. Nothing is split into an explicit correction, so there is no condition on the wall slope.

import numpy as np

import conatus as cn

box = cn.MappedBox(x=cn.periodic(2 * np.pi),
                   z=cn.between(bottom=lambda x: 0.1 * np.cos(x), top=1.0))
m = cn.Model(box, name="Poisson")
m.scalar("u")
m.given("f", lambda x, z: np.sin(x) * z)
m.equations("lap(u) = f")
m.at("walls", u=0)
print(box.explain())
print(m.discretize((24, 16)).explain())

Geometry: MappedBox(x=periodic(6.28319), z=between(function(x), 1)) (mapped box (terrain-following))
  bases (fixed per geometry, sizes given at discretization):
    x: periodic on [0, 6.28319) -> RealFourier (cos/sin)(N_x)
    z: between the surfaces function(x) and 1 -> Chebyshev T(xi)(N_z) in the terrain-following coordinate xi in [-1, 1]
    map: z = b(x) + (t(x) - b(x)) (1 + xi)/2, xi in [-1, 1] (b: function(x); t: 1)
    metric (chain rule, coefficients of the operators): dz = J dxi, J = 2/(t - b); dx|_z = dx|_xi + (A_x + B_x xi) dxi, A_x = -J (b_x + t_x)/2, B_x = -J (t_x - b_x)/2; dV = H dxi dx, H = (t - b)/2
    all Fourier modes coupled by the metric: one sparse system (banded along xi); traces 'bottom'/'top' are at xi = -1/+1; dz(f) in a condition is the physical vertical derivative, not the normal one
  named boundaries: bottom = z left; top = z right; walls = z left, z right
TerrainPlan of LBVP 'Poisson' on a mapped domain
  map: z = b(x) + (t(x) - b(x)) (1 + xi)/2, xi in [-1, 1] (b: function(x); t: 1)
  chain rule: dz = J dxi, J = 2/(t - b); dx|_z = dx|_xi + (A_x + B_x xi) dxi, A_x = -J (b_x + t_x)/2, B_x = -J (t_x - b_x)/2; dV = H dxi dx, H = (t - b)/2
  metric terms implicit (one coupled sparse system, all Fourier modes; banded along xi: metric polynomial of degree <= 1 in xi, exact)
  Fourier modes |k| <= N/2 - 1 per direction (23 complex modes of real fields), 16 Chebyshev modes along xi
  system: 414 x 414, nnz(L) = 37397, nnz(M) = 0; direct sparse LU (SuperLU), refused above 100000 unknowns or nnz 1e+07
  metric J: 23 of 45 Fourier coefficients kept (Fourier bandwidth), sampled tail 5.0e-15 on the grid (accepted up to ncc_tail_tol = 1e-08)
  metric H: 3 of 45 Fourier coefficients kept (Fourier bandwidth), sampled tail 5.3e-18 on the grid (accepted up to ncc_tail_tol = 1e-08)
  metric A_x: 24 of 45 Fourier coefficients kept (Fourier bandwidth), sampled tail 9.7e-14 on the grid (accepted up to ncc_tail_tol = 1e-08)
  metric B_x: 24 of 45 Fourier coefficients kept (Fourier bandwidth), sampled tail 9.7e-14 on the grid (accepted up to ncc_tail_tol = 1e-08)
  sampled Chebyshev tail of every non-constant coefficient <= ncc_tail_tol = 1e-08 (declared default); above: refused, never truncated silently
  data and explicit terms: pointwise on the dealiased physical grid (x_j, z(x_j, xi_l))

The Fourier coefficients of each metric function are computed once on a fine grid. Their tail on the grid of the unknowns must be below ncc_tail_tol (default 1e-8), like any non-constant coefficient: a wall too steep for the horizontal resolution is refused, with the resolution that resolves it, never used truncated (Errors and what they mean). The same holds for a known field m.given(...), which on a mapped domain is a function of \(x\) and of the physical \(z\).

Convergence

A manufactured solution: the error falls geometrically with the resolution down to rounding.

import numpy as np

import conatus as cn

eps = 0.05
h = lambda x: eps * np.cos(x)                               # noqa: E731
ue = lambda x, z: np.exp(np.sin(x)) * np.sin(z) + z**2      # noqa: E731  exact solution
fe = lambda x, z: ((np.cos(x)**2 - np.sin(x) - 1) * np.exp(np.sin(x)) * np.sin(z)  # noqa: E731
                   + 2)                                     # its Laplacian
for N in (16, 20, 24, 28):
    m = cn.Model(cn.MappedBox(x=cn.periodic(2 * np.pi), z=cn.between(bottom=h, top=1.0)))
    m.scalar("u")
    m.given("f", fe)
    m.equations("lap(u) = f")
    m.at("bottom", u=lambda x: ue(x, h(x)))
    m.at("top", u=lambda x: ue(x, 1.0))
    s = m.solve((N, N))
    err = np.abs(s.grid("u") - ue(s.coordinate("x"), s.coordinate("z"))).max()
    print(f"N = {N:2d}: max error {err:.1e}, residual {s.max_residual:.1e}")

N = 16: max error 1.7e-07, residual 2.3e-16
N = 20: max error 4.6e-10, residual 1.9e-16
N = 24: max error 8.8e-13, residual 6.9e-16
N = 28: max error 2.7e-15, residual 2.1e-16

Conditions on a curved wall

A value (T=..., u=0) is imposed on the surface itself. A derivative in a condition is the physical one along the coordinate: dz(T) at bottom is the vertical derivative, not the derivative along the normal. To impose a flux, write the normal as known fields of \(x\):

slope = lambda x: -eps * np.sin(x)                 # noqa: E731  b'(x)
norm = lambda x: np.sqrt(1 + slope(x)**2)          # noqa: E731
m = cn.Model(box, name="heat flux through a wavy wall")
m.scalar("T")
m.given("nx", lambda x: slope(x) / norm(x))       # outward unit normal of z = b(x)
m.given("nz", lambda x: -1 / norm(x))
m.equations("lap(T) = 0")
m.at("bottom", "nx*dx(T) + nz*dz(T) = 1")          # unit flux into the fluid
m.at("top", T=0)
sol = m.solve((24, 16))
print(sol.explain())

cn.free_slip and cn.navier_slip are refused on a curved surface: their operators use the normal of a flat wall. On a flat surface of a MappedBox (top=1.0, say) they are accepted and are the conditions of cn.Box.

Problems, schemes and limits

  • m.solve(N) (linear boundary-value problems), m.eigenvalues(N, ...) (one coupled spectrum for all modes) and m.simulate(N, scheme=..., dt=...) with SBDF1–SBDF4, RK111, RK222, RK443. The implicit part is the coupled operator; nonlinear terms on the right are evaluated on the dealiased physical grid. Boundary data that depend on time are imposed at each stage time (bc_time="stage"; the default "auto" resolves to it and explain() says why; "rk" and "step" are refused). RK443 is then of order 2: the same errors as cn.Box with bc_time="stage", and the same order as its default "rk" on this manufactured diffusion test.

  • Cost: one sparse LU of the coupled system (all Fourier modes). Measured on one CPU core: 2D Stokes at 128 × 32 (12 701 unknowns) about 90 s and 1.8 GB; 3D Poisson at 24 × 24 × 16 (9 522 unknowns) about 16 s and 1.1 GB. Refused above 100 000 unknowns or 10⁷ nonzeros in the operators (MAX_COUPLED_UNKNOWNS, MAX_COUPLED_NNZ of conatus.discretize.terrain_plan, declared in cn.tolerances.explain()).

  • A given field that multiplies an unknown must be resolved on the grid: a tail above ncc_tail_tol along \(\xi\) or in the kept Fourier modes, or energy at the Nyquist wavenumber (aliasing), is refused with the fix.

  • 3D: cn.MappedBox(x=cn.periodic(Lx), y=cn.periodic(Ly), z=cn.between(bottom=lambda x, y: ..., top=...)). The coupled system grows as \((N_x N_y)^2\) in the worst case: moderate resolutions.

  • With flat surfaces the system is the one of cn.Box (identical solutions to rounding, no coupling between Fourier modes).

  • Refused with a message: traces at other heights, horizontal averages at fixed \(z\), periodic(..., complex=True), symmetry subspaces, GPU backends.