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

Introduce Diff type for euler

This commit is contained in:
Magnus Ulimoen 2021-03-30 18:46:28 +02:00
parent 27abaebc2c
commit 5f7d38dd55
2 changed files with 67 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
euler/src/lib.rs
View File

@ -16,7 +16,7 @@ pub const GAMMA: Float = 1.4;
#[derive(Debug)]
pub struct System<SBP: SbpOperator2d> {
sys: (Field, Field),
k: [Field; 4],
k: [Diff; 4],
wb: WorkBuffers,
grid: (Grid, Metrics),
op: SBP,
@ -34,10 +34,10 @@ impl<SBP: SbpOperator2d> System<SBP> {
sys: (Field::new(ny, nx), Field::new(ny, nx)),
grid: (grid, metrics),
k: [
Field::new(ny, nx),
Field::new(ny, nx),
Field::new(ny, nx),
Field::new(ny, nx),
Diff::zeros((ny, nx)),
Diff::zeros((ny, nx)),
Diff::zeros((ny, nx)),
Diff::zeros((ny, nx)),
],
wb: WorkBuffers::new(ny, nx),
op,
@ -55,7 +55,7 @@ impl<SBP: SbpOperator2d> System<SBP> {
let wb = &mut self.wb.0;
let grid = &self.grid.0;
let metrics = &self.grid.1;
let rhs_trad = |k: &mut Field, y: &Field, _time: Float| {
let rhs_trad = |k: &mut Diff, y: &Field, _time: Float| {
let boundaries = boundary_extractor(y, grid, &bc);
RHS_trad(op, k, y, metrics, &boundaries, wb)
};
@ -128,7 +128,7 @@ impl<UO: UpwindOperator2d + SbpOperator2d> System<UO> {
let op = &self.op;
let grid = &self.grid;
let wb = &mut self.wb.0;
let rhs_upwind = |k: &mut Field, y: &Field, _time: Float| {
let rhs_upwind = |k: &mut Diff, y: &Field, _time: Float| {
let (grid, metrics) = grid;
let boundaries = boundary_extractor(y, grid, &bc);
RHS_upwind(op, k, y, metrics, &boundaries, wb)
@ -153,7 +153,7 @@ impl<UO: UpwindOperator2d + SbpOperator2d> System<UO> {
let op = &self.op;
let grid = &self.grid;
let wb = &mut self.wb.0;
let mut rhs_upwind = |k: &mut Field, y: &Field, _time: Float| {
let mut rhs_upwind = |k: &mut Diff, y: &Field, _time: Float| {
let (grid, metrics) = grid;
let boundaries = boundary_extractor(y, grid, &bc);
RHS_upwind(op, k, y, metrics, &boundaries, wb)
@ -186,9 +186,13 @@ impl<UO: UpwindOperator2d + SbpOperator2d> System<UO> {
/// A 4 x ny x nx array
pub struct Field(pub(crate) Array3<Float>);
#[derive(Clone, Debug)]
/// A 4 x ny x nx array
pub struct Diff(pub(crate) Array3<Float>);
impl integrate::Integrable for Field {
type State = Field;
type Diff = Field;
type Diff = Diff;
fn assign(s: &mut Self::State, o: &Self::State) {
s.0.assign(&o.0);
@ -288,17 +292,21 @@ impl Field {
pub fn west(&self) -> ArrayView2<Float> {
self.slice(s![.., .., 0])
}
#[allow(unused)]
fn north_mut(&mut self) -> ArrayViewMut2<Float> {
let ny = self.ny();
self.slice_mut(s![.., ny - 1, ..])
}
#[allow(unused)]
fn south_mut(&mut self) -> ArrayViewMut2<Float> {
self.slice_mut(s![.., 0, ..])
}
#[allow(unused)]
fn east_mut(&mut self) -> ArrayViewMut2<Float> {
let nx = self.nx();
self.slice_mut(s![.., .., nx - 1])
}
#[allow(unused)]
fn west_mut(&mut self) -> ArrayViewMut2<Float> {
self.slice_mut(s![.., .., 0])
}
@ -448,6 +456,34 @@ fn h2_diff() {
assert!((field0.h2_err(&field1, &SBP8).powi(2) - 4.0).abs() < 1e-3);
}
// Lazy inheritance...
impl core::ops::Deref for Diff {
type Target = Array3<Float>;
fn deref(&self) -> &Self::Target {
&self.0
}
}
impl Diff {
pub fn zeros((ny, nx): (usize, usize)) -> Self {
Self(Array3::zeros((4, ny, nx)))
}
fn north_mut(&mut self) -> ArrayViewMut2<Float> {
let ny = self.shape()[1];
self.0.slice_mut(s![.., ny - 1, ..])
}
fn south_mut(&mut self) -> ArrayViewMut2<Float> {
self.0.slice_mut(s![.., 0, ..])
}
fn east_mut(&mut self) -> ArrayViewMut2<Float> {
let nx = self.shape()[2];
self.0.slice_mut(s![.., .., nx - 1])
}
fn west_mut(&mut self) -> ArrayViewMut2<Float> {
self.0.slice_mut(s![.., .., 0])
}
}
fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Float {
(gamma - 1.0) * (e - (rhou * rhou + rhov * rhov) / (2.0 * rho))
}
@ -455,7 +491,7 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo
#[allow(non_snake_case)]
pub fn RHS_trad(
op: &dyn SbpOperator2d,
k: &mut Field,
k: &mut Diff,
y: &Field,
metrics: &Metrics,
boundaries: &BoundaryTerms,
@ -490,7 +526,7 @@ pub fn RHS_trad(
#[allow(non_snake_case)]
pub fn RHS_upwind(
op: &dyn SbpOperator2d,
k: &mut Field,
k: &mut Diff,
y: &Field,
metrics: &Metrics,
boundaries: &BoundaryTerms,
@ -880,17 +916,20 @@ fn vortexify(
/// Boundary conditions (SAT)
fn SAT_characteristics(
op: &dyn SbpOperator2d,
k: &mut Field,
k: &mut Diff,
y: &Field,
metrics: &Metrics,
boundaries: &BoundaryTerms,
) {
let ny = y.ny();
let nx = y.nx();
// North boundary
{
let hi = if op.is_h2eta() {
(k.ny() - 2) as Float / op.heta()[0]
(ny - 2) as Float / op.heta()[0]
} else {
(k.ny() - 1) as Float / op.heta()[0]
(ny - 1) as Float / op.heta()[0]
};
let sign = -1.0;
let tau = 1.0;
@ -910,9 +949,9 @@ fn SAT_characteristics(
// South boundary
{
let hi = if op.is_h2eta() {
(k.ny() - 2) as Float / op.heta()[0]
(ny - 2) as Float / op.heta()[0]
} else {
(k.ny() - 1) as Float / op.heta()[0]
(ny - 1) as Float / op.heta()[0]
};
let sign = 1.0;
let tau = -1.0;
@ -932,9 +971,9 @@ fn SAT_characteristics(
// West Boundary
{
let hi = if op.is_h2xi() {
(k.nx() - 2) as Float / op.hxi()[0]
(nx - 2) as Float / op.hxi()[0]
} else {
(k.nx() - 1) as Float / op.hxi()[0]
(nx - 1) as Float / op.hxi()[0]
};
let sign = 1.0;
let tau = -1.0;
@ -954,9 +993,9 @@ fn SAT_characteristics(
// East Boundary
{
let hi = if op.is_h2xi() {
(k.nx() - 2) as Float / op.hxi()[0]
(nx - 2) as Float / op.hxi()[0]
} else {
(k.nx() - 1) as Float / op.hxi()[0]
(nx - 1) as Float / op.hxi()[0]
};
let sign = -1.0;
let tau = 1.0;

12
multigrid/src/main.rs
View File

@ -12,7 +12,7 @@ struct System {
fnow: Vec<euler::Field>,
fnext: Vec<euler::Field>,
wb: Vec<euler::WorkBuffers>,
k: [Vec<euler::Field>; 4],
k: [Vec<euler::Diff>; 4],
grids: Vec<grid::Grid>,
metrics: Vec<grid::Metrics>,
bt: Vec<euler::BoundaryCharacteristics>,
@ -26,7 +26,7 @@ pub(crate) static MULTITHREAD: AtomicBool = AtomicBool::new(false);
impl integrate::Integrable for System {
type State = Vec<euler::Field>;
type Diff = Vec<euler::Field>;
type Diff = Vec<euler::Diff>;
fn assign(s: &mut Self::State, o: &Self::State) {
if MULTITHREAD.load(Ordering::Acquire) {
s.par_iter_mut()
@ -66,7 +66,11 @@ impl System {
.iter()
.map(|g| euler::WorkBuffers::new(g.ny(), g.nx()))
.collect();
let k = [fnow.clone(), fnow.clone(), fnow.clone(), fnow.clone()];
let k = grids
.iter()
.map(|g| euler::Diff::zeros((g.ny(), g.nx())))
.collect::<Vec<_>>();
let k = [k.clone(), k.clone(), k.clone(), k];
let metrics = grids
.iter()
.zip(&operators)
@ -107,7 +111,7 @@ impl System {
let eb = &mut self.eb;
let operators = &self.operators;
let rhs = move |fut: &mut Vec<euler::Field>, prev: &Vec<euler::Field>, time: Float| {
let rhs = move |fut: &mut Vec<euler::Diff>, prev: &Vec<euler::Field>, time: Float| {
let prev_all = &prev;
if MULTITHREAD.load(Ordering::Acquire) {
rayon::scope(|s| {