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

revisit SBP Traits

This commit is contained in:
Magnus Ulimoen 2021-01-17 15:37:45 +01:00
parent 6f3a810cd3
commit 1f15bcc056
13 changed files with 229 additions and 219 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
euler/src/lib.rs
View File

@ -115,7 +115,7 @@ impl<SBP: SbpOperator2d> System<SBP> {
}
}
impl<UO: UpwindOperator2d> System<UO> {
impl<UO: UpwindOperator2d + SbpOperator2d> System<UO> {
pub fn advance_upwind(&mut self, dt: Float) {
let bc = BoundaryCharacteristics {
north: BoundaryCharacteristic::This,
@ -558,7 +558,7 @@ pub fn RHS_trad(
#[allow(non_snake_case)]
pub fn RHS_upwind(
op: &dyn UpwindOperator2d,
op: &dyn SbpOperator2d,
k: &mut Field,
y: &Field,
metrics: &Metrics,
@ -583,7 +583,14 @@ pub fn RHS_upwind(
let ad_xi = &mut tmp.4;
let ad_eta = &mut tmp.5;
upwind_dissipation(op, (ad_xi, ad_eta), y, metrics, (&mut tmp.0, &mut tmp.1));
let diss_op = op.upwind().expect("This is not an upwind operator");
upwind_dissipation(
&*diss_op,
(ad_xi, ad_eta),
y,
metrics,
(&mut tmp.0, &mut tmp.1),
);
azip!((out in &mut k.0,
eflux in &dE.0,
@ -594,7 +601,7 @@ pub fn RHS_upwind(
*out = (-eflux - fflux + ad_xi + ad_eta)/detj
});
SAT_characteristics(op.as_sbp(), k, y, metrics, boundaries);
SAT_characteristics(op, k, y, metrics, boundaries);
}
#[allow(clippy::many_single_char_names)]

4
maxwell/src/lib.rs
View File

@ -196,7 +196,7 @@ impl<SBP: SbpOperator2d> System<SBP> {
}
}
impl<UO: UpwindOperator2d> System<UO> {
impl<UO: SbpOperator2d + UpwindOperator2d> System<UO> {
/// Using artificial dissipation with the upwind operator
pub fn advance_upwind(&mut self, dt: Float) {
let op = &self.op;
@ -271,7 +271,7 @@ fn RHS<SBP: SbpOperator2d>(
}
#[allow(non_snake_case)]
fn RHS_upwind<UO: UpwindOperator2d>(
fn RHS_upwind<UO: SbpOperator2d + UpwindOperator2d>(
op: &UO,
k: &mut Field,
y: &Field,

1
multigrid/Cargo.toml
View File

@ -13,7 +13,6 @@ rayon = "1.3.0"
indicatif = "0.15.0"
structopt = "0.3.14"
ndarray = { version = "0.13.1", features = ["serde"] }
either = "1.5.3"
serde = { version = "1.0.115", features = ["derive"] }
json5 = "0.2.8"
indexmap = { version = "1.5.2", features = ["serde-1"] }

38
multigrid/src/main.rs
View File

@ -1,15 +1,12 @@
use either::*;
use structopt::StructOpt;
use sbp::operators::{SbpOperator2d, UpwindOperator2d};
use sbp::operators::SbpOperator2d;
use sbp::*;
mod file;
mod parsing;
use file::*;
pub(crate) type DiffOp = Either<Box<dyn SbpOperator2d>, Box<dyn UpwindOperator2d>>;
struct System {
fnow: Vec<euler::Field>,
fnext: Vec<euler::Field>,
@ -20,14 +17,14 @@ struct System {
bt: Vec<euler::BoundaryCharacteristics>,
eb: Vec<euler::BoundaryStorage>,
time: Float,
operators: Vec<DiffOp>,
operators: Vec<Box<dyn SbpOperator2d>>,
}
impl System {
fn new(
grids: Vec<grid::Grid>,
bt: Vec<euler::BoundaryCharacteristics>,
operators: Vec<DiffOp>,
operators: Vec<Box<dyn SbpOperator2d>>,
) -> Self {
let fnow = grids
.iter()
@ -42,10 +39,7 @@ impl System {
let metrics = grids
.iter()
.zip(&operators)
.map(|(g, op)| {
let sbpop: &dyn SbpOperator2d = op.as_ref().either(|op| &**op, |uo| uo.as_sbp());
g.metrics(sbpop).unwrap()
})
.map(|(g, op)| g.metrics(&**op).unwrap())
.collect::<Vec<_>>();
let eb = bt
@ -97,13 +91,10 @@ impl System {
{
s.spawn(move |_| {
let bc = euler::boundary_extracts(prev_all, bt, prev, grid, eb, time);
match op.as_ref() {
Left(sbp) => {
euler::RHS_trad(&**sbp, fut, prev, metrics, &bc, &mut wb.0);
}
Right(uo) => {
euler::RHS_upwind(&**uo, fut, prev, metrics, &bc, &mut wb.0);
}
if op.upwind().is_some() {
euler::RHS_upwind(&**op, fut, prev, metrics, &bc, &mut wb.0);
} else {
euler::RHS_trad(&**op, fut, prev, metrics, &bc, &mut wb.0);
}
})
}
@ -130,12 +121,10 @@ impl System {
/// Suggested maximum dt for this problem
fn max_dt(&self) -> Float {
let is_h2 = self.operators.iter().any(|op| {
op.as_ref().either(
|op| op.is_h2xi() || op.is_h2eta(),
|op| op.is_h2xi() || op.is_h2eta(),
)
});
let is_h2 = self
.operators
.iter()
.any(|op| op.is_h2xi() || op.is_h2eta());
let c_max = if is_h2 { 0.5 } else { 1.0 };
let mut max_dt: Float = Float::INFINITY;
@ -283,8 +272,7 @@ fn main() {
for ((fmod, grid), op) in sys.fnow.iter().zip(&sys.grids).zip(&sys.operators) {
let mut fvort = fmod.clone();
fvort.vortex(grid.x(), grid.y(), time, &vortexparams);
let sbpop: &dyn SbpOperator2d = op.as_ref().either(|op| &**op, |uo| uo.as_sbp());
e += fmod.h2_err(&fvort, sbpop);
e += fmod.h2_err(&fvort, &**op);
}
println!("Total error: {:e}", e);
}

42
multigrid/src/parsing.rs
View File

@ -1,5 +1,4 @@
use super::DiffOp;
use either::*;
use sbp::operators::SbpOperator2d;
use sbp::utils::h2linspace;
use sbp::Float;
@ -148,7 +147,7 @@ pub struct RuntimeConfiguration {
pub names: Vec<String>,
pub grids: Vec<sbp::grid::Grid>,
pub bc: Vec<euler::BoundaryCharacteristics>,
pub op: Vec<DiffOp>,
pub op: Vec<Box<dyn SbpOperator2d>>,
pub integration_time: Float,
pub vortex: euler::VortexParameters,
}
@ -223,32 +222,19 @@ impl Configuration {
use sbp::operators::*;
use Operator as op;
match (eta, xi) {
(op::Upwind4, op::Upwind4) => {
Right(Box::new(Upwind4) as Box<dyn UpwindOperator2d>)
let matcher = |op| -> Box<dyn SbpOperator2d> {
match op {
op::Upwind4 => Box::new(Upwind4),
op::Upwind4h2 => Box::new(Upwind4h2),
op::Upwind9 => Box::new(Upwind9),
op::Upwind9h2 => Box::new(Upwind9h2),
op::Sbp4 => Box::new(SBP4),
op::Sbp8 => Box::new(SBP8),
}
(op::Upwind4h2, op::Upwind4h2) => {
Right(Box::new(Upwind4h2) as Box<dyn UpwindOperator2d>)
}
(op::Upwind9, op::Upwind9) => {
Right(Box::new(Upwind9) as Box<dyn UpwindOperator2d>)
}
(op::Upwind9h2, op::Upwind9h2) => {
Right(Box::new(Upwind9h2) as Box<dyn UpwindOperator2d>)
}
(op::Upwind4, op::Upwind4h2) => {
Right(Box::new((&Upwind4, &Upwind4h2)) as Box<dyn UpwindOperator2d>)
}
(op::Upwind9, op::Upwind9h2) => {
Right(Box::new((&Upwind9, &Upwind9h2)) as Box<dyn UpwindOperator2d>)
}
(op::Upwind9h2, op::Upwind9) => {
Right(Box::new((&Upwind9h2, &Upwind9)) as Box<dyn UpwindOperator2d>)
}
(op::Sbp4, op::Sbp4) => Left(Box::new(SBP4) as Box<dyn SbpOperator2d>),
(op::Sbp8, op::Sbp8) => Left(Box::new(SBP8) as Box<dyn SbpOperator2d>),
_ => todo!("Combination {:?}, {:?} not implemented", eta, xi),
}
};
Box::new((matcher(eta), matcher(xi))) as Box<dyn SbpOperator2d>
})
.collect();
let bc = self

195
sbp/src/operators.rs
View File

@ -40,80 +40,6 @@ pub trait SbpOperator1d2: SbpOperator1d {
fn d1_vec(&self, n: usize, front: bool) -> sprs::CsMat<Float>;
}
pub trait SbpOperator2d: Send + Sync {
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;
fn op_xi(&self) -> &dyn SbpOperator1d;
fn op_eta(&self) -> &dyn SbpOperator1d;
}
impl<SBPeta: SbpOperator1d, SBPxi: SbpOperator1d> SbpOperator2d for (&SBPeta, &SBPxi) {
default 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()
}
fn op_xi(&self) -> &dyn SbpOperator1d {
self.1
}
fn op_eta(&self) -> &dyn SbpOperator1d {
self.0
}
}
impl<SBP: SbpOperator1d + Copy> SbpOperator2d for SBP {
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
<(&SBP, &SBP) as SbpOperator2d>::diffxi(&(self, self), prev, fut)
}
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
<(&SBP, &SBP) as SbpOperator2d>::diffeta(&(self, self), prev, fut)
}
fn hxi(&self) -> &'static [Float] {
<(&SBP, &SBP) as SbpOperator2d>::hxi(&(self, self))
}
fn heta(&self) -> &'static [Float] {
<(&SBP, &SBP) as SbpOperator2d>::heta(&(self, self))
}
fn is_h2xi(&self) -> bool {
<(&SBP, &SBP) as SbpOperator2d>::is_h2xi(&(self, self))
}
fn is_h2eta(&self) -> bool {
<(&SBP, &SBP) as SbpOperator2d>::is_h2eta(&(self, self))
}
fn op_xi(&self) -> &dyn SbpOperator1d {
self
}
fn op_eta(&self) -> &dyn SbpOperator1d {
self
}
}
pub trait UpwindOperator1d: SbpOperator1d + Send + Sync {
/// Dissipation operator
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>);
@ -123,57 +49,55 @@ pub trait UpwindOperator1d: SbpOperator1d + Send + Sync {
fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float>;
}
pub trait UpwindOperator2d: SbpOperator2d + Send + Sync {
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>);
fn as_sbp(&self) -> &dyn SbpOperator2d;
pub trait SbpOperator2d: Send + Sync {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
for (p, f) in prev.outer_iter().zip(fut.outer_iter_mut()) {
self.op_xi().diff(p, f)
}
}
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.diffxi(prev.reversed_axes(), fut.reversed_axes())
}
fn hxi(&self) -> &'static [Float] {
self.op_xi().h()
}
fn heta(&self) -> &'static [Float] {
self.op_eta().h()
}
fn is_h2xi(&self) -> bool {
self.op_xi().is_h2()
}
fn is_h2eta(&self) -> bool {
self.op_eta().is_h2()
}
fn op_xi(&self) -> &dyn SbpOperator1d;
fn op_eta(&self) -> &dyn SbpOperator1d;
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
None
}
}
pub trait UpwindOperator2d: Send + Sync {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
for (p, f) in prev.outer_iter().zip(fut.outer_iter_mut()) {
UpwindOperator2d::op_xi(self).diss(p, f)
}
}
// Assuming operator is symmetrical x/y
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.dissxi(prev.reversed_axes(), fut.reversed_axes())
}
fn op_xi(&self) -> &dyn UpwindOperator1d;
fn op_eta(&self) -> &dyn UpwindOperator1d;
}
impl<UOeta: UpwindOperator1d, UOxi: UpwindOperator1d> UpwindOperator2d for (&UOeta, &UOxi) {
default 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.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())
}
fn as_sbp(&self) -> &dyn SbpOperator2d {
self
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
self.1
}
fn op_eta(&self) -> &dyn UpwindOperator1d {
self.0
}
}
impl<UO: UpwindOperator1d + Copy> UpwindOperator2d for UO {
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
<(&UO, &UO) as UpwindOperator2d>::dissxi(&(self, self), prev, fut)
}
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
<(&UO, &UO) as UpwindOperator2d>::disseta(&(self, self), prev, fut)
}
fn as_sbp(&self) -> &dyn SbpOperator2d {
self
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
self
}
fn op_eta(&self) -> &dyn UpwindOperator1d {
self
}
}
pub trait InterpolationOperator: Send + Sync {
/// Interpolation from a grid with twice resolution
fn fine2coarse(&self, fine: ArrayView1<Float>, coarse: ArrayViewMut1<Float>);
@ -181,6 +105,37 @@ pub trait InterpolationOperator: Send + Sync {
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>);
}
impl SbpOperator2d for (Box<dyn SbpOperator2d>, Box<dyn SbpOperator2d>) {
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.1.diffxi(prev, fut)
}
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.0.diffeta(prev, fut)
}
fn op_xi(&self) -> &dyn SbpOperator1d {
self.1.op_xi()
}
fn op_eta(&self) -> &dyn SbpOperator1d {
self.0.op_eta()
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
match (self.0.upwind(), self.1.upwind()) {
(Some(u), Some(v)) => Some(Box::new((u, v))),
_ => None,
}
}
}
impl UpwindOperator2d for (Box<dyn UpwindOperator2d>, Box<dyn UpwindOperator2d>) {
fn op_xi(&self) -> &dyn UpwindOperator1d {
self.1.op_xi()
}
fn op_eta(&self) -> &dyn UpwindOperator1d {
self.0.op_eta()
}
}
#[inline(always)]
fn diff_op_1d(
block: &[&[Float]],

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

@ -71,7 +71,7 @@ impl SbpOperator1d for SBP4 {
}
}
impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &SBP4) {
impl SbpOperator2d for SBP4 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * SBP4::BLOCK.len());
@ -95,6 +95,12 @@ impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &SBP4) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
}
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
}
impl super::SbpOperator1d2 for SBP4 {

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

@ -59,7 +59,7 @@ impl SbpOperator1d for SBP8 {
}
}
impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &SBP8) {
impl SbpOperator2d for SBP8 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * SBP8::BLOCK.len());
@ -83,6 +83,12 @@ impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &SBP8) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
}
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
}
#[test]

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

