Swift–Hohenberg: onset of rolls, weakly nonlinear amplitude, roll branch by continuation, 2D rolls

Family: patterns · Script: examples/patterns/swift_hohenberg.py

Obtained

threshold (r, q) on a single pencil (q a parameter, RealFourier); A(r = 0.1) = 0.3652433 (A² − 4r/3 → r²/144); periodic NLBVP (declared phase) = odd subspace to 1e-12; 2D IVP on the branch to 1e-9

Reference

Swift & Hohenberg 1977; Cross & Hohenberg 1993

 
du/dt = r u - (1 + lap)^2 u - u^3                         (Swift & Hohenberg 1977)

(a) Linear threshold. A normal mode exp(i q x + s t) of the state u = 0 grows at
    s = r - (1 - q^2)^2, so the neutral curve is r = (1 - q^2)^2 and the onset is r_c = 0 at
    q_c = 1 (Cross & Hohenberg 1993, Rev. Mod. Phys. 65, 851, section on the SH model). Computed
    here by linearize() of the IVP around u = 0 and one parametric pencil in (r, q): the Fourier
    wavenumber q is a declared pencil parameter (``wavenumbers={"x": q}``), so critical() minimizes
    over q without rediscretizing (one pencil per q before). RealFourier: each group
    holds the (cos, sin) pair, a double eigenvalue (translation invariance); track() follows a
    declared member of the eigenspace (it raised LinAlgError before).

(b) Roll branch. Steady rolls u = A sin x + B sin 3x + ... obey, projecting on sin x and sin 3x,
    B = A^3 / (4 (64 - r)) and r = 3 A^2 / 4 - 3 A B / 4 + ..., i.e.

        A^2 = 4 r / 3 + r^2 / 144 + O(r^3)

    (leading term: Cross & Hohenberg 1993; the r^2/144 correction is the third-harmonic feedback
    derived above). The rolls are solved as a periodic NLBVP (PeriodicNewtonSolver: dense coupled
    Jacobian, translation fixed by a declared phase condition, the -sin coefficient of mode 1 set
    to zero, with its multiplier reported), and, independently, in the odd subspace (the banded
    NewtonSolver refuses the periodic form: "varies along a Fourier direction"), *exactly
    equivalent* for the periodic rolls and followed in r by continuation: on x in [0, pi] with
    u = u'' = 0 at both ends, every solution extends to an odd 2 pi-periodic one (sin(n x) modes;
    (1 + d_xx)^2 and u^3 preserve oddness), and conversely. Chebyshev in x, four BCs for the
    fourth-order operator. The amplitude A = (2/pi) int_0^pi u sin x dx is an extra unknown of the
    NLBVP, so it is carried by Newton and by the continuation.

(c) 2D rolls. The same equation as an IVP on a doubly periodic 2 pi x 2 pi box (RealFourier x
    RealFourier, dealias = 2 declared for the cubic term), started from a mixture of x- and
    y-rolls: the y-rolls decay, the x-rolls saturate at the amplitude of the branch (b).

Run:  python examples/patterns/swift_hohenberg.py

