Coming from Dedalus¶

If you know Dedalus, you already know how to write a Conatus model; what changes is what you no longer write, and what you can do with the model afterwards.

The same problem, side by side¶

Two-dimensional Boussinesq convection between a hot and a cold wall, periodic sideways (14. Coming from Dedalus runs both and compares them). On the left, the problem part of a Dedalus 3 script written for that tutorial (rb_cell_dedalus.py); on the right, the same cell as a Conatus model:

Dedalus

plane = d3.CartesianCoordinates('x', 'z')
mesh = d3.Distributor(plane, dtype=np.float64)
horiz = d3.RealFourier(plane['x'], size=N_x, bounds=(0, period), dealias=1.5)
vert = d3.ChebyshevT(plane['z'], size=N_z, bounds=(0, 1), dealias=1.5)
u = mesh.VectorField(plane, name='u', bases=(horiz, vert))
p = mesh.Field(name='p', bases=(horiz, vert))
th = mesh.Field(name='th', bases=(horiz, vert))
u_bot = mesh.VectorField(plane, name='u_bot', bases=horiz)    # tau terms, one per wall condition
u_top = mesh.VectorField(plane, name='u_top', bases=horiz)
th_bot = mesh.Field(name='th_bot', bases=horiz)
th_top = mesh.Field(name='th_top', bases=horiz)
gauge = mesh.Field(name='gauge')                               # constant tau of the pressure

# Tau terms live in the first-derivative basis of z; the bottom ones enter through the gradients
e_x, e_z = plane.unit_vector_fields(mesh)
dz_basis = vert.derivative_basis(1)
tau_z = lambda f: d3.Lift(f, dz_basis, -1)  # noqa: E731 (a Dedalus idiom)
Gu = d3.grad(u) + e_z*tau_z(u_bot)
Gth = d3.grad(th) + e_z*tau_z(th_bot)

onset = d3.IVP([u, p, th, u_bot, u_top, th_bot, th_top, gauge], namespace=locals())
onset.add_equation("dt(u) - Pr*div(Gu) + grad(p) - Ra*Pr*th*e_z + tau_z(u_top) = -u@grad(u)")
onset.add_equation("trace(Gu) + gauge = 0")
onset.add_equation("dt(th) - div(Gth) - u@e_z + tau_z(th_top) = -u@grad(th)")
onset.add_equation("u(z=0) = 0")
onset.add_equation("th(z=0) = 0")
onset.add_equation("u(z=1) = 0")
onset.add_equation("th(z=1) = 0")
onset.add_equation("integ(p) = 0")

Conatus

# Onset of convection in diffusive units: th is the temperature minus the conduction profile 1 - z,
# so it vanishes on both walls; u.ez carries the conduction profile into the heat equation.
RA, PR, PERIOD, NX, DT = 3000.0, 7.0, 2.0, 24, 5e-4


def start(x, z):
    """One roll of the temperature deviation."""
    return 1e-3 * np.cos(2 * np.pi * x / PERIOD) * np.sin(np.pi * z)


def native(nz):
    """The cell as a Conatus model: unknowns by role, equations as written on paper, conditions by
    wall name. No tau field and no lift: the native layer places one tau per condition."""
    m = cn.Model(cn.Box(x=cn.periodic(PERIOD), z=cn.walls(0, 1)), name="convection onset",
                 Ra=RA, Pr=PR)
    m.velocity("u"), m.pressure("p"), m.scalar("th")
    m.equations("""
        momentum:  dt(u) - Pr*lap(u) + grad(p) - Ra*Pr*th*ez = -u·grad(u)
        mass:      div(u) = 0
        heat:      dt(th) - lap(th) - u·ez = -u·grad(th)
    """)
    m.at("bottom", u=0, th=0)
    m.at("top", u=0, th=0)
    sim = m.simulate(N=(NX, nz), scheme="SBDF2", dt=DT)
    sim.initial(th=start)
    return sim.run(steps=40)


cell = native(16)
print(f"native model: t = {cell.t:.3f}, max|u| = {np.abs(cell.grid('u')).max():.3e}")

What you stop writing¶

Tau fields and lifts. Each condition is placed by a fixed rule, printed by explain() (Boundary conditions); the pressure gauge is recognised from the structure.
Bases with sizes. The geometry fixes the bases; the number of modes is given when you solve.
The split into left and right sides. Linear terms are implicit, nonlinear terms explicit, wherever you wrote them (Equations).
The time loop. sim.run(until=t) with a fixed dt; Conatus never adapts dt for you.