@ -229,7 +229,7 @@ impl SbpOperator1d for Upwind4 {
}
}
impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &Upwind4) {
impl SbpOperator2d for Upwind4 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind4::BLOCK.len());
@ -259,6 +259,15 @@ impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &Upwind4) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
}
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
Some(Box::new(Self))
}
}
#[test]
@ -377,7 +386,7 @@ impl UpwindOperator1d for Upwind4 {
}
}
impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind4) {
impl UpwindOperator2d for Upwind4 {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind4::BLOCK.len());
@ -407,6 +416,13 @@ impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind4) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self
}
fn op_eta(&self) -> &dyn UpwindOperator1d {
&Self
}
}
#[test]

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

@ -77,7 +77,7 @@ impl SbpOperator1d for Upwind4h2 {
}
}
impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &Upwind4h2) {
impl SbpOperator2d for Upwind4h2 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind4h2::BLOCK.len());
@ -101,9 +101,18 @@ impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &Upwind4h2) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
}
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
Some(Box::new(Self))
}
}
impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind4h2) {
impl UpwindOperator2d for Upwind4h2 {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind4h2::BLOCK.len());
@ -137,6 +146,12 @@ impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind4h2) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self
}
fn op_eta(&self) -> &dyn UpwindOperator1d {
&Self
}
}
#[test]

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