Run it: python examples/patterns/swift_hohenberg.py (in the conatus environment).

  1"""Swift–Hohenberg: onset of rolls, weakly nonlinear amplitude, roll branch by continuation, 2D rolls.
  2
  3    du/dt = r u - (1 + lap)^2 u - u^3                         (Swift & Hohenberg 1977)
  4
  5(a) Linear threshold. A normal mode exp(i q x + s t) of the state u = 0 grows at
  6    s = r - (1 - q^2)^2, so the neutral curve is r = (1 - q^2)^2 and the onset is r_c = 0 at
  7    q_c = 1 (Cross & Hohenberg 1993, Rev. Mod. Phys. 65, 851, section on the SH model). Computed
  8    here by linearize() of the IVP around u = 0 and one parametric pencil in (r, q): the Fourier
  9    wavenumber q is a declared pencil parameter (``wavenumbers={"x": q}``), so critical() minimizes
 10    over q without rediscretizing (one pencil per q before). RealFourier: each group
 11    holds the (cos, sin) pair, a double eigenvalue (translation invariance); track() follows a
 12    declared member of the eigenspace (it raised LinAlgError before).
 13
 14(b) Roll branch. Steady rolls u = A sin x + B sin 3x + ... obey, projecting on sin x and sin 3x,
 15    B = A^3 / (4 (64 - r)) and r = 3 A^2 / 4 - 3 A B / 4 + ..., i.e.
 16
 17        A^2 = 4 r / 3 + r^2 / 144 + O(r^3)
 18
 19    (leading term: Cross & Hohenberg 1993; the r^2/144 correction is the third-harmonic feedback
 20    derived above). The rolls are solved as a periodic NLBVP (PeriodicNewtonSolver: dense coupled
 21    Jacobian, translation fixed by a declared phase condition, the -sin coefficient of mode 1 set
 22    to zero, with its multiplier reported), and, independently, in the odd subspace (the banded
 23    NewtonSolver refuses the periodic form: "varies along a Fourier direction"), *exactly
 24    equivalent* for the periodic rolls and followed in r by continuation: on x in [0, pi] with
 25    u = u'' = 0 at both ends, every solution extends to an odd 2 pi-periodic one (sin(n x) modes;
 26    (1 + d_xx)^2 and u^3 preserve oddness), and conversely. Chebyshev in x, four BCs for the
 27    fourth-order operator. The amplitude A = (2/pi) int_0^pi u sin x dx is an extra unknown of the
 28    NLBVP, so it is carried by Newton and by the continuation.
 29
 30(c) 2D rolls. The same equation as an IVP on a doubly periodic 2 pi x 2 pi box (RealFourier x
 31    RealFourier, dealias = 2 declared for the cubic term), started from a mixture of x- and
 32    y-rolls: the y-rolls decay, the x-rolls saturate at the amplitude of the branch (b).
 33
 34Run:  python examples/patterns/swift_hohenberg.py
 35"""
 36import numpy as np
 37
 38from conatus.expert import (IVP, NLBVP, CartesianCoordinates, Chebyshev, ContinuationSettings,
 39                            Domain, Field, IVPSolver, NewtonSolver, NLBVPContinuation, Parameter,
 40                            ParametricEVP, PeriodicNewtonSolver, RealFourier, critical, diff,
 41                            discretize, dt, integ, interp, lap, linearize)
 42from conatus.ir import Zero
 43from conatus.solvers import PhaseCondition
 44
 45
 46def swift_hohenberg_ivp(domain, lap_):
 47    """du/dt - r u + (1 + lap)^2 u = -u^3 on ``domain`` (lap_ = the Laplacian of that domain)."""
 48    u, r = Field("u", domain), Parameter("r", 0.1)
 49
 50    def one_plus_lap(f):
 51        return f + lap_(f)
 52
 53    ivp = IVP([u], name="Swift-Hohenberg")
 54    ivp.add_equation((dt(u) - r * u + one_plus_lap(one_plus_lap(u)), -u * u * u))
 55    return ivp, u, r
 56
 57
 58# ------------------------------------------------------------------ (a) linear threshold
 59def growth_pencil():
 60    """Pencil in (r, q) of the perturbation of u = 0; q = wavenumber of x (a pencil parameter,
 61    so the period of the domain plays no role)."""
 62    cs = CartesianCoordinates("x")
 63    x = cs["x"]
 64    D = Domain(cs, [RealFourier(x, 8, (0, 2 * np.pi), dealias=1)])
 65    ivp, u, r = swift_hohenberg_ivp(D, lambda f: diff(diff(f, x), x))
 66    evp = linearize(ivp, around={u: Zero(())}, eigenvalue=Parameter("s", dtype=complex))
 67    return ParametricEVP(evp, [r], wavenumbers={"x": Parameter("q", 1.0)}), evp
 68
 69
 70def onset():
 71    """(q_c, r_c): neutral r of the least stable mode, minimized over q on one pencil."""
 72    cp = critical(growth_pencil()[0], "r", bracket=(-0.5, 1.0), over={"q": (0.5, 1.5)}, xtol=1e-10)
 73    return cp.over["q"], cp.value
 74
 75
 76# ------------------------------------------------------------------ (b) roll branch (NLBVP)
 77def rolls(N: int):
 78    """Steady rolls in the odd subspace: x in [0, pi], u = u'' = 0; unknowns u and amplitude A."""
 79    cs = CartesianCoordinates("x")
 80    x = cs["x"]
 81    D = Domain(cs, [Chebyshev(x, N, (0, np.pi), dealias=1)])
 82    u, A, sinx, r = Field("u", D), Field("A", D.without(x)), Field("sinx", D), Parameter("r", 0.1)
 83
 84    def one_plus_lap(f):
 85        return f + diff(diff(f, x), x)
 86
 87    pb = NLBVP([u, A], name="Swift-Hohenberg rolls (odd subspace)")
 88    pb.add_equation((r * u - one_plus_lap(one_plus_lap(u)) - u * u * u, 0))
 89    pb.add_equation((A - (2 / np.pi) * integ(sinx * u), 0))      # A = first sine coefficient
 90    uxx = diff(diff(u, x), x)
 91    pb.add_bc(interp(u, x, 0), interp(u, x, np.pi), interp(uxx, x, 0), interp(uxx, x, np.pi), into=0)
 92    return pb, u, {sinx: np.sin}
 93
 94
 95def rolls_periodic(N: int, r: float):
 96    """Steady rolls on the 2 pi-periodic line (RealFourier, dealias 2 for the cubic term), phase
 97    fixed by msin_1 = 0; returns the solution (first harmonic via .harmonic("u", x, 1))."""
 98    cs = CartesianCoordinates("x")
 99    x = cs["x"]
