use workspaces
This commit is contained in:
Magnus Ulimoen
695
src/euler.rs
695
src/euler.rs
@@ -1,695 +0,0 @@
|
||||
use super::integrate;
|
||||
use super::operators::{SbpOperator, UpwindOperator};
|
||||
use super::Grid;
|
||||
use ndarray::azip;
|
||||
use ndarray::prelude::*;
|
||||
|
||||
pub const GAMMA: f32 = 1.4;
|
||||
|
||||
// A collection of buffers that allows one to efficiently
|
||||
// move to the next state
|
||||
pub struct System<SBP: SbpOperator> {
|
||||
sys: (Field, Field),
|
||||
wb: WorkBuffers,
|
||||
grid: Grid<SBP>,
|
||||
}
|
||||
|
||||
impl<SBP: SbpOperator> System<SBP> {
|
||||
pub fn new(x: ndarray::Array2<f32>, y: ndarray::Array2<f32>) -> Self {
|
||||
let grid = Grid::new(x, y).expect(
|
||||
"Could not create grid. Different number of elements compared to width*height?",
|
||||
);
|
||||
let nx = grid.nx();
|
||||
let ny = grid.ny();
|
||||
Self {
|
||||
sys: (Field::new(ny, nx), Field::new(ny, nx)),
|
||||
grid,
|
||||
wb: WorkBuffers::new(ny, nx),
|
||||
}
|
||||
}
|
||||
|
||||
pub fn advance(&mut self, dt: f32) {
|
||||
integrate::rk4(
|
||||
RHS_trad,
|
||||
&self.sys.0,
|
||||
&mut self.sys.1,
|
||||
dt,
|
||||
&self.grid,
|
||||
&mut self.wb.k,
|
||||
&mut self.wb.tmp,
|
||||
);
|
||||
std::mem::swap(&mut self.sys.0, &mut self.sys.1);
|
||||
}
|
||||
|
||||
#[allow(clippy::many_single_char_names)]
|
||||
pub fn init_with_vortex(&mut self, x0: f32, y0: f32) {
|
||||
// Should parametrise such that we have radius, drop in pressure at center, etc
|
||||
let rstar = 1.0;
|
||||
let eps = 3.0;
|
||||
#[allow(non_snake_case)]
|
||||
let M = 0.5;
|
||||
|
||||
let p_inf = 1.0 / (GAMMA * M * M);
|
||||
let t = 0.0;
|
||||
|
||||
let nx = self.grid.nx();
|
||||
let ny = self.grid.ny();
|
||||
|
||||
for j in 0..ny {
|
||||
for i in 0..nx {
|
||||
let x = self.grid.x[(j, i)];
|
||||
let y = self.grid.y[(j, i)];
|
||||
|
||||
let dx = (x - x0) - t;
|
||||
let dy = y - y0;
|
||||
let f = (1.0 - (dx * dx + dy * dy)) / (rstar * rstar);
|
||||
|
||||
use std::f32::consts::PI;
|
||||
let u =
|
||||
1.0 - eps * dy / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
|
||||
let v =
|
||||
0.0 + eps * dx / (2.0 * PI * p_inf.sqrt() * rstar * rstar) * (f / 2.0).exp();
|
||||
let rho = f32::powf(
|
||||
1.0 - eps * eps * (GAMMA - 1.0) * M * M
|
||||
/ (8.0 * PI * PI * p_inf * rstar * rstar)
|
||||
* f.exp(),
|
||||
1.0 / (GAMMA - 1.0),
|
||||
);
|
||||
assert!(rho > 0.0);
|
||||
let p = p_inf * rho.powf(GAMMA);
|
||||
assert!(p > 0.0);
|
||||
let e = p / (GAMMA - 1.0) + rho * (u * u + v * v) / 2.0;
|
||||
assert!(e > 0.0);
|
||||
|
||||
self.sys.0[(0, j, i)] = rho;
|
||||
self.sys.0[(1, j, i)] = rho * u;
|
||||
self.sys.0[(2, j, i)] = rho * v;
|
||||
self.sys.0[(3, j, i)] = e;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub fn field(&self) -> &Field {
|
||||
&self.sys.0
|
||||
}
|
||||
}
|
||||
|
||||
impl<UO: UpwindOperator> System<UO> {
|
||||
pub fn advance_upwind(&mut self, dt: f32) {
|
||||
integrate::rk4(
|
||||
RHS_upwind,
|
||||
&self.sys.0,
|
||||
&mut self.sys.1,
|
||||
dt,
|
||||
&self.grid,
|
||||
&mut self.wb.k,
|
||||
&mut self.wb.tmp,
|
||||
);
|
||||
std::mem::swap(&mut self.sys.0, &mut self.sys.1);
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
/// A 4 x ny x nx array
|
||||
pub struct Field(pub(crate) Array3<f32>);
|
||||
|
||||
impl std::ops::Deref for Field {
|
||||
type Target = Array3<f32>;
|
||||
fn deref(&self) -> &Self::Target {
|
||||
&self.0
|
||||
}
|
||||
}
|
||||
|
||||
impl std::ops::DerefMut for Field {
|
||||
fn deref_mut(&mut self) -> &mut Self::Target {
|
||||
&mut self.0
|
||||
}
|
||||
}
|
||||
|
||||
impl Field {
|
||||
pub fn new(ny: usize, nx: usize) -> Self {
|
||||
let field = Array3::zeros((4, ny, nx));
|
||||
|
||||
Self(field)
|
||||
}
|
||||
|
||||
pub fn nx(&self) -> usize {
|
||||
self.0.shape()[2]
|
||||
}
|
||||
pub fn ny(&self) -> usize {
|
||||
self.0.shape()[1]
|
||||
}
|
||||
|
||||
pub fn rho(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![0, .., ..])
|
||||
}
|
||||
pub fn rhou(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![1, .., ..])
|
||||
}
|
||||
pub fn rhov(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![2, .., ..])
|
||||
}
|
||||
pub fn e(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![3, .., ..])
|
||||
}
|
||||
|
||||
pub fn rho_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![0, .., ..])
|
||||
}
|
||||
pub fn rhou_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![1, .., ..])
|
||||
}
|
||||
pub fn rhov_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![2, .., ..])
|
||||
}
|
||||
pub fn e_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![3, .., ..])
|
||||
}
|
||||
|
||||
#[allow(unused)]
|
||||
pub fn components(
|
||||
&self,
|
||||
) -> (
|
||||
ArrayView2<f32>,
|
||||
ArrayView2<f32>,
|
||||
ArrayView2<f32>,
|
||||
ArrayView2<f32>,
|
||||
) {
|
||||
(self.rho(), self.rhou(), self.rhov(), self.e())
|
||||
}
|
||||
#[allow(unused)]
|
||||
pub fn components_mut(
|
||||
&mut self,
|
||||
) -> (
|
||||
ArrayViewMut2<f32>,
|
||||
ArrayViewMut2<f32>,
|
||||
ArrayViewMut2<f32>,
|
||||
ArrayViewMut2<f32>,
|
||||
) {
|
||||
let mut iter = self.0.outer_iter_mut();
|
||||
|
||||
let rho = iter.next().unwrap();
|
||||
let rhou = iter.next().unwrap();
|
||||
let rhov = iter.next().unwrap();
|
||||
let e = iter.next().unwrap();
|
||||
assert_eq!(iter.next(), None);
|
||||
|
||||
(rho, rhou, rhov, e)
|
||||
}
|
||||
|
||||
fn north(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![.., self.ny() - 1, ..])
|
||||
}
|
||||
fn south(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![.., 0, ..])
|
||||
}
|
||||
fn east(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![.., .., self.nx() - 1])
|
||||
}
|
||||
fn west(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![.., .., 0])
|
||||
}
|
||||
fn north_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
let ny = self.ny();
|
||||
self.slice_mut(s![.., ny - 1, ..])
|
||||
}
|
||||
fn south_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![.., 0, ..])
|
||||
}
|
||||
fn east_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
let nx = self.nx();
|
||||
self.slice_mut(s![.., .., nx - 1])
|
||||
}
|
||||
fn west_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![.., .., 0])
|
||||
}
|
||||
}
|
||||
|
||||
fn pressure(gamma: f32, rho: f32, rhou: f32, rhov: f32, e: f32) -> f32 {
|
||||
(gamma - 1.0) * (e - (rhou * rhou + rhov * rhov) / (2.0 * rho))
|
||||
}
|
||||
|
||||
#[allow(non_snake_case)]
|
||||
pub(crate) fn RHS_trad<SBP: SbpOperator>(
|
||||
k: &mut Field,
|
||||
y: &Field,
|
||||
grid: &Grid<SBP>,
|
||||
tmp: &mut (Field, Field, Field, Field, Field, Field),
|
||||
) {
|
||||
let ehat = &mut tmp.0;
|
||||
let fhat = &mut tmp.1;
|
||||
fluxes((ehat, fhat), y, grid);
|
||||
let dE = &mut tmp.2;
|
||||
let dF = &mut tmp.3;
|
||||
|
||||
SBP::diffxi(ehat.rho(), dE.rho_mut());
|
||||
SBP::diffxi(ehat.rhou(), dE.rhou_mut());
|
||||
SBP::diffxi(ehat.rhov(), dE.rhov_mut());
|
||||
SBP::diffxi(ehat.e(), dE.e_mut());
|
||||
|
||||
SBP::diffeta(fhat.rho(), dF.rho_mut());
|
||||
SBP::diffeta(fhat.rhou(), dF.rhou_mut());
|
||||
SBP::diffeta(fhat.rhov(), dF.rhov_mut());
|
||||
SBP::diffeta(fhat.e(), dF.e_mut());
|
||||
|
||||
azip!((out in &mut k.0,
|
||||
eflux in &dE.0,
|
||||
fflux in &dF.0,
|
||||
detj in &grid.detj.broadcast((4, y.ny(), y.nx())).unwrap()) {
|
||||
*out = (-eflux - fflux)/detj
|
||||
});
|
||||
|
||||
let boundaries = BoundaryTerms {
|
||||
north: Boundary::This,
|
||||
south: Boundary::This,
|
||||
west: Boundary::This,
|
||||
east: Boundary::This,
|
||||
};
|
||||
|
||||
SAT_characteristics(k, y, grid, &boundaries);
|
||||
}
|
||||
|
||||
#[allow(non_snake_case)]
|
||||
pub(crate) fn RHS_upwind<UO: UpwindOperator>(
|
||||
k: &mut Field,
|
||||
y: &Field,
|
||||
grid: &Grid<UO>,
|
||||
tmp: &mut (Field, Field, Field, Field, Field, Field),
|
||||
) {
|
||||
let ehat = &mut tmp.0;
|
||||
let fhat = &mut tmp.1;
|
||||
fluxes((ehat, fhat), y, grid);
|
||||
let dE = &mut tmp.2;
|
||||
let dF = &mut tmp.3;
|
||||
|
||||
UO::diffxi(ehat.rho(), dE.rho_mut());
|
||||
UO::diffxi(ehat.rhou(), dE.rhou_mut());
|
||||
UO::diffxi(ehat.rhov(), dE.rhov_mut());
|
||||
UO::diffxi(ehat.e(), dE.e_mut());
|
||||
|
||||
UO::diffeta(fhat.rho(), dF.rho_mut());
|
||||
UO::diffeta(fhat.rhou(), dF.rhou_mut());
|
||||
UO::diffeta(fhat.rhov(), dF.rhov_mut());
|
||||
UO::diffeta(fhat.e(), dF.e_mut());
|
||||
|
||||
let ad_xi = &mut tmp.4;
|
||||
let ad_eta = &mut tmp.5;
|
||||
upwind_dissipation((ad_xi, ad_eta), y, grid, (&mut tmp.0, &mut tmp.1));
|
||||
|
||||
azip!((out in &mut k.0,
|
||||
eflux in &dE.0,
|
||||
fflux in &dF.0,
|
||||
ad_xi in &ad_xi.0,
|
||||
ad_eta in &ad_eta.0,
|
||||
detj in &grid.detj.broadcast((4, y.ny(), y.nx())).unwrap()) {
|
||||
*out = (-eflux - fflux + ad_xi + ad_eta)/detj
|
||||
});
|
||||
|
||||
let boundaries = BoundaryTerms {
|
||||
north: Boundary::This,
|
||||
south: Boundary::This,
|
||||
west: Boundary::This,
|
||||
east: Boundary::This,
|
||||
};
|
||||
SAT_characteristics(k, y, grid, &boundaries);
|
||||
}
|
||||
|
||||
#[allow(clippy::many_single_char_names)]
|
||||
fn upwind_dissipation<UO: UpwindOperator>(
|
||||
k: (&mut Field, &mut Field),
|
||||
y: &Field,
|
||||
grid: &Grid<UO>,
|
||||
tmp: (&mut Field, &mut Field),
|
||||
) {
|
||||
let n = y.nx() * y.ny();
|
||||
let yview = y.view().into_shape((4, n)).unwrap();
|
||||
let mut tmp0 = tmp.0.view_mut().into_shape((4, n)).unwrap();
|
||||
let mut tmp1 = tmp.1.view_mut().into_shape((4, n)).unwrap();
|
||||
|
||||
for (
|
||||
((((((y, mut tmp0), mut tmp1), detj), detj_dxi_dx), detj_dxi_dy), detj_deta_dx),
|
||||
detj_deta_dy,
|
||||
) in yview
|
||||
.axis_iter(ndarray::Axis(1))
|
||||
.zip(tmp0.axis_iter_mut(ndarray::Axis(1)))
|
||||
.zip(tmp1.axis_iter_mut(ndarray::Axis(1)))
|
||||
.zip(grid.detj.iter())
|
||||
.zip(grid.detj_dxi_dx.iter())
|
||||
.zip(grid.detj_dxi_dy.iter())
|
||||
.zip(grid.detj_deta_dx.iter())
|
||||
.zip(grid.detj_deta_dy.iter())
|
||||
{
|
||||
let rho = y[0];
|
||||
assert!(rho > 0.0);
|
||||
let rhou = y[1];
|
||||
let rhov = y[2];
|
||||
let e = y[3];
|
||||
|
||||
let u = rhou / rho;
|
||||
let v = rhov / rho;
|
||||
|
||||
let uhat = detj_dxi_dx / detj * u + detj_dxi_dy / detj * v;
|
||||
let vhat = detj_deta_dx / detj * u + detj_deta_dy / detj * v;
|
||||
|
||||
let p = pressure(GAMMA, rho, rhou, rhov, e);
|
||||
assert!(p > 0.0);
|
||||
let c = (GAMMA * p / rho).sqrt();
|
||||
|
||||
let alpha_u = uhat.abs() + c;
|
||||
let alpha_v = vhat.abs() + c;
|
||||
|
||||
tmp0[0] = alpha_u * rho * detj;
|
||||
tmp1[0] = alpha_v * rho * detj;
|
||||
|
||||
tmp0[1] = alpha_u * rhou * detj;
|
||||
tmp1[1] = alpha_v * rhou * detj;
|
||||
|
||||
tmp0[2] = alpha_u * rhov * detj;
|
||||
tmp1[2] = alpha_v * rhov * detj;
|
||||
|
||||
tmp0[3] = alpha_u * e * detj;
|
||||
tmp1[3] = alpha_v * e * detj;
|
||||
}
|
||||
|
||||
UO::dissxi(tmp.0.rho(), k.0.rho_mut());
|
||||
UO::dissxi(tmp.0.rhou(), k.0.rhou_mut());
|
||||
UO::dissxi(tmp.0.rhov(), k.0.rhov_mut());
|
||||
UO::dissxi(tmp.0.e(), k.0.e_mut());
|
||||
|
||||
UO::disseta(tmp.1.rho(), k.1.rho_mut());
|
||||
UO::disseta(tmp.1.rhou(), k.1.rhou_mut());
|
||||
UO::disseta(tmp.1.rhov(), k.1.rhov_mut());
|
||||
UO::disseta(tmp.1.e(), k.1.e_mut());
|
||||
}
|
||||
|
||||
fn fluxes<SBP: SbpOperator>(k: (&mut Field, &mut Field), y: &Field, grid: &Grid<SBP>) {
|
||||
let j_dxi_dx = grid.detj_dxi_dx.view();
|
||||
let j_dxi_dy = grid.detj_dxi_dy.view();
|
||||
let j_deta_dx = grid.detj_deta_dx.view();
|
||||
let j_deta_dy = grid.detj_deta_dy.view();
|
||||
|
||||
let rho = y.rho();
|
||||
let rhou = y.rhou();
|
||||
let rhov = y.rhov();
|
||||
let e = y.e();
|
||||
|
||||
for j in 0..y.ny() {
|
||||
for i in 0..y.nx() {
|
||||
let rho = rho[(j, i)];
|
||||
assert!(rho > 0.0);
|
||||
let rhou = rhou[(j, i)];
|
||||
let rhov = rhov[(j, i)];
|
||||
let e = e[(j, i)];
|
||||
let p = pressure(GAMMA, rho, rhou, rhov, e);
|
||||
|
||||
assert!(p > 0.0);
|
||||
|
||||
let ef = [
|
||||
rhou,
|
||||
rhou * rhou / rho + p,
|
||||
rhou * rhov / rho,
|
||||
rhou * (e + p) / rho,
|
||||
];
|
||||
let ff = [
|
||||
rhov,
|
||||
rhou * rhov / rho,
|
||||
rhov * rhov / rho + p,
|
||||
rhov * (e + p) / rho,
|
||||
];
|
||||
|
||||
for comp in 0..4 {
|
||||
let eflux = j_dxi_dx[(j, i)] * ef[comp] + j_dxi_dy[(j, i)] * ff[comp];
|
||||
let fflux = j_deta_dx[(j, i)] * ef[comp] + j_deta_dy[(j, i)] * ff[comp];
|
||||
|
||||
k.0[(comp, j, i)] = eflux;
|
||||
k.1[(comp, j, i)] = fflux;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub enum Boundary {
|
||||
This,
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct BoundaryTerms {
|
||||
pub north: Boundary,
|
||||
pub south: Boundary,
|
||||
pub east: Boundary,
|
||||
pub west: Boundary,
|
||||
}
|
||||
|
||||
#[allow(non_snake_case)]
|
||||
/// Boundary conditions (SAT)
|
||||
fn SAT_characteristics<SBP: SbpOperator>(
|
||||
k: &mut Field,
|
||||
y: &Field,
|
||||
grid: &Grid<SBP>,
|
||||
_boundaries: &BoundaryTerms,
|
||||
) {
|
||||
/* // Whean using infinite boundaries, use this...
|
||||
let steady_v = [1.0, 1.0, 0.0, {
|
||||
let M = 0.1;
|
||||
let p_inf = 1.0 / (GAMMA * M * M);
|
||||
p_inf / (GAMMA - 1.0) + 0.5
|
||||
}];
|
||||
let steady_a = ndarray::Array1::from(steady_v.to_vec());
|
||||
let steady = steady_a.broadcast((k.nx(), 4)).unwrap().reversed_axes();
|
||||
assert_eq!(steady.shape(), [4, k.nx()]);
|
||||
*/
|
||||
// North boundary
|
||||
{
|
||||
let hi = (k.ny() - 1) as f32 * SBP::h()[0];
|
||||
let sign = -1.0;
|
||||
let tau = 1.0;
|
||||
let slice = s![y.ny() - 1, ..];
|
||||
SAT_characteristic(
|
||||
k.north_mut(),
|
||||
y.north(),
|
||||
y.south(), // Self South
|
||||
//steady.view(),
|
||||
hi,
|
||||
sign,
|
||||
tau,
|
||||
grid.detj.slice(slice),
|
||||
grid.detj_deta_dx.slice(slice),
|
||||
grid.detj_deta_dy.slice(slice),
|
||||
);
|
||||
}
|
||||
// South boundary
|
||||
{
|
||||
let hi = (k.ny() - 1) as f32 * SBP::h()[0];
|
||||
let sign = 1.0;
|
||||
let tau = -1.0;
|
||||
let slice = s![0, ..];
|
||||
SAT_characteristic(
|
||||
k.south_mut(),
|
||||
y.south(),
|
||||
y.north(), // Self North
|
||||
//steady.view(),
|
||||
hi,
|
||||
sign,
|
||||
tau,
|
||||
grid.detj.slice(slice),
|
||||
grid.detj_deta_dx.slice(slice),
|
||||
grid.detj_deta_dy.slice(slice),
|
||||
);
|
||||
}
|
||||
/*let steady = ndarray::Array2::from_shape_fn((4, k.ny()), |(k, _)| match k {
|
||||
0 => 1.0,
|
||||
1 => 1.0,
|
||||
2 => 0.0,
|
||||
3 => {
|
||||
let M = 0.1;
|
||||
let p_inf = 1.0 / (GAMMA * M * M);
|
||||
p_inf / (GAMMA - 1.0) + 0.5
|
||||
}
|
||||
_ => unreachable!(),
|
||||
});*/
|
||||
// West Boundary
|
||||
{
|
||||
let hi = (k.nx() - 1) as f32 * SBP::h()[0];
|
||||
let sign = 1.0;
|
||||
let tau = -1.0;
|
||||
let slice = s![.., 0];
|
||||
SAT_characteristic(
|
||||
k.west_mut(),
|
||||
y.west(),
|
||||
y.east(), // Self East
|
||||
//steady.view(),
|
||||
hi,
|
||||
sign,
|
||||
tau,
|
||||
grid.detj.slice(slice),
|
||||
grid.detj_dxi_dx.slice(slice),
|
||||
grid.detj_dxi_dy.slice(slice),
|
||||
);
|
||||
}
|
||||
// East Boundary
|
||||
{
|
||||
let hi = (k.nx() - 1) as f32 * SBP::h()[0];
|
||||
let sign = -1.0;
|
||||
let tau = 1.0;
|
||||
let slice = s![.., y.nx() - 1];
|
||||
SAT_characteristic(
|
||||
k.east_mut(),
|
||||
y.east(),
|
||||
y.west(), // Self West
|
||||
//steady.view(),
|
||||
hi,
|
||||
sign,
|
||||
tau,
|
||||
grid.detj.slice(slice),
|
||||
grid.detj_dxi_dx.slice(slice),
|
||||
grid.detj_dxi_dy.slice(slice),
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[allow(non_snake_case)]
|
||||
#[allow(clippy::many_single_char_names)]
|
||||
#[allow(clippy::too_many_arguments)]
|
||||
/// Boundary conditions (SAT)
|
||||
fn SAT_characteristic(
|
||||
mut k: ArrayViewMut2<f32>,
|
||||
y: ArrayView2<f32>,
|
||||
z: ArrayView2<f32>, // Size 4 x n (all components in line)
|
||||
hi: f32,
|
||||
sign: f32,
|
||||
tau: f32,
|
||||
detj: ArrayView1<f32>,
|
||||
detj_d_dx: ArrayView1<f32>,
|
||||
detj_d_dy: ArrayView1<f32>,
|
||||
) {
|
||||
assert_eq!(detj.shape(), detj_d_dx.shape());
|
||||
assert_eq!(detj.shape(), detj_d_dy.shape());
|
||||
assert_eq!(y.shape(), z.shape());
|
||||
assert_eq!(y.shape()[0], 4);
|
||||
assert_eq!(y.shape()[1], detj.shape()[0]);
|
||||
|
||||
for (((((mut k, y), z), detj), detj_d_dx), detj_d_dy) in k
|
||||
.axis_iter_mut(ndarray::Axis(1))
|
||||
.zip(y.axis_iter(ndarray::Axis(1)))
|
||||
.zip(z.axis_iter(ndarray::Axis(1)))
|
||||
.zip(detj.iter())
|
||||
.zip(detj_d_dx.iter())
|
||||
.zip(detj_d_dy.iter())
|
||||
{
|
||||
let rho = y[0];
|
||||
let rhou = y[1];
|
||||
let rhov = y[2];
|
||||
let e = y[3];
|
||||
|
||||
let kx_ = detj_d_dx / detj;
|
||||
let ky_ = detj_d_dy / detj;
|
||||
|
||||
let (kx, ky) = {
|
||||
let r = f32::hypot(kx_, ky_);
|
||||
(kx_ / r, ky_ / r)
|
||||
};
|
||||
|
||||
let u = rhou / rho;
|
||||
let v = rhov / rho;
|
||||
|
||||
let theta = kx * u + ky * v;
|
||||
|
||||
let p = pressure(GAMMA, rho, rhou, rhov, e);
|
||||
let c = (GAMMA * p / rho).sqrt();
|
||||
let phi2 = (GAMMA - 1.0) * (u * u + v * v) / 2.0;
|
||||
|
||||
let phi2_c2 = (phi2 + c * c) / (GAMMA - 1.0);
|
||||
|
||||
let T = [
|
||||
[1.0, 0.0, 1.0, 1.0],
|
||||
[u, ky, u + kx * c, u - kx * c],
|
||||
[v, -kx, v + ky * c, v - ky * c],
|
||||
[
|
||||
phi2 / (GAMMA - 1.0),
|
||||
ky * u - kx * v,
|
||||
phi2_c2 + c * theta,
|
||||
phi2_c2 - c * theta,
|
||||
],
|
||||
];
|
||||
let U = kx_ * u + ky_ * v;
|
||||
let L = [
|
||||
U,
|
||||
U,
|
||||
U + c * f32::hypot(kx_, ky_),
|
||||
U - c * f32::hypot(kx_, ky_),
|
||||
];
|
||||
let beta = 1.0 / (2.0 * c * c);
|
||||
let TI = [
|
||||
[
|
||||
1.0 - phi2 / (c * c),
|
||||
(GAMMA - 1.0) * u / (c * c),
|
||||
(GAMMA - 1.0) * v / (c * c),
|
||||
-(GAMMA - 1.0) / (c * c),
|
||||
],
|
||||
[-(ky * u - kx * v), ky, -kx, 0.0],
|
||||
[
|
||||
beta * (phi2 - c * theta),
|
||||
beta * (kx * c - (GAMMA - 1.0) * u),
|
||||
beta * (ky * c - (GAMMA - 1.0) * v),
|
||||
beta * (GAMMA - 1.0),
|
||||
],
|
||||
[
|
||||
beta * (phi2 + c * theta),
|
||||
-beta * (kx * c + (GAMMA - 1.0) * u),
|
||||
-beta * (ky * c + (GAMMA - 1.0) * v),
|
||||
beta * (GAMMA - 1.0),
|
||||
],
|
||||
];
|
||||
|
||||
let res = [rho - z[0], rhou - z[1], rhov - z[2], e - z[3]];
|
||||
let mut TIres = [0.0; 4];
|
||||
#[allow(clippy::needless_range_loop)]
|
||||
for row in 0..4 {
|
||||
for col in 0..4 {
|
||||
TIres[row] += TI[row][col] * res[col];
|
||||
}
|
||||
}
|
||||
|
||||
// L + sign(abs(L)) * TIres
|
||||
let mut LTIres = [0.0; 4];
|
||||
for row in 0..4 {
|
||||
LTIres[row] = (L[row] + sign * L[row].abs()) * TIres[row];
|
||||
}
|
||||
|
||||
// T*LTIres
|
||||
let mut TLTIres = [0.0; 4];
|
||||
#[allow(clippy::needless_range_loop)]
|
||||
for row in 0..4 {
|
||||
for col in 0..4 {
|
||||
TLTIres[row] += T[row][col] * LTIres[col];
|
||||
}
|
||||
}
|
||||
|
||||
for comp in 0..4 {
|
||||
k[comp] += hi * tau * TLTIres[comp];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub struct WorkBuffers {
|
||||
k: [Field; 4],
|
||||
tmp: (Field, Field, Field, Field, Field, Field),
|
||||
}
|
||||
|
||||
impl WorkBuffers {
|
||||
pub fn new(nx: usize, ny: usize) -> Self {
|
||||
let arr3 = Field::new(nx, ny);
|
||||
Self {
|
||||
k: [arr3.clone(), arr3.clone(), arr3.clone(), arr3.clone()],
|
||||
tmp: (
|
||||
arr3.clone(),
|
||||
arr3.clone(),
|
||||
arr3.clone(),
|
||||
arr3.clone(),
|
||||
arr3.clone(),
|
||||
arr3,
|
||||
),
|
||||
}
|
||||
}
|
||||
}
|
||||
75
src/grid.rs
75
src/grid.rs
@@ -1,75 +0,0 @@
|
||||
use ndarray::Array2;
|
||||
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct Grid<SBP>
|
||||
where
|
||||
SBP: super::operators::SbpOperator,
|
||||
{
|
||||
pub(crate) x: Array2<f32>,
|
||||
pub(crate) y: Array2<f32>,
|
||||
|
||||
pub(crate) detj: Array2<f32>,
|
||||
pub(crate) detj_dxi_dx: Array2<f32>,
|
||||
pub(crate) detj_dxi_dy: Array2<f32>,
|
||||
pub(crate) detj_deta_dx: Array2<f32>,
|
||||
pub(crate) detj_deta_dy: Array2<f32>,
|
||||
|
||||
operator: std::marker::PhantomData<SBP>,
|
||||
}
|
||||
|
||||
impl<SBP: super::operators::SbpOperator> Grid<SBP> {
|
||||
pub fn new(x: Array2<f32>, y: Array2<f32>) -> Result<Self, ndarray::ShapeError> {
|
||||
assert_eq!(x.shape(), y.shape());
|
||||
let ny = x.shape()[0];
|
||||
let nx = x.shape()[1];
|
||||
|
||||
let mut dx_dxi = Array2::zeros((ny, nx));
|
||||
SBP::diffxi(x.view(), dx_dxi.view_mut());
|
||||
let mut dx_deta = Array2::zeros((ny, nx));
|
||||
SBP::diffeta(x.view(), dx_deta.view_mut());
|
||||
let mut dy_dxi = Array2::zeros((ny, nx));
|
||||
SBP::diffxi(y.view(), dy_dxi.view_mut());
|
||||
let mut dy_deta = Array2::zeros((ny, nx));
|
||||
SBP::diffeta(y.view(), dy_deta.view_mut());
|
||||
|
||||
let mut detj = Array2::zeros((ny, nx));
|
||||
ndarray::azip!((detj in &mut detj,
|
||||
&dx_dxi in &dx_dxi,
|
||||
&dx_deta in &dx_deta,
|
||||
&dy_dxi in &dy_dxi,
|
||||
&dy_deta in &dy_deta) {
|
||||
*detj = dx_dxi * dy_deta - dx_deta * dy_dxi;
|
||||
assert!(*detj > 0.0);
|
||||
});
|
||||
|
||||
let detj_dxi_dx = dy_deta;
|
||||
let detj_dxi_dy = {
|
||||
let mut dx_deta = dx_deta;
|
||||
dx_deta.mapv_inplace(|v| -v);
|
||||
dx_deta
|
||||
};
|
||||
let detj_deta_dx = {
|
||||
let mut dy_dxi = dy_dxi;
|
||||
dy_dxi.mapv_inplace(|v| -v);
|
||||
dy_dxi
|
||||
};
|
||||
let detj_deta_dy = dx_dxi;
|
||||
|
||||
Ok(Self {
|
||||
x,
|
||||
y,
|
||||
detj,
|
||||
detj_dxi_dx,
|
||||
detj_dxi_dy,
|
||||
detj_deta_dx,
|
||||
detj_deta_dy,
|
||||
operator: std::marker::PhantomData,
|
||||
})
|
||||
}
|
||||
pub fn nx(&self) -> usize {
|
||||
self.x.shape()[1]
|
||||
}
|
||||
pub fn ny(&self) -> usize {
|
||||
self.x.shape()[0]
|
||||
}
|
||||
}
|
||||
@@ -1,46 +0,0 @@
|
||||
use super::grid::Grid;
|
||||
use super::operators::SbpOperator;
|
||||
use ndarray::{Array3, Zip};
|
||||
|
||||
pub(crate) fn rk4<'a, F: 'a, SBP, RHS, WB>(
|
||||
rhs: RHS,
|
||||
prev: &F,
|
||||
fut: &mut F,
|
||||
dt: f32,
|
||||
grid: &Grid<SBP>,
|
||||
k: &mut [F; 4],
|
||||
mut wb: &mut WB,
|
||||
) where
|
||||
F: std::ops::Deref<Target = Array3<f32>> + std::ops::DerefMut<Target = Array3<f32>>,
|
||||
SBP: SbpOperator,
|
||||
RHS: Fn(&mut F, &F, &Grid<SBP>, &mut WB),
|
||||
{
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
|
||||
for i in 0..4 {
|
||||
// y = y0 + c*kn
|
||||
fut.assign(prev);
|
||||
match i {
|
||||
0 => {}
|
||||
1 | 2 => {
|
||||
fut.scaled_add(1.0 / 2.0 * dt, &k[i - 1]);
|
||||
}
|
||||
3 => {
|
||||
fut.scaled_add(dt, &k[i - 1]);
|
||||
}
|
||||
_ => {
|
||||
unreachable!();
|
||||
}
|
||||
};
|
||||
|
||||
rhs(&mut k[i], &fut, grid, &mut wb);
|
||||
}
|
||||
|
||||
Zip::from(&mut **fut)
|
||||
.and(&**prev)
|
||||
.and(&*k[0])
|
||||
.and(&*k[1])
|
||||
.and(&*k[2])
|
||||
.and(&*k[3])
|
||||
.apply(|y1, &y0, &k1, &k2, &k3, &k4| *y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4));
|
||||
}
|
||||
125
src/lib.rs
125
src/lib.rs
@@ -1,125 +0,0 @@
|
||||
use wasm_bindgen::prelude::*;
|
||||
|
||||
pub mod euler;
|
||||
mod grid;
|
||||
pub(crate) mod integrate;
|
||||
pub mod maxwell;
|
||||
pub mod operators;
|
||||
pub(crate) use grid::Grid;
|
||||
|
||||
#[cfg(feature = "wee_alloc")]
|
||||
#[global_allocator]
|
||||
static ALLOC: wee_alloc::WeeAlloc = wee_alloc::WeeAlloc::INIT;
|
||||
|
||||
#[wasm_bindgen]
|
||||
pub fn set_panic_hook() {
|
||||
#[cfg(feature = "console_error_panic_hook")]
|
||||
console_error_panic_hook::set_once();
|
||||
}
|
||||
|
||||
#[wasm_bindgen]
|
||||
pub struct MaxwellUniverse(maxwell::System<operators::Upwind4>);
|
||||
|
||||
#[wasm_bindgen]
|
||||
impl MaxwellUniverse {
|
||||
#[wasm_bindgen(constructor)]
|
||||
pub fn new(height: usize, width: usize, x: &[f32], y: &[f32]) -> Self {
|
||||
let x = ndarray::Array2::from_shape_vec((height, width), x.to_vec()).unwrap();
|
||||
let y = ndarray::Array2::from_shape_vec((height, width), y.to_vec()).unwrap();
|
||||
Self(maxwell::System::new(x, y))
|
||||
}
|
||||
|
||||
pub fn init(&mut self, x0: f32, y0: f32) {
|
||||
self.0.set_gaussian(x0, y0);
|
||||
}
|
||||
|
||||
pub fn advance(&mut self, dt: f32) {
|
||||
self.0.advance(dt)
|
||||
}
|
||||
|
||||
pub fn advance_upwind(&mut self, dt: f32) {
|
||||
self.0.advance_upwind(dt)
|
||||
}
|
||||
|
||||
pub fn get_ex_ptr(&self) -> *const u8 {
|
||||
self.0.field().ex().as_ptr() as *const u8
|
||||
}
|
||||
|
||||
pub fn get_ey_ptr(&self) -> *const u8 {
|
||||
self.0.field().ey().as_ptr() as *const u8
|
||||
}
|
||||
|
||||
pub fn get_hz_ptr(&self) -> *const u8 {
|
||||
self.0.field().hz().as_ptr() as *const u8
|
||||
}
|
||||
}
|
||||
|
||||
#[wasm_bindgen]
|
||||
pub struct EulerUniverse(euler::System<operators::Upwind4>);
|
||||
|
||||
impl EulerUniverse {
|
||||
pub fn new(x: ndarray::Array2<f32>, y: ndarray::Array2<f32>) -> Self {
|
||||
Self(euler::System::new(x, y))
|
||||
}
|
||||
}
|
||||
|
||||
#[wasm_bindgen]
|
||||
impl EulerUniverse {
|
||||
#[wasm_bindgen(constructor)]
|
||||
pub fn new_with_slice(height: usize, width: usize, x: &[f32], y: &[f32]) -> Self {
|
||||
let x = ndarray::Array2::from_shape_vec((height, width), x.to_vec()).unwrap();
|
||||
let y = ndarray::Array2::from_shape_vec((height, width), y.to_vec()).unwrap();
|
||||
Self(euler::System::new(x, y))
|
||||
}
|
||||
|
||||
pub fn init(&mut self, x0: f32, y0: f32) {
|
||||
self.0.init_with_vortex(x0, y0)
|
||||
}
|
||||
|
||||
pub fn advance(&mut self, dt: f32) {
|
||||
self.0.advance(dt)
|
||||
}
|
||||
|
||||
pub fn advance_upwind(&mut self, dt: f32) {
|
||||
self.0.advance_upwind(dt)
|
||||
}
|
||||
|
||||
pub fn get_rho_ptr(&self) -> *const u8 {
|
||||
self.0.field().rho().as_ptr() as *const u8
|
||||
}
|
||||
pub fn get_rhou_ptr(&self) -> *const u8 {
|
||||
self.0.field().rhou().as_ptr() as *const u8
|
||||
}
|
||||
pub fn get_rhov_ptr(&self) -> *const u8 {
|
||||
self.0.field().rhov().as_ptr() as *const u8
|
||||
}
|
||||
pub fn get_e_ptr(&self) -> *const u8 {
|
||||
self.0.field().e().as_ptr() as *const u8
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn start_and_advance_euler() {
|
||||
let x = ndarray::Array2::from_shape_fn((20, 20), |(_j, i)| {
|
||||
5.0 * 2.0 * ((i as f32 / (20 - 1) as f32) - 0.5)
|
||||
});
|
||||
let y = ndarray::Array2::from_shape_fn((20, 20), |(j, _i)| {
|
||||
5.0 * 2.0 * ((j as f32 / (20 - 1) as f32) - 0.5)
|
||||
});
|
||||
let mut universe = EulerUniverse::new(x, y);
|
||||
universe.init(-1.0, 0.0);
|
||||
for _ in 0..50 {
|
||||
universe.advance(0.01);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn start_and_advance_upwind_euler() {
|
||||
let x = ndarray::Array2::from_shape_fn((20, 10), |(_j, i)| i as f32 / (10 - 1) as f32);
|
||||
let y = ndarray::Array2::from_shape_fn((20, 10), |(j, _i)| j as f32 / (20 - 1) as f32);
|
||||
let mut universe = EulerUniverse::new(x, y);
|
||||
universe.init(0.5, 0.5);
|
||||
for _ in 0..50 {
|
||||
universe.advance_upwind(0.01);
|
||||
}
|
||||
}
|
||||
|
||||
575
src/maxwell.rs
575
src/maxwell.rs
@@ -1,575 +0,0 @@
|
||||
use super::integrate;
|
||||
use super::operators::{SbpOperator, UpwindOperator};
|
||||
use super::Grid;
|
||||
use ndarray::azip;
|
||||
use ndarray::prelude::*;
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct Field(pub(crate) Array3<f32>);
|
||||
|
||||
impl std::ops::Deref for Field {
|
||||
type Target = Array3<f32>;
|
||||
fn deref(&self) -> &Self::Target {
|
||||
&self.0
|
||||
}
|
||||
}
|
||||
|
||||
impl std::ops::DerefMut for Field {
|
||||
fn deref_mut(&mut self) -> &mut Self::Target {
|
||||
&mut self.0
|
||||
}
|
||||
}
|
||||
|
||||
impl Field {
|
||||
pub fn new(height: usize, width: usize) -> Self {
|
||||
let field = Array3::zeros((3, height, width));
|
||||
|
||||
Self(field)
|
||||
}
|
||||
|
||||
pub fn nx(&self) -> usize {
|
||||
self.0.shape()[2]
|
||||
}
|
||||
pub fn ny(&self) -> usize {
|
||||
self.0.shape()[1]
|
||||
}
|
||||
|
||||
pub fn ex(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![0, .., ..])
|
||||
}
|
||||
pub fn hz(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![1, .., ..])
|
||||
}
|
||||
pub fn ey(&self) -> ArrayView2<f32> {
|
||||
self.slice(s![2, .., ..])
|
||||
}
|
||||
|
||||
pub fn ex_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![0, .., ..])
|
||||
}
|
||||
pub fn hz_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![1, .., ..])
|
||||
}
|
||||
pub fn ey_mut(&mut self) -> ArrayViewMut2<f32> {
|
||||
self.slice_mut(s![2, .., ..])
|
||||
}
|
||||
|
||||
pub fn components_mut(
|
||||
&mut self,
|
||||
) -> (ArrayViewMut2<f32>, ArrayViewMut2<f32>, ArrayViewMut2<f32>) {
|
||||
let mut iter = self.0.outer_iter_mut();
|
||||
|
||||
let ex = iter.next().unwrap();
|
||||
let hz = iter.next().unwrap();
|
||||
let ey = iter.next().unwrap();
|
||||
assert_eq!(iter.next(), None);
|
||||
|
||||
(ex, hz, ey)
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct System<SBP: SbpOperator> {
|
||||
sys: (Field, Field),
|
||||
wb: WorkBuffers,
|
||||
grid: Grid<SBP>,
|
||||
}
|
||||
|
||||
impl<SBP: SbpOperator> System<SBP> {
|
||||
pub fn new(x: Array2<f32>, y: Array2<f32>) -> Self {
|
||||
assert_eq!(x.shape(), y.shape());
|
||||
let ny = x.shape()[0];
|
||||
let nx = x.shape()[1];
|
||||
|
||||
let grid = Grid::new(x, y).unwrap();
|
||||
|
||||
Self {
|
||||
sys: (Field::new(ny, nx), Field::new(ny, nx)),
|
||||
grid,
|
||||
wb: WorkBuffers::new(ny, nx),
|
||||
}
|
||||
}
|
||||
|
||||
pub fn field(&self) -> &Field {
|
||||
&self.sys.0
|
||||
}
|
||||
|
||||
pub fn set_gaussian(&mut self, x0: f32, y0: f32) {
|
||||
let (ex, hz, ey) = self.sys.0.components_mut();
|
||||
ndarray::azip!(
|
||||
(ex in ex, hz in hz, ey in ey,
|
||||
&x in &self.grid.x, &y in &self.grid.y)
|
||||
{
|
||||
*ex = 0.0;
|
||||
*ey = 0.0;
|
||||
*hz = gaussian(x, x0, y, y0)/32.0;
|
||||
});
|
||||
}
|
||||
|
||||
pub fn advance(&mut self, dt: f32) {
|
||||
integrate::rk4(
|
||||
RHS,
|
||||
&self.sys.0,
|
||||
&mut self.sys.1,
|
||||
dt,
|
||||
&self.grid,
|
||||
&mut self.wb.k,
|
||||
&mut self.wb.tmp,
|
||||
);
|
||||
std::mem::swap(&mut self.sys.0, &mut self.sys.1);
|
||||
}
|
||||
}
|
||||
|
||||
impl<UO: UpwindOperator> System<UO> {
|
||||
/// Using artificial dissipation with the upwind operator
|
||||
pub fn advance_upwind(&mut self, dt: f32) {
|
||||
integrate::rk4(
|
||||
RHS_upwind,
|
||||
&self.sys.0,
|
||||
&mut self.sys.1,
|
||||
dt,
|
||||
&self.grid,
|
||||
&mut self.wb.k,
|
||||
&mut self.wb.tmp,
|
||||
);
|
||||
std::mem::swap(&mut self.sys.0, &mut self.sys.1);
|
||||
}
|
||||
}
|
||||
|
||||
fn gaussian(x: f32, x0: f32, y: f32, y0: f32) -> f32 {
|
||||
use std::f32;
|
||||
let x = x - x0;
|
||||
let y = y - y0;
|
||||
|
||||
let sigma = 0.05;
|
||||
|
||||
1.0 / (2.0 * f32::consts::PI * sigma * sigma) * (-(x * x + y * y) / (2.0 * sigma * sigma)).exp()
|
||||
}
|
||||
|
||||
#[allow(non_snake_case)]
|
||||
/// Solving (Au)_x + (Bu)_y
|
||||
/// with:
|
||||
/// A B
|
||||
/// [ 0, 0, 0] [ 0, 1, 0]
|
||||
/// [ 0, 0, -1] [ 1, 0, 0]
|
||||
/// [ 0, -1, 0] [ 0, 0, 0]
|
||||
///
|
||||
/// This flux is rotated by the grid metrics
|
||||
/// (Au)_x + (Bu)_y = 1/J [
|
||||
/// (J xi_x Au)_xi + (J eta_x Au)_eta
|
||||
/// (J xi_y Bu)_xi + (J eta_y Bu)_eta
|
||||
/// ]
|
||||
/// where J is the grid determinant
|
||||
///
|
||||
/// This is used both in fluxes and SAT terms
|
||||
fn RHS<SBP: SbpOperator>(
|
||||
k: &mut Field,
|
||||
y: &Field,
|
||||
grid: &Grid<SBP>,
|
||||
tmp: &mut (Array2<f32>, Array2<f32>, Array2<f32>, Array2<f32>),
|
||||
) {
|
||||
fluxes(k, y, grid, tmp);
|
||||
|
||||
let boundaries = BoundaryTerms {
|
||||
north: Boundary::This,
|
||||
south: Boundary::This,
|
||||
west: Boundary::This,
|
||||
east: Boundary::This,
|
||||
};
|
||||
SAT_characteristics(k, y, grid, &boundaries);
|
||||
|
||||
azip!((k in &mut k.0,
|
||||
&detj in &grid.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
|
||||
*k /= detj;
|
||||
});
|
||||
}
|
||||
|
||||
#[allow(non_snake_case)]
|
||||
fn RHS_upwind<UO: UpwindOperator>(
|
||||
k: &mut Field,
|
||||
y: &Field,
|
||||
grid: &Grid<UO>,
|
||||
tmp: &mut (Array2<f32>, Array2<f32>, Array2<f32>, Array2<f32>),
|
||||
) {
|
||||
fluxes(k, y, grid, tmp);
|
||||
dissipation(k, y, grid, tmp);
|
||||
|
||||
let boundaries = BoundaryTerms {
|
||||
north: Boundary::This,
|
||||
south: Boundary::This,
|
||||
west: Boundary::This,
|
||||
east: Boundary::This,
|
||||
};
|
||||
SAT_characteristics(k, y, grid, &boundaries);
|
||||
|
||||
azip!((k in &mut k.0,
|
||||
&detj in &grid.detj.broadcast((3, y.ny(), y.nx())).unwrap()) {
|
||||
*k /= detj;
|
||||
});
|
||||
}
|
||||
|
||||
fn fluxes<SBP: SbpOperator>(
|
||||
k: &mut Field,
|
||||
y: &Field,
|
||||
grid: &Grid<SBP>,
|
||||
tmp: &mut (Array2<f32>, Array2<f32>, Array2<f32>, Array2<f32>),
|
||||
) {
|
||||
// ex = hz_y
|
||||
{
|
||||
ndarray::azip!((a in &mut tmp.0,
|
||||
&dxi_dy in &grid.detj_dxi_dy,
|
||||
&hz in &y.hz())
|
||||
*a = dxi_dy * hz
|
||||
);
|
||||
SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
|
||||
|
||||
ndarray::azip!((b in &mut tmp.2,
|
||||
&deta_dy in &grid.detj_deta_dy,
|
||||
&hz in &y.hz())
|
||||
*b = deta_dy * hz
|
||||
);
|
||||
SBP::diffeta(tmp.2.view(), tmp.3.view_mut());
|
||||
|
||||
ndarray::azip!((flux in &mut k.ex_mut(), &ax in &tmp.1, &by in &tmp.3)
|
||||
*flux = ax + by
|
||||
);
|
||||
}
|
||||
|
||||
{
|
||||
// hz = -ey_x + ex_y
|
||||
ndarray::azip!((a in &mut tmp.0,
|
||||
&dxi_dx in &grid.detj_dxi_dx,
|
||||
&dxi_dy in &grid.detj_dxi_dy,
|
||||
&ex in &y.ex(),
|
||||
&ey in &y.ey())
|
||||
*a = dxi_dx * -ey + dxi_dy * ex
|
||||
);
|
||||
SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
|
||||
|
||||
ndarray::azip!((b in &mut tmp.2,
|
||||
&deta_dx in &grid.detj_deta_dx,
|
||||
&deta_dy in &grid.detj_deta_dy,
|
||||
&ex in &y.ex(),
|
||||
&ey in &y.ey())
|
||||
*b = deta_dx * -ey + deta_dy * ex
|
||||
);
|
||||
SBP::diffeta(tmp.2.view(), tmp.3.view_mut());
|
||||
|
||||
ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3)
|
||||
*flux = ax + by
|
||||
);
|
||||
}
|
||||
|
||||
// ey = -hz_x
|
||||
{
|
||||
ndarray::azip!((a in &mut tmp.0,
|
||||
&dxi_dx in &grid.detj_dxi_dx,
|
||||
&hz in &y.hz())
|
||||
*a = dxi_dx * -hz
|
||||
);
|
||||
SBP::diffxi(tmp.0.view(), tmp.1.view_mut());
|
||||
|
||||
azip!((b in &mut tmp.2,
|
||||
&deta_dx in &grid.detj_deta_dx,
|
||||
&hz in &y.hz())
|
||||
*b = deta_dx * -hz
|
||||
);
|
||||
SBP::diffeta(tmp.2.view(), tmp.3.view_mut());
|
||||
|
||||
azip!((flux in &mut k.ey_mut(), &ax in &tmp.1, &by in &tmp.3)
|
||||
*flux = ax + by
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
fn dissipation<UO: UpwindOperator>(
|
||||
k: &mut Field,
|
||||
y: &Field,
|
||||
grid: &Grid<UO>,
|
||||
tmp: &mut (Array2<f32>, Array2<f32>, Array2<f32>, Array2<f32>),
|
||||
) {
|
||||
// ex component
|
||||
{
|
||||
ndarray::azip!((a in &mut tmp.0,
|
||||
&kx in &grid.detj_dxi_dx,
|
||||
&ky in &grid.detj_dxi_dy,
|
||||
&ex in &y.ex(),
|
||||
&ey in &y.ey()) {
|
||||
let r = f32::hypot(kx, ky);
|
||||
*a = ky*ky/r * ex + -kx*ky/r*ey;
|
||||
});
|
||||
UO::dissxi(tmp.0.view(), tmp.1.view_mut());
|
||||
|
||||
ndarray::azip!((b in &mut tmp.2,
|
||||
&kx in &grid.detj_deta_dx,
|
||||
&ky in &grid.detj_deta_dy,
|
||||
&ex in &y.ex(),
|
||||
&ey in &y.ey()) {
|
||||
let r = f32::hypot(kx, ky);
|
||||
*b = ky*ky/r * ex + -kx*ky/r*ey;
|
||||
});
|
||||
UO::disseta(tmp.2.view(), tmp.3.view_mut());
|
||||
|
||||
ndarray::azip!((flux in &mut k.ex_mut(), &ax in &tmp.1, &by in &tmp.3)
|
||||
*flux += ax + by
|
||||
);
|
||||
}
|
||||
|
||||
// hz component
|
||||
{
|
||||
ndarray::azip!((a in &mut tmp.0,
|
||||
&kx in &grid.detj_dxi_dx,
|
||||
&ky in &grid.detj_dxi_dy,
|
||||
&hz in &y.hz()) {
|
||||
let r = f32::hypot(kx, ky);
|
||||
*a = r * hz;
|
||||
});
|
||||
UO::dissxi(tmp.0.view(), tmp.1.view_mut());
|
||||
|
||||
ndarray::azip!((b in &mut tmp.2,
|
||||
&kx in &grid.detj_deta_dx,
|
||||
&ky in &grid.detj_deta_dy,
|
||||
&hz in &y.hz()) {
|
||||
let r = f32::hypot(kx, ky);
|
||||
*b = r * hz;
|
||||
});
|
||||
UO::disseta(tmp.2.view(), tmp.3.view_mut());
|
||||
|
||||
ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3)
|
||||
*flux += ax + by
|
||||
);
|
||||
}
|
||||
|
||||
// ey
|
||||
{
|
||||
ndarray::azip!((a in &mut tmp.0,
|
||||
&kx in &grid.detj_dxi_dx,
|
||||
&ky in &grid.detj_dxi_dy,
|
||||
&ex in &y.ex(),
|
||||
&ey in &y.ey()) {
|
||||
let r = f32::hypot(kx, ky);
|
||||
*a = -kx*ky/r * ex + kx*kx/r*ey;
|
||||
});
|
||||
UO::dissxi(tmp.0.view(), tmp.1.view_mut());
|
||||
|
||||
ndarray::azip!((b in &mut tmp.2,
|
||||
&kx in &grid.detj_deta_dx,
|
||||
&ky in &grid.detj_deta_dy,
|
||||
&ex in &y.ex(),
|
||||
&ey in &y.ey()) {
|
||||
let r = f32::hypot(kx, ky);
|
||||
*b = -kx*ky/r * ex + kx*kx/r*ey;
|
||||
});
|
||||
UO::disseta(tmp.2.view(), tmp.3.view_mut());
|
||||
|
||||
ndarray::azip!((flux in &mut k.hz_mut(), &ax in &tmp.1, &by in &tmp.3)
|
||||
*flux += ax + by
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub enum Boundary {
|
||||
This,
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct BoundaryTerms {
|
||||
pub north: Boundary,
|
||||
pub south: Boundary,
|
||||
pub east: Boundary,
|
||||
pub west: Boundary,
|
||||
}
|
||||
|
||||
#[allow(non_snake_case)]
|
||||
/// Boundary conditions (SAT)
|
||||
fn SAT_characteristics<SBP: SbpOperator>(
|
||||
k: &mut Field,
|
||||
y: &Field,
|
||||
grid: &Grid<SBP>,
|
||||
boundaries: &BoundaryTerms,
|
||||
) {
|
||||
let ny = y.ny();
|
||||
let nx = y.nx();
|
||||
|
||||
fn positive_flux(kx: f32, ky: f32) -> [[f32; 3]; 3] {
|
||||
let r = (kx * kx + ky * ky).sqrt();
|
||||
[
|
||||
[ky * ky / r / 2.0, ky / 2.0, -kx * ky / r / 2.0],
|
||||
[ky / 2.0, r / 2.0, -kx / 2.0],
|
||||
[-kx * ky / r / 2.0, -kx / 2.0, kx * kx / r / 2.0],
|
||||
]
|
||||
}
|
||||
fn negative_flux(kx: f32, ky: f32) -> [[f32; 3]; 3] {
|
||||
let r = (kx * kx + ky * ky).sqrt();
|
||||
[
|
||||
[-ky * ky / r / 2.0, ky / 2.0, kx * ky / r / 2.0],
|
||||
[ky / 2.0, -r / 2.0, -kx / 2.0],
|
||||
[kx * ky / r / 2.0, -kx / 2.0, -kx * kx / r / 2.0],
|
||||
]
|
||||
}
|
||||
|
||||
{
|
||||
let g = match boundaries.east {
|
||||
Boundary::This => y.slice(s![.., .., 0]),
|
||||
};
|
||||
// East boundary
|
||||
let hinv = 1.0 / (SBP::h()[0] / (nx - 1) as f32);
|
||||
for ((((mut k, v), g), &kx), &ky) in k
|
||||
.slice_mut(s![.., .., nx - 1])
|
||||
.gencolumns_mut()
|
||||
.into_iter()
|
||||
.zip(y.slice(s![.., .., nx - 1]).gencolumns())
|
||||
.zip(g.gencolumns())
|
||||
.zip(grid.detj_dxi_dx.slice(s![.., nx - 1]))
|
||||
.zip(grid.detj_dxi_dy.slice(s![.., nx - 1]))
|
||||
{
|
||||
// East boundary, positive flux
|
||||
let tau = -1.0;
|
||||
|
||||
let v = (v[0], v[1], v[2]);
|
||||
let g = (g[0], g[1], g[2]);
|
||||
|
||||
let plus = positive_flux(kx, ky);
|
||||
|
||||
k[0] += tau
|
||||
* hinv
|
||||
* (plus[0][0] * (v.0 - g.0) + plus[0][1] * (v.1 - g.1) + plus[0][2] * (v.2 - g.2));
|
||||
k[1] += tau
|
||||
* hinv
|
||||
* (plus[1][0] * (v.0 - g.0) + plus[1][1] * (v.1 - g.1) + plus[1][2] * (v.2 - g.2));
|
||||
k[2] += tau
|
||||
* hinv
|
||||
* (plus[2][0] * (v.0 - g.0) + plus[2][1] * (v.1 - g.1) + plus[2][2] * (v.2 - g.2));
|
||||
}
|
||||
}
|
||||
{
|
||||
// West boundary, negative flux
|
||||
let g = match boundaries.east {
|
||||
Boundary::This => y.slice(s![.., .., nx - 1]),
|
||||
};
|
||||
let hinv = 1.0 / (SBP::h()[0] / (nx - 1) as f32);
|
||||
for ((((mut k, v), g), &kx), &ky) in k
|
||||
.slice_mut(s![.., .., 0])
|
||||
.gencolumns_mut()
|
||||
.into_iter()
|
||||
.zip(y.slice(s![.., .., 0]).gencolumns())
|
||||
.zip(g.gencolumns())
|
||||
.zip(grid.detj_dxi_dx.slice(s![.., 0]))
|
||||
.zip(grid.detj_dxi_dy.slice(s![.., 0]))
|
||||
{
|
||||
let tau = 1.0;
|
||||
|
||||
let v = (v[0], v[1], v[2]);
|
||||
let g = (g[0], g[1], g[2]);
|
||||
|
||||
let minus = negative_flux(kx, ky);
|
||||
|
||||
k[0] += tau
|
||||
* hinv
|
||||
* (minus[0][0] * (v.0 - g.0)
|
||||
+ minus[0][1] * (v.1 - g.1)
|
||||
+ minus[0][2] * (v.2 - g.2));
|
||||
k[1] += tau
|
||||
* hinv
|
||||
* (minus[1][0] * (v.0 - g.0)
|
||||
+ minus[1][1] * (v.1 - g.1)
|
||||
+ minus[1][2] * (v.2 - g.2));
|
||||
k[2] += tau
|
||||
* hinv
|
||||
* (minus[2][0] * (v.0 - g.0)
|
||||
+ minus[2][1] * (v.1 - g.1)
|
||||
+ minus[2][2] * (v.2 - g.2));
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
let g = match boundaries.north {
|
||||
Boundary::This => y.slice(s![.., 0, ..]),
|
||||
};
|
||||
let hinv = 1.0 / (SBP::h()[0] / (ny - 1) as f32);
|
||||
for ((((mut k, v), g), &kx), &ky) in k
|
||||
.slice_mut(s![.., ny - 1, ..])
|
||||
.gencolumns_mut()
|
||||
.into_iter()
|
||||
.zip(y.slice(s![.., ny - 1, ..]).gencolumns())
|
||||
.zip(g.gencolumns())
|
||||
.zip(grid.detj_deta_dx.slice(s![ny - 1, ..]))
|
||||
.zip(grid.detj_deta_dy.slice(s![ny - 1, ..]))
|
||||
{
|
||||
// North boundary, positive flux
|
||||
let tau = -1.0;
|
||||
let v = (v[0], v[1], v[2]);
|
||||
let g = (g[0], g[1], g[2]);
|
||||
|
||||
let plus = positive_flux(kx, ky);
|
||||
|
||||
k[0] += tau
|
||||
* hinv
|
||||
* (plus[0][0] * (v.0 - g.0) + plus[0][1] * (v.1 - g.1) + plus[0][2] * (v.2 - g.2));
|
||||
k[1] += tau
|
||||
* hinv
|
||||
* (plus[1][0] * (v.0 - g.0) + plus[1][1] * (v.1 - g.1) + plus[1][2] * (v.2 - g.2));
|
||||
k[2] += tau
|
||||
* hinv
|
||||
* (plus[2][0] * (v.0 - g.0) + plus[2][1] * (v.1 - g.1) + plus[2][2] * (v.2 - g.2));
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
let g = match boundaries.south {
|
||||
Boundary::This => y.slice(s![.., ny - 1, ..]),
|
||||
};
|
||||
let hinv = 1.0 / (SBP::h()[0] / (ny - 1) as f32);
|
||||
for ((((mut k, v), g), &kx), &ky) in k
|
||||
.slice_mut(s![.., 0, ..])
|
||||
.gencolumns_mut()
|
||||
.into_iter()
|
||||
.zip(y.slice(s![.., 0, ..]).gencolumns())
|
||||
.zip(g.gencolumns())
|
||||
.zip(grid.detj_deta_dx.slice(s![0, ..]))
|
||||
.zip(grid.detj_deta_dy.slice(s![0, ..]))
|
||||
{
|
||||
// South boundary, negative flux
|
||||
|
||||
let tau = 1.0;
|
||||
let v = (v[0], v[1], v[2]);
|
||||
let g = (g[0], g[1], g[2]);
|
||||
|
||||
let minus = negative_flux(kx, ky);
|
||||
|
||||
k[0] += tau
|
||||
* hinv
|
||||
* (minus[0][0] * (v.0 - g.0)
|
||||
+ minus[0][1] * (v.1 - g.1)
|
||||
+ minus[0][2] * (v.2 - g.2));
|
||||
k[1] += tau
|
||||
* hinv
|
||||
* (minus[1][0] * (v.0 - g.0)
|
||||
+ minus[1][1] * (v.1 - g.1)
|
||||
+ minus[1][2] * (v.2 - g.2));
|
||||
k[2] += tau
|
||||
* hinv
|
||||
* (minus[2][0] * (v.0 - g.0)
|
||||
+ minus[2][1] * (v.1 - g.1)
|
||||
+ minus[2][2] * (v.2 - g.2));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct WorkBuffers {
|
||||
k: [Field; 4],
|
||||
tmp: (Array2<f32>, Array2<f32>, Array2<f32>, Array2<f32>),
|
||||
}
|
||||
|
||||
impl WorkBuffers {
|
||||
pub fn new(ny: usize, nx: usize) -> Self {
|
||||
let arr2 = Array2::zeros((ny, nx));
|
||||
let arr3 = Field::new(ny, nx);
|
||||
Self {
|
||||
k: [arr3.clone(), arr3.clone(), arr3.clone(), arr3],
|
||||
tmp: (arr2.clone(), arr2.clone(), arr2.clone(), arr2),
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -1,78 +0,0 @@
|
||||
#![allow(clippy::excessive_precision)]
|
||||
#![allow(clippy::unreadable_literal)]
|
||||
|
||||
use ndarray::{ArrayView2, ArrayViewMut2};
|
||||
|
||||
pub trait SbpOperator {
|
||||
fn diffxi(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>);
|
||||
fn diffeta(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>);
|
||||
fn h() -> &'static [f32];
|
||||
}
|
||||
|
||||
pub trait UpwindOperator: SbpOperator {
|
||||
fn dissxi(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>);
|
||||
fn disseta(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>);
|
||||
}
|
||||
|
||||
#[macro_export]
|
||||
macro_rules! diff_op_1d {
|
||||
($self: ty, $name: ident, $BLOCK: expr, $DIAG: expr, $symmetric: expr) => {
|
||||
impl $self {
|
||||
fn $name(prev: ArrayView1<f32>, mut fut: ArrayViewMut1<f32>) {
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
let nx = prev.shape()[0];
|
||||
assert!(nx >= 2 * $BLOCK.len());
|
||||
|
||||
let dx = 1.0 / (nx - 1) as f32;
|
||||
let idx = 1.0 / dx;
|
||||
|
||||
let block = ::ndarray::arr2($BLOCK);
|
||||
let diag = ::ndarray::arr1($DIAG);
|
||||
|
||||
|
||||
let first_elems = prev.slice(::ndarray::s!(..block.len_of(::ndarray::Axis(1))));
|
||||
for (bl, f) in block.outer_iter().zip(&mut fut) {
|
||||
let diff = first_elems.dot(&bl);
|
||||
*f = diff * idx;
|
||||
}
|
||||
|
||||
// The window needs to be aligned to the diagonal elements,
|
||||
// based on the block size
|
||||
let window_elems_to_skip =
|
||||
block.len_of(::ndarray::Axis(0)) - ((diag.len() - 1) / 2);
|
||||
|
||||
for (window, f) in prev
|
||||
.windows(diag.len())
|
||||
.into_iter()
|
||||
.skip(window_elems_to_skip)
|
||||
.zip(fut.iter_mut().skip(block.len_of(::ndarray::Axis(0))))
|
||||
.take(nx - 2 * block.len_of(::ndarray::Axis(0)))
|
||||
{
|
||||
let diff = diag.dot(&window);
|
||||
*f = diff * idx;
|
||||
}
|
||||
|
||||
let last_elems = prev.slice(::ndarray::s!(nx - block.len_of(::ndarray::Axis(1))..;-1));
|
||||
for (bl, f) in block.outer_iter()
|
||||
.zip(&mut fut.slice_mut(s![nx - block.len_of(::ndarray::Axis(0))..;-1]))
|
||||
{
|
||||
let diff = if $symmetric {
|
||||
bl.dot(&last_elems)
|
||||
} else {
|
||||
-bl.dot(&last_elems)
|
||||
};
|
||||
*f = diff * idx;
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
mod upwind4;
|
||||
pub use upwind4::Upwind4;
|
||||
mod upwind9;
|
||||
pub use upwind9::Upwind9;
|
||||
mod traditional4;
|
||||
pub use traditional4::SBP4;
|
||||
mod traditional8;
|
||||
pub use traditional8::SBP8;
|
||||
@@ -1,45 +0,0 @@
|
||||
use super::SbpOperator;
|
||||
use crate::diff_op_1d;
|
||||
use ndarray::{s, ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
|
||||
|
||||
pub struct SBP4 {}
|
||||
|
||||
diff_op_1d!(SBP4, diff_1d, SBP4::BLOCK, SBP4::DIAG, false);
|
||||
|
||||
impl SBP4 {
|
||||
#[rustfmt::skip]
|
||||
const HBLOCK: &'static [f32] = &[
|
||||
17.0 / 48.0, 59.0 / 48.0, 43.0 / 48.0, 49.0 / 48.0,
|
||||
];
|
||||
#[rustfmt::skip]
|
||||
const DIAG: &'static [f32] = &[
|
||||
1.0 / 12.0, -2.0 / 3.0, 0.0, 2.0 / 3.0, -1.0 / 12.0,
|
||||
];
|
||||
#[rustfmt::skip]
|
||||
const BLOCK: &'static [[f32; 6]] = &[
|
||||
[-1.41176470588235e+00, 1.73529411764706e+00, -2.35294117647059e-01, -8.82352941176471e-02, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-5.00000000000000e-01, 0.00000000000000e+00, 5.00000000000000e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[9.30232558139535e-02, -6.86046511627907e-01, 0.00000000000000e+00, 6.86046511627907e-01, -9.30232558139535e-02, 0.00000000000000e+00],
|
||||
[3.06122448979592e-02, 0.00000000000000e+00, -6.02040816326531e-01, 0.00000000000000e+00, 6.53061224489796e-01, -8.16326530612245e-02],
|
||||
];
|
||||
}
|
||||
|
||||
impl SbpOperator for SBP4 {
|
||||
fn diffxi(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
|
||||
|
||||
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
||||
Self::diff_1d(r0, r1);
|
||||
}
|
||||
}
|
||||
|
||||
fn diffeta(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>) {
|
||||
// transpose then use diffxi
|
||||
Self::diffxi(prev.reversed_axes(), fut.reversed_axes());
|
||||
}
|
||||
|
||||
fn h() -> &'static [f32] {
|
||||
Self::HBLOCK
|
||||
}
|
||||
}
|
||||
@@ -1,49 +0,0 @@
|
||||
use super::SbpOperator;
|
||||
use crate::diff_op_1d;
|
||||
use ndarray::{s, ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
|
||||
|
||||
pub struct SBP8 {}
|
||||
|
||||
diff_op_1d!(SBP8, diff_1d, SBP8::BLOCK, SBP8::DIAG, false);
|
||||
|
||||
impl SBP8 {
|
||||
#[rustfmt::skip]
|
||||
const HBLOCK: &'static [f32] = &[
|
||||
2.94890676177879e-01, 1.52572062389771e+00, 2.57452876984127e-01, 1.79811370149912e+00, 4.12708057760141e-01, 1.27848462301587e+00, 9.23295579805997e-01, 1.00933386085916e+00
|
||||
];
|
||||
#[rustfmt::skip]
|
||||
const DIAG: &'static [f32] = &[
|
||||
3.57142857142857e-03, -3.80952380952381e-02, 2.00000000000000e-01, -8.00000000000000e-01, -0.00000000000000e+00, 8.00000000000000e-01, -2.00000000000000e-01, 3.80952380952381e-02, -3.57142857142857e-03
|
||||
];
|
||||
#[rustfmt::skip]
|
||||
const BLOCK: &'static [[f32; 12]] = &[
|
||||
[-1.69554360443190e+00, 2.24741246341404e+00, -3.38931922601500e-02, -7.81028168126749e-01, 2.54881486107905e-02, 3.43865227388873e-01, -8.62858162633335e-02, -2.00150583315761e-02, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-4.34378988266985e-01, 0.00000000000000e+00, 9.18511925072956e-02, 4.94008626807984e-01, -2.46151762937235e-02, -1.86759403432935e-01, 5.27267838475813e-02, 7.16696483080115e-03, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[3.88218088704253e-02, -5.44329744454984e-01, 0.00000000000000e+00, 3.89516189693211e-01, 1.36433486528546e-01, 1.03290582800845e-01, -1.79720579323281e-01, 5.59882558852296e-02, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[1.28088632226564e-01, -4.19172130036008e-01, -5.57707021445779e-02, 0.00000000000000e+00, 1.24714160903055e-01, 2.81285212519100e-01, -3.94470423942641e-02, -1.96981310738430e-02, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-1.82119472519009e-02, 9.09986646154550e-02, -8.51090570277506e-02, -5.43362886365301e-01, 0.00000000000000e+00, 6.37392455438558e-01, -1.02950081118829e-01, 2.98964956216039e-02, -8.65364391190110e-03, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-7.93147196245203e-02, 2.22875323171502e-01, -2.07999824391436e-02, -3.95611167748401e-01, -2.05756876210586e-01, 0.00000000000000e+00, 5.45876519966127e-01, -9.42727926638298e-02, 2.97971812952850e-02, -2.79348574643297e-03, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[2.75587615266177e-02, -8.71295642560637e-02, 5.01135077563584e-02, 7.68229253600969e-02, 4.60181213406519e-02, -7.55873581663580e-01, 0.00000000000000e+00, 8.21713248844682e-01, -2.16615355227872e-01, 4.12600676624518e-02, -3.86813134335486e-03, 0.00000000000000e+00],
|
||||
[5.84767272160451e-03, -1.08336661209337e-02, -1.42810403117803e-02, 3.50919361287023e-02, -1.22244235731112e-02, 1.19411743193552e-01, -7.51668243727123e-01, 0.00000000000000e+00, 7.92601963555477e-01, -1.98150490888869e-01, 3.77429506454989e-02, -3.53840162301552e-03],
|
||||
];
|
||||
}
|
||||
|
||||
impl SbpOperator for SBP8 {
|
||||
fn diffxi(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
|
||||
|
||||
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
||||
Self::diff_1d(r0, r1);
|
||||
}
|
||||
}
|
||||
|
||||
fn diffeta(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>) {
|
||||
// transpose then use diffxi
|
||||
Self::diffxi(prev.reversed_axes(), fut.reversed_axes());
|
||||
}
|
||||
|
||||
fn h() -> &'static [f32] {
|
||||
Self::HBLOCK
|
||||
}
|
||||
}
|
||||
@@ -1,455 +0,0 @@
|
||||
use super::{SbpOperator, UpwindOperator};
|
||||
use crate::diff_op_1d;
|
||||
use ndarray::{s, ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Axis};
|
||||
|
||||
pub struct Upwind4 {}
|
||||
|
||||
/// Simdtype used in diff_simd_col
|
||||
type SimdT = packed_simd::f32x8;
|
||||
|
||||
diff_op_1d!(Upwind4, diff_1d, Upwind4::BLOCK, Upwind4::DIAG, false);
|
||||
diff_op_1d!(
|
||||
Upwind4,
|
||||
diss_1d,
|
||||
Upwind4::DISS_BLOCK,
|
||||
Upwind4::DISS_DIAG,
|
||||
true
|
||||
);
|
||||
|
||||
macro_rules! diff_simd_row_7_47 {
|
||||
($self: ident, $name: ident, $BLOCK: expr, $DIAG: expr, $symmetric: expr) => {
|
||||
impl $self {
|
||||
#[inline(never)]
|
||||
fn $name(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
|
||||
use packed_simd::{f32x8, u32x8};
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
assert!(prev.len_of(Axis(1)) >= 2 * $BLOCK.len());
|
||||
assert!(prev.len() >= f32x8::lanes());
|
||||
// The prev and fut array must have contigous last dimension
|
||||
assert_eq!(prev.strides()[1], 1);
|
||||
assert_eq!(fut.strides()[1], 1);
|
||||
|
||||
let nx = prev.len_of(Axis(1));
|
||||
let dx = 1.0 / (nx - 1) as f32;
|
||||
let idx = 1.0 / dx;
|
||||
|
||||
for j in 0..prev.len_of(Axis(0)) {
|
||||
use std::slice;
|
||||
let prev =
|
||||
unsafe { slice::from_raw_parts(prev.uget((j, 0)) as *const f32, nx) };
|
||||
let fut =
|
||||
unsafe { slice::from_raw_parts_mut(fut.uget_mut((j, 0)) as *mut f32, nx) };
|
||||
//let mut fut = fut.slice_mut(s![j, ..]);
|
||||
|
||||
let first_elems = unsafe { f32x8::from_slice_unaligned_unchecked(prev) };
|
||||
let block = {
|
||||
let bl = $BLOCK;
|
||||
[
|
||||
f32x8::new(
|
||||
bl[0][0], bl[0][1], bl[0][2], bl[0][3], bl[0][4], bl[0][5],
|
||||
bl[0][6], 0.0,
|
||||
),
|
||||
f32x8::new(
|
||||
bl[1][0], bl[1][1], bl[1][2], bl[1][3], bl[1][4], bl[1][5],
|
||||
bl[1][6], 0.0,
|
||||
),
|
||||
f32x8::new(
|
||||
bl[2][0], bl[2][1], bl[2][2], bl[2][3], bl[2][4], bl[2][5],
|
||||
bl[2][6], 0.0,
|
||||
),
|
||||
f32x8::new(
|
||||
bl[3][0], bl[3][1], bl[3][2], bl[3][3], bl[3][4], bl[3][5],
|
||||
bl[3][6], 0.0,
|
||||
),
|
||||
]
|
||||
};
|
||||
fut[0] = idx * (block[0] * first_elems).sum();
|
||||
fut[1] = idx * (block[1] * first_elems).sum();
|
||||
fut[2] = idx * (block[2] * first_elems).sum();
|
||||
fut[3] = idx * (block[3] * first_elems).sum();
|
||||
|
||||
let diag = {
|
||||
let diag = $DIAG;
|
||||
f32x8::new(
|
||||
diag[0], diag[1], diag[2], diag[3], diag[4], diag[5], diag[6], 0.0,
|
||||
)
|
||||
};
|
||||
for (f, p) in fut
|
||||
.iter_mut()
|
||||
.skip(block.len())
|
||||
.zip(
|
||||
prev.windows(f32x8::lanes())
|
||||
.map(f32x8::from_slice_unaligned)
|
||||
.skip(1),
|
||||
)
|
||||
.take(nx - 2 * block.len())
|
||||
{
|
||||
*f = idx * (p * diag).sum();
|
||||
}
|
||||
|
||||
let last_elems =
|
||||
unsafe { f32x8::from_slice_unaligned_unchecked(&prev[nx - 8..]) }
|
||||
.shuffle1_dyn(u32x8::new(7, 6, 5, 4, 3, 2, 1, 0));
|
||||
if $symmetric {
|
||||
fut[nx - 4] = idx * (block[3] * last_elems).sum();
|
||||
fut[nx - 3] = idx * (block[2] * last_elems).sum();
|
||||
fut[nx - 2] = idx * (block[1] * last_elems).sum();
|
||||
fut[nx - 1] = idx * (block[0] * last_elems).sum();
|
||||
} else {
|
||||
fut[nx - 4] = -idx * (block[3] * last_elems).sum();
|
||||
fut[nx - 3] = -idx * (block[2] * last_elems).sum();
|
||||
fut[nx - 2] = -idx * (block[1] * last_elems).sum();
|
||||
fut[nx - 1] = -idx * (block[0] * last_elems).sum();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
diff_simd_row_7_47!(Upwind4, diff_simd_row, Upwind4::BLOCK, Upwind4::DIAG, false);
|
||||
diff_simd_row_7_47!(
|
||||
Upwind4,
|
||||
diss_simd_row,
|
||||
Upwind4::DISS_BLOCK,
|
||||
Upwind4::DISS_DIAG,
|
||||
true
|
||||
);
|
||||
|
||||
macro_rules! diff_simd_col_7_47 {
|
||||
($self: ident, $name: ident, $BLOCK: expr, $DIAG: expr, $symmetric: expr) => {
|
||||
impl $self {
|
||||
#[inline(never)]
|
||||
fn $name(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
|
||||
use std::slice;
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
assert_eq!(prev.stride_of(Axis(0)), 1);
|
||||
assert_eq!(prev.stride_of(Axis(0)), 1);
|
||||
let ny = prev.len_of(Axis(0));
|
||||
let nx = prev.len_of(Axis(1));
|
||||
assert!(nx >= 2 * $BLOCK.len());
|
||||
assert!(ny >= SimdT::lanes());
|
||||
assert!(ny % SimdT::lanes() == 0);
|
||||
|
||||
let dx = 1.0 / (nx - 1) as f32;
|
||||
let idx = 1.0 / dx;
|
||||
|
||||
for j in (0..ny).step_by(SimdT::lanes()) {
|
||||
let a = unsafe {
|
||||
[
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, 0)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, 1)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, 2)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, 3)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, 4)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, 5)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, 6)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
]
|
||||
};
|
||||
|
||||
for (i, bl) in $BLOCK.iter().enumerate() {
|
||||
let b = idx
|
||||
* (a[0] * bl[0]
|
||||
+ a[1] * bl[1]
|
||||
+ a[2] * bl[2]
|
||||
+ a[3] * bl[3]
|
||||
+ a[4] * bl[4]
|
||||
+ a[5] * bl[5]
|
||||
+ a[6] * bl[6]);
|
||||
unsafe {
|
||||
b.write_to_slice_unaligned(slice::from_raw_parts_mut(
|
||||
fut.uget_mut((j, i)) as *mut f32,
|
||||
SimdT::lanes(),
|
||||
));
|
||||
}
|
||||
}
|
||||
|
||||
let mut a = a;
|
||||
for i in $BLOCK.len()..nx - $BLOCK.len() {
|
||||
// Push a onto circular buffer
|
||||
a = [a[1], a[2], a[3], a[4], a[5], a[6], unsafe {
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, i + 3)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
))
|
||||
}];
|
||||
let b = idx
|
||||
* (a[0] * $DIAG[0]
|
||||
+ a[1] * $DIAG[1]
|
||||
+ a[2] * $DIAG[2]
|
||||
+ a[3] * $DIAG[3]
|
||||
+ a[4] * $DIAG[4]
|
||||
+ a[5] * $DIAG[5]
|
||||
+ a[6] * $DIAG[6]);
|
||||
unsafe {
|
||||
b.write_to_slice_unaligned(slice::from_raw_parts_mut(
|
||||
fut.uget_mut((j, i)) as *mut f32,
|
||||
SimdT::lanes(),
|
||||
));
|
||||
}
|
||||
}
|
||||
|
||||
let a = unsafe {
|
||||
[
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, nx - 1)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, nx - 2)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, nx - 3)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, nx - 4)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, nx - 5)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, nx - 6)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
SimdT::from_slice_unaligned(slice::from_raw_parts(
|
||||
prev.uget((j, nx - 7)) as *const f32,
|
||||
SimdT::lanes(),
|
||||
)),
|
||||
]
|
||||
};
|
||||
|
||||
for (i, bl) in $BLOCK.iter().enumerate() {
|
||||
let idx = if $symmetric { idx } else { -idx };
|
||||
let b = idx
|
||||
* (a[0] * bl[0]
|
||||
+ a[1] * bl[1]
|
||||
+ a[2] * bl[2]
|
||||
+ a[3] * bl[3]
|
||||
+ a[4] * bl[4]
|
||||
+ a[5] * bl[5]
|
||||
+ a[6] * bl[6]);
|
||||
unsafe {
|
||||
b.write_to_slice_unaligned(slice::from_raw_parts_mut(
|
||||
fut.uget_mut((j, nx - 1 - i)) as *mut f32,
|
||||
SimdT::lanes(),
|
||||
));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
diff_simd_col_7_47!(Upwind4, diff_simd_col, Upwind4::BLOCK, Upwind4::DIAG, false);
|
||||
diff_simd_col_7_47!(
|
||||
Upwind4,
|
||||
diss_simd_col,
|
||||
Upwind4::DISS_BLOCK,
|
||||
Upwind4::DISS_DIAG,
|
||||
true
|
||||
);
|
||||
|
||||
impl Upwind4 {
|
||||
#[rustfmt::skip]
|
||||
const HBLOCK: &'static [f32] = &[
|
||||
49.0 / 144.0, 61.0 / 48.0, 41.0 / 48.0, 149.0 / 144.0
|
||||
];
|
||||
#[rustfmt::skip]
|
||||
const DIAG: &'static [f32] = &[
|
||||
-1.0 / 24.0, 1.0 / 4.0, -7.0 / 8.0, 0.0, 7.0 / 8.0, -1.0 / 4.0, 1.0 / 24.0
|
||||
];
|
||||
#[rustfmt::skip]
|
||||
const BLOCK: &'static [[f32; 7]] = &[
|
||||
[ -72.0 / 49.0, 187.0 / 98.0, -20.0 / 49.0, -3.0 / 98.0, 0.0, 0.0, 0.0],
|
||||
[-187.0 / 366.0, 0.0, 69.0 / 122.0, -16.0 / 183.0, 2.0 / 61.0, 0.0, 0.0],
|
||||
[ 20.0 / 123.0, -69.0 / 82.0, 0.0, 227.0 / 246.0, -12.0 / 41.0, 2.0 / 41.0, 0.0],
|
||||
[ 3.0 / 298.0, 16.0 / 149.0, -227.0 / 298.0, 0.0, 126.0 / 149.0, -36.0 / 149.0, 6.0 / 149.0],
|
||||
];
|
||||
|
||||
#[rustfmt::skip]
|
||||
const DISS_BLOCK: &'static [[f32; 7]; 4] = &[
|
||||
[-3.0 / 49.0, 9.0 / 49.0, -9.0 / 49.0, 3.0 / 49.0, 0.0, 0.0, 0.0],
|
||||
[ 3.0 / 61.0, -11.0 / 61.0, 15.0 / 61.0, -9.0 / 61.0, 2.0 / 61.0, 0.0, 0.0],
|
||||
[-3.0 / 41.0, 15.0 / 41.0, -29.0 / 41.0, 27.0 / 41.0, -12.0 / 41.0, 2.0 / 41.0, 0.0],
|
||||
[3.0 / 149.0, -27.0 / 149.0, 81.0 / 149.0, -117.0 / 149.0, 90.0 / 149.0, -36.0 / 149.0, 6.0 / 149.0],
|
||||
];
|
||||
|
||||
#[rustfmt::skip]
|
||||
const DISS_DIAG: &'static [f32; 7] = &[
|
||||
1.0 / 24.0, -1.0 / 4.0, 5.0 / 8.0, -5.0 / 6.0, 5.0 / 8.0, -1.0 / 4.0, 1.0 / 24.0
|
||||
];
|
||||
}
|
||||
|
||||
impl SbpOperator for Upwind4 {
|
||||
fn diffxi(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
|
||||
|
||||
match (prev.strides(), fut.strides()) {
|
||||
([_, 1], [_, 1]) => {
|
||||
Self::diff_simd_row(prev, fut);
|
||||
}
|
||||
([1, _], [1, _]) if prev.len_of(Axis(0)) % SimdT::lanes() == 0 => {
|
||||
Self::diff_simd_col(prev, fut);
|
||||
}
|
||||
([_, _], [_, _]) => {
|
||||
// Fallback, work row by row
|
||||
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
||||
Self::diff_1d(r0, r1);
|
||||
}
|
||||
}
|
||||
_ => unreachable!("Should only be two elements in the strides vectors"),
|
||||
}
|
||||
}
|
||||
|
||||
fn diffeta(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>) {
|
||||
// transpose then use diffxi
|
||||
Self::diffxi(prev.reversed_axes(), fut.reversed_axes());
|
||||
}
|
||||
|
||||
fn h() -> &'static [f32] {
|
||||
Self::HBLOCK
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn upwind4_test() {
|
||||
use ndarray::prelude::*;
|
||||
let nx = 20;
|
||||
let dx = 1.0 / (nx - 1) as f32;
|
||||
let mut source: ndarray::Array1<f32> = ndarray::Array1::zeros(nx);
|
||||
let mut res = ndarray::Array1::zeros(nx);
|
||||
let mut target = ndarray::Array1::zeros(nx);
|
||||
|
||||
for i in 0..nx {
|
||||
source[i] = i as f32 * dx;
|
||||
target[i] = 1.0;
|
||||
}
|
||||
res.fill(0.0);
|
||||
Upwind4::diff_1d(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4);
|
||||
{
|
||||
let source = source.to_owned().insert_axis(ndarray::Axis(0));
|
||||
let mut res = res.to_owned().insert_axis(ndarray::Axis(0));
|
||||
let target = target.to_owned().insert_axis(ndarray::Axis(0));
|
||||
res.fill(0.0);
|
||||
Upwind4::diffxi(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
|
||||
}
|
||||
|
||||
{
|
||||
let source = Array2::from_shape_fn((nx, 8), |(i, _)| source[i]);
|
||||
let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]);
|
||||
let mut res = Array2::zeros((nx, 8));
|
||||
res.fill(0.0);
|
||||
Upwind4::diffeta(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2);
|
||||
}
|
||||
|
||||
for i in 0..nx {
|
||||
let x = i as f32 * dx;
|
||||
source[i] = x * x;
|
||||
target[i] = 2.0 * x;
|
||||
}
|
||||
res.fill(0.0);
|
||||
Upwind4::diff_1d(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-4);
|
||||
{
|
||||
let source = source.to_owned().insert_axis(ndarray::Axis(0));
|
||||
let mut res = res.to_owned().insert_axis(ndarray::Axis(0));
|
||||
let target = target.to_owned().insert_axis(ndarray::Axis(0));
|
||||
res.fill(0.0);
|
||||
Upwind4::diffxi(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
|
||||
}
|
||||
|
||||
{
|
||||
let source = Array2::from_shape_fn((nx, 8), |(i, _)| source[i]);
|
||||
let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]);
|
||||
let mut res = Array2::zeros((nx, 8));
|
||||
res.fill(0.0);
|
||||
Upwind4::diffeta(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2);
|
||||
}
|
||||
|
||||
for i in 0..nx {
|
||||
let x = i as f32 * dx;
|
||||
source[i] = x * x * x;
|
||||
target[i] = 3.0 * x * x;
|
||||
}
|
||||
res.fill(0.0);
|
||||
Upwind4::diff_1d(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
|
||||
|
||||
{
|
||||
let source = source.to_owned().insert_axis(ndarray::Axis(0));
|
||||
let mut res = res.to_owned().insert_axis(ndarray::Axis(0));
|
||||
let target = target.to_owned().insert_axis(ndarray::Axis(0));
|
||||
res.fill(0.0);
|
||||
Upwind4::diffxi(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res, &target, epsilon = 1e-2);
|
||||
}
|
||||
|
||||
{
|
||||
let source = Array2::from_shape_fn((nx, 8), |(i, _)| source[i]);
|
||||
let target = Array2::from_shape_fn((nx, 8), |(i, _)| target[i]);
|
||||
let mut res = Array2::zeros((nx, 8));
|
||||
res.fill(0.0);
|
||||
Upwind4::diffeta(source.view(), res.view_mut());
|
||||
approx::assert_abs_diff_eq!(&res.to_owned(), &target.to_owned(), epsilon = 1e-2);
|
||||
}
|
||||
}
|
||||
|
||||
impl UpwindOperator for Upwind4 {
|
||||
fn dissxi(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
|
||||
|
||||
match (prev.strides(), fut.strides()) {
|
||||
([_, 1], [_, 1]) => {
|
||||
Self::diss_simd_row(prev, fut);
|
||||
}
|
||||
([1, _], [1, _]) if prev.len_of(Axis(0)) % SimdT::lanes() == 0 => {
|
||||
Self::diss_simd_col(prev, fut);
|
||||
}
|
||||
([_, _], [_, _]) => {
|
||||
// Fallback, work row by row
|
||||
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
||||
Self::diss_1d(r0, r1);
|
||||
}
|
||||
}
|
||||
_ => unreachable!("Should only be two elements in the strides vectors"),
|
||||
}
|
||||
}
|
||||
|
||||
fn disseta(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>) {
|
||||
// diffeta = transpose then use dissxi
|
||||
Self::dissxi(prev.reversed_axes(), fut.reversed_axes());
|
||||
}
|
||||
}
|
||||
@@ -1,89 +0,0 @@
|
||||
use super::{SbpOperator, UpwindOperator};
|
||||
use crate::diff_op_1d;
|
||||
use ndarray::{s, ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
|
||||
|
||||
pub struct Upwind9 {}
|
||||
|
||||
diff_op_1d!(Upwind9, diff_1d, Upwind9::BLOCK, Upwind9::DIAG, false);
|
||||
diff_op_1d!(
|
||||
Upwind9,
|
||||
diss_1d,
|
||||
Upwind9::DISS_BLOCK,
|
||||
Upwind9::DISS_DIAG,
|
||||
true
|
||||
);
|
||||
|
||||
impl Upwind9 {
|
||||
#[rustfmt::skip]
|
||||
const HBLOCK: &'static [f32] = &[
|
||||
1070017.0/3628800.0, 5537111.0/3628800.0, 103613.0/403200.0, 261115.0/145152.0, 298951.0/725760.0, 515677.0/403200.0, 3349879.0/3628800.0, 3662753.0/3628800.0
|
||||
];
|
||||
#[rustfmt::skip]
|
||||
const DIAG: &'static [f32] = &[
|
||||
-1.0/1260.0, 5.0/504.0, -5.0/84.0, 5.0/21.0, -5.0/6.0, 0.0, 5.0/6.0, -5.0/21.0, 5.0/84.0, -5.0/504.0, 1.0/1260.0,
|
||||
];
|
||||
#[rustfmt::skip]
|
||||
const BLOCK: &'static [[f32; 13]] = &[
|
||||
[-1.69567399396458e+00, 2.29023358159400e+00, -2.16473500425698e-01, -5.05879766354449e-01, -1.01161106778154e-01, 2.59147072064383e-01, 1.93922119400659e-02, -4.95844980755642e-02, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-4.42575354959737e-01, 0.00000000000000e+00, 1.91582959381899e-01, 2.82222626681305e-01, 1.12083989713257e-01, -1.51334868892111e-01, -2.23600502721044e-02, 3.03806983474913e-02, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[2.48392603571843e-01, -1.13758367065272e+00, 0.00000000000000e+00, 1.95334726810969e+00, -1.58879011773212e+00, 3.93797129320378e-01, 2.52140821030291e-01, -1.21304033647356e-01, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[8.29213105268236e-02, -2.39388470313226e-01, -2.79038666398460e-01, 0.00000000000000e+00, 3.43018053395471e-01, 1.10370852514749e-01, 1.72029988649808e-03, -2.00445645303789e-02, 4.41184918522490e-04, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[7.24159504343116e-02, -4.15199475743626e-01, 9.91181694804303e-01, -1.49802407438608e+00, 0.00000000000000e+00, 1.30188867830442e+00, -6.03535071819214e-01, 1.73429775718218e-01, -2.40842144699299e-02, 1.92673715759439e-03, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-5.97470838462221e-02, 1.80551858630298e-01, -7.91241454636765e-02, -1.55240829877729e-01, -4.19298775383066e-01, 0.00000000000000e+00, 6.42287612546289e-01, -1.48833147569152e-01, 4.65407609802260e-02, -7.75679349670433e-03, 6.20543479736347e-04, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-6.19425252179959e-03, 3.69595678895333e-02, -7.01892820620398e-02, -3.35233082197107e-03, 2.69304373763091e-01, -8.89857974743355e-01, 0.00000000000000e+00, 8.66656645522330e-01, -2.57919763669076e-01, 6.44799409172690e-02, -1.07466568195448e-02, 8.59732545563586e-04, 0.00000000000000e+00],
|
||||
[1.44853491014330e-02, -4.59275574977554e-02, 3.08833474560615e-02, 3.57240610228828e-02, -7.07760049349999e-02, 1.88587240076292e-01, -7.92626447113877e-01, 0.00000000000000e+00, 8.25608497215073e-01, -2.35888142061449e-01, 5.89720355153623e-02, -9.82867258589373e-03, 7.86293806871498e-04],
|
||||
];
|
||||
|
||||
#[rustfmt::skip]
|
||||
const DISS_BLOCK: &'static [[f32; 13]] = &[
|
||||
[-3.99020778658945e-04, 2.05394169917502e-03, -4.24493243399805e-03, 4.38126393542801e-03, -2.18883813216888e-03, 2.98565988131608e-04, 1.38484104084115e-04, -3.94643819928825e-05, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[3.96913216138553e-04, -2.28230530115522e-03, 5.43069719436758e-03, -6.81086901935894e-03, 4.69064759201504e-03, -1.61429862514855e-03, 1.62083873811316e-04, 2.71310693302277e-05, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-4.87084939816571e-03, 3.22464611075207e-02, -9.06094757860846e-02, 1.39830191253413e-01, -1.27675500367419e-01, 6.87310321912961e-02, -2.00917702215270e-02, 2.43991122096699e-03, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[7.18155125886276e-04, -5.77715378536685e-03, 1.99749582302141e-02, -3.87940986951101e-02, 4.62756436981388e-02, -3.46770570075288e-02, 1.59058082995305e-02, -4.06744078428648e-03, 4.41184918522490e-04, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[-1.56687484682703e-03, 1.73758484693946e-02, -7.96515646886111e-02, 2.02094401829054e-01, -3.16098733124618e-01, 3.17999240131250e-01, -2.06522928911140e-01, 8.37112455598470e-02, -1.92673715759439e-02, 1.92673715759439e-03, 0.00000000000000e+00, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[6.88352254356072e-05, -1.92595810396278e-03, 1.38098624496279e-02, -4.87746083763075e-02, 1.02417890394006e-01, -1.38292226669620e-01, 1.23829022892659e-01, -7.34723830823462e-02, 2.79244565881356e-02, -6.20543479736347e-03, 6.20543479736347e-04, 0.00000000000000e+00, 0.00000000000000e+00],
|
||||
[4.42345367100640e-05, 2.67913080025652e-04, -5.59301314813691e-03, 3.09954862110834e-02, -9.21529346596015e-02, 1.71559035817103e-01, -2.12738289547735e-01, 1.79835101537893e-01, -1.03167905467630e-01, 3.86879645503614e-02, -8.59732545563586e-03, 8.59732545563586e-04, 0.00000000000000e+00],
|
||||
[-1.15289127131636e-05, 4.10149803795578e-05, 6.21188131452618e-04, -7.24912245235322e-03, 3.41622279353287e-02, -9.30972311856124e-02, 1.64473506705108e-01, -1.98013074867399e-01, 1.65121699443015e-01, -9.43552568245798e-02, 3.53832213092174e-02, -7.86293806871498e-03, 7.86293806871498e-04]
|
||||
];
|
||||
|
||||
#[rustfmt::skip]
|
||||
const DISS_DIAG: &'static [f32] = &[
|
||||
1.0/1260.0, -1.0/126.0, 1.0/28.0, -2.0/21.0, 1.0/6.0, -1.0/5.0, 1.0/6.0, -2.0/21.0, 1.0/28.0, -1.0/126.0, 1.0/1260.0,
|
||||
];
|
||||
}
|
||||
|
||||
impl SbpOperator for Upwind9 {
|
||||
fn diffxi(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
|
||||
|
||||
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
||||
Self::diff_1d(r0, r1);
|
||||
}
|
||||
}
|
||||
|
||||
fn diffeta(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>) {
|
||||
// transpose then use diffxi
|
||||
Self::diffxi(prev.reversed_axes(), fut.reversed_axes());
|
||||
}
|
||||
|
||||
fn h() -> &'static [f32] {
|
||||
Self::HBLOCK
|
||||
}
|
||||
}
|
||||
|
||||
impl UpwindOperator for Upwind9 {
|
||||
fn dissxi(prev: ArrayView2<f32>, mut fut: ArrayViewMut2<f32>) {
|
||||
assert_eq!(prev.shape(), fut.shape());
|
||||
assert!(prev.shape()[1] >= 2 * Self::BLOCK.len());
|
||||
|
||||
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
||||
Self::diss_1d(r0, r1);
|
||||
}
|
||||
}
|
||||
|
||||
fn disseta(prev: ArrayView2<f32>, fut: ArrayViewMut2<f32>) {
|
||||
// diffeta = transpose then use dissxi
|
||||
Self::dissxi(prev.reversed_axes(), fut.reversed_axes());
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user