What stays the same: the spectral method itself (Fourier, Chebyshev and curvilinear bases with sparse ultraspherical operators, the generalised tau method, IMEX schemes with Dedalus’s coefficients), the string syntax of the equations, and the dealiasing convention. On the problems where both codes apply the same method, they agree to round-off.

What is added: the GPU backend, and the dynamics toolkit on the same model: eigenvalues with residuals, adjoints, Newton–Krylov searches for invariant solutions, continuation, Floquet and Lyapunov analysis (Dynamics).

Bases and conditions you may know from Dedalus 2 are declarations of the axis in the native layer: sine and cosine series between stress-free or fixed-value walls (cn.walls(0, 1, basis="sincos"), the parity of each field read from its conditions), Hermite functions on the whole line (cn.line(scale=...)), Laguerre functions on the half line (cn.half_line(scale=...)) and compound Chebyshev intervals (cn.walls(a, c, b), continuity at the cuts added for you); a stress-free wall is u=cn.free_slip, on plane and curved walls (Bases and domains, Boundary conditions). All thirteen time steppers of Dedalus 3 are there, with the same coefficients (Simulations and time steppers).

Importing a Dedalus script¶

cn.from_dedalus(script) translates a Dedalus 3 script into Conatus. You get back:

the Conatus problem;
the declared time stepping;
the initial and given data;
a line-by-line report.

from pathlib import Path

import conatus as cn

Path("kdv_dedalus.py").write_text("""
import numpy as np
import dedalus.public as d3
period, Nx, damping, dispersion = 8, 64, 3e-4, 5e-4
s = d3.Coordinate('x')
mesh = d3.Distributor(s, dtype=np.float64)
ring = d3.RealFourier(s, size=Nx, bounds=(0, period), dealias=1.5)
h = mesh.Field(name='h', bases=ring)
D = lambda F: d3.Differentiate(F, s)
wave = d3.IVP([h], namespace=locals())
wave.add_equation("dt(h) - damping*D(D(h)) - dispersion*D(D(D(h))) = -h*D(h)")
xs = mesh.local_grid(ring)
h['g'] = np.exp(-2*(xs - 4)**2)
run = wave.build_solver(d3.SBDF3)
run.stop_sim_time = 0.025
dt_fixed = 2.5e-3
while run.proceed:
    run.step(dt_fixed)
""")                                                 # a Dedalus 3 script (KdV-Burgers)

tr = cn.from_dedalus("kdv_dedalus.py")            # a path, or the source text
print(tr.explain())                               # counts + one line per statement
tr.problem                                        # an ordinary Conatus Problem
solver = tr.ivp_solver(startup="lower_order")     # scheme and fixed dt from the script
solver.run(tr.steps())                            # stop_sim_time / dt steps, no Python loop
print(tr.to_script())                             # a runnable program, report as comments

In a notebook, tr displays the report as a colour-coded table. 14. Coming from Dedalus goes through a complete example and compares the result with Dedalus.

What happens to your script¶

The script is read, not run. Python’s ast module parses it, and a restricted interpreter walks it statement by statement. It understands:

literals, arithmetic, tuples and lambdas;
the dedalus.public names it knows;
a whitelist of NumPy/SciPy functions, used for constants and initial data;
the equation strings, which go through the Equations parser.

No statement of the script is executed by Python, and Dedalus is never imported.

Each statement gets one of four statuses:

Status	Meaning
translated	same mathematics and same method; the note lists any declared rewrite (`trace(grad(u)) -> div(u)`) or Dedalus default made explicit (`dealias=1`)
approximated	translated, but the method differs in a way that can change digits; the note says how
unsupported	not translated, with the reason and what to write instead; everything that depends on it is reported as well
skipped	output only: logging, plots, reading results after the solve

tr.complete is true only when a problem was translated with no line unsupported and none approximated. A time loop whose body changes the simulation (writes field data, calls a mutating method, stops on a condition) is unsupported: solver.run(n) would not do it. **options in a call are expanded, and a problem’s own namespace={...} alone resolves the names of its equations, as in Dedalus. tr.to_script() is a runnable program: it translates the embedded script again, refuses to run if a line is unsupported, then runs the IVP or solves the LBVP (an EVP or NLBVP stops at the discretized plan: the group or the initial guess is your choice).

The translated subset¶

