Documentation in SBP

This commit is contained in:
Magnus Ulimoen 2020-09-19 15:46:02 +02:00
parent e3a15b3844
commit 4eb653da23
4 changed files with 60 additions and 0 deletions

View File

@ -1,8 +1,22 @@
//! Integration of explicit PDEs using different Butcher Tableaus
//!
//! Integration can be performed on all systems that can be represented
//! as using a transform into an [`ndarray::ArrayView`] for both the state
//! and the state difference.
//!
//! The integration functions are memory efficient, and relies
//! on the `k` parameter to hold the system state differences.
//! This parameter is tied to the Butcher Tableau
use super::Float; use super::Float;
use ndarray::{ArrayView, ArrayViewMut}; use ndarray::{ArrayView, ArrayViewMut};
/// The Butcher Tableau, with the state transitions described as
/// [on wikipedia](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Explicit_Runge%E2%80%93Kutta_methods).
pub trait ButcherTableau { pub trait ButcherTableau {
/// This bound should not be overridden
const S: usize = Self::B.len(); const S: usize = Self::B.len();
/// Only the lower triangle will be used (explicit integration)
const A: &'static [&'static [Float]]; const A: &'static [&'static [Float]];
const B: &'static [Float]; const B: &'static [Float];
const C: &'static [Float]; const C: &'static [Float];
@ -137,6 +151,16 @@ impl EmbeddedButcherTableau for BogackiShampine {
} }
#[allow(clippy::too_many_arguments)] #[allow(clippy::too_many_arguments)]
/// Integrates using the [`ButcherTableau`] specified. `rhs` should be the result
/// of the right hand side of $u_t = rhs$
///
/// rhs takes the old state and the current time, and outputs the state difference
/// in the first parameter
///
/// Should be called as
/// ```rust,ignore
/// integrate::<Rk4, _, _, _, _>(...)
/// ```
pub fn integrate<BTableau: ButcherTableau, F, RHS, D>( pub fn integrate<BTableau: ButcherTableau, F, RHS, D>(
mut rhs: RHS, mut rhs: RHS,
prev: &F, prev: &F,
@ -191,6 +215,10 @@ pub fn integrate<BTableau: ButcherTableau, F, RHS, D>(
} }
#[allow(clippy::too_many_arguments)] #[allow(clippy::too_many_arguments)]
/// Integrate using an [`EmbeddedButcherTableau`], else similar to [`integrate`]
///
/// This produces two results, the most accurate result in `fut`, and the less accurate
/// result in `fut2`. This can be used for convergence testing and adaptive timesteps.
pub fn integrate_embedded_rk<BTableau: EmbeddedButcherTableau, F, RHS, D>( pub fn integrate_embedded_rk<BTableau: EmbeddedButcherTableau, F, RHS, D>(
rhs: RHS, rhs: RHS,
prev: &F, prev: &F,
@ -217,6 +245,16 @@ pub fn integrate_embedded_rk<BTableau: EmbeddedButcherTableau, F, RHS, D>(
#[cfg(feature = "rayon")] #[cfg(feature = "rayon")]
#[allow(clippy::too_many_arguments)] #[allow(clippy::too_many_arguments)]
/// Integrates a multigrid problem, much the same as [`integrate`],
/// using a `rayon` threadpool for parallelisation.
///
/// note that `rhs` accepts the full system state, and is responsible
/// for computing the full state difference.
/// `rhs` can be a mutable closure, so buffers can be used
/// and mutated inside the closure.
///
/// This function requires the `rayon` feature, and is not callable in
/// a `wasm` context.
pub fn integrate_multigrid<BTableau: ButcherTableau, F, RHS, D>( pub fn integrate_multigrid<BTableau: ButcherTableau, F, RHS, D>(
mut rhs: RHS, mut rhs: RHS,
prev: &[F], prev: &[F],

View File

@ -1,11 +1,14 @@
#![feature(min_specialization)] #![feature(min_specialization)]
#![feature(core_intrinsics)] #![feature(core_intrinsics)]
/// Type used for floats, configure with the `f32` feature
#[cfg(feature = "f32")] #[cfg(feature = "f32")]
pub type Float = f32; pub type Float = f32;
#[cfg(not(feature = "f32"))] #[cfg(not(feature = "f32"))]
/// Type used for floats, configure with the `f32` feature
pub type Float = f64; pub type Float = f64;
/// Associated constants for [`Float`]
pub mod consts { pub mod consts {
#[cfg(feature = "f32")] #[cfg(feature = "f32")]
pub use std::f32::consts::*; pub use std::f32::consts::*;
@ -13,7 +16,11 @@ pub mod consts {
pub use std::f64::consts::*; pub use std::f64::consts::*;
} }
/// Grid and grid metrics
pub mod grid; pub mod grid;
/// RK operators and methods for implicit integration
pub mod integrate; pub mod integrate;
/// SBP and interpolation operators
pub mod operators; pub mod operators;
/// General utilities
pub mod utils; pub mod utils;

View File

@ -4,10 +4,14 @@
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
/// One-dimensional Summation By Parts operator
pub trait SbpOperator1d: Send + Sync { pub trait SbpOperator1d: Send + Sync {
/// Differentiate on unit grid
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>); fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
/// block component of `H`, with diagonal component 1
fn h(&self) -> &'static [Float]; fn h(&self) -> &'static [Float];
/// Whether the operator acts on a `h2` spaced computational grid
fn is_h2(&self) -> bool { fn is_h2(&self) -> bool {
false false
} }
@ -18,6 +22,7 @@ pub trait SbpOperator1d: Send + Sync {
} }
pub trait SbpOperator1d2: SbpOperator1d { pub trait SbpOperator1d2: SbpOperator1d {
/// Double differentiation on unit grid
fn diff2(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>); fn diff2(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
/// Result of H^-1 * d1, without the 1/h scaling /// Result of H^-1 * d1, without the 1/h scaling
fn d1(&self) -> &[Float]; fn d1(&self) -> &[Float];
@ -103,6 +108,7 @@ impl<SBP: SbpOperator1d + Copy> SbpOperator2d for SBP {
} }
pub trait UpwindOperator1d: SbpOperator1d + Send + Sync { pub trait UpwindOperator1d: SbpOperator1d + Send + Sync {
/// Dissipation operator
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>); fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
fn as_sbp(&self) -> &dyn SbpOperator1d; fn as_sbp(&self) -> &dyn SbpOperator1d;
@ -162,7 +168,9 @@ impl<UO: UpwindOperator1d + Copy> UpwindOperator2d for UO {
} }
pub trait InterpolationOperator: Send + Sync { pub trait InterpolationOperator: Send + Sync {
/// Interpolation from a grid with twice resolution
fn fine2coarse(&self, fine: ArrayView1<Float>, coarse: ArrayViewMut1<Float>); fn fine2coarse(&self, fine: ArrayView1<Float>, coarse: ArrayViewMut1<Float>);
/// Interpolation from a grid with half resolution
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>); fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>);
} }

View File

@ -13,6 +13,7 @@ use serde::{Deserialize, Serialize};
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
#[derive(Copy, Clone, Debug, Default)] #[derive(Copy, Clone, Debug, Default)]
/// struct to hold output for four directions
pub struct Direction<T> { pub struct Direction<T> {
pub north: T, pub north: T,
pub south: T, pub south: T,
@ -45,6 +46,7 @@ impl<T> Direction<T> {
east: f(self.east), east: f(self.east),
} }
} }
/// Combines two [`Direction`] into one
pub fn zip<U>(self, other: Direction<U>) -> Direction<(T, U)> { pub fn zip<U>(self, other: Direction<U>) -> Direction<(T, U)> {
Direction { Direction {
north: (self.north, other.north), north: (self.north, other.north),
@ -56,6 +58,8 @@ impl<T> Direction<T> {
} }
impl<T> Direction<Option<T>> { impl<T> Direction<Option<T>> {
/// Partially unwrap individual components in `self` or replace with
/// `other`
pub fn unwrap_or(self, other: Direction<T>) -> Direction<T> { pub fn unwrap_or(self, other: Direction<T>) -> Direction<T> {
Direction { Direction {
north: self.north.unwrap_or(other.north), north: self.north.unwrap_or(other.north),
@ -66,6 +70,7 @@ impl<T> Direction<Option<T>> {
} }
} }
/// Methods for borrowing
impl<T> Direction<T> { impl<T> Direction<T> {
pub fn north(&self) -> &T { pub fn north(&self) -> &T {
&self.north &self.north
@ -93,6 +98,8 @@ impl<T> Direction<T> {
} }
} }
/// Linearly spaced parameters, apart from the boundaries which
/// only have a distance of `h/2` from the boundary
pub fn h2linspace(start: Float, end: Float, n: usize) -> ndarray::Array1<Float> { pub fn h2linspace(start: Float, end: Float, n: usize) -> ndarray::Array1<Float> {
let h = (end - start) / (n - 2) as Float; let h = (end - start) / (n - 2) as Float;
ndarray::Array1::from_shape_fn(n, |i| match i { ndarray::Array1::from_shape_fn(n, |i| match i {