mulimoen/SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

change sbp operators for dimensionality

Browse Source
This commit is contained in:
Magnus Ulimoen
2020-04-15 00:12:54 +02:00
parent 85f3c46430
commit 1667eaaca0
12 changed files with 230 additions and 162 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
sbp/benches/euler.rs
View File

@@ -1,14 +1,14 @@
use criterion::{black_box, criterion_group, criterion_main, Criterion}; use criterion::{black_box, criterion_group, criterion_main, Criterion};
use sbp::euler::System; use sbp::euler::System;
use sbp::operators::{SbpOperator, Upwind4, UpwindOperator, SBP4}; use sbp::operators::{SbpOperator2d, Upwind4, UpwindOperator2d, SBP4};
fn advance_system<SBP: SbpOperator>(universe: &mut System<SBP>, n: usize) { fn advance_system<SBP: SbpOperator2d>(universe: &mut System<SBP>, n: usize) {
for _ in 0..n { for _ in 0..n {
universe.advance(1.0 / 40.0 * 0.2); universe.advance(1.0 / 40.0 * 0.2);
} }
} }
fn advance_system_upwind<UO: UpwindOperator>(universe: &mut System<UO>, n: usize) { fn advance_system_upwind<UO: UpwindOperator2d>(universe: &mut System<UO>, n: usize) {
for _ in 0..n { for _ in 0..n {
universe.advance_upwind(1.0 / 40.0 * 0.2); universe.advance_upwind(1.0 / 40.0 * 0.2);
} }

6
sbp/benches/maxwell.rs
View File

@@ -1,15 +1,15 @@
use criterion::{black_box, criterion_group, criterion_main, Criterion}; use criterion::{black_box, criterion_group, criterion_main, Criterion};
use sbp::maxwell::System; use sbp::maxwell::System;
use sbp::operators::{SbpOperator, Upwind4, UpwindOperator, SBP4}; use sbp::operators::{SbpOperator2d, Upwind4, UpwindOperator2d, SBP4};
use sbp::Float; use sbp::Float;
fn advance_system<SBP: SbpOperator>(universe: &mut System<SBP>, n: usize) { fn advance_system<SBP: SbpOperator2d>(universe: &mut System<SBP>, n: usize) {
for _ in 0..n { for _ in 0..n {
universe.advance(0.01); universe.advance(0.01);
} }
} }
fn advance_system_upwind<UO: UpwindOperator>(universe: &mut System<UO>, n: usize) { fn advance_system_upwind<UO: UpwindOperator2d>(universe: &mut System<UO>, n: usize) {
for _ in 0..n { for _ in 0..n {
universe.advance_upwind(0.01); universe.advance_upwind(0.01);
} }

128
sbp/src/euler.rs
View File

@@ -1,6 +1,6 @@
use super::grid::{Grid, Metrics}; use super::grid::{Grid, Metrics};
use super::integrate; use super::integrate;
use super::operators::{InterpolationOperator, SbpOperator, UpwindOperator}; use super::operators::{InterpolationOperator, SbpOperator2d, UpwindOperator2d};
use super::Float; use super::Float;
use ndarray::azip; use ndarray::azip;
use ndarray::prelude::*; use ndarray::prelude::*;
@@ -10,7 +10,7 @@ pub const GAMMA: Float = 1.4;
// A collection of buffers that allows one to efficiently // A collection of buffers that allows one to efficiently
// move to the next state // move to the next state
#[derive(Debug)] #[derive(Debug)]
pub struct System<SBP: SbpOperator> { pub struct System<SBP: SbpOperator2d> {
sys: (Field, Field), sys: (Field, Field),
k: [Field; 4], k: [Field; 4],
wb: WorkBuffers, wb: WorkBuffers,
@@ -18,12 +18,12 @@ pub struct System<SBP: SbpOperator> {
op: SBP, op: SBP,
} }
impl<SBP: SbpOperator> System<SBP> { impl<SBP: SbpOperator2d> System<SBP> {
pub fn new(x: ndarray::Array2<Float>, y: ndarray::Array2<Float>, op: SBP) -> Self { pub fn new(x: ndarray::Array2<Float>, y: ndarray::Array2<Float>, op: SBP) -> Self {
let grid = Grid::new(x, y).expect( let grid = Grid::new(x, y).expect(
"Could not create grid. Different number of elements compared to width*height?", "Could not create grid. Different number of elements compared to width*height?",
); );
let metrics = grid.metrics(op, op).unwrap(); let metrics = grid.metrics(op).unwrap();
let nx = grid.nx(); let nx = grid.nx();
let ny = grid.ny(); let ny = grid.ny();
Self { Self {
@@ -51,7 +51,7 @@ impl<SBP: SbpOperator> System<SBP> {
let rhs_trad = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| { let rhs_trad = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| {
let (grid, metrics) = gm; let (grid, metrics) = gm;
let boundaries = boundary_extractor(y, grid, &bc); let boundaries = boundary_extractor(y, grid, &bc);
RHS_trad((op, op), k, y, metrics, &boundaries, wb) RHS_trad(op, k, y, metrics, &boundaries, wb)
}; };
integrate::integrate::<integrate::Rk4, _, _, _, _>( integrate::integrate::<integrate::Rk4, _, _, _, _>(
rhs_trad, rhs_trad,
@@ -100,7 +100,7 @@ impl<SBP: SbpOperator> System<SBP> {
} }
} }
impl<UO: UpwindOperator> System<UO> { impl<UO: UpwindOperator2d> System<UO> {
pub fn advance_upwind(&mut self, dt: Float) { pub fn advance_upwind(&mut self, dt: Float) {
let bc = BoundaryCharacteristics { let bc = BoundaryCharacteristics {
north: BoundaryCharacteristic::This, north: BoundaryCharacteristic::This,
@@ -112,7 +112,7 @@ impl<UO: UpwindOperator> System<UO> {
let rhs_upwind = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| { let rhs_upwind = |k: &mut Field, y: &Field, _time: Float, gm: &(_, _), wb: &mut _| {
let (grid, metrics) = gm; let (grid, metrics) = gm;
let boundaries = boundary_extractor(y, grid, &bc); let boundaries = boundary_extractor(y, grid, &bc);
RHS_upwind((op, op), k, y, metrics, &boundaries, wb) RHS_upwind(op, k, y, metrics, &boundaries, wb)
}; };
integrate::integrate::<integrate::Rk4, _, _, _, _>( integrate::integrate::<integrate::Rk4, _, _, _, _>(
rhs_upwind, rhs_upwind,
@@ -272,11 +272,7 @@ impl Field {
impl Field { impl Field {
/// sqrt((self-other)^T*H*(self-other)) /// sqrt((self-other)^T*H*(self-other))
pub fn h2_err<SBPeta: SbpOperator, SBPxi: SbpOperator>( pub fn h2_err<SBP: SbpOperator2d>(&self, op: SBP, other: &Self) -> Float {
&self,
(opeta, opxi): (SBPeta, SBPxi),
other: &Self,
) -> Float {
assert_eq!(self.nx(), other.nx()); assert_eq!(self.nx(), other.nx());
assert_eq!(self.ny(), other.ny()); assert_eq!(self.ny(), other.ny());
@@ -297,9 +293,9 @@ impl Field {
}; };
let hxiterator = itermaker( let hxiterator = itermaker(
opxi.h(), op.hxi(),
self.nx(), self.nx(),
if opxi.is_h2() { if op.is_h2xi() {
1.0 / (self.nx() - 2) as Float 1.0 / (self.nx() - 2) as Float
} else { } else {
1.0 / (self.nx() - 1) as Float 1.0 / (self.nx() - 1) as Float
@@ -310,9 +306,9 @@ impl Field {
let hxiterator = hxiterator.cycle().take(self.nx() * self.ny()); let hxiterator = hxiterator.cycle().take(self.nx() * self.ny());
let hyiterator = itermaker( let hyiterator = itermaker(
opeta.h(), op.heta(),
self.ny(), self.ny(),
1.0 / if opeta.is_h2() { 1.0 / if op.is_h2eta() {
(self.ny() - 2) as Float (self.ny() - 2) as Float
} else { } else {
(self.ny() - 1) as Float (self.ny() - 1) as Float
@@ -410,8 +406,8 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo
} }
#[allow(non_snake_case)] #[allow(non_snake_case)]
pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>( pub fn RHS_trad<SBP: SbpOperator2d>(
(opeta, opxi): (SBPeta, SBPxi), op: SBP,
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics, metrics: &Metrics,
@@ -424,15 +420,15 @@ pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>(
let dE = &mut tmp.2; let dE = &mut tmp.2;
let dF = &mut tmp.3; let dF = &mut tmp.3;
opxi.diffxi(ehat.rho(), dE.rho_mut()); op.diffxi(ehat.rho(), dE.rho_mut());
opxi.diffxi(ehat.rhou(), dE.rhou_mut()); op.diffxi(ehat.rhou(), dE.rhou_mut());
opxi.diffxi(ehat.rhov(), dE.rhov_mut()); op.diffxi(ehat.rhov(), dE.rhov_mut());
opxi.diffxi(ehat.e(), dE.e_mut()); op.diffxi(ehat.e(), dE.e_mut());
opeta.diffeta(fhat.rho(), dF.rho_mut()); op.diffeta(fhat.rho(), dF.rho_mut());
opeta.diffeta(fhat.rhou(), dF.rhou_mut()); op.diffeta(fhat.rhou(), dF.rhou_mut());
opeta.diffeta(fhat.rhov(), dF.rhov_mut()); op.diffeta(fhat.rhov(), dF.rhov_mut());
opeta.diffeta(fhat.e(), dF.e_mut()); op.diffeta(fhat.e(), dF.e_mut());
azip!((out in &mut k.0, azip!((out in &mut k.0,
eflux in &dE.0, eflux in &dE.0,
@@ -441,12 +437,12 @@ pub fn RHS_trad<SBPeta: SbpOperator, SBPxi: SbpOperator>(
*out = (-eflux - fflux)/detj *out = (-eflux - fflux)/detj
}); });
SAT_characteristics((opeta, opxi), k, y, metrics, boundaries); SAT_characteristics(op, k, y, metrics, boundaries);
} }
#[allow(non_snake_case)] #[allow(non_snake_case)]
pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>( pub fn RHS_upwind<UO: UpwindOperator2d>(
(opeta, opxi): (UOeta, UOxi), op: UO,
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics, metrics: &Metrics,
@@ -459,25 +455,19 @@ pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>(
let dE = &mut tmp.2; let dE = &mut tmp.2;
let dF = &mut tmp.3; let dF = &mut tmp.3;
opxi.diffxi(ehat.rho(), dE.rho_mut()); op.diffxi(ehat.rho(), dE.rho_mut());
opxi.diffxi(ehat.rhou(), dE.rhou_mut()); op.diffxi(ehat.rhou(), dE.rhou_mut());
opxi.diffxi(ehat.rhov(), dE.rhov_mut()); op.diffxi(ehat.rhov(), dE.rhov_mut());
opxi.diffxi(ehat.e(), dE.e_mut()); op.diffxi(ehat.e(), dE.e_mut());
opeta.diffeta(fhat.rho(), dF.rho_mut()); op.diffeta(fhat.rho(), dF.rho_mut());
opeta.diffeta(fhat.rhou(), dF.rhou_mut()); op.diffeta(fhat.rhou(), dF.rhou_mut());
opeta.diffeta(fhat.rhov(), dF.rhov_mut()); op.diffeta(fhat.rhov(), dF.rhov_mut());
opeta.diffeta(fhat.e(), dF.e_mut()); op.diffeta(fhat.e(), dF.e_mut());
let ad_xi = &mut tmp.4; let ad_xi = &mut tmp.4;
let ad_eta = &mut tmp.5; let ad_eta = &mut tmp.5;
upwind_dissipation( upwind_dissipation(op, (ad_xi, ad_eta), y, metrics, (&mut tmp.0, &mut tmp.1));
(opeta, opxi),
(ad_xi, ad_eta),
y,
metrics,
(&mut tmp.0, &mut tmp.1),
);
azip!((out in &mut k.0, azip!((out in &mut k.0,
eflux in &dE.0, eflux in &dE.0,
@@ -488,12 +478,12 @@ pub fn RHS_upwind<UOeta: UpwindOperator, UOxi: UpwindOperator>(
*out = (-eflux - fflux + ad_xi + ad_eta)/detj *out = (-eflux - fflux + ad_xi + ad_eta)/detj
}); });
SAT_characteristics((opeta, opxi), k, y, metrics, boundaries); SAT_characteristics(op, k, y, metrics, boundaries);
} }
#[allow(clippy::many_single_char_names)] #[allow(clippy::many_single_char_names)]
fn upwind_dissipation<UOeta: UpwindOperator, UOxi: UpwindOperator>( fn upwind_dissipation<UO: UpwindOperator2d>(
(opeta, opxi): (UOeta, UOxi), op: UO,
k: (&mut Field, &mut Field), k: (&mut Field, &mut Field),
y: &Field, y: &Field,
metrics: &Metrics, metrics: &Metrics,
@@ -549,15 +539,15 @@ fn upwind_dissipation<UOeta: UpwindOperator, UOxi: UpwindOperator>(
tmp1[3] = alpha_v * e * detj; tmp1[3] = alpha_v * e * detj;
} }
opxi.dissxi(tmp.0.rho(), k.0.rho_mut()); op.dissxi(tmp.0.rho(), k.0.rho_mut());
opxi.dissxi(tmp.0.rhou(), k.0.rhou_mut()); op.dissxi(tmp.0.rhou(), k.0.rhou_mut());
opxi.dissxi(tmp.0.rhov(), k.0.rhov_mut()); op.dissxi(tmp.0.rhov(), k.0.rhov_mut());
opxi.dissxi(tmp.0.e(), k.0.e_mut()); op.dissxi(tmp.0.e(), k.0.e_mut());
opeta.disseta(tmp.1.rho(), k.1.rho_mut()); op.disseta(tmp.1.rho(), k.1.rho_mut());
opeta.disseta(tmp.1.rhou(), k.1.rhou_mut()); op.disseta(tmp.1.rhou(), k.1.rhou_mut());
opeta.disseta(tmp.1.rhov(), k.1.rhov_mut()); op.disseta(tmp.1.rhov(), k.1.rhov_mut());
opeta.disseta(tmp.1.e(), k.1.e_mut()); op.disseta(tmp.1.e(), k.1.e_mut());
} }
fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics) { fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics) {
@@ -853,8 +843,8 @@ fn vortexify(
#[allow(non_snake_case)] #[allow(non_snake_case)]
/// Boundary conditions (SAT) /// Boundary conditions (SAT)
fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>( fn SAT_characteristics<SBP: SbpOperator2d>(
(opeta, opxi): (SBPeta, SBPxi), op: SBP,
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics, metrics: &Metrics,
@@ -862,10 +852,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
) { ) {
// North boundary // North boundary
{ {
let hi = if opeta.is_h2() { let hi = if op.is_h2eta() {
(k.ny() - 2) as Float / opeta.h()[0] (k.ny() - 2) as Float / op.heta()[0]
} else { } else {
(k.ny() - 1) as Float / opeta.h()[0] (k.ny() - 1) as Float / op.heta()[0]
}; };
let sign = -1.0; let sign = -1.0;
let tau = 1.0; let tau = 1.0;
@@ -884,10 +874,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
} }
// South boundary // South boundary
{ {
let hi = if opeta.is_h2() { let hi = if op.is_h2eta() {
(k.ny() - 2) as Float / opeta.h()[0] (k.ny() - 2) as Float / op.heta()[0]
} else { } else {
(k.ny() - 1) as Float / opeta.h()[0] (k.ny() - 1) as Float / op.heta()[0]
}; };
let sign = 1.0; let sign = 1.0;
let tau = -1.0; let tau = -1.0;
@@ -906,10 +896,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
} }
// West Boundary // West Boundary
{ {
let hi = if opxi.is_h2() { let hi = if op.is_h2xi() {
(k.nx() - 2) as Float / opxi.h()[0] (k.nx() - 2) as Float / op.hxi()[0]
} else { } else {
(k.nx() - 1) as Float / opxi.h()[0] (k.nx() - 1) as Float / op.hxi()[0]
}; };
let sign = 1.0; let sign = 1.0;
let tau = -1.0; let tau = -1.0;
@@ -928,10 +918,10 @@ fn SAT_characteristics<SBPeta: SbpOperator, SBPxi: SbpOperator>(
} }
// East Boundary // East Boundary
{ {
let hi = if opxi.is_h2() { let hi = if op.is_h2xi() {
(k.nx() - 2) as Float / opxi.h()[0] (k.nx() - 2) as Float / op.hxi()[0]
} else { } else {
(k.nx() - 1) as Float / opxi.h()[0] (k.nx() - 1) as Float / op.hxi()[0]
}; };
let sign = -1.0; let sign = -1.0;
let tau = 1.0; let tau = 1.0;

20
sbp/src/grid.rs
View File

@@ -36,12 +36,11 @@ impl Grid {
self.y.view() self.y.view()
} }
pub fn metrics<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>( pub fn metrics<SBP: super::operators::SbpOperator2d>(
&self, &self,
opeta: SBPeta, op: SBP,
opxi: SBPxi,
) -> Result<Metrics, ndarray::ShapeError> { ) -> Result<Metrics, ndarray::ShapeError> {
Metrics::new(self, opeta, opxi) Metrics::new(self, op)
} }
pub fn north(&self) -> (ndarray::ArrayView1<Float>, ndarray::ArrayView1<Float>) { pub fn north(&self) -> (ndarray::ArrayView1<Float>, ndarray::ArrayView1<Float>) {
@@ -71,10 +70,9 @@ impl Grid {
} }
impl Metrics { impl Metrics {
fn new<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>( fn new<SBP: super::operators::SbpOperator2d>(
grid: &Grid, grid: &Grid,
opeta: SBPeta, op: SBP,
opxi: SBPxi,
) -> Result<Self, ndarray::ShapeError> { ) -> Result<Self, ndarray::ShapeError> {
let ny = grid.ny(); let ny = grid.ny();
let nx = grid.nx(); let nx = grid.nx();
@@ -82,13 +80,13 @@ impl Metrics {
let y = &grid.y; let y = &grid.y;
let mut dx_dxi = Array2::zeros((ny, nx)); let mut dx_dxi = Array2::zeros((ny, nx));
opxi.diffxi(x.view(), dx_dxi.view_mut()); op.diffxi(x.view(), dx_dxi.view_mut());
let mut dx_deta = Array2::zeros((ny, nx)); let mut dx_deta = Array2::zeros((ny, nx));
opeta.diffeta(x.view(), dx_deta.view_mut()); op.diffeta(x.view(), dx_deta.view_mut());
let mut dy_dxi = Array2::zeros((ny, nx)); let mut dy_dxi = Array2::zeros((ny, nx));
opxi.diffxi(y.view(), dy_dxi.view_mut()); op.diffxi(y.view(), dy_dxi.view_mut());
let mut dy_deta = Array2::zeros((ny, nx)); let mut dy_deta = Array2::zeros((ny, nx));
opeta.diffeta(y.view(), dy_deta.view_mut()); op.diffeta(y.view(), dy_deta.view_mut());
let mut detj = Array2::zeros((ny, nx)); let mut detj = Array2::zeros((ny, nx));
ndarray::azip!((detj in &mut detj, ndarray::azip!((detj in &mut detj,

44
sbp/src/maxwell.rs
View File

@@ -1,6 +1,6 @@
use super::grid::{Grid, Metrics}; use super::grid::{Grid, Metrics};
use super::integrate; use super::integrate;
use super::operators::{SbpOperator, UpwindOperator}; use super::operators::{SbpOperator2d, UpwindOperator2d};
use crate::Float; use crate::Float;
use ndarray::azip; use ndarray::azip;
use ndarray::prelude::*; use ndarray::prelude::*;
@@ -74,7 +74,7 @@ impl Field {
} }
#[derive(Debug, Clone)] #[derive(Debug, Clone)]
pub struct System<SBP: SbpOperator> { pub struct System<SBP: SbpOperator2d> {
sys: (Field, Field), sys: (Field, Field),
wb: WorkBuffers, wb: WorkBuffers,
grid: Grid, grid: Grid,
@@ -82,14 +82,14 @@ pub struct System<SBP: SbpOperator> {
op: SBP, op: SBP,
} }
impl<SBP: SbpOperator> System<SBP> { impl<SBP: SbpOperator2d> System<SBP> {
pub fn new(x: Array2<Float>, y: Array2<Float>, op: SBP) -> Self { pub fn new(x: Array2<Float>, y: Array2<Float>, op: SBP) -> Self {
assert_eq!(x.shape(), y.shape()); assert_eq!(x.shape(), y.shape());
let ny = x.shape()[0]; let ny = x.shape()[0];
let nx = x.shape()[1]; let nx = x.shape()[1];
let grid = Grid::new(x, y).unwrap(); let grid = Grid::new(x, y).unwrap();
let metrics = grid.metrics(op, op).unwrap(); let metrics = grid.metrics(op).unwrap();
Self { Self {
op, op,
@@ -146,7 +146,7 @@ impl<SBP: SbpOperator> System<SBP> {
} }
} }
impl<UO: UpwindOperator> System<UO> { impl<UO: UpwindOperator2d> System<UO> {
/// Using artificial dissipation with the upwind operator /// Using artificial dissipation with the upwind operator
pub fn advance_upwind(&mut self, dt: Float) { pub fn advance_upwind(&mut self, dt: Float) {
let op = self.op; let op = self.op;
@@ -205,7 +205,7 @@ fn gaussian(x: Float, x0: Float, y: Float, y0: Float) -> Float {
/// where J is the grid determinant /// where J is the grid determinant
/// ///
/// This is used both in fluxes and SAT terms /// This is used both in fluxes and SAT terms
fn RHS<SBP: SbpOperator>( fn RHS<SBP: SbpOperator2d>(
op: SBP, op: SBP,
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
@@ -230,7 +230,7 @@ fn RHS<SBP: SbpOperator>(
} }
#[allow(non_snake_case)] #[allow(non_snake_case)]
fn RHS_upwind<UO: UpwindOperator>( fn RHS_upwind<UO: UpwindOperator2d>(
op: UO, op: UO,
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
@@ -255,7 +255,7 @@ fn RHS_upwind<UO: UpwindOperator>(
}); });
} }
fn fluxes<SBP: SbpOperator>( fn fluxes<SBP: super::operators::SbpOperator2d>(
op: SBP, op: SBP,
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
@@ -330,7 +330,7 @@ fn fluxes<SBP: SbpOperator>(
} }
} }
fn dissipation<UO: UpwindOperator>( fn dissipation<UO: UpwindOperator2d>(
op: UO, op: UO,
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
@@ -432,7 +432,7 @@ pub struct BoundaryTerms {
#[allow(non_snake_case)] #[allow(non_snake_case)]
/// Boundary conditions (SAT) /// Boundary conditions (SAT)
fn SAT_characteristics<SBP: SbpOperator>( fn SAT_characteristics<SBP: SbpOperator2d>(
op: SBP, op: SBP,
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
@@ -464,10 +464,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
Boundary::This => y.slice(s![.., .., 0]), Boundary::This => y.slice(s![.., .., 0]),
}; };
// East boundary // East boundary
let hinv = if op.is_h2() { let hinv = if op.is_h2xi() {
(nx - 2) as Float / op.h()[0] (nx - 2) as Float / op.hxi()[0]
} else { } else {
(nx - 1) as Float / op.h()[0] (nx - 1) as Float / op.hxi()[0]
}; };
for ((((mut k, v), g), &kx), &ky) in k for ((((mut k, v), g), &kx), &ky) in k
.slice_mut(s![.., .., nx - 1]) .slice_mut(s![.., .., nx - 1])
@@ -502,10 +502,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
let g = match boundaries.east { let g = match boundaries.east {
Boundary::This => y.slice(s![.., .., nx - 1]), Boundary::This => y.slice(s![.., .., nx - 1]),
}; };
let hinv = if op.is_h2() { let hinv = if op.is_h2xi() {
(nx - 2) as Float / op.h()[0] (nx - 2) as Float / op.hxi()[0]
} else { } else {
(nx - 1) as Float / op.h()[0] (nx - 1) as Float / op.hxi()[0]
}; };
for ((((mut k, v), g), &kx), &ky) in k for ((((mut k, v), g), &kx), &ky) in k
.slice_mut(s![.., .., 0]) .slice_mut(s![.., .., 0])
@@ -545,10 +545,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
let g = match boundaries.north { let g = match boundaries.north {
Boundary::This => y.slice(s![.., 0, ..]), Boundary::This => y.slice(s![.., 0, ..]),
}; };
let hinv = if op.is_h2() { let hinv = if op.is_h2eta() {
(ny - 2) as Float / op.h()[0] (ny - 2) as Float / op.heta()[0]
} else { } else {
(ny - 1) as Float / op.h()[0] (ny - 1) as Float / op.heta()[0]
}; };
for ((((mut k, v), g), &kx), &ky) in k for ((((mut k, v), g), &kx), &ky) in k
.slice_mut(s![.., ny - 1, ..]) .slice_mut(s![.., ny - 1, ..])
@@ -582,10 +582,10 @@ fn SAT_characteristics<SBP: SbpOperator>(
let g = match boundaries.south { let g = match boundaries.south {
Boundary::This => y.slice(s![.., ny - 1, ..]), Boundary::This => y.slice(s![.., ny - 1, ..]),
}; };
let hinv = if op.is_h2() { let hinv = if op.is_h2eta() {
(ny - 2) as Float / op.h()[0] (ny - 2) as Float / op.heta()[0]
} else { } else {
(ny - 1) as Float / op.h()[0] (ny - 1) as Float / op.heta()[0]
}; };
for ((((mut k, v), g), &kx), &ky) in k for ((((mut k, v), g), &kx), &ky) in k
.slice_mut(s![.., 0, ..]) .slice_mut(s![.., 0, ..])

104
sbp/src/operators.rs
View File

@@ -2,32 +2,104 @@
#![allow(clippy::unreadable_literal)] #![allow(clippy::unreadable_literal)]
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
pub trait SbpOperator: Send + Sync + Copy { pub trait SbpOperator1d: Copy + Clone + core::fmt::Debug {
fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>); fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
self.diff1d(r0, r1);
}
}
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.diffxi(prev.reversed_axes(), fut.reversed_axes())
}
fn h(&self) -> &'static [Float]; fn h(&self) -> &'static [Float];
fn is_h2(&self) -> bool { fn is_h2(&self) -> bool {
false false
} }
} }
pub trait UpwindOperator: SbpOperator { pub trait SbpOperator2d: Copy + Clone {
fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>); fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
fn hxi(&self) -> &'static [Float];
fn heta(&self) -> &'static [Float];
fn is_h2xi(&self) -> bool;
fn is_h2eta(&self) -> bool;
}
impl<SBPeta: SbpOperator1d, SBPxi: SbpOperator1d> SbpOperator2d for (SBPeta, SBPxi) {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
self.1.diff(r0, r1)
}
}
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
let ba = (self.1, self.0);
ba.diffxi(prev.reversed_axes(), fut.reversed_axes())
}
fn hxi(&self) -> &'static [Float] {
self.1.h()
}
fn heta(&self) -> &'static [Float] {
self.0.h()
}
fn is_h2xi(&self) -> bool {
self.1.is_h2()
}
fn is_h2eta(&self) -> bool {
self.0.is_h2()
}
}
impl<SBP: SbpOperator1d> SbpOperator2d for SBP {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
self.diff(r0, r1)
}
}
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.diffxi(prev.reversed_axes(), fut.reversed_axes())
}
fn hxi(&self) -> &'static [Float] {
self.h()
}
fn heta(&self) -> &'static [Float] {
self.h()
}
fn is_h2xi(&self) -> bool {
self.is_h2()
}
fn is_h2eta(&self) -> bool {
self.is_h2()
}
}
pub trait UpwindOperator1d: SbpOperator1d + Copy + Clone {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
}
pub trait UpwindOperator2d: SbpOperator2d + Copy + Clone {
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
}
impl<UOeta: UpwindOperator1d, UOxi: UpwindOperator1d> UpwindOperator2d for (UOeta, UOxi) {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) { for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
self.diss1d(r0, r1); self.1.diss(r0, r1);
}
}
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
let ba = (self.1, self.0);
ba.dissxi(prev.reversed_axes(), fut.reversed_axes())
}
}
impl<UO: UpwindOperator1d> UpwindOperator2d for UO {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
self.diss(r0, r1);
} }
} }
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) { fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
@@ -140,7 +212,7 @@ pub(crate) mod testing {
dfdy: FY, dfdy: FY,
eps: Float, eps: Float,
) where ) where
SBP: SbpOperator, SBP: SbpOperator2d,
F: Fn(Float, Float) -> Float, F: Fn(Float, Float) -> Float,
FX: Fn(Float, Float) -> Float, FX: Fn(Float, Float) -> Float,
FY: Fn(Float, Float) -> Float, FY: Fn(Float, Float) -> Float,

6
sbp/src/operators/traditional4.rs
View File

@@ -1,4 +1,4 @@
use super::SbpOperator; use super::*;
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayViewMut1}; use ndarray::{ArrayView1, ArrayViewMut1};
@@ -23,8 +23,8 @@ impl SBP4 {
]; ];
} }
impl SbpOperator for SBP4 { impl SbpOperator1d for SBP4 {
fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::BLOCK).view(), ndarray::arr2(Self::BLOCK).view(),
ndarray::arr1(Self::DIAG).view(), ndarray::arr1(Self::DIAG).view(),

6
sbp/src/operators/traditional8.rs
View File

@@ -1,4 +1,4 @@
use super::SbpOperator; use super::*;
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayViewMut1}; use ndarray::{ArrayView1, ArrayViewMut1};
@@ -27,8 +27,8 @@ impl SBP8 {
]; ];
} }
impl SbpOperator for SBP8 { impl SbpOperator1d for SBP8 {
fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::BLOCK).view(), ndarray::arr2(Self::BLOCK).view(),
ndarray::arr1(Self::DIAG).view(), ndarray::arr1(Self::DIAG).view(),

30
sbp/src/operators/upwind4.rs
View File

@@ -1,4 +1,4 @@
use super::{SbpOperator, UpwindOperator}; use super::*;
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Axis}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Axis};
@@ -275,8 +275,8 @@ impl Upwind4 {
]; ];
} }
impl SbpOperator for Upwind4 { impl SbpOperator1d for Upwind4 {
fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::BLOCK).view(), ndarray::arr2(Self::BLOCK).view(),
ndarray::arr1(Self::DIAG).view(), ndarray::arr1(Self::DIAG).view(),
@@ -286,6 +286,12 @@ impl SbpOperator for Upwind4 {
fut, fut,
) )
} }
fn h(&self) -> &'static [Float] {
Self::HBLOCK
}
}
/*
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Self::BLOCK.len()); assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
@@ -307,10 +313,8 @@ impl SbpOperator for Upwind4 {
} }
} }
fn h(&self) -> &'static [Float] {
Self::HBLOCK
}
} }
*/
#[test] #[test]
fn upwind4_test() { fn upwind4_test() {
@@ -326,7 +330,7 @@ fn upwind4_test() {
target[i] = 1.0; target[i] = 1.0;
} }
res.fill(0.0); res.fill(0.0);
Upwind4.diff1d(source.view(), res.view_mut()); Upwind4.diff(source.view(), res.view_mut());
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4);
{ {
let source = source.to_owned().insert_axis(ndarray::Axis(0)); let source = source.to_owned().insert_axis(ndarray::Axis(0));
@@ -352,7 +356,7 @@ fn upwind4_test() {
target[i] = 2.0 * x; target[i] = 2.0 * x;
} }
res.fill(0.0); res.fill(0.0);
Upwind4.diff1d(source.view(), res.view_mut()); Upwind4.diff(source.view(), res.view_mut());
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4);
{ {
let source = source.to_owned().insert_axis(ndarray::Axis(0)); let source = source.to_owned().insert_axis(ndarray::Axis(0));
@@ -378,7 +382,7 @@ fn upwind4_test() {
target[i] = 3.0 * x * x; target[i] = 3.0 * x * x;
} }
res.fill(0.0); res.fill(0.0);
Upwind4.diff1d(source.view(), res.view_mut()); Upwind4.diff(source.view(), res.view_mut());
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2); approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
{ {
@@ -400,8 +404,8 @@ fn upwind4_test() {
} }
} }
impl UpwindOperator for Upwind4 { impl UpwindOperator1d for Upwind4 {
fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr2(Self::DISS_BLOCK).view(),
ndarray::arr1(Self::DISS_DIAG).view(), ndarray::arr1(Self::DISS_DIAG).view(),
@@ -411,6 +415,9 @@ impl UpwindOperator for Upwind4 {
fut, fut,
) )
} }
}
/*
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Self::BLOCK.len()); assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
@@ -432,6 +439,7 @@ impl UpwindOperator for Upwind4 {
} }
} }
} }
*/
#[test] #[test]
fn upwind4_test2() { fn upwind4_test2() {

16
sbp/src/operators/upwind4h2.rs
View File

@@ -1,4 +1,4 @@
use super::{SbpOperator, UpwindOperator}; use super::*;
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayViewMut1}; use ndarray::{ArrayView1, ArrayViewMut1};
@@ -36,8 +36,8 @@ impl Upwind4h2 {
]; ];
} }
impl SbpOperator for Upwind4h2 { impl SbpOperator1d for Upwind4h2 {
fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::BLOCK).view(), ndarray::arr2(Self::BLOCK).view(),
ndarray::arr1(Self::DIAG).view(), ndarray::arr1(Self::DIAG).view(),
@@ -64,25 +64,25 @@ fn upwind4h2_test() {
let mut res = ndarray::Array1::zeros(nx); let mut res = ndarray::Array1::zeros(nx);
Upwind4h2.diff1d(x.view(), res.view_mut()); Upwind4h2.diff(x.view(), res.view_mut());
let ans = &x * 0.0 + 1.0; let ans = &x * 0.0 + 1.0;
approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4); approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
res.fill(0.0); res.fill(0.0);
let y = &x * &x / 2.0; let y = &x * &x / 2.0;
Upwind4h2.diff1d(y.view(), res.view_mut()); Upwind4h2.diff(y.view(), res.view_mut());
let ans = &x; let ans = &x;
approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4); approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
res.fill(0.0); res.fill(0.0);
let y = &x * &x * &x / 3.0; let y = &x * &x * &x / 3.0;
Upwind4h2.diff1d(y.view(), res.view_mut()); Upwind4h2.diff(y.view(), res.view_mut());
let ans = &x * &x; let ans = &x * &x;
approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2); approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2);
} }
impl UpwindOperator for Upwind4h2 { impl UpwindOperator1d for Upwind4h2 {
fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr2(Self::DISS_BLOCK).view(),
ndarray::arr1(Self::DISS_DIAG).view(), ndarray::arr1(Self::DISS_DIAG).view(),

10
sbp/src/operators/upwind9.rs
View File

@@ -1,4 +1,4 @@
use super::{SbpOperator, UpwindOperator}; use super::*;
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayViewMut1}; use ndarray::{ArrayView1, ArrayViewMut1};
@@ -44,8 +44,8 @@ impl Upwind9 {
]; ];
} }
impl SbpOperator for Upwind9 { impl SbpOperator1d for Upwind9 {
fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::BLOCK).view(), ndarray::arr2(Self::BLOCK).view(),
ndarray::arr1(Self::DIAG).view(), ndarray::arr1(Self::DIAG).view(),
@@ -61,8 +61,8 @@ impl SbpOperator for Upwind9 {
} }
} }
impl UpwindOperator for Upwind9 { impl UpwindOperator1d for Upwind9 {
fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr2(Self::DISS_BLOCK).view(),
ndarray::arr1(Self::DISS_DIAG).view(), ndarray::arr1(Self::DISS_DIAG).view(),

16
sbp/src/operators/upwind9h2.rs
View File

@@ -1,4 +1,4 @@
use super::{SbpOperator, UpwindOperator}; use super::*;
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayViewMut1}; use ndarray::{ArrayView1, ArrayViewMut1};
@@ -44,8 +44,8 @@ impl Upwind9h2 {
]; ];
} }
impl SbpOperator for Upwind9h2 { impl SbpOperator1d for Upwind9h2 {
fn diff1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::BLOCK).view(), ndarray::arr2(Self::BLOCK).view(),
ndarray::arr1(Self::DIAG).view(), ndarray::arr1(Self::DIAG).view(),
@@ -72,25 +72,25 @@ fn upwind9h2_test() {
let mut res = ndarray::Array1::zeros(nx); let mut res = ndarray::Array1::zeros(nx);
Upwind9h2.diff1d(x.view(), res.view_mut()); Upwind9h2.diff(x.view(), res.view_mut());
let ans = &x * 0.0 + 1.0; let ans = &x * 0.0 + 1.0;
approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4); approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
res.fill(0.0); res.fill(0.0);
let y = &x * &x / 2.0; let y = &x * &x / 2.0;
Upwind9h2.diff1d(y.view(), res.view_mut()); Upwind9h2.diff(y.view(), res.view_mut());
let ans = &x; let ans = &x;
approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4); approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-4);
res.fill(0.0); res.fill(0.0);
let y = &x * &x * &x / 3.0; let y = &x * &x * &x / 3.0;
Upwind9h2.diff1d(y.view(), res.view_mut()); Upwind9h2.diff(y.view(), res.view_mut());
let ans = &x * &x; let ans = &x * &x;
approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2); approx::assert_abs_diff_eq!(&res, &ans, epsilon = 1e-2);
} }
impl UpwindOperator for Upwind9h2 { impl UpwindOperator1d for Upwind9h2 {
fn diss1d(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(
ndarray::arr2(Self::DISS_BLOCK).view(), ndarray::arr2(Self::DISS_BLOCK).view(),
ndarray::arr1(Self::DISS_DIAG).view(), ndarray::arr1(Self::DISS_DIAG).view(),