russell_pde 2.8.0

Essential tools to solve partial differential equations; not a full-fledged PDE solver
Documentation

Russell PDE - Essential tools to solve partial differential equations; not a full-fledged PDE solver

documentation

This crate is part of Russell - Rust Scientific Library

Contents

Introduction

This library implements essential tools to solve partial differential equations (PDEs). It does not implement full-fledge PDE solvers for general problems and, hence, this library is quite limited.

A goal is to provide tools to test other crates such as russell_ode and russell_nonlin because they employ PDE problems as testing platforms.

Currently, simple finite differences operators are implemented, in addition to spectral collocation methods in 1D and 2D. The library also implements the transfinite mapping method to generate meshes on non-rectangular domains.

The linear systems are solved using the System Partitioning Strategy (SPS) or the Lagrange Multipliers Method (LMM).

Documentation

References

  1. Kopriva DA (2009) Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists and Engineers, Springer, 394p
  2. Trefethen LN (2000) Spectral Methods in MATLAB, SIAM, 165p
  3. Pozrikidis C (2014) Introduction to Finite and Spectral Element Methods Using MATLAB, 2nd ed, Chapman and Hall/CRC, 806p

Installation

This crate depends on russell_lab and russell_sparse, which require non-Rust high-performance libraries. See the main README file for the steps to install these dependencies.

Setting Cargo.toml

Crates.io

👆 Check the crate version and update your Cargo.toml accordingly:

[dependencies]
russell_pde = "*"

Optional features

The following (Rust) features are available:

  • intel_mkl: Use Intel MKL instead of OpenBLAS
  • local_sparse: Use locally compiled SuiteSparse and MUMPS
  • cudss: Enable the NVIDIA cuDSS GPU solver

Note that the main README file presents the steps to compile the required libraries according to each feature.

🌟 Examples

This section illustrates how to use russell_pde. See also:

Example 1: Solving 1D Poisson equation with Finite Differences

This example solves the Poisson equation in 1D using the Finite Difference Method (FDM):

  ∂²ϕ
- ——— = x    on [0, 1]
  ∂x²

With boundary conditions: ϕ(0) = 0, ϕ(1) = 0

The analytical solution is: ϕ(x) = (x - x³) / 6

use russell_lab::approx_eq;
use russell_pde::{Fdm1d, Grid1d, EssentialBcs1d, NaturalBcs1d, StrError};

fn main() -> Result<(), StrError> {
    // Define the problem domain and diffusion coefficient
    let (xmin, xmax) = (0.0, 1.0);
    let kx = 1.0;

    // Set up boundary conditions
    let mut ebcs = EssentialBcs1d::new();
    ebcs.set_homogeneous(); // ϕ(0) = 0, ϕ(1) = 0

    let nbcs = NaturalBcs1d::new();

    // Create uniform grid with 10 subdivisions
    let nx = 10;
    let grid = Grid1d::new_uniform(xmin, xmax, nx)?;

    // Create the solver
    let fdm = Fdm1d::new(grid, ebcs, nbcs, kx)?;

    // Define the source term f(x) = x
    let source = |x: f64| x;

    // Solve using System Partitioning Strategy
    let solution = fdm.solve_sps(0.0, source)?;

    // Verify against analytical solution
    let analytical = |x: f64| (x - x.powi(3)) / 6.0;
    fdm.for_each_coord(|m, x| {
        let error = f64::abs(solution[m] - analytical(x));
        approx_eq(solution[m], analytical(x), 1e-15);
        println!("x = {:.3}, ϕ = {:.6}, error = {:.e}", x, solution[m], error);
    });

    Ok(())
}

The output looks like this:

x = 0.000, ϕ = 0.000000, error = 0e0
x = 0.111, ϕ = 0.018290, error = 0e0
x = 0.222, ϕ = 0.035208, error = 0e0
x = 0.333, ϕ = 0.049383, error = 6.938893903907228e-18
x = 0.444, ϕ = 0.059442, error = 6.938893903907228e-18
x = 0.556, ϕ = 0.064015, error = 0e0
x = 0.667, ϕ = 0.061728, error = 0e0
x = 0.778, ϕ = 0.051212, error = 2.7755575615628914e-17
x = 0.889, ϕ = 0.031093, error = 3.469446951953614e-18
x = 1.000, ϕ = 0.000000, error = 0e0

Example 2: Solving 1D problems with Spectral Collocation

This example solves the Poisson equation in 1D using Spectral Collocation:

  ∂²ϕ
- ——— = x    on [0, 1]
  ∂x²

With boundary conditions: ϕ(0) = 0, ϕ(1) = 0

The analytical solution is: ϕ(x) = (x - x³) / 6

use russell_lab::approx_eq;
use russell_pde::{Spc1d, EssentialBcs1d, NaturalBcs1d, StrError};

fn main() -> Result<(), StrError> {
    let (xmin, xmax) = (0.0, 1.0);
    let kx = 1.0;

    // Set up boundary conditions
    let mut ebcs = EssentialBcs1d::new();
    ebcs.set_homogeneous();
    let nbcs = NaturalBcs1d::new();

    // Create spectral collocation solver with N=8 polynomial degree
    let nx = 8;
    let spc = Spc1d::new(xmin, xmax, nx, ebcs, nbcs, kx)?;

    // Solve the problem
    let source = |x: f64| x;
    let solution = spc.solve_sps(0.0, source)?;

    // Verify against analytical solution
    let analytical = |x: f64| (x - x.powi(3)) / 6.0;
    spc.for_each_coord(|m, x| {
        let error = f64::abs(solution[m] - analytical(x));
        approx_eq(solution[m], analytical(x), 1e-15);
        println!("x = {:.3}, ϕ = {:.6}, error = {:.e}", x, solution[m], error);
    });

    Ok(())
}

