move boundaryextractor

This commit is contained in:
Magnus Ulimoen 2020-04-10 12:30:18 +02:00
parent 032c1d9436
commit e7e40c759b
2 changed files with 110 additions and 105 deletions

View File

@ -17,7 +17,7 @@ struct System<T: operators::UpwindOperator> {
grids: Vec<grid::Grid>,
metrics: Vec<grid::Metrics<T>>,
bt: Vec<euler::BoundaryCharacteristics>,
eb: Vec<MaybeBoundary>,
eb: Vec<euler::BoundaryStorage>,
time: Float,
}
@ -40,32 +40,7 @@ impl<T: operators::UpwindOperator> System<T> {
let eb = bt
.iter()
.zip(&grids)
.map(|(bt, grid)| MaybeBoundary {
n: match bt.north {
euler::BoundaryCharacteristic::Vortex(_) => {
Some(ndarray::Array2::zeros((4, grid.nx())))
}
_ => None,
},
s: match bt.south {
euler::BoundaryCharacteristic::Vortex(_) => {
Some(ndarray::Array2::zeros((4, grid.nx())))
}
_ => None,
},
e: match bt.east {
euler::BoundaryCharacteristic::Vortex(_) => {
Some(ndarray::Array2::zeros((4, grid.ny())))
}
_ => None,
},
w: match bt.west {
euler::BoundaryCharacteristic::Vortex(_) => {
Some(ndarray::Array2::zeros((4, grid.ny())))
}
_ => None,
},
})
.map(|(bt, grid)| euler::BoundaryStorage::new(bt, grid))
.collect();
Self {
@ -168,7 +143,13 @@ impl<T: operators::UpwindOperator> System<T> {
s.scope(|s| {
let fields = &self.fnext;
let bt = extract_boundaries(&fields, &mut self.bt, &mut self.eb, &self.grids, time);
let bt = euler::extract_boundaries(
&fields,
&mut self.bt,
&mut self.eb,
&self.grids,
time,
);
for ((((prev, fut), metrics), wb), bt) in fields
.iter()
.zip(&mut self.k[i])
@ -183,83 +164,6 @@ impl<T: operators::UpwindOperator> System<T> {
}
}
fn extract_boundaries<'a>(
fields: &'a [euler::Field],
bt: &[euler::BoundaryCharacteristics],
eb: &'a mut [MaybeBoundary],
grids: &[grid::Grid],
time: Float,
) -> Vec<euler::BoundaryTerms<'a>> {
bt.iter()
.zip(eb)
.zip(grids)
.zip(fields)
.map(|(((bt, eb), grid), field)| euler::BoundaryTerms {
north: match bt.north {
euler::BoundaryCharacteristic::This => field.south(),
euler::BoundaryCharacteristic::Grid(g) => fields[g].south(),
euler::BoundaryCharacteristic::Vortex(v) => {
let field = eb.n.as_mut().unwrap();
vortexify(field.view_mut(), grid.north(), v, time);
field.view()
}
},
south: match bt.south {
euler::BoundaryCharacteristic::This => field.north(),
euler::BoundaryCharacteristic::Grid(g) => fields[g].north(),
euler::BoundaryCharacteristic::Vortex(v) => {
let field = eb.s.as_mut().unwrap();
vortexify(field.view_mut(), grid.south(), v, time);
field.view()
}
},
west: match bt.west {
euler::BoundaryCharacteristic::This => field.east(),
euler::BoundaryCharacteristic::Grid(g) => fields[g].east(),
euler::BoundaryCharacteristic::Vortex(v) => {
let field = eb.w.as_mut().unwrap();
vortexify(field.view_mut(), grid.west(), v, time);
field.view()
}
},
east: match bt.east {
euler::BoundaryCharacteristic::This => field.west(),
euler::BoundaryCharacteristic::Grid(g) => fields[g].west(),
euler::BoundaryCharacteristic::Vortex(v) => {
let field = eb.e.as_mut().unwrap();
vortexify(field.view_mut(), grid.east(), v, time);
field.view()
}
},
})
.collect()
}
#[derive(Debug, Clone)]
struct MaybeBoundary {
n: Option<ndarray::Array2<Float>>,
s: Option<ndarray::Array2<Float>>,
e: Option<ndarray::Array2<Float>>,
w: Option<ndarray::Array2<Float>>,
}
fn vortexify(
mut field: ndarray::ArrayViewMut2<Float>,
yx: (ndarray::ArrayView1<Float>, ndarray::ArrayView1<Float>),
v: euler::VortexParameters,
t: Float,
) {
let mut fiter = field.outer_iter_mut();
let (rho, rhou, rhov, e) = (
fiter.next().unwrap(),
fiter.next().unwrap(),
fiter.next().unwrap(),
fiter.next().unwrap(),
);
let (y, x) = yx;
euler::vortex(rho, rhou, rhov, e, x, y, t, v);
}
#[derive(Debug, StructOpt)]
struct Options {
json: std::path::PathBuf,

View File

@ -620,6 +620,107 @@ fn boundary_extractor<'a>(
}
}
pub fn extract_boundaries<'a>(
fields: &'a [Field],
bt: &[BoundaryCharacteristics],
eb: &'a mut [BoundaryStorage],
grids: &[Grid],
time: Float,
) -> Vec<BoundaryTerms<'a>> {
bt.iter()
.zip(eb)
.zip(grids)
.zip(fields)
.map(|(((bt, eb), grid), field)| BoundaryTerms {
north: match bt.north {
BoundaryCharacteristic::This => field.south(),
BoundaryCharacteristic::Grid(g) => fields[g].south(),
BoundaryCharacteristic::Vortex(v) => {
let field = eb.n.as_mut().unwrap();
vortexify(field.view_mut(), grid.north(), v, time);
field.view()
}
},
south: match bt.south {
BoundaryCharacteristic::This => field.north(),
BoundaryCharacteristic::Grid(g) => fields[g].north(),
BoundaryCharacteristic::Vortex(v) => {
let field = eb.s.as_mut().unwrap();
vortexify(field.view_mut(), grid.south(), v, time);
field.view()
}
},
west: match bt.west {
BoundaryCharacteristic::This => field.east(),
BoundaryCharacteristic::Grid(g) => fields[g].east(),
BoundaryCharacteristic::Vortex(v) => {
let field = eb.w.as_mut().unwrap();
vortexify(field.view_mut(), grid.west(), v, time);
field.view()
}
},
east: match bt.east {
BoundaryCharacteristic::This => field.west(),
BoundaryCharacteristic::Grid(g) => fields[g].west(),
BoundaryCharacteristic::Vortex(v) => {
let field = eb.e.as_mut().unwrap();
vortexify(field.view_mut(), grid.east(), v, time);
field.view()
}
},
})
.collect()
}
#[derive(Debug, Clone)]
/// Used for storing boundary elements
pub struct BoundaryStorage {
pub n: Option<ndarray::Array2<Float>>,
pub s: Option<ndarray::Array2<Float>>,
pub e: Option<ndarray::Array2<Float>>,
pub w: Option<ndarray::Array2<Float>>,
}
impl BoundaryStorage {
pub fn new(bt: &BoundaryCharacteristics, grid: &Grid) -> Self {
Self {
n: match bt.north {
BoundaryCharacteristic::Vortex(_) => Some(ndarray::Array2::zeros((4, grid.nx()))),
_ => None,
},
s: match bt.south {
BoundaryCharacteristic::Vortex(_) => Some(ndarray::Array2::zeros((4, grid.nx()))),
_ => None,
},
e: match bt.east {
BoundaryCharacteristic::Vortex(_) => Some(ndarray::Array2::zeros((4, grid.ny()))),
_ => None,
},
w: match bt.west {
BoundaryCharacteristic::Vortex(_) => Some(ndarray::Array2::zeros((4, grid.ny()))),
_ => None,
},
}
}
}
fn vortexify(
mut field: ndarray::ArrayViewMut2<Float>,
yx: (ndarray::ArrayView1<Float>, ndarray::ArrayView1<Float>),
v: VortexParameters,
t: Float,
) {
let mut fiter = field.outer_iter_mut();
let (rho, rhou, rhov, e) = (
fiter.next().unwrap(),
fiter.next().unwrap(),
fiter.next().unwrap(),
fiter.next().unwrap(),
);
let (y, x) = yx;
vortex(rho, rhou, rhov, e, x, y, t, v);
}
#[allow(non_snake_case)]
/// Boundary conditions (SAT)
fn SAT_characteristics<SBP: SbpOperator>(