volterra_solver/
runner_dec.rs

1//! Dry active nematic runner on 2D DEC meshes.
2//!
3//! Solves the Beris-Edwards equation (no flow coupling) on a triangulated
4//! 2-manifold using RK4 time integration.
5//!
6//! ## Governing equation
7//!
8//! ```text
9//! dQ/dt = gamma_r * H
10//! ```
11//!
12//! where H is the molecular field from [`volterra_dec::molecular_field_dec`].
13//! Flow coupling (Stokes solver) is deferred to Sub-project B.
14
15use cartan_core::Manifold;
16use cartan_dec::Operators;
17use volterra_core::ActiveNematicParams;
18use volterra_dec::{molecular_field_dec, QFieldDec};
19
20/// Per-snapshot statistics from a DEC simulation.
21#[derive(Debug, Clone)]
22pub struct SnapStatsDec {
23    /// Simulation time.
24    pub time: f64,
25    /// Mean scalar order parameter.
26    pub mean_s: f64,
27    /// Number of vertices in the mesh.
28    pub n_vertices: usize,
29}
30
31/// Beris-Edwards RHS for the dry active nematic (no flow).
32///
33/// dQ/dt = gamma_r * H
34fn beris_edwards_rhs_dec<M: Manifold>(
35    q: &QFieldDec,
36    params: &ActiveNematicParams,
37    ops: &Operators<M, 3, 2>,
38    curvature_correction: Option<&dyn Fn(usize) -> [[f64; 3]; 3]>,
39) -> QFieldDec {
40    let h = molecular_field_dec(q, params, ops, curvature_correction);
41    h.scale(params.gamma_r)
42}
43
44/// One RK4 step for the Q-tensor field on a DEC mesh.
45fn rk4_step<M: Manifold>(
46    q: &QFieldDec,
47    dt: f64,
48    params: &ActiveNematicParams,
49    ops: &Operators<M, 3, 2>,
50    curvature_correction: Option<&dyn Fn(usize) -> [[f64; 3]; 3]>,
51) -> QFieldDec {
52    let k1 = beris_edwards_rhs_dec(q, params, ops, curvature_correction);
53
54    let q2 = q.add(&k1.scale(0.5 * dt));
55    let k2 = beris_edwards_rhs_dec(&q2, params, ops, curvature_correction);
56
57    let q3 = q.add(&k2.scale(0.5 * dt));
58    let k3 = beris_edwards_rhs_dec(&q3, params, ops, curvature_correction);
59
60    let q4 = q.add(&k3.scale(dt));
61    let k4 = beris_edwards_rhs_dec(&q4, params, ops, curvature_correction);
62
63    // q_next = q + (dt/6)(k1 + 2*k2 + 2*k3 + k4)
64    let rhs = k1
65        .add(&k2.scale(2.0))
66        .add(&k3.scale(2.0))
67        .add(&k4);
68    q.add(&rhs.scale(dt / 6.0))
69}
70
71/// Run the dry active nematic on a 2D DEC mesh.
72///
73/// Evolves Q via RK4 for `n_steps` time steps with no hydrodynamic flow.
74/// Records statistics every `snap_every` steps.
75///
76/// # Arguments
77///
78/// * `q_init` - Initial Q-tensor field on the mesh.
79/// * `params` - Physical and numerical parameters (uses `k_r`, `gamma_r`,
80///   `a_landau`, `zeta_eff`, `c_landau`, `dt`).
81/// * `ops` - Precomputed DEC operators (Laplace-Beltrami, Hodge stars).
82/// * `curvature_correction` - Optional Weitzenboeck endomorphism for curved
83///   surfaces. Pass `None` for flat meshes.
84/// * `n_steps` - Total number of time steps.
85/// * `snap_every` - Steps between snapshot statistics.
86///
87/// # Returns
88///
89/// `(q_final, stats)`: final Q-field and one [`SnapStatsDec`] per snapshot.
90pub fn run_dry_active_nematic_dec<M: Manifold>(
91    q_init: &QFieldDec,
92    params: &ActiveNematicParams,
93    ops: &Operators<M, 3, 2>,
94    curvature_correction: Option<&dyn Fn(usize) -> [[f64; 3]; 3]>,
95    n_steps: usize,
96    snap_every: usize,
97) -> (QFieldDec, Vec<SnapStatsDec>) {
98    let mut q = q_init.clone();
99    let mut stats = Vec::new();
100
101    for step in 0..=n_steps {
102        if step % snap_every == 0 {
103            stats.push(SnapStatsDec {
104                time: step as f64 * params.dt,
105                mean_s: q.mean_order_param(),
106                n_vertices: q.n_vertices,
107            });
108        }
109        if step < n_steps {
110            q = rk4_step(&q, params.dt, params, ops, curvature_correction);
111        }
112    }
113
114    (q, stats)
115}
volterra_solver/runner_dec.rs

volterra_solver/
runner_dec.rs
