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

introduce two diff SBP operators on same domain

This commit is contained in:
Magnus Ulimoen 2020-04-13 00:56:53 +02:00
parent a58892044d
commit 9c917abb16
4 changed files with 49 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
multigrid/src/main.rs
View File

@ -16,11 +16,11 @@ struct System {
} }
enum Metrics { enum Metrics {
Upwind4(grid::Metrics<operators::Upwind4>), Upwind4(grid::Metrics<operators::Upwind4, operators::Upwind4>),
Upwind9(grid::Metrics<operators::Upwind9>), Upwind9(grid::Metrics<operators::Upwind9, operators::Upwind9>),
Upwind4h2(grid::Metrics<operators::Upwind4h2>), Upwind4h2(grid::Metrics<operators::Upwind4h2, operators::Upwind4h2>),
Trad4(grid::Metrics<operators::SBP4>), Trad4(grid::Metrics<operators::SBP4, operators::SBP4>),
Trad8(grid::Metrics<operators::SBP8>), Trad8(grid::Metrics<operators::SBP8, operators::SBP8>),
} }
impl System { impl System {
@ -42,11 +42,20 @@ impl System {
let metrics = grids let metrics = grids
.iter() .iter()
.map(|g| match operator { .map(|g| match operator {
"upwind4" => Metrics::Upwind4(g.metrics::<operators::Upwind4>().unwrap()), "upwind4" => Metrics::Upwind4(
"upwind9" => Metrics::Upwind9(g.metrics::<operators::Upwind9>().unwrap()), g.metrics::<operators::Upwind4, operators::Upwind4>()
"upwind4h2" => Metrics::Upwind4h2(g.metrics::<operators::Upwind4h2>().unwrap()), .unwrap(),
"trad4" => Metrics::Trad4(g.metrics::<operators::SBP4>().unwrap()), ),
"trad8" => Metrics::Trad8(g.metrics::<operators::SBP8>().unwrap()), "upwind9" => Metrics::Upwind9(
g.metrics::<operators::Upwind9, operators::Upwind9>()
.unwrap(),
),
"upwind4h2" => Metrics::Upwind4h2(
g.metrics::<operators::Upwind4h2, operators::Upwind4h2>()
.unwrap(),
),
"trad4" => Metrics::Trad4(g.metrics::<operators::SBP4, operators::SBP4>().unwrap()),
"trad8" => Metrics::Trad8(g.metrics::<operators::SBP8, operators::SBP8>().unwrap()),
op => panic!("operator {} not known", op), op => panic!("operator {} not known", op),
}) })
.collect::<Vec<_>>(); .collect::<Vec<_>>();

12
sbp/src/euler.rs
View File

@ -14,7 +14,7 @@ pub struct System<SBP: SbpOperator> {
sys: (Field, Field), sys: (Field, Field),
k: [Field; 4], k: [Field; 4],
wb: WorkBuffers, wb: WorkBuffers,
grid: (Grid, Metrics<SBP>), grid: (Grid, Metrics<SBP, SBP>),
} }
impl<SBP: SbpOperator> System<SBP> { impl<SBP: SbpOperator> System<SBP> {
@ -403,7 +403,7 @@ fn pressure(gamma: Float, rho: Float, rhou: Float, rhov: Float, e: Float) -> Flo
pub fn RHS_trad<SBP: SbpOperator>( pub fn RHS_trad<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics<SBP>, metrics: &Metrics<SBP, SBP>,
boundaries: &BoundaryTerms, boundaries: &BoundaryTerms,
tmp: &mut (Field, Field, Field, Field, Field, Field), tmp: &mut (Field, Field, Field, Field, Field, Field),
) { ) {
@ -437,7 +437,7 @@ pub fn RHS_trad<SBP: SbpOperator>(
pub fn RHS_upwind<UO: UpwindOperator>( pub fn RHS_upwind<UO: UpwindOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics<UO>, metrics: &Metrics<UO, UO>,
boundaries: &BoundaryTerms, boundaries: &BoundaryTerms,
tmp: &mut (Field, Field, Field, Field, Field, Field), tmp: &mut (Field, Field, Field, Field, Field, Field),
) { ) {
@ -477,7 +477,7 @@ pub fn RHS_upwind<UO: UpwindOperator>(
fn upwind_dissipation<UO: UpwindOperator>( fn upwind_dissipation<UO: UpwindOperator>(
k: (&mut Field, &mut Field), k: (&mut Field, &mut Field),
y: &Field, y: &Field,
metrics: &Metrics<UO>, metrics: &Metrics<UO, UO>,
tmp: (&mut Field, &mut Field), tmp: (&mut Field, &mut Field),
) { ) {
let n = y.nx() * y.ny(); let n = y.nx() * y.ny();
@ -541,7 +541,7 @@ fn upwind_dissipation<UO: UpwindOperator>(
UO::disseta(tmp.1.e(), k.1.e_mut()); UO::disseta(tmp.1.e(), k.1.e_mut());
} }
fn fluxes<SBP: SbpOperator>(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics<SBP>) { fn fluxes<SBP: SbpOperator>(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics<SBP, SBP>) {
let j_dxi_dx = metrics.detj_dxi_dx.view(); let j_dxi_dx = metrics.detj_dxi_dx.view();
let j_dxi_dy = metrics.detj_dxi_dy.view(); let j_dxi_dy = metrics.detj_dxi_dy.view();
let j_deta_dx = metrics.detj_deta_dx.view(); let j_deta_dx = metrics.detj_deta_dx.view();
@ -815,7 +815,7 @@ fn vortexify(
fn SAT_characteristics<SBP: SbpOperator>( fn SAT_characteristics<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics<SBP>, metrics: &Metrics<SBP, SBP>,
boundaries: &BoundaryTerms, boundaries: &BoundaryTerms,
) { ) {
// North boundary // North boundary

27
sbp/src/grid.rs
View File

@ -8,9 +8,10 @@ pub struct Grid {
} }
#[derive(Debug, Clone)] #[derive(Debug, Clone)]
pub struct Metrics<SBP> pub struct Metrics<SBPeta, SBPxi>
where where
SBP: super::operators::SbpOperator, SBPeta: super::operators::SbpOperator,
SBPxi: super::operators::SbpOperator,
{ {
pub(crate) detj: Array2<Float>, pub(crate) detj: Array2<Float>,
pub(crate) detj_dxi_dx: Array2<Float>, pub(crate) detj_dxi_dx: Array2<Float>,
@ -18,7 +19,8 @@ where
pub(crate) detj_deta_dx: Array2<Float>, pub(crate) detj_deta_dx: Array2<Float>,
pub(crate) detj_deta_dy: Array2<Float>, pub(crate) detj_deta_dy: Array2<Float>,
operator: std::marker::PhantomData<SBP>, operatoreta: std::marker::PhantomData<SBPeta>,
operatorxi: std::marker::PhantomData<SBPxi>,
} }
impl Grid { impl Grid {
@ -41,9 +43,9 @@ impl Grid {
self.y.view() self.y.view()
} }
pub fn metrics<SBP: super::operators::SbpOperator>( pub fn metrics<SBPx: super::operators::SbpOperator, SBPy: super::operators::SbpOperator>(
&self, &self,
) -> Result<Metrics<SBP>, ndarray::ShapeError> { ) -> Result<Metrics<SBPx, SBPy>, ndarray::ShapeError> {
Metrics::new(self) Metrics::new(self)
} }
@ -73,7 +75,9 @@ impl Grid {
} }
} }
impl<SBP: super::operators::SbpOperator> Metrics<SBP> { impl<SBPeta: super::operators::SbpOperator, SBPxi: super::operators::SbpOperator>
Metrics<SBPeta, SBPxi>
{
fn new(grid: &Grid) -> Result<Self, ndarray::ShapeError> { fn new(grid: &Grid) -> Result<Self, ndarray::ShapeError> {
let ny = grid.ny(); let ny = grid.ny();
let nx = grid.nx(); let nx = grid.nx();
@ -81,13 +85,13 @@ impl<SBP: super::operators::SbpOperator> Metrics<SBP> {
let y = &grid.y; let y = &grid.y;
let mut dx_dxi = Array2::zeros((ny, nx)); let mut dx_dxi = Array2::zeros((ny, nx));
SBP::diffxi(x.view(), dx_dxi.view_mut()); SBPxi::diffxi(x.view(), dx_dxi.view_mut());
let mut dx_deta = Array2::zeros((ny, nx)); let mut dx_deta = Array2::zeros((ny, nx));
SBP::diffeta(x.view(), dx_deta.view_mut()); SBPeta::diffeta(x.view(), dx_deta.view_mut());
let mut dy_dxi = Array2::zeros((ny, nx)); let mut dy_dxi = Array2::zeros((ny, nx));
SBP::diffxi(y.view(), dy_dxi.view_mut()); SBPxi::diffxi(y.view(), dy_dxi.view_mut());
let mut dy_deta = Array2::zeros((ny, nx)); let mut dy_deta = Array2::zeros((ny, nx));
SBP::diffeta(y.view(), dy_deta.view_mut()); SBPeta::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,
@ -118,7 +122,8 @@ impl<SBP: super::operators::SbpOperator> Metrics<SBP> {
detj_dxi_dy, detj_dxi_dy,
detj_deta_dx, detj_deta_dx,
detj_deta_dy, detj_deta_dy,
operator: std::marker::PhantomData, operatorxi: std::marker::PhantomData,
operatoreta: std::marker::PhantomData,
}) })
} }
} }

16
sbp/src/maxwell.rs
View File

@ -78,7 +78,7 @@ pub struct System<SBP: SbpOperator> {
sys: (Field, Field), sys: (Field, Field),
wb: WorkBuffers, wb: WorkBuffers,
grid: Grid, grid: Grid,
metrics: Metrics<SBP>, metrics: Metrics<SBP, SBP>,
} }
impl<SBP: SbpOperator> System<SBP> { impl<SBP: SbpOperator> System<SBP> {
@ -119,7 +119,7 @@ impl<SBP: SbpOperator> System<SBP> {
fut: &mut Field, fut: &mut Field,
prev: &Field, prev: &Field,
_time: Float, _time: Float,
c: &(&Grid, &Metrics<SBP>), c: &(&Grid, &Metrics<SBP, SBP>),
m: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), m: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
let (grid, metrics) = c; let (grid, metrics) = c;
@ -147,7 +147,7 @@ impl<UO: UpwindOperator> System<UO> {
fut: &mut Field, fut: &mut Field,
prev: &Field, prev: &Field,
_time: Float, _time: Float,
c: &(&Grid, &Metrics<UO>), c: &(&Grid, &Metrics<UO, UO>),
m: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), m: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
let (grid, metrics) = c; let (grid, metrics) = c;
@ -199,7 +199,7 @@ fn RHS<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
_grid: &Grid, _grid: &Grid,
metrics: &Metrics<SBP>, metrics: &Metrics<SBP, SBP>,
tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
fluxes(k, y, metrics, tmp); fluxes(k, y, metrics, tmp);
@ -223,7 +223,7 @@ fn RHS_upwind<UO: UpwindOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
_grid: &Grid, _grid: &Grid,
metrics: &Metrics<UO>, metrics: &Metrics<UO, UO>,
tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
fluxes(k, y, metrics, tmp); fluxes(k, y, metrics, tmp);
@ -246,7 +246,7 @@ fn RHS_upwind<UO: UpwindOperator>(
fn fluxes<SBP: SbpOperator>( fn fluxes<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics<SBP>, metrics: &Metrics<SBP, SBP>,
tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
// ex = hz_y // ex = hz_y
@ -320,7 +320,7 @@ fn fluxes<SBP: SbpOperator>(
fn dissipation<UO: UpwindOperator>( fn dissipation<UO: UpwindOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics<UO>, metrics: &Metrics<UO, UO>,
tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>), tmp: &mut (Array2<Float>, Array2<Float>, Array2<Float>, Array2<Float>),
) { ) {
// ex component // ex component
@ -421,7 +421,7 @@ pub struct BoundaryTerms {
fn SAT_characteristics<SBP: SbpOperator>( fn SAT_characteristics<SBP: SbpOperator>(
k: &mut Field, k: &mut Field,
y: &Field, y: &Field,
metrics: &Metrics<SBP>, metrics: &Metrics<SBP, SBP>,
boundaries: &BoundaryTerms, boundaries: &BoundaryTerms,
) { ) {
let ny = y.ny(); let ny = y.ny();