@ -82,7 +82,7 @@ impl SbpOperator1d for Upwind9 {
}
}
impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &Upwind9) {
impl SbpOperator2d for Upwind9 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind9::BLOCK.len());
@ -106,6 +106,15 @@ impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &Upwind9) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
}
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
Some(Box::new(Self))
}
}
impl UpwindOperator1d for Upwind9 {
@ -136,7 +145,7 @@ impl UpwindOperator1d for Upwind9 {
}
}
impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind9) {
impl UpwindOperator2d for Upwind9 {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind9::BLOCK.len());
@ -160,6 +169,12 @@ impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind9) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self
}
fn op_eta(&self) -> &dyn UpwindOperator1d {
&Self
}
}
#[test]

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

@ -85,7 +85,7 @@ impl SbpOperator1d for Upwind9h2 {
}
}
impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &Upwind9h2) {
impl SbpOperator2d for Upwind9h2 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind9h2::BLOCK.len());
@ -109,6 +109,15 @@ impl<SBP: SbpOperator1d> SbpOperator2d for (&SBP, &Upwind9h2) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
}
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
Some(Box::new(Self))
}
}
#[test]
@ -163,7 +172,7 @@ impl UpwindOperator1d for Upwind9h2 {
}
}
impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind9h2) {
impl UpwindOperator2d for Upwind9h2 {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind9h2::BLOCK.len());
@ -197,4 +206,10 @@ impl<UO: UpwindOperator1d> UpwindOperator2d for (&UO, &Upwind9h2) {
_ => unreachable!("Should only be two elements in the strides vectors"),
}
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self
}
fn op_eta(&self) -> &dyn UpwindOperator1d {
&Self
}
}