Output:

x = 0.000, ϕ = 0.000000, error = 0e0
x = 0.050, ϕ = 0.008232, error = 0e0
x = 0.188, ϕ = 0.030264, error = 1.0408340855860843e-17
x = 0.389, ϕ = 0.054999, error = 4.163336342344337e-17
x = 0.611, ϕ = 0.063812, error = 4.163336342344337e-17
x = 0.812, ϕ = 0.046144, error = 3.469446951953614e-17
x = 0.950, ϕ = 0.015300, error = 2.949029909160572e-17
x = 1.000, ϕ = 0.000000, error = 0e0

Example 3: Spectral collocation in 2D with transfinite mapping

Example: Solving a 2D Poisson equation on a rotated square domain

This example demonstrates the use of SpcMap2d (spectral collocation with transfinite mapping) to solve the Poisson equation:

  -k · ∇²u = f    on a unit square rotated by angle α
  u = g           on the boundary (Dirichlet conditions)

The analytical solution used for verification is:

  u(x,y) = sin(π·x·cos(α) + π·y·sin(α)) · exp(π·y·cos(α) - π·x·sin(α))

The domain is mapped from the reference square (r,s) ∈ [-1,1]×[-1,1] to the physical rotated square via transfinite interpolation.

See the code

The output looks like this:

N = 20 max(err) = 3.90799e-14

And the plot looks like this:

Solution

Test Problems

The figures below are generated by the tests in the tests directory. Each test verifies the numerical solution against an analytical solution with a given tolerance.

Problem 01 — 2D Poisson with mixed BCs

On [0,1]², two cases: (a) homogeneous Dirichlet all sides, (b) mixed Neumann/Dirichlet. test code (FDM) · test code (SPC)

FDM SPC contour (case a) SPC contour (case b)
prob01 FDM prob01 SPC a prob01 SPC b

Problem 02 — 2D Poisson ϕ = y·sin(πx)

On [0,1]² with Dirichlet BCs. Source: s = −π²y·sin(πx). test code

prob02 FDM

Problem 03 — 2D Helmholtz/Poisson with mixed BCs

On [0,1]² with five combinations of Dirichlet/Neumann BCs (DDDD, NNDD, NDND, DNND, DDNN), tested with and without Helmholtz term. test code (FDM) · test code (SPC)

FDM (case a) FDM (case b)
prob03 FDM a prob03 FDM b
SPC DDDD contour SPC DDDD surface
prob03 SPC DDDD hz prob03 SPC DDDD hz surface
SPC Map DDDD contour SPC Map DDDD surface
prob03 SPC Map DDDD hz prob03 SPC Map DDDD hz surface

Problem 04 — 2D Poisson on [-1,1]² (benchmark)

Laplace with source s = 1 and homogeneous Dirichlet BCs. test code (FDM) · test code (SPC)

FDM SPC
prob04 FDM prob04 SPC

Problem 05 — 2D Poisson polynomial

On [-1,1]², source s = −6x, analytical ϕ = 1 + x³. Dirichlet on x-edges, Neumann on y-edges. test code (FDM) · test code (SPC)

FDM grid (a) FDM solution (b) SPC grid (a) SPC solution (b)
prob05 FDM a prob05 FDM b prob05 SPC a prob05 SPC b

Problem 06 — 2D Poisson ϕ = tanh(1−x+y)

On [-1,1]² with mixed Dirichlet/Neumann BCs. Source is the Laplacian of the analytical solution. test code (FDM) · test code (SPC)

FDM grid (a) FDM solution (b) SPC grid (a) SPC solution (b)
prob06 FDM a prob06 FDM b prob06 SPC a prob06 SPC b

Problem 08 — 2D Poisson on curvilinear domains

Kopriva benchmark 7.1.4: ∇²ϕ = −16·ln(r)/r²·sin(4θ) on quarter-annulus or perforated lozenge. test code

Ring domain Lozenge domain
prob08 ring prob08 lozenge

Problem 09 — 2D Laplace (potential flow)

Kopriva benchmark 7.1.5: Potential flow around a cylinder on a half-ring domain. test code

prob09 SPC Map

1D Problem 02 — Helmholtz heat conduction-convection

−k·∂²ϕ/∂x² + α·ϕ = α·ϕ∞ on [0,0.05] with fixed left temperature (320°C) and insulated right. test code (FDM) · test code (SPC)

FDM SPC
1d prob02 FDM 1d prob02 SPC

1D Problem 03 — Helmholtz with flux BC

−∂²ϕ/∂x² + ϕ = x² on [0,1] with fixed left (2°C) and flux input (−3 W) at right. test code (FDM) · test code (SPC)

FDM SPC
1d prob03 FDM 1d prob03 SPC

1D Problem 04 — Trefethen programs 13 and 33

Program 13: homogeneous Dirichlet; Program 33: Neumann left + Dirichlet right. PDE: ∂²ϕ/∂x² = exp(4x) on [−1,1]. test code

Program 13 (Dirichlet) Program 33 (Neumann+Dirichlet)
1d prob04a SPC 1d prob04b SPC

1D Problem 05 — Pozrikidis Helmholtz

∂²ϕ/∂x² + β²·ϕ = 0 on [0,L] with Neumann left (g0=1) and Dirichlet right (ϕL=0.2). test code (FDM) · test code (SPC)

FDM SPC
1d prob05 FDM 1d prob05 SPC

Metrics — Curvilinear coordinates

Covariant and contravariant basis vectors on a curved quadrilateral domain. test code

Metrics

For developers

  • This crate is pure Rust with no C dependencies
  • Run the examples with cargo run --example <name>