Weakly nonlinear theory derived from the IR: Landau and Ginzburg–Landau coefficients¶
Family: bifurcation · Script: examples/bifurcation/weakly_nonlinear.py
(a) Swift–Hohenberg, du/dt = r u − (1 + ∂x²)² u − u³ (Swift & Hohenberg 1977). Rolls u = A cos x:
the pitchfork normal form on the even subspace gives dA/dt = r A − (3/4) A³, i.e. A² = 4r/3
(Cross & Hohenberg 1993, RMP 65, §IV); the envelope equation has τ0 = 1, ξ0² = 4, g0 = 3.
(b) Rayleigh–Bénard between stress-free isothermal walls (ψ–θ form, thermal-diffusion units,
z ∈ (0, 1)), onset Ra_c = 27π⁴/4 at k_c = π/√2 (Rayleigh 1916). The Landau equation and the
second-order mean temperature give the heat flux of the rolls: Nu − 1 = 2 (Ra − Ra_c)/Ra_c + ...
for every Pr (Malkus & Veronis 1958; the same number follows from the Lorenz 1963 truncation,
exact at this order since the O(a²) correction is the sin 2πz mean mode only).
Every coefficient comes from ``normal_form``: symbolic 2nd/3rd Fréchet derivatives of the IR
right-hand side, bordered solves on the discrete Jacobian. Nothing is derived by hand here.
Run: python examples/bifurcation/weakly_nonlinear.py
Run it: python examples/bifurcation/weakly_nonlinear.py (in the conatus environment).
1"""Weakly nonlinear theory derived from the IR: Landau and Ginzburg–Landau coefficients.
2
3(a) Swift–Hohenberg, du/dt = r u − (1 + ∂x²)² u − u³ (Swift & Hohenberg 1977). Rolls u = A cos x:
4 the pitchfork normal form on the even subspace gives dA/dt = r A − (3/4) A³, i.e. A² = 4r/3
5 (Cross & Hohenberg 1993, RMP 65, §IV); the envelope equation has τ0 = 1, ξ0² = 4, g0 = 3.
6(b) Rayleigh–Bénard between stress-free isothermal walls (ψ–θ form, thermal-diffusion units,
7 z ∈ (0, 1)), onset Ra_c = 27π⁴/4 at k_c = π/√2 (Rayleigh 1916). The Landau equation and the
8 second-order mean temperature give the heat flux of the rolls: Nu − 1 = 2 (Ra − Ra_c)/Ra_c + ...
9 for every Pr (Malkus & Veronis 1958; the same number follows from the Lorenz 1963 truncation,
10 exact at this order since the O(a²) correction is the sin 2πz mean mode only).
11
12Every coefficient comes from ``normal_form``: symbolic 2nd/3rd Fréchet derivatives of the IR
13right-hand side, bordered solves on the discrete Jacobian. Nothing is derived by hand here.
14
15Run: python examples/bifurcation/weakly_nonlinear.py
16"""
17import numpy as np
18
19import conatus.expert as C
20from conatus.bifurcation import (amplitude_equation, cos_mode_observable, ginzburg_landau, ir_system,
21 normal_form)
22
23
24# ------------------------------------------------------------------ (a) Swift–Hohenberg
25def swift_hohenberg(k: float, N: int = 16):
26 """SH on one period 2π/k, even subspace (kills the translation mode), r free."""
27 cs = C.CartesianCoordinates("x")
28 x = cs["x"]
29 D = C.Domain(cs, [C.RealFourier(x, N, (0, 2 * np.pi / k), dealias=2)]) # cubic: dealias 2
30 u, r = C.Field("u", D), C.Parameter("r", 0.0)
31 pb = C.IVP([u], name="Swift-Hohenberg", symmetry=C.SymmetrySubspace(C.Reflection(x, about=0), u=C.EVEN))
32 L = lambda f: f + C.diff(C.diff(f, x), x) # noqa: E731
33 pb.add_equation((C.dt(u) - r * u + L(L(u)), -u * u * u))
34 s = ir_system(pb, ["r"])
35 return s, np.zeros(s.n)
36
37
38def sh_coefficients():
39 s, x0 = swift_hohenberg(1.0)
40 le = amplitude_equation(s, x0, "r", kind="stationary", observable=cos_mode_observable(s, "u"))
41 gl = ginzburg_landau(swift_hohenberg, 1.0, "r", kind="stationary", observable=cos_mode_observable(s, "u"), h=1e-2)
42 return le, gl
43
44
45# ------------------------------------------------------------------ (b) stress-free convection
46def rb_free_slip(Nx: int, Nz: int, Pr: float):
47 """ψ–θ Boussinesq rolls between stress-free walls, x-period 2π/k_c, reflection subspace
48 (ψ odd, θ even in x). Ra free, starting at Ra_c = 27π⁴/4."""
49 kc = np.pi / np.sqrt(2)
50 cs = C.CartesianCoordinates("x", "z")
51 x, z = cs.coords
52 bz = C.Chebyshev(z, Nz, (0, 1), dealias=3 / 2)
53 D = C.Domain(cs, [C.RealFourier(x, Nx, (0, 2 * np.pi / kc), dealias=3 / 2), bz])
54 Dx = D.without(z)
55 psi, th = C.Field("psi", D), C.Field("theta", D)
56 tps = [C.Field(f"tpsi{i}", Dx) for i in range(1, 5)]
57 tth = [C.Field(f"tth{i}", Dx) for i in range(1, 3)]
58 Ra, Pr_ = C.Parameter("Ra", 27 * np.pi ** 4 / 4), C.Parameter("Pr", Pr)
59 G = C.SymmetrySubspace(C.Reflection(x, about=0), psi=C.ODD, theta=C.EVEN,
60 **{t.name: C.ODD for t in tps}, **{t.name: C.EVEN for t in tth})
61 pb = C.IVP([psi, th, *tps, *tth], name="RB stress-free (psi-theta)", symmetry=G)
62 u, w, om = C.diff(psi, z), -C.diff(psi, x), C.lap(psi)
63 pb.add_equation((C.dt(om) - Pr_ * C.lap(om) + Ra * Pr_ * C.diff(th, x)
64 + sum(C.lift(t, bz, -i) for i, t in enumerate(tps, 1)),
65 -(u * C.diff(om, x) + w * C.diff(om, z))))
66 pb.add_equation((C.dt(th) - C.lap(th) + C.diff(psi, x)
67 + sum(C.lift(t, bz, -i) for i, t in enumerate(tth, 1)),
68 -(u * C.diff(th, x) + w * C.diff(th, z))))
69 for side in ("left", "right"):
70 pb.add_equation((C.interp(psi, z, side), 0))
71 pb.add_equation((C.interp(C.diff(C.diff(psi, z), z), z, side), 0))
72 pb.add_equation((C.interp(th, z, side), 0))
73 s = ir_system(pb, ["Ra"])
74 return s, np.zeros(s.n)
75
76
77def nusselt_slope(Nx: int, Nz: int, Pr: float):
78 """d Nu / dε at onset: Nu − 1 = −(a²/2) ∂zθ̄_{h2}(0) with a² = −μ ΔRa / a3 (normal form)."""
79 s, x0 = rb_free_slip(Nx, Nz, Pr)
80 nf = normal_form(s, x0, kind="zero", parameter="Ra", kernel_tol=1e-6)
81 theta = s.solver.system.scatter(nf.h2)[1] # (x modes, Chebyshev); row 0 = x-mean
82 n = np.arange(theta.shape[1])
83 dz0 = np.sum(theta[0] * 2 * (-1.0) ** (n + 1) * n ** 2) # d/dz T_n(2z − 1) at z = 0
84 Rac = 27 * np.pi ** 4 / 4
85 return Rac * (-nf.mu / nf.a3) * (-0.5 * dz0), nf
86
87
88if __name__ == "__main__":
89 le, gl = sh_coefficients()
90 print(le.explain(), gl.explain(), sep="\n")
91 for Pr in (0.1, 1.0, 7.0):
92 slope, nf = nusselt_slope(8, 16, Pr)
93 print(f"stress-free RB, Pr = {Pr}: dNu/dε = {slope:.10f} (Malkus & Veronis 1958: 2), "
94 f"Ra_c = {nf.params['Ra']:.10f}, |λ0| = {nf.checks['|λ_0|']:.1e}")
95 # the script checks its references: Rayleigh (1916) Ra_c = 27 pi^4 / 4, Malkus & Veronis (1958) slope 2
96 assert abs(nf.params["Ra"] - 27 * np.pi ** 4 / 4) < 1e-9 and abs(slope - 2) < 1e-8, (Pr, slope)
97 assert abs(le.amplitude(0.1) ** 2 - 4 * 0.1 / 3) < 1e-13 # Cross & Hohenberg (1993): A^2 = 4r/3