Crate cons_laws 

Deterministic particle methods for 1D conservation laws.

This crate implements the deterministic particle schemes described in the article Entropy solutions of non-local scalar conservation laws with congestion via deterministic particle method, E. Radici, F. Stra, (2021), https://arxiv.org/abs/2107.10760.

You can cite the article as

@online{RadiciStra2021,
    title={Entropy solutions of non-local scalar conservation laws with congestion via deterministic particle method},
    author={Emanuela Radici and Federico Stra},
    year={2021},
    eprint={2107.10760},
    archivePrefix={arXiv},
    primaryClass={math.AP},
    url={https://arxiv.org/abs/2107.10760}
}

This is a reimplementation in Rust of the Julia package ConservationLawsParticles.jl.

What’s the goal?

The conservation law to be solved is

$$ \partial_t\rho(t,x) + \mathop{\mathrm{div}}_x\bigl[\rho(t,x)v\bigl(\rho(t,x)\bigr)\bigl(V(t,x)-(\partial_xW*\rho)(t,x)\bigr)\bigr] = 0 , $$

where

• $V:[0,\infty)\times\mathbb{R}\to\mathbb{R}$ is the external velocity,
• $W:[0,\infty)\to\mathbb{R}$ is the interaction potential,
• $v:[0,\infty)\to[0,\infty)$ is a non increasing function to model congestion,
• the convolution is performed in space only.

The approach to solve the conservation law is by tracking the motion of $N+1$ particles $X=(x_0,x_1,\dots,x_N)$ such that between each consecutive pair of particles there is a fraction $1/N$ of the total mass.

Two deterministic particle schemes are described: one with integrated interaction $(\mathrm{ODE}_I)$ and one with sampled interaction $(\mathrm{ODE}_S)$:

\begin{aligned}
(\mathrm{ODE}_I): & \left\{\begin{aligned}
x_i'(t) &= v_i(t) \bar U_i(t), \\
\bar U_i(t) &= V\bigl(t,x_i(t)\bigr) - (W*\partial_x\bar\rho)\bigl(t,x_i(t)\bigr)
    = V\bigl(t,x_i(t)\bigr) - \sum_{j=0}^N(\rho_{j+1}(t) - \rho_j(t)) W\bigl(t,x_i(t)-x_j(t)\bigr), \\
v_i(t) &= \begin{cases}
    v\bigl(\rho_i(t)\bigr), & \text{if } \bar U_i(t) < 0, \\
    v\bigl(\rho_{i+1}(t)\bigr), & \text{if } \bar U_i(t) \geq 0, \end{cases}
\end{aligned}\right. \\
(\mathrm{ODE}_S): & \left\{\begin{aligned}
x_i'(t) &= v_i(t) \dot U_i(t), \\
\dot U_i(t) &= V\bigl(t,x_i(t)\bigr) - (\partial_xW*\dot\rho)\bigl(t,x_i(t)\bigr)
    = V\bigl(t,x_i(t)\bigr) - \frac1N\sum_{j=0}^NW'\bigl(t,x_i(t)-x_j(t)\bigr), \\
v_i(t) &= \begin{cases}
    v\bigl(\rho_i(t)\bigr), & \text{if } \dot U_i(t) < 0, \\
    v\bigl(\rho_{i+1}(t)\bigr), & \text{if } \dot U_i(t) \geq 0. \end{cases}
\end{aligned}\right.
\end{aligned}

Here $\rho_i$ denotes the quantity $\frac1{N(x_i-x_{i-1})}$, $\bar\rho$ is the piecewise constant density $\sum_{i=1}^N \rho_i 1_{[x_{i-1},x_i]}$ and $\dot\rho$ is the atomic measure $\frac1N \sum_{i=0}^N \delta_{x_i}$. Notice that $\dot\rho$ is not a probability.

Example

This example requires the ode_solver feature.

Consider the following external velocity, interaction potential and mobility:

\begin{aligned}
V(t,x) &= \mathop{\mathrm{sign}}(7.5-t)
    [2 + 0.5\sin(2x) + 0.3\sin(2\sqrt2x) + 0.2\cos(2\sqrt7x)] , \\
W(t,x) &= (1 + \sin(t)) (x^2 - |x|) , \\
v(t,\rho) &= (1-\rho)_+ .
\end{aligned}

The initial datum $\rho_0 = 1_{[-1.5,1.5]}$ is discretized $N=201$ with equally spaced particles.

use cons_laws::*;
use ode_solvers::{dopri5::Dopri5, SVector};

#[inline]
fn V(t: f64, x: f64) -> f64 {
    let sign = if t < 7.5 { 1.0 } else { -1.0 };
    let magn = 2.0 + 0.5 * (2.0 * x).sin() + 0.3 * (2.0 * 2.0f64.sqrt() * x).sin()
        + 0.2 * (2.0 * 7.0f64.sqrt() * x).cos();
    sign * magn
}
#[inline]
fn Wprime(t: f64, x: f64) -> f64 { (x + x - x.signum()) * (1.0 + t.sin()) }
#[inline]
fn W(t: f64, x: f64) -> f64 { (x * x - x.abs()) * (1.0 + t.sin()) }
#[inline]
fn v(t: f64, rho: f64) -> f64 { 0.0f64.max(1.0 - rho) }

// define the model, using the sampled or integrated interaction
let sampl_inter = SampledInteraction::new(Wprime);
let integ_inter = SampledInteraction::new(W);
let model = SingleSpeciesModel::new(V, sampl_inter, v);

// initial condition
const N: usize = 201;
let x0 = SVector::<f64, N>::from_iterator((0..N).map(|i| 3.0 * i as f64 / (N - 1) as f64 - 1.5));

// configure the solver
let t_start = 0.0;
let t_end = 15.0;
let dt = 0.001;
let rtol = 1e-10;
let atol = 1e-10;
let mut solver = Dopri5::new(model, t_start, t_end, dt, x0, rtol, atol);

// solve the system of ODEs
solver.integrate();
let t = solver.x_out(); // times
let x = solver.y_out(); // particles positions
let n_steps = t.len();

Plotting (some of) the trajectories with plotters produces the following picture:

Structs

IntegratedInteraction

OneMobility

SampledInteraction

SingleSpeciesModel

Velocity

ZeroInteraction

ZeroVelocity

Traits

ExternalVelocity

A time-dependent velocity field that affects the particles.

Interaction

A time-dependent interaction between the particles.

Mobility

A time-dependent mobility.