Bases and domains¶
In the expert layer a domain is a coordinate system with one spectral basis per coordinate, each with its size and dealiasing factor.
A Domain is a coordinate system and one basis per coordinate. Fourier directions are separable:
the discrete problem splits into independent groups (one per wavenumber, or per azimuthal
number \(m\) in polar geometry), each solved with banded Olver–Townsend operators.
Dealiasing is always explicit¶
Every basis takes a dealias factor, and there is no default:
from conatus.expert import CartesianCoordinates, Chebyshev
z = CartesianCoordinates("z")["z"]
try:
Chebyshev(z, 32, (-1, 1))
except ValueError as err:
print(err)
# Chebyshev: dealias must be given explicitly
# why: the dealiasing factor sets the grid of the nonlinear products; Conatus never picks it for you
# fix: dealias=1 (no padding, linear problems) or dealias=3/2 (quadratic nonlinearities)
try:
Chebyshev(z, 32, (-1, 1), dealias=True)
except TypeError as err:
print(err)
# Chebyshev: dealias=True/False is ambiguous; give the factor (1, 3/2, 2, ...)
Linear problems (EVP, LBVP) use dealias=1; quadratic nonlinearities need 3/2. The grid size is
floor(N * dealias) (the Dedalus convention). With grid="smooth7" the grid is raised to the
smallest size whose FFT length has no prime factor larger than 7, so dealiasing stays exact and
the FFTs avoid large-prime lengths:
import numpy as np
from conatus.expert import CartesianCoordinates, RealFourier
x = CartesianCoordinates("x")["x"]
print(RealFourier(x, 170, (0, 2*np.pi), dealias=3/2).grid_size) # 255
print(RealFourier(x, 170, (0, 2*np.pi), dealias=3/2, grid="smooth7").grid_size) # 256
The choice is printed by basis.explain() and recorded in the provenance and Methods text.
Cartesian bases¶
Basis |
Use |
Notes |
|---|---|---|
|
real periodic data |
cos/sin pairs |
|
complex periodic data, normal modes |
\(e^{ikx}\); unknowns declared with |
|
bounded interval, walls |
\(T_n\); tau conditions (Boundary conditions) |
|
other Jacobi families |
sparse operators (Olver–Townsend); taus lifted with |
|
(semi-)infinite intervals |
declared map |
from conatus.expert import CartesianCoordinates, Chebyshev, MappedChebyshev, TanhMap
z = CartesianCoordinates("z")["z"]
y = CartesianCoordinates("y")["y"]
print(Chebyshev(z, 32, (-1, 1), dealias=3/2).explain())
# Chebyshev(z, size=32, bounds=(-1.0, 1.0), dealias=3/2)
# family : Chebyshev T (Jacobi -1/2,-1/2)
# coefficients : T_n(t), n = 0..31, t = x
# dealias : 3/2 (explicit; grid size 48, floor(N*d))
print(MappedChebyshev(y, 96, TanhMap(L=4.0), dealias=1).explain())
# MappedChebyshev(y, size=96, map=TanhMap(L=4), dealias=1)
# family : Chebyshev T in a mapped coordinate
# coefficients : T_n(xi), n = 0..95, TanhMap(L=4): y = 0 + 4 artanh(xi), dxi/dy = (1 - xi^2)/4, ...
Legendre, Jacobi and ultraspherical axes. A wall-bounded axis can be expanded in any Jacobi
family \(P_n^{(a,b)}\) (\(a, b > -1\); Legendre is \(a = b = 0\), \(C_n^{(\lambda)}\) the ultraspherical
family, and ChebyshevU, ChebyshevV are \(C^{(1)}\), \(C^{(2)}\) as in Dedalus). An equation of
differential order \(q\) is written in the family \(P^{(a+q,\,b+q)}\), so its \(q\) tau terms are lifted
there: lift(tau, basis, -1, lam=q). Every operator is sparse and exact on the retained modes, the
grid is the Gauss–Jacobi grid of the family, and a non-constant coefficient reports its spectral
tail as on a Chebyshev axis. The checks: Poisson with an exact solution to round-off in each family,
the eigenvalues of a string to \(10^{-10}\), and a viscous Burgers run equal to the Chebyshev \(T\) run
to \(5\times10^{-13}\). In the native layer cn.walls is Chebyshev \(T\) by default, or a sine/cosine series with
cn.walls(a, b, basis="sincos") (Geometries); these Jacobi families are
for the expert layer and for Dedalus scripts translated by cn.from_dedalus (Coming from Dedalus). Their GPU path
(backend="cuda") has run on an A100 only (GPU).
Mapped bases put the resolution where a free shear layer or a critical layer is, without
truncating the line to \([-L, L]\): the Taylor–Goldstein example of the gallery reaches
\(J_c = 1/4\) to \(2.5\times10^{-12}\) this way. TanhMap(L, contour=...) also offers a shifted
complex path. The map and its metric are declared and appear in explain().
The helper channel(N, base="poiseuille" | "couette" | None) builds the wall-normal Chebyshev
domain of a plane channel with its base flow as an NCC; everything it chooses is visible in
ch.explain() (dealias=1 by default, since linear problems have no products).
Curvilinear geometries¶
The same IR, the same discretize → plan → eigenvalues / solve / IVPSolver path; the bases are
regular at the axis or centre by construction (no \(1/r \to 1/(r+\varepsilon)\), no condition at
\(r = 0\)).
Geometry |
Coordinates |
Bases |
Notes |
|---|---|---|---|
Disk |
|
|
Zernike-type \(r^{\lvert m\rvert} w(2r^2-1)\) (Vasil et al. 2016) |
Annulus, Taylor–Couette |
|
× |
Chebyshev in \(r\), equations multiplied by \(r^p\) |
Pipe, periodic cylinder |
|
Fourier × |
Hagen–Poiseuille = Meseguer & Trefethen (2003) |
Sphere \(S^2\) |
|
|
spin-weighted harmonics (Vasil et al. 2019) |
Shell, ball |
|
× |
regular at the centre (ball) |
Taylor–Couette and pipe flow are written with these bases in the expert layer (see
8. Polar geometry: a drum, a pipe, two rotating cylinders); the native geometries cn.Annulus and cn.Pipe choose them for
you (Geometries).
Symmetry subspaces¶
A SymmetrySubspace(Reflection(z), T=EVEN, ...) passed to a problem (symmetry=) restricts every
field to a declared parity, and each system shrinks to the size of the subspace. Parities are
checked on every equation; an inconsistent declaration is refused. See
3. Symmetries and symmetry subspaces and Dynamics.
Sine/cosine walls (native layer)¶
cn.walls(a, b, basis="sincos") discretizes a wall axis with sine and cosine series,
\(\sin(n\pi(z-a)/L)\) and \(\cos(n\pi(z-a)/L)\), \(n < N_z\), \(L = b - a\). It is exact for
stress-free walls, fixed-temperature or insulating walls, hinged plates: every unknown is a pure
sine or cosine series, so its wall conditions hold identically and no tau term is needed. The
parity of each unknown is read off the wall conditions you already write (never guessed):
Condition (at both walls) |
Series |
|---|---|
|
sine |
|
cosine |
|
sine |
|
tangential components cosine, normal component sine |
no condition (pressure) |
the unique parity that makes every equation equivariant |
The same model text runs with basis="chebyshev" (taus) and basis="sincos". Free-slip
Rayleigh–Bénard onset, at \(k_c = \pi/\sqrt{2}\), on both sides of \(Ra_c = 27\pi^4/4\):
import numpy as np
import conatus as cn
kc = np.pi / np.sqrt(2) # critical wavenumber of stress-free RB
Ra_c = 27 * np.pi**4 / 4 # 657.511...
def model(basis, Ra):
box = cn.Box(x=cn.periodic(2 * np.pi / kc), z=cn.walls(0, 1, basis=basis))
m = cn.Model(box, name="RB stress-free", Ra=Ra, Pr=1)
u, T = m.velocity("u"), m.scalar("T")
m.pressure("p")
m.equations("""
momentum: dt(u) + grad(p) - Pr*lap(u) - Pr*Ra*T*ez = -u·grad(u)
continuity: div(u) = 0
heat: dt(T) - lap(T) - u·ez = -u·grad(T)
""")
m.at("walls", u=cn.free_slip, T=0) # stress-free, fixed temperature (perturbation)
return m, u, T
for Ra in (0.99 * Ra_c, 1.01 * Ra_c):
m, u, T = model("sincos", Ra)
s = m.linearize(about={u: 0, T: 0}).eigenvalues(8, k=kc, n=1)
print(f"sincos Ra/Ra_c = {Ra / Ra_c:.2f}: s = {s.values.real.max():+.15f}")
m, u, T = model("chebyshev", Ra)
s = m.linearize(about={u: 0, T: 0}).eigenvalues(16, k=kc)
print(f"chebyshev Ra/Ra_c = {Ra / Ra_c:.2f}: s = {s.values.real.max():+.15f}")
sincos Ra/Ra_c = 0.99: s = -0.074208019189858
chebyshev Ra/Ra_c = 0.99: s = -0.074208019189862
sincos Ra/Ra_c = 1.01: s = +0.073837897458271
chebyshev Ra/Ra_c = 1.01: s = +0.073837897458291
The closed form \((s/Pr + K^2)(s + K^2)K^2 = Ra\,k^2\), \(K^2 = k^2 + \pi^2\), gives
\(-0.07420801918986\) and \(+0.07383789745827\): the sine/cosine result is exact to round-off (each
group holds one vertical wavenumber, n=1 selects it; without n= all of them are gathered),
the Chebyshev one is spectrally accurate. m.explain() declares the choices:
m, u, T = model("sincos", Ra_c)
print("\n".join(l.strip() for l in m.explain().splitlines() if l.lstrip().startswith(("z:", "u:", "p:", "T:"))))
m, u, T = model("chebyshev", Ra_c)
print("\n".join(l.strip() for l in m.explain().splitlines() if ("free-slip" in l and "tau" in l) or l.strip().startswith("divergence tau")))
z: walls at 0, 1 -> sine/cosine bases (N_z wavenumbers n pi/L, n < N_z, L = 1), realised as the reflection subspace of RealFourier(2 N_z) on [0, 2); no tau
u: cosine/sine by component (x: cosine, z: sine) (free-slip (u·n = 0, dz(u_t) = 0) at walls)
T: sine series (T = 0 at walls)
p: cosine series (deduced from the equivariance of momentum)
How it is computed: the axis is the reflection-symmetric subspace (about \(z = a\)) of a
RealFourier(2 N_z) series on the doubled period \([a, 2b - a)\), declared on the IR problem as
SymmetrySubspace(Reflection(z, about=a), ...); plans, eigenvalues and IVPSolver run in that
subspace (reduced=True, set explicitly; reduced=False is refused). The dealiasing factor
applies to the doubled grid. sim.initial(...) evaluates your function on \([a, b]\) only and
extends it by the parity (data of a sine field that do not vanish at the walls are refused, never
truncated); sim.grid(name) and sim.coordinate("z") return the \(N_z + 1\) points of \([a, b]\);
sim.coeffs(name) and sim.solver (expert) keep the doubled-period layout.
What is refused, with what to write instead: a non-homogeneous condition (T=1 at a wall: write
the perturbation about the conduction profile, or use Chebyshev walls), a Robin or coupled
condition, conditions asking for both a sine and a cosine, no-slip velocity (normal and tangential
components would need opposite parities: use cn.free_slip or Chebyshev walls), a coefficient or
base state without parity (m.given(U=lambda z: 1 - z): give it as a cosine/sine series or use
Chebyshev walls), a source term of the wrong parity, and drift= (no tau, no spurious mode).
Whole line, half line, compound intervals¶
Three more axes for cn.Box, each discretized by a fixed basis printed by explain()
(IR: Hermite, Laguerre, CompoundChebyshev):
Axis |
Basis |
Conditions |
Resolution |
|---|---|---|---|
|
Hermite functions \(\psi_n((x-x_0)/s)\), Gaussian decay built in |
none (no boundary) |
|
|
Laguerre functions \(L_n(y)e^{-y/2}\), \(y=(x-a)/s\) |
one wall: |
|
|
Chebyshev \(T_n\) on each piece \([a, c], \dots\) |
walls at \(a\) and \(b\); continuity at each \(c\) is automatic |
one size per piece, |
scale is required: it is the width of the basis functions, a numerical choice like \(N\), never
picked for you (choose it of the order of the width of the solution). On the half line an
equation of order \(q\) is written in the derivative basis \(L^{(q)}_n e^{-y/2}\) (bidiagonal
derivative and conversions, the Laguerre analogue of the ultraspherical method) and takes \(q/2\)
conditions at its wall. On a compound interval, for each unknown of order \(q\) the value and the
derivatives up to \(q-1\) are continuous at every interface; each piece receives \(q\) tau terms in
its last modes, and explain() lists every one of them.
import numpy as np
import conatus as cn
# whole line: quantum harmonic oscillator, E_n = 2n + 1
m = cn.Model(cn.Box(x=cn.line(scale=1.0)), name="harmonic oscillator")
m.scalar("u")
m.eigenvalue("E", growth="real")
m.given(V=lambda x: x**2) # polynomial coefficient: banded, exact
m.equations("E*u + lap(u) - V*u = 0")
E = np.sort(m.eigenvalues(N=40).values.real)
print("E_0..4 =", np.round(E[:5], 12), " max error (n < 20):",
f"{np.max(abs(E[:20] - (2*np.arange(20) + 1))):.1e}")
# half line: heat equation with u(0) = 0, CNAB2
m = cn.Model(cn.Box(y=cn.half_line(scale=0.5)), name="heat on y >= 0")
m.scalar("u")
m.equations("dt(u) - lap(u) = 0")
m.at("wall", u=0)
exact = lambda y, t: y / (1 + t)**1.5 * np.exp(-y**2 / (4*(1 + t))) # noqa: E731
sim = m.simulate(N=40, scheme="CNAB2", dt=1e-3)
sim.initial(u=lambda y: exact(y, 0))
sim.run(until=0.5)
err = np.max(abs(sim.grid("u") - exact(sim.coordinate("y"), sim.t)))
print(f"half line, t = 0.5: max error {err:.1e}")
# compound interval: an interior layer tanh(20 x), same total N = 192
k = 20.0
for label, axis, N in (("walls(-1, 1)", cn.walls(-1, 1), 192),
("walls(-1, -0.15, 0.15, 1)", cn.walls(-1, -0.15, 0.15, 1), (64, 64, 64))):
m = cn.Model(cn.Box(x=axis), name="layer")
m.scalar("u")
m.given(f=lambda x: -2*k*k*np.tanh(k*x) / np.cosh(k*x)**2)
m.equations("lap(u) = f")
m.at("bottom", u=np.tanh(-k))
m.at("top", u=np.tanh(k))
sol = m.solve(N=N) # linear boundary-value model: direct solve
print(f"{label:>26}: max error {np.max(abs(sol.grid('u') - np.tanh(k*sol.coordinate('x')))):.1e}")
E_0..4 = [1. 3. 5. 7. 9.] max error (n < 20): 2.8e-14
half line, t = 0.5: max error 7.1e-08
walls(-1, 1): max error 1.3e-07
walls(-1, -0.15, 0.15, 1): max error 2.5e-14
The half-line error is the time error of CNAB2 at dt = 1e-3 (second order); the spatial error
is below \(10^{-9}\) at \(N = 40\). What m.explain(N=(64, 64, 64)) says about the compound layer:
print("\n".join(l.strip() for l in m.explain(N=(64, 64, 64)).splitlines()
if l.strip().startswith(("x:", "u = ", "u continuous", "dx(u) continuous"))))
x: compound interval [-1, -0.15], [-0.15, 0.15], [0.15, 1] -> Chebyshev T per piece (N_x = one size per piece), continuity of value and derivatives at x = -0.15, x = 0.15 (automatic)
u = -1 at bottom: tau_u1 lifted into mode 63 (N_x,0-1) of piece 0 of [0] eq0 -- u has its highest x-order (2) there; Chebyshev lift into the last modes of piece 0
u continuous at interface x=-0.15: tau_u2 lifted into mode 62 (N_x,0-2) of piece 0 of [0] eq0 -- continuity of u across the interface (order 2 of u); Chebyshev lift into the last modes of piece 0
dx(u) continuous at interface x=-0.15: tau_u3 lifted into mode 63 (N_x,1-1) of piece 1 of [0] eq0 -- continuity of dx(u) across the interface (order 2 of u); Chebyshev lift into the last modes of piece 1
u continuous at interface x=0.15: tau_u4 lifted into mode 62 (N_x,1-2) of piece 1 of [0] eq0 -- continuity of u across the interface (order 2 of u); Chebyshev lift into the last modes of piece 1
dx(u) continuous at interface x=0.15: tau_u5 lifted into mode 63 (N_x,2-1) of piece 2 of [0] eq0 -- continuity of dx(u) across the interface (order 2 of u); Chebyshev lift into the last modes of piece 2
u = 1 at top: tau_u6 lifted into mode 62 (N_x,2-2) of piece 2 of [0] eq0 -- u has its highest x-order (2) there; Chebyshev lift into the last modes of piece 2
When to use which: a compound interval pays off for a localized feature (interior layer, boundary layer, front); for a smooth solution one interval with the same total \(N\) is more accurate (each piece spends \(q\) modes on its taus: \(e^x\) on \([-1, 1]\) is \(10^{-8}\) with \((8, 8)\) against \(10^{-15}\) with 16). A mapped Chebyshev basis (above) remains the choice for a bounded but non-decaying solution on an infinite interval; Hermite and Laguerre for decaying ones.
What is refused, with what to write instead: a constant along an unbounded axis (a constant
source, tau or gauge: write the problem for the decaying deviation from the far field), a
non-polynomial coefficient on an unbounded axis (tanh(x): make the term explicit or use a mapped
basis), a condition on cn.line (no boundary), one integer for a compound axis
(N_x=(24, 24): one size per piece), incompressible flow on these axes (the gauge/divergence tau
recipe is derived for walls(a, b) only: use the IR), symmetry subspaces along them, and the
fused GPU kernel (f_backend="fused": the evaluator runs on the GPU instead, transforms by dense
matrices on the device).