100    u, rr = Field("u", Domain(cs, [RealFourier(x, N, (0, 2 * np.pi), dealias=2)])), Parameter("r", r)
101
102    def one_plus_lap(f):
103        return f + diff(diff(f, x), x)
104
105    pb = NLBVP([u], name="Swift-Hohenberg rolls (periodic)")
106    pb.add_equation((rr * u - one_plus_lap(one_plus_lap(u)) - u * u * u, 0))
107    solver = PeriodicNewtonSolver(pb, tol=1e-13, damping=1.0, max_iterations=20,
108                                  phase=PhaseCondition(u, x, kind="sin1"))
109    return solver.solve({u: lambda X: np.sqrt(4 * r / 3) * np.cos(X)}), x
110
111
112def roll_solver(N: int, r: float):
113    """Newton for the rolls at r, started from the leading-order amplitude sqrt(4 r / 3)."""
114    pb, u, ncc = rolls(N)
115    solver = NewtonSolver(pb, tol=1e-13, damping=1.0, max_iterations=30, params={"r": r}, ncc=ncc)
116    return solver, solver.solve({u: lambda X: np.sqrt(4 * r / 3) * np.sin(X)})
117
118
119def roll_branch(N: int, r0: float, r1: float, ds_max: float):
120    """Pseudo-arclength continuation of the rolls in r from r0 to r1: [(r, A, residual)]."""
121    solver, sol = roll_solver(N, r0)
122    cfg = ContinuationSettings(ds=0.02, ds_min=1e-6, ds_max=ds_max, adapt=1.5, fast_iterations=3,
123                               tolerance=1e-12, max_newton=15, max_points=40, dp_fd=1e-6)
124    cont = NLBVPContinuation(solver, "r", cfg)
125    branch = cont.run(sol.coeffs, r0, +1, p_bounds=(r0, r1))
126    pts = [(q.p, float(np.real(cont.coeffs_of(q.x)["A"]).ravel()[0]), q.residual) for q in branch.points]
127    return pts, branch
128
129
130# ------------------------------------------------------------------ (c) 2D rolls (IVP)
131def rolls_2d(N: int, r: float, timestep: float, steps: int):
132    """IVP on the 2 pi-periodic square; returns (A_x, A_y, |du/dt| at the end, solver)."""
133    cs = CartesianCoordinates("x", "y")
134    x, y = cs.coords
135    D = Domain(cs, [RealFourier(x, N, (0, 2 * np.pi), dealias=2), RealFourier(y, N, (0, 2 * np.pi), dealias=2)])
136    ivp, u, _ = swift_hohenberg_ivp(D, lap)
137    solver = IVPSolver(discretize(ivp, params={"r": r}), "SBDF2", timestep)
138    solver.set_state({u: lambda X, Y: 0.05 * np.cos(X) + 0.04 * np.cos(Y) + 0.01 * np.sin(2 * X + Y)})
139    solver.run(steps - 1)
140    before = solver.grid("u")
141    solver.run(1)
142    g = solver.grid("u")
143    e = np.exp(-1j * 2 * np.pi * np.arange(N) / N)
144    A_x = 2 * abs(np.mean(g * e[:, None]))           # |cos x, sin x| content (phase free)
145    A_y = 2 * abs(np.mean(g * e[None, :]))
146    return A_x, A_y, float(np.abs(g - before).max() / timestep), solver
147
148
149if __name__ == "__main__":
150    q_c, r_c = onset()
151    print(f"(a) onset: q_c = {q_c:.8f}, r_c = {r_c:.2e}  (theory: q_c = 1, r_c = 0)")
152    assert abs(q_c - 1) < 1e-4 and abs(r_c) < 1e-8, (q_c, r_c)    # Swift & Hohenberg (1977)
153    P, evp = growth_pencil()
154    print(evp.explain())
155    print(P.explain())
156    print("(b) rolls at r = 0.1, convergence in N (A = first sine coefficient):")
157    for N in (16, 20, 24, 32):
158        _, sol = roll_solver(N, 0.1)
159        print(f"  N = {N}: A = {sol.value('A'):.15f}, Newton residual {sol.residual:.1e}")
160    A = sol.value("A")                                             # N = 32, converged in N
161    sol, x = rolls_periodic(16, 0.1)
162    print(f"  periodic NLBVP (N = 16, phase msin_1 = 0): A = {sol.harmonic('u', x, 1):.15f}")
163    assert abs(sol.harmonic("u", x, 1) - A) < 1e-12               # periodic NLBVP = odd subspace
164    print(sol.explain())
165    pts, branch = roll_branch(24, 0.02, 0.5, 0.05)
166    print(branch.explain())
167    print("  r         A              A^2 - 4r/3     (A^2 - 4r/3)/r^2  (theory 1/144 = 0.006944)")
168    for r, A_r, res in pts:
169        gap = A_r * A_r - 4 * r / 3
170        print(f"  {r:.5f}  {A_r:.12f}  {gap:.3e}      {gap / r ** 2:.6f}   residual {res:.1e}")
171        assert abs(gap) < 0.01 * r ** 2, (r, gap)                 # A^2 = 4r/3 + O(r^2) (Cross & Hohenberg 1993)
172    A_x, A_y, rate, s = rolls_2d(16, 0.1, 0.5, 800)
173    print(f"(c) 2D IVP at r = 0.1, t = {s.t}: A_x = {A_x:.12f}, A_y = {A_y:.1e}, |du/dt| = {rate:.1e}")
174    assert rate < 1e-8 and A_y < 1e-6 and abs(A_x - A) < 1e-9, (rate, A_y, A_x, A)   # 2D IVP = 1D rolls
175    print(s.explain())