60
shallow_water/src/lib.rs
View File

@ -62,7 +62,7 @@ pub struct System {
fnext: Field,
x: (Float, Float, usize),
y: (Float, Float, usize),
op: Box<dyn operators::UpwindOperator2d>,
op: Box<dyn operators::SbpOperator2d>,
k: [Field; 4],
}
@ -175,35 +175,37 @@ impl System {
// Upwind dissipation
if false {
let mut temp_dx = temp_dy;
azip!((dest in &mut temp, eta in now.eta(), etau in now.etau()) {
*dest = -(eta.powf(3.0/2.0)*G.sqrt() + etau.abs())/eta
});
op.dissxi(temp.view(), temp_dx.view_mut());
azip!((dest in &mut next_eta, eta in now.eta(), diss in &temp_dx) {
*dest -= eta*diss;
});
azip!((dest in &mut next_etau, etau in now.etau(), diss in &temp_dx) {
*dest -= etau*diss;
});
azip!((dest in &mut next_etav, etav in now.etav(), diss in &temp_dx) {
*dest -= etav*diss;
});
if let Some(op) = op.upwind() {
let mut temp_dx = temp_dy;
azip!((dest in &mut temp, eta in now.eta(), etau in now.etau()) {
*dest = -(eta.powf(3.0/2.0)*G.sqrt() + etau.abs())/eta
});
op.dissxi(temp.view(), temp_dx.view_mut());
azip!((dest in &mut next_eta, eta in now.eta(), diss in &temp_dx) {
*dest -= eta*diss;
});
azip!((dest in &mut next_etau, etau in now.etau(), diss in &temp_dx) {
*dest -= etau*diss;
});
azip!((dest in &mut next_etav, etav in now.etav(), diss in &temp_dx) {
*dest -= etav*diss;
});
let mut temp_dy = temp_dx;
azip!((dest in &mut temp, eta in now.eta(), etav in now.etav()) {
*dest = -(eta.powf(3.0/2.0)*G.sqrt() + etav.abs())/eta
});
op.disseta(temp.view(), temp_dy.view_mut());
azip!((dest in &mut next_eta, eta in now.eta(), diss in &temp_dy) {
*dest -= eta*diss;
});
azip!((dest in &mut next_etau, etau in now.etau(), diss in &temp_dy) {
*dest -= etau*diss;
});
azip!((dest in &mut next_etav, etav in now.etav(), diss in &temp_dy) {
*dest -= etav*diss;
});
let mut temp_dy = temp_dx;
azip!((dest in &mut temp, eta in now.eta(), etav in now.etav()) {
*dest = -(eta.powf(3.0/2.0)*G.sqrt() + etav.abs())/eta
});
op.disseta(temp.view(), temp_dy.view_mut());
azip!((dest in &mut next_eta, eta in now.eta(), diss in &temp_dy) {
*dest -= eta*diss;
});
azip!((dest in &mut next_etau, etau in now.etau(), diss in &temp_dy) {
*dest -= etau*diss;
});
azip!((dest in &mut next_etav, etav in now.etav(), diss in &temp_dy) {
*dest -= etav*diss;
});
}
}
// SAT boundaries