Dedalus	Conatus	Note
`Coordinate`, `CartesianCoordinates`, `PolarCoordinates`, `Distributor(dtype=...)`	coordinate systems; field dtype	process mesh and `comm` are ignored (parallel layout is decided at run time)
`RealFourier`, `ChebyshevT`/`Chebyshev`, `Legendre`, `Jacobi(a, b)`	same bases	`dealias` omitted means Dedalus’s default of 1, written explicitly
`ComplexFourier(N)`	`ComplexFourier(N)`	approximated: Dedalus keeps the Nyquist mode \(k=-N/2\), Conatus drops it; that group is absent
`DiskBasis(coords, shape=(Nphi, Nr))`	`RealFourier`/`ComplexFourier(phi)` × `DiskBasis(r)`	same polynomial space for each \(m\); the pipe-flow eigenvalues equal Dedalus’s to round-off
`AnnulusBasis`	× `AnnulusBasis(r)`	approximated: Conatus multiplies the equations by a declared power of \(r\)
`SphereBasis`, `ShellBasis`, `BallBasis`, `S2Coordinates`, `SphericalCoordinates`	—	unsupported in the translator: the geometries exist in Conatus (Geometries), but the shape-to-size mapping is not automatic
`dist.Field`, `VectorField` (named or not)	`Field`, `VectorField`	an unnamed field takes its variable’s name (reported)
`IVP`, `LBVP`, `EVP(eigenvalue=s)`, `NLBVP`, `IVP(time=t)`	same classes	a 0-dimensional eigenvalue or time field becomes a `Parameter`
`add_equation("...")`, `add_equation((lhs, rhs))`	`add_equation`	`u*v` of two tensors is `outer(u, v)` (Dedalus convention, reported)
`Lift(A, basis.derivative_basis(q), n)`	`lift(A, basis, n, lam=q)`	Chebyshev, Legendre, Jacobi, `ChebyshevU/V`, `Ultraspherical` and disk bases (data given in coefficient space on `ChebyshevU/V` is refused: their coefficients are normalised differently)
`trace`, `div(ez*lift(tau))` of first-order formulations	`div`, `dot`, `dz`	exact identities, listed in the note
`SBDF1-4`, `CNAB1`, `CNAB2`, `MCNAB2`, `CNLF2`, `RK111`, `RK222`, `RK443`, `RKSMR`, `RKGFY`	same schemes, same start-up (lower order first)	results agree with Dedalus to round-off; time-dependent boundary data: `tr.ivp_solver()` keeps the Dedalus convention (`bc_time="step"`, said by `explain()`), while a solver built by hand imposes them so that each Runge–Kutta scheme keeps its order (`bc_time="auto"`, which is `"rk"`: the implicit tableau applied to \(g'\); `"stage"` imposes \(g\) at the stage times; pass either to `ivp_solver`; Simulations and time steppers); `SBDF2_R`, `RK443_IMEX` are unsupported
`solver.step(dt)` in `while solver.proceed`	`solver.run(n)` with fixed `dt`	the time loop becomes one call
`CFL(...)`, `compute_timestep()`	—	unsupported: Conatus never changes `dt` by itself; give `dt`, `conatus.expert.CFL` reports the CFL step
`f['g'] = ...`, `f['g'][i] += ...` with grids from `local_grids`	grid functions `tr.data`, used by `set_state` / `ncc`	NumPy/SciPy whitelist
`fill_random`, `low_pass_filter`, `f['c'] = ...`	—	unsupported: Dedalus’s random layout cannot be reproduced; write an explicit initial condition
file handlers, `GlobalFlowProperty`	—	unsupported: use `conatus.runs` and the Studio
`solve_dense` / `solve_sparse(subproblems[i])`	`eigenvalues(plan, group)`	approximated: a Dedalus subproblem index is not a Conatus group; select the group by wavenumber

Scripts that build their problem inside a function (for example rayleigh_benard_evp.py) are translated with from_dedalus(script, entry=("max_growth_rate", dict(Rayleigh=1710, ...))). overrides={"Nx": 64} replaces top-level assignments (also a name bound by unpacking, Nx, Nz = 256, 64), and each override is listed in the report; an override that matches no top-level assignment is refused with the names that can be overridden.

The official Dedalus examples¶

Every official Dedalus 3 example was translated unmodified, which tests the coverage of the translator. Where the comparison is possible, both codes then ran the same script (smaller sizes, a short time, a deterministic initial condition instead of fill_random, a fixed dt instead of the CFL): KdV–Burgers, Rayleigh–Bénard, the 2D shear flow, the Poisson problem, the Rayleigh–Bénard eigenvalue problem and the pipe-flow eigenvalues agree to round-off. The spherical examples (sphere, shell, ball) use geometries that Conatus has, but their translation is not automatic yet: write them as native models (Geometries).