Numerical Solver

Numerical Solver#

The soliton_solver implements a theory-agnostic numerical solver for 2D PDEs describing topological solitons. This document describes the finite-difference method, the arrested Newton flow minimization algorithm, and solver configuration.

Overview#

The numerical core solves the relaxation equation:

\[ \ddot{\phi}(t) = -\nabla_\phi E_h[\phi] \]

where \(E_h[\phi]\) is a discretized energy functional and \(\phi\) is the field configuration. The goal is to find local minima of \(E_h\).

Key features:

  • Second-order relaxation — Velocity helps escape shallow regions of the energy landscape

  • Flow arrest — Velocity resets if energy increases, preventing oscillations

  • Theory-agnostic — Core solver doesn’t know specific physics

  • GPU-accelerated — All kernels run on NVIDIA GPUs via Numba CUDA

Finite Difference Method#

Fields are discretized on a 2D Cartesian grid with uniform spacing. Spatial derivatives are computed via fourth-order central finite differences.

Grid setup#

The computational domain is discretized into \(n_x \times n_y\) lattice points:

  • Domain: \([0, L_x] \times [0, L_y]\) in physical space

  • Grid points: \(n_x \times n_y\) lattice sites

  • Spacing: \(\Delta x = L_x / n_x\), \(\Delta y = L_y / n_y\)

  • Halo: Additional boundary points for stencil computations

Parameters:

params = theory.params.default_params(
    xlen=256,        # Number of lattice points x
    ylen=256,        # Number of lattice points y
    xsize=80.0,      # Physical domain size x
    ysize=80.0,      # Physical domain size y
    halo=2,          # Halo width (default: 2)
)
rp = params.resolved()
print(f"Lattice spacing: dx={rp.lsx}, dy={rp.lsy}")
print(f"Grid volume: {rp.grid_volume}")

Halo width — Boundary region used by finite-difference stencils. For fourth-order stencils with ±2 offset, halo=2 is standard.

Stencils#

First derivatives#

Fourth-order central difference for \(\partial_x f\):

\[ \frac{\partial f}{\partial x}\bigg|_{i,j} \approx \frac{1}{12\Delta x}(f_{i-2,j} - 8f_{i-1,j} + 8f_{i+1,j} - f_{i+2,j}) \]

Implementation (from derivatives.py):

@cuda.jit(device=True)
def compute_derivative_first(d1fd1x, Field, a, x, y, p_i, p_f):
    """Compute first derivatives at (x, y) for field component a."""
    xlen = p_i[0]
    ylen = p_i[1]
    halo = p_i[2]
    lsx = p_f[2]  # Lattice spacing
    
    # Check interior (not in halo)
    if x > halo - 1 and x < xlen - halo:
        d1fd1x[idx_d1(0, a, x, y, p_i)] = (
            (1.0/12.0) * Field[idx_field(a, x-2, y, p_i)]
            - (2.0/3.0) * Field[idx_field(a, x-1, y, p_i)]
            + (2.0/3.0) * Field[idx_field(a, x+1, y, p_i)]
            - (1.0/12.0) * Field[idx_field(a, x+2, y, p_i)]
        ) / lsx
    else:
        d1fd1x[idx_d1(0, a, x, y, p_i)] = 0.0

Accuracy — Fourth-order: error \(\mathcal{O}(\Delta x^4)\)

Second derivatives#

Fourth-order Laplacian \(\partial_{xx} f\):

\[ \frac{\partial^2 f}{\partial x^2}\bigg|_{i,j} \approx \frac{1}{(\Delta x)^2}\left(-\frac{1}{12}f_{i-2,j} + \frac{4}{3}f_{i-1,j} - \frac{5}{2}f_{i,j} + \frac{4}{3}f_{i+1,j} - \frac{1}{12}f_{i+2,j}\right) \]

Similarly for \(\partial_{yy}\) and mixed derivatives \(\partial_{xy}\) (via chain rule from diagonals).

Boundary conditions#

Fields satisfy Dirichlet boundary conditions on the computational domain. The halo region fixes boundary values.

For a field with halo=2:

  • Interior points: 2 ≤ x < nx - 2

  • Boundary points: 0 ≤ x < 2 or nx - 2 ≤ x < nx

Stability and accuracy#

Stability criterion — The Courant condition for explicit time-stepping:

\[ \Delta t \leq C \cdot \min(\Delta x, \Delta y)^2 \]

where \(C\) is a theory-dependent constant (default Courant number ≈ 0.5).

Accuracy — Overall fourth-order in space, second-order in time (RK4).

Arrested Newton Flow Algorithm#

ANF is a second-order relaxation scheme combining fast convergence with stability.

Algorithm description#

The continuous ANF equation is:

\[ \ddot{\phi} = -\nabla_\phi E_h[\phi] \]

This is a system of ODEs in velocity \(v = \dot{\phi}\) and position \(\phi\):

\[\begin{split} \begin{aligned} \dot{\phi} &= v \\ \dot{v} &= -\nabla_\phi E_h[\phi] \end{aligned} \end{split}\]

Arrest rule — If energy increases, velocity is reset:

\[ v^{n+1} \leftarrow 0 \quad \text{if} \quad E_h[\phi^{n+1}] > E_h[\phi^n] \]

Convergence properties#

  • Second-order dynamics — Velocity accelerates motion; useful for escaping shallow regions

  • Energy stable — Arrest prevents oscillations; energy is non-increasing (except at resets)

  • Faster than gradient descent — Particularly for multi-soliton configurations

  • Convergence criterion — Stop when \(\|\nabla E\|_\infty < \epsilon\) (convergence tolerance)

GPU implementation: RK4 integration#

One full ANF step consists of:

  1. Compute gradient — \(\nabla E\) at all grid points via finite differences

  2. RK4 step — Advance \((φ, v)\) using fourth-order Runge-Kutta

  3. Arrest check — Reset velocity if energy increases

  4. Convergence check — Test \(\|\nabla E\|_\infty\) against tolerance

RK4 stages#

The RK4 method for \(\dot{\mathbf{y}} = \mathbf{f}(t, \mathbf{y})\) is:

\[\begin{split} \begin{aligned} \mathbf{k}_1 &= \mathbf{f}(t, \mathbf{y}) \\ \mathbf{k}_2 &= \mathbf{f}(t + \Delta t/2, \mathbf{y} + \Delta t \mathbf{k}_1/2) \\ \mathbf{k}_3 &= \mathbf{f}(t + \Delta t/2, \mathbf{y} + \Delta t \mathbf{k}_2/2) \\ \mathbf{k}_4 &= \mathbf{f}(t + \Delta t, \mathbf{y} + \Delta t \mathbf{k}_3) \\ \mathbf{y}^{n+1} &= \mathbf{y}^n + \frac{\Delta t}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4) \end{aligned} \end{split}\]

For the ANF system, \(\mathbf{y} = (\phi, v)\) and:

\[\begin{split} \mathbf{f}(t, \mathbf{y}) = \begin{pmatrix} v \\ -\nabla_\phi E_h[\phi] \end{pmatrix} \end{split}\]

The implementation launches four gradient computation kernels (one per RK4 stage) and updates field and velocity buffers.

Solver Configuration#

The solver is configured via the Params and ResolvedParams classes.

Core parameters#

Grid and domain:

Params(
    xlen=256,              # Grid points in x
    ylen=256,              # Grid points in y
    xsize=80.0,            # Physical size in x
    ysize=80.0,            # Physical size in y
    halo=2,                # Halo width for stencils
)

Time stepping:

Params(
    time_step=0.01,        # Explicit time step (if set)
    courant=0.5,           # Courant number for auto time stepping
)

If time_step is None, it is computed from the Courant number:

\[ \Delta t = \text{courant} \cdot \min(\Delta x, \Delta y)^2 \]

Arrested Newton flow:

Params(
    newtonflow=True,       # Enable Newton flow (vs. gradient descent)
    killkinen=True,        # Reset velocity on energy increase
)

Constraints:

Params(
    unit_magnetization=False,  # Enforce unit norm on magnetization (some theories)
)

Resolved parameters#

After calling .resolved(), derived parameters are computed:

rp = params.resolved()

print(rp.xlen, rp.ylen)              # Grid dimensions
print(rp.lsx, rp.lsy)                # Lattice spacings
print(rp.grid_volume)                # Cell area
print(rp.time_step)                  # Actual time step
print(rp.dim_grid)                   # Total grid points
print(rp.dim_fields)                 # Total field values

Example configuration#

from soliton_solver.theories import load_theory

theory = load_theory("Chiral magnet")

# Create base parameters
params = theory.params.default_params(
    xlen=512, ylen=512,      # Finer grid
    xsize=20.0, ysize=20.0,  # Larger domain
    courant=0.3,             # Conservative time step
    newtonflow=True,
    killkinen=True,
)

# Resolve and create simulation
sim = Simulation(params, theory)
sim.initialize({"mode": "ground"})

# Relax via arrested Newton flow
from soliton_solver.core.simulation import Simulation
energy = sim.observables_mod.compute_energy(
    sim.Field, sim.d1fd1x, sim.en, sim.entmp, 
    sim.gridsum_partial, sim.p_i_d, sim.p_f_d, 
    sim.p_i_h, sim.p_f_h
)

# Advance one step
new_energy, err = sim.step(prev_energy=energy)
print(f"Energy: {new_energy}, Gradient norm: {err}")

Performance considerations#

Time step selection#

  • Too small — Slow convergence, many unnecessary steps

  • Too large — Instability, divergence, energy oscillations

The Courant condition provides automatic stability; manual override via time_step should be conservative.

Grid resolution#

  • Coarse grids — Fast but inaccurate; may miss fine soliton details

  • Fine grids — Accurate but slow; scales as \(O(n_x n_y)\) per step

Recommended: start coarse, refine around features of interest.

Convergence tolerance#

Typically \(10^{-4}\) to \(10^{-6}\) depending on application:

# Run until convergence
energy = initial_energy
for step in range(max_steps):
    new_energy, err = sim.step(energy)
    energy = new_energy
    
    if err < 1e-4:
        print(f"Converged in {step} steps")
        break
    
    if step % 100 == 0:
        print(f"Step {step}: E={energy:.6e}, ||∇E||_∞={err:.6e}")

Memory usage#

Field arrays are stored in GPU device memory:

\[ \text{Memory per field} = n_x \cdot n_y \cdot 8 \text{ bytes} \quad (\text{float64}) \]

For multiple fields and buffers (Field, Velocity, derivatives, RK4 stages):

\[ \text{Total} \approx 10 \times n_x \cdot n_y \cdot 8 \text{ bytes} \]

Example: 1024×1024 grid with 4 fields ≈ 80 MB on GPU.

Typical runtime#

A single RK4 step on NVIDIA A100:

  • 512×512, 2–4 fields: ~1–2 ms

  • 1024×1024, 4–8 fields: ~5–10 ms

Convergence typically requires 1,000–10,000 steps depending on configuration complexity.

Contents