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

use matrices everywhere

This commit is contained in:
Magnus Ulimoen 2021-02-01 23:08:55 +01:00
parent 31ac46e386
commit 6f7268bf33
9 changed files with 405 additions and 1289 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
sbp/src/operators.rs
View File

@ -78,7 +78,9 @@ pub trait SbpOperator2d: Send + Sync {
} }
fn op_xi(&self) -> &dyn SbpOperator1d; fn op_xi(&self) -> &dyn SbpOperator1d;
fn op_eta(&self) -> &dyn SbpOperator1d; fn op_eta(&self) -> &dyn SbpOperator1d {
self.op_xi()
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> { fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
None None
@ -98,7 +100,9 @@ pub trait UpwindOperator2d: Send + Sync {
} }
fn op_xi(&self) -> &dyn UpwindOperator1d; fn op_xi(&self) -> &dyn UpwindOperator1d;
fn op_eta(&self) -> &dyn UpwindOperator1d; fn op_eta(&self) -> &dyn UpwindOperator1d {
self.op_xi()
}
} }
pub trait InterpolationOperator: Send + Sync { pub trait InterpolationOperator: Send + Sync {

616
sbp/src/operators/algos.rs
View File

@ -1,4 +1,5 @@
use super::*; use super::*;
use ndarray::s;
pub(crate) mod constmatrix; pub(crate) mod constmatrix;
pub(crate) use constmatrix::{flip_lr, flip_sign, flip_ud, ColVector, Matrix, RowVector}; pub(crate) use constmatrix::{flip_lr, flip_sign, flip_ud, ColVector, Matrix, RowVector};
@ -8,11 +9,55 @@ mod fastfloat;
#[cfg(feature = "fast-float")] #[cfg(feature = "fast-float")]
use fastfloat::FastFloat; use fastfloat::FastFloat;
#[derive(Clone, Debug, PartialEq)]
pub(crate) struct DiagonalMatrix<const B: usize> {
pub start: [Float; B],
pub diag: Float,
pub end: [Float; B],
}
impl<const B: usize> DiagonalMatrix<B> {
pub const fn new(block: [Float; B]) -> Self {
let start = block;
let diag = 1.0;
let mut end = block;
let mut i = 0;
while i < B {
end[i] = block[B - 1 - i];
i += 1;
}
Self { start, diag, end }
}
}
#[derive(Clone, Debug, PartialEq)]
pub(crate) struct BlockMatrix<const M: usize, const N: usize, const D: usize> {
pub start: Matrix<Float, M, N>,
pub diag: RowVector<Float, D>,
pub end: Matrix<Float, M, N>,
}
impl<const M: usize, const N: usize, const D: usize> BlockMatrix<M, N, D> {
pub const fn new(
start: Matrix<Float, M, N>,
diag: RowVector<Float, D>,
end: Matrix<Float, M, N>,
) -> Self {
Self { start, diag, end }
}
}
#[derive(PartialEq, Copy, Clone)]
pub(crate) enum OperatorType {
Normal,
H2,
// TODO: D2
}
#[inline(always)] #[inline(always)]
pub(crate) fn diff_op_1d_matrix<const M: usize, const N: usize, const D: usize>( /// Works on all 1d vectors
block: &Matrix<Float, M, N>, pub(crate) fn diff_op_1d_fallback<const M: usize, const N: usize, const D: usize>(
blockend: &Matrix<Float, M, N>, matrix: &BlockMatrix<M, N, D>,
diag: &RowVector<Float, D>,
optype: OperatorType, optype: OperatorType,
prev: ArrayView1<Float>, prev: ArrayView1<Float>,
mut fut: ArrayViewMut1<Float>, mut fut: ArrayViewMut1<Float>,
@ -29,12 +74,11 @@ pub(crate) fn diff_op_1d_matrix<const M: usize, const N: usize, const D: usize>(
}; };
let idx = 1.0 / dx; let idx = 1.0 / dx;
for (bl, f) in block.iter_rows().zip(&mut fut) { let (futstart, futmid, futend) =
let diff = bl fut.multi_slice_mut((s![..M], s![M..nx - 2 * M], s![nx - 2 * M..]));
.iter()
.zip(prev.iter()) for (bl, f) in matrix.start.iter_rows().zip(futstart) {
.map(|(x, y)| x * y) let diff = dotproduct(bl.iter(), prev.iter());
.sum::<Float>();
*f = diff * idx; *f = diff * idx;
} }
@ -46,42 +90,33 @@ pub(crate) fn diff_op_1d_matrix<const M: usize, const N: usize, const D: usize>(
.windows(D) .windows(D)
.into_iter() .into_iter()
.skip(window_elems_to_skip) .skip(window_elems_to_skip)
.zip(fut.iter_mut().skip(M)) .zip(futmid)
.take(nx - 2 * M)
{ {
let diff = diag.iter().zip(&window).map(|(x, y)| x * y).sum::<Float>(); let diff = dotproduct(matrix.diag.row(0), window);
*f = diff * idx; *f = diff * idx;
} }
let prev = prev.slice(ndarray::s![nx - N..]); let prev = prev.slice(ndarray::s![nx - N..]);
for (bl, f) in blockend.iter_rows().zip(fut.iter_mut().rev().take(M).rev()) { for (bl, f) in matrix.end.iter_rows().zip(futend) {
let diff = bl let diff = dotproduct(bl, prev);
.iter()
.zip(prev.iter())
.map(|(x, y)| x * y)
.sum::<Float>();
*f = diff * idx; *f = diff * idx;
} }
} }
#[inline(always)] #[inline(always)]
pub(crate) fn diff_op_1d_slice_matrix<const M: usize, const N: usize, const D: usize>( /// diff op in 1d for slices
block: &Matrix<Float, M, N>, pub(crate) fn diff_op_1d_slice<const M: usize, const N: usize, const D: usize>(
endblock: &Matrix<Float, M, N>, matrix: &BlockMatrix<M, N, D>,
diag: &RowVector<Float, D>,
optype: OperatorType, optype: OperatorType,
prev: &[Float], prev: &[Float],
fut: &mut [Float], fut: &mut [Float],
) { ) {
#[cfg(feature = "fast-float")] #[cfg(feature = "fast-float")]
let (block, endblock, diag, prev, fut) = { let (matrix, prev, fut) = {
use std::mem::transmute; use std::mem::transmute;
unsafe { unsafe {
( (
transmute::<_, &Matrix<FastFloat, M, N>>(block), transmute::<_, &BlockMatrix<FastFloat, M, N, D>>(matrix),
transmute::<_, &Matrix<FastFloat, M, N>>(endblock),
transmute::<_, &RowVector<FastFloat, D>>(diag),
transmute::<_, &[FastFloat]>(prev), transmute::<_, &[FastFloat]>(prev),
transmute::<_, &mut [FastFloat]>(fut), transmute::<_, &mut [FastFloat]>(fut),
) )
@ -113,7 +148,7 @@ pub(crate) fn diff_op_1d_slice_matrix<const M: usize, const N: usize, const D: u
let prev = ColVector::<_, N>::map_to_col(prev.array_windows::<N>().next().unwrap()); let prev = ColVector::<_, N>::map_to_col(prev.array_windows::<N>().next().unwrap());
let fut = ColVector::<_, M>::map_to_col_mut(futb1.try_into().unwrap()); let fut = ColVector::<_, M>::map_to_col_mut(futb1.try_into().unwrap());
fut.matmul_into(block, prev); fut.matmul_into(&matrix.start, prev);
*fut *= idx; *fut *= idx;
} }
@ -129,7 +164,7 @@ pub(crate) fn diff_op_1d_slice_matrix<const M: usize, const N: usize, const D: u
let fut = ColVector::<_, 1>::map_to_col_mut(f); let fut = ColVector::<_, 1>::map_to_col_mut(f);
let prev = ColVector::<_, D>::map_to_col(window); let prev = ColVector::<_, D>::map_to_col(window);
fut.matmul_into(diag, prev); fut.matmul_into(&matrix.diag, prev);
*fut *= idx; *fut *= idx;
} }
@ -138,154 +173,35 @@ pub(crate) fn diff_op_1d_slice_matrix<const M: usize, const N: usize, const D: u
let prev = ColVector::<_, N>::map_to_col(prev); let prev = ColVector::<_, N>::map_to_col(prev);
let fut = ColVector::<_, M>::map_to_col_mut(futb2.try_into().unwrap()); let fut = ColVector::<_, M>::map_to_col_mut(futb2.try_into().unwrap());
fut.matmul_into(endblock, prev); fut.matmul_into(&matrix.end, prev);
*fut *= idx; *fut *= idx;
} }
} }
#[inline(always)] #[inline(always)]
pub(crate) fn diff_op_1d( /// Will always work on 1d, delegated based on slicedness
block: &[&[Float]], pub(crate) fn diff_op_1d<const M: usize, const N: usize, const D: usize>(
diag: &[Float], matrix: &BlockMatrix<M, N, D>,
symmetry: Symmetry,
optype: OperatorType, optype: OperatorType,
prev: ArrayView1<Float>, prev: ArrayView1<Float>,
mut fut: ArrayViewMut1<Float>, mut fut: ArrayViewMut1<Float>,
) { ) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
let nx = prev.shape()[0]; let nx = prev.shape()[0];
assert!(nx >= 2 * block.len()); assert!(nx >= 2 * M);
let dx = if optype == OperatorType::H2 { if let Some((prev, fut)) = prev.as_slice().zip(fut.as_slice_mut()) {
1.0 / (nx - 2) as Float diff_op_1d_slice(matrix, optype, prev, fut)
} else { } else {
1.0 / (nx - 1) as Float diff_op_1d_fallback(matrix, optype, prev, fut)
};
let idx = 1.0 / dx;
for (bl, f) in block.iter().zip(&mut fut) {
let diff = bl
.iter()
.zip(prev.iter())
.map(|(x, y)| x * y)
.sum::<Float>();
*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() - ((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()))
.take(nx - 2 * block.len())
{
let diff = diag.iter().zip(&window).map(|(x, y)| x * y).sum::<Float>();
*f = diff * idx;
}
for (bl, f) in block.iter().zip(fut.iter_mut().rev()) {
let diff = bl
.iter()
.zip(prev.iter().rev())
.map(|(x, y)| x * y)
.sum::<Float>();
*f = idx
* if symmetry == Symmetry::Symmetric {
diff
} else {
-diff
};
}
}
#[derive(PartialEq, Copy, Clone)]
pub(crate) enum Symmetry {
Symmetric,
AntiSymmetric,
}
#[derive(PartialEq, Copy, Clone)]
pub(crate) enum OperatorType {
Normal,
H2,
}
#[inline(always)]
#[allow(unused)]
pub(crate) fn diff_op_col_naive(
block: &'static [&'static [Float]],
diag: &'static [Float],
symmetry: Symmetry,
optype: OperatorType,
) -> impl Fn(ArrayView2<Float>, ArrayViewMut2<Float>) {
#[inline(always)]
move |prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>| {
assert_eq!(prev.shape(), fut.shape());
let nx = prev.shape()[1];
assert!(nx >= 2 * block.len());
assert_eq!(prev.strides()[0], 1);
assert_eq!(fut.strides()[0], 1);
let dx = if optype == OperatorType::H2 {
1.0 / (nx - 2) as Float
} else {
1.0 / (nx - 1) as Float
};
let idx = 1.0 / dx;
fut.fill(0.0);
// First block
for (bl, mut fut) in block.iter().zip(fut.axis_iter_mut(ndarray::Axis(1))) {
debug_assert_eq!(fut.len(), prev.shape()[0]);
for (&bl, prev) in bl.iter().zip(prev.axis_iter(ndarray::Axis(1))) {
debug_assert_eq!(prev.len(), fut.len());
fut.scaled_add(idx * bl, &prev);
}
}
let half_diag_width = (diag.len() - 1) / 2;
assert!(half_diag_width <= block.len());
// Diagonal entries
for (ifut, mut fut) in fut
.axis_iter_mut(ndarray::Axis(1))
.enumerate()
.skip(block.len())
.take(nx - 2 * block.len())
{
for (id, &d) in diag.iter().enumerate() {
let offset = ifut - half_diag_width + id;
fut.scaled_add(idx * d, &prev.slice(ndarray::s![.., offset]))
}
}
// End block
for (bl, mut fut) in block.iter().zip(fut.axis_iter_mut(ndarray::Axis(1)).rev()) {
fut.fill(0.0);
for (&bl, prev) in bl.iter().zip(prev.axis_iter(ndarray::Axis(1)).rev()) {
if symmetry == Symmetry::Symmetric {
fut.scaled_add(idx * bl, &prev);
} else {
fut.scaled_add(-idx * bl, &prev);
}
}
}
} }
} }
#[inline(always)] #[inline(always)]
#[allow(unused)] #[allow(unused)]
pub(crate) fn diff_op_col_naive_matrix<const M: usize, const N: usize, const D: usize>( /// 2D diff fallback for when matrices are not slicable
block: &Matrix<Float, M, N>, pub(crate) fn diff_op_2d_fallback<const M: usize, const N: usize, const D: usize>(
blockend: &Matrix<Float, M, N>, matrix: &BlockMatrix<M, N, D>,
diag: &RowVector<Float, D>,
optype: OperatorType, optype: OperatorType,
prev: ArrayView2<Float>, prev: ArrayView2<Float>,
mut fut: ArrayViewMut2<Float>, mut fut: ArrayViewMut2<Float>,
@ -295,9 +211,6 @@ pub(crate) fn diff_op_col_naive_matrix<const M: usize, const N: usize, const D:
let ny = prev.shape()[0]; let ny = prev.shape()[0];
assert!(nx >= 2 * M); assert!(nx >= 2 * M);
assert_eq!(prev.strides()[0], 1);
assert_eq!(fut.strides()[0], 1);
let dx = if optype == OperatorType::H2 { let dx = if optype == OperatorType::H2 {
1.0 / (nx - 2) as Float 1.0 / (nx - 2) as Float
} else { } else {
@ -314,7 +227,11 @@ pub(crate) fn diff_op_col_naive_matrix<const M: usize, const N: usize, const D:
)); ));
// First block // First block
for (bl, mut fut) in block.iter_rows().zip(fut0.axis_iter_mut(ndarray::Axis(1))) { for (bl, mut fut) in matrix
.start
.iter_rows()
.zip(fut0.axis_iter_mut(ndarray::Axis(1)))
{
debug_assert_eq!(fut.len(), prev.shape()[0]); debug_assert_eq!(fut.len(), prev.shape()[0]);
for (&bl, prev) in bl.iter().zip(prev.axis_iter(ndarray::Axis(1))) { for (&bl, prev) in bl.iter().zip(prev.axis_iter(ndarray::Axis(1))) {
if bl == 0.0 { if bl == 0.0 {
@ -332,7 +249,7 @@ pub(crate) fn diff_op_col_naive_matrix<const M: usize, const N: usize, const D:
.axis_iter_mut(ndarray::Axis(1)) .axis_iter_mut(ndarray::Axis(1))
.zip(prev.windows((ny, D)).into_iter().skip(window_elems_to_skip)) .zip(prev.windows((ny, D)).into_iter().skip(window_elems_to_skip))
{ {
for (&d, id) in diag.iter().zip(id.axis_iter(ndarray::Axis(1))) { for (&d, id) in matrix.diag.iter().zip(id.axis_iter(ndarray::Axis(1))) {
if d == 0.0 { if d == 0.0 {
continue; continue;
} }
@ -342,7 +259,8 @@ pub(crate) fn diff_op_col_naive_matrix<const M: usize, const N: usize, const D:
// End block // End block
let prev = prev.slice(ndarray::s!(.., nx - N..)); let prev = prev.slice(ndarray::s!(.., nx - N..));
for (bl, mut fut) in blockend for (bl, mut fut) in matrix
.end
.iter_rows() .iter_rows()
.zip(futn.axis_iter_mut(ndarray::Axis(1))) .zip(futn.axis_iter_mut(ndarray::Axis(1)))
{ {
@ -356,215 +274,49 @@ pub(crate) fn diff_op_col_naive_matrix<const M: usize, const N: usize, const D:
} }
} }
#[inline(always)] pub(crate) fn diff_op_2d_sliceable<const M: usize, const N: usize, const D: usize>(
pub(crate) fn diff_op_col( matrix: &BlockMatrix<M, N, D>,
block: &'static [&'static [Float]],
diag: &'static [Float],
symmetry: Symmetry,
optype: OperatorType, optype: OperatorType,
) -> impl Fn(ArrayView2<Float>, ArrayViewMut2<Float>) { prev: ArrayView2<Float>,
diff_op_col_simd(block, diag, symmetry, optype) mut fut: ArrayViewMut2<Float>,
} ) {
#[inline(always)]
pub(crate) fn diff_op_col_simd(
block: &'static [&'static [Float]],
diag: &'static [Float],
symmetry: Symmetry,
optype: OperatorType,
) -> impl Fn(ArrayView2<Float>, ArrayViewMut2<Float>) {
#[inline(always)]
move |prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>| {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
let nx = prev.shape()[1]; let nx = prev.shape()[1];
assert!(nx >= 2 * block.len()); for (prev, mut fut) in prev.outer_iter().zip(fut.outer_iter_mut()) {
let prev = &prev.as_slice().unwrap()[..nx];
assert_eq!(prev.strides()[0], 1); let fut = &mut fut.as_slice_mut().unwrap()[..nx];
assert_eq!(fut.strides()[0], 1); diff_op_1d_slice(matrix, optype, prev, fut)
let dx = if optype == OperatorType::H2 {
1.0 / (nx - 2) as Float
} else {
1.0 / (nx - 1) as Float
};
let idx = 1.0 / dx;
#[cfg(not(feature = "f32"))]
type SimdT = packed_simd::f64x8;
#[cfg(feature = "f32")]
type SimdT = packed_simd::f32x16;
let ny = prev.shape()[0];
// How many elements that can be simdified
let simdified = SimdT::lanes() * (ny / SimdT::lanes());
let half_diag_width = (diag.len() - 1) / 2;
assert!(half_diag_width <= block.len());
let fut_base_ptr = fut.as_mut_ptr();
let fut_stride = fut.strides()[1];
let fut_ptr = |j, i| {
debug_assert!(j < ny && i < nx);
unsafe { fut_base_ptr.offset(fut_stride * i as isize + j as isize) }
};
let prev_base_ptr = prev.as_ptr();
let prev_stride = prev.strides()[1];
let prev_ptr = |j, i| {
debug_assert!(j < ny && i < nx);
unsafe { prev_base_ptr.offset(prev_stride * i as isize + j as isize) }
};
// Not algo necessary, but gives performance increase
assert_eq!(fut_stride, prev_stride);
// First block
{
for (ifut, &bl) in block.iter().enumerate() {
for j in (0..simdified).step_by(SimdT::lanes()) {
let index_to_simd = |i| unsafe {
// j never moves past end of slice due to step_by and
// rounding down
SimdT::from_slice_unaligned(std::slice::from_raw_parts(
prev_ptr(j, i),
SimdT::lanes(),
))
};
let mut f = SimdT::splat(0.0);
for (iprev, &bl) in bl.iter().enumerate() {
f = index_to_simd(iprev).mul_adde(SimdT::splat(bl), f);
}
f *= idx;
unsafe {
f.write_to_slice_unaligned(std::slice::from_raw_parts_mut(
fut_ptr(j, ifut),
SimdT::lanes(),
));
}
}
for j in simdified..ny {
unsafe {
let mut f = 0.0;
for (iprev, bl) in bl.iter().enumerate() {
f += bl * *prev_ptr(j, iprev);
}
*fut_ptr(j, ifut) = f * idx;
}
}
}
}
// Diagonal elements
{
for ifut in block.len()..nx - block.len() {
for j in (0..simdified).step_by(SimdT::lanes()) {
let index_to_simd = |i| unsafe {
// j never moves past end of slice due to step_by and
// rounding down
SimdT::from_slice_unaligned(std::slice::from_raw_parts(
prev_ptr(j, i),
SimdT::lanes(),
))
};
let mut f = SimdT::splat(0.0);
for (id, &d) in diag.iter().enumerate() {
let offset = ifut - half_diag_width + id;
f = index_to_simd(offset).mul_adde(SimdT::splat(d), f);
}
f *= idx;
unsafe {
// puts simd along stride 1, j never goes past end of slice
f.write_to_slice_unaligned(std::slice::from_raw_parts_mut(
fut_ptr(j, ifut),
SimdT::lanes(),
));
}
}
for j in simdified..ny {
let mut f = 0.0;
for (id, &d) in diag.iter().enumerate() {
let offset = ifut - half_diag_width + id;
unsafe {
f += d * *prev_ptr(j, offset);
}
}
unsafe {
*fut_ptr(j, ifut) = idx * f;
}
}
}
}
// End block
{
// Get blocks and corresponding offsets
// (rev to iterate in ifut increasing order)
for (bl, ifut) in block.iter().zip((0..nx).rev()) {
for j in (0..simdified).step_by(SimdT::lanes()) {
let index_to_simd = |i| unsafe {
// j never moves past end of slice due to step_by and
// rounding down
SimdT::from_slice_unaligned(std::slice::from_raw_parts(
prev_ptr(j, i),
SimdT::lanes(),
))
};
let mut f = SimdT::splat(0.0);
for (&bl, iprev) in bl.iter().zip((0..nx).rev()) {
f = index_to_simd(iprev).mul_adde(SimdT::splat(bl), f);
}
f = if symmetry == Symmetry::Symmetric {
f * idx
} else {
-f * idx
};
unsafe {
f.write_to_slice_unaligned(std::slice::from_raw_parts_mut(
fut_ptr(j, ifut),
SimdT::lanes(),
));
}
}
for j in simdified..ny {
unsafe {
let mut f = 0.0;
for (&bl, iprev) in bl.iter().zip((0..nx).rev()).rev() {
f += bl * *prev_ptr(j, iprev);
}
*fut_ptr(j, ifut) = if symmetry == Symmetry::Symmetric {
f * idx
} else {
-f * idx
};
}
}
}
}
} }
} }
#[inline(always)] #[inline(always)]
fn dotproduct<'a>(u: impl Iterator<Item = &'a Float>, v: impl Iterator<Item = &'a Float>) -> Float { /// Dispatch based on strides
u.zip(v).fold(0.0, |acc, (&u, &v)| { pub(crate) fn diff_op_2d<const M: usize, const N: usize, const D: usize>(
#[cfg(feature = "fast-float")] matrix: &BlockMatrix<M, N, D>,
{ optype: OperatorType,
// We do not care about the order of multiplication nor addition prev: ArrayView2<Float>,
(FastFloat::from(acc) + FastFloat::from(u) * FastFloat::from(v)).into() fut: ArrayViewMut2<Float>,
) {
assert_eq!(prev.shape(), fut.shape());
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => diff_op_2d_sliceable(matrix, optype, prev, fut),
_ => diff_op_2d_fallback(matrix, optype, prev, fut),
} }
#[cfg(not(feature = "fast-float"))]
{
acc + u * v
}
})
} }
/*
#[inline(always)] #[inline(always)]
/// Way to too much overhead with SIMD:
/// output SIMD oriented:
/// |S | = |P0 P1| |P0 P1|
/// |S | = a1|P0 P1| + b1|P0 P1|
/// |S | = |P0 P1| |P0 P1|
///
/// | S | = |P0 P1| |P0 P1|
/// | S | = a2|P0 P1| + b1|P0 P1|
/// | S | = |P0 P1| |P0 P1|
pub(crate) fn diff_op_col_matrix<const M: usize, const N: usize, const D: usize>( pub(crate) fn diff_op_col_matrix<const M: usize, const N: usize, const D: usize>(
block: &Matrix<Float, M, N>, matrix: &BlockMatrix<M, N, D>,
block2: &Matrix<Float, M, N>,
diag: &RowVector<Float, D>,
optype: OperatorType, optype: OperatorType,
prev: ArrayView2<Float>, prev: ArrayView2<Float>,
fut: ArrayViewMut2<Float>, fut: ArrayViewMut2<Float>,
@ -615,7 +367,7 @@ pub(crate) fn diff_op_col_matrix<const M: usize, const N: usize, const D: usize>
{ {
let prev = prev.slice(ndarray::s![.., ..N]); let prev = prev.slice(ndarray::s![.., ..N]);
let (prevb, prevl) = prev.split_at(Axis(0), simdified); let (prevb, prevl) = prev.split_at(Axis(0), simdified);
for (mut fut, &bl) in fut1.axis_iter_mut(Axis(1)).zip(block.iter_rows()) { for (mut fut, &bl) in fut1.axis_iter_mut(Axis(1)).zip(matrix.start.iter_rows()) {
let fut = fut.as_slice_mut().unwrap(); let fut = fut.as_slice_mut().unwrap();
let fut = &mut fut[..ny]; let fut = &mut fut[..ny];
@ -662,7 +414,7 @@ pub(crate) fn diff_op_col_matrix<const M: usize, const N: usize, const D: usize>
for (fut, prev) in fut.by_ref().zip(prev) { for (fut, prev) in fut.by_ref().zip(prev) {
let mut f = SimdT::splat(0.0); let mut f = SimdT::splat(0.0);
for (&d, prev) in diag.iter().zip(prev.axis_iter(Axis(1))) { for (&d, prev) in matrix.diag.iter().zip(prev.axis_iter(Axis(1))) {
let prev = prev.to_slice().unwrap(); let prev = prev.to_slice().unwrap();
let prev = SimdT::from_slice_unaligned(prev); let prev = SimdT::from_slice_unaligned(prev);
f = prev.mul_adde(SimdT::splat(d), f); f = prev.mul_adde(SimdT::splat(d), f);
@ -677,7 +429,7 @@ pub(crate) fn diff_op_col_matrix<const M: usize, const N: usize, const D: usize>
.zip(prevl.axis_iter(Axis(0))) .zip(prevl.axis_iter(Axis(0)))
{ {
let mut f = 0.0; let mut f = 0.0;
for (&d, prev) in diag.iter().zip(prev) { for (&d, prev) in matrix.diag.iter().zip(prev) {
f += d * prev; f += d * prev;
} }
*fut = idx * f; *fut = idx * f;
@ -687,7 +439,7 @@ pub(crate) fn diff_op_col_matrix<const M: usize, const N: usize, const D: usize>
// End block // End block
{ {
for (mut fut, &bl) in fut2.axis_iter_mut(Axis(1)).zip(block2.iter_rows()) { for (mut fut, &bl) in fut2.axis_iter_mut(Axis(1)).zip(matrix.end.iter_rows()) {
let fut = fut.as_slice_mut().unwrap(); let fut = fut.as_slice_mut().unwrap();
let fut = &mut fut[..ny]; let fut = &mut fut[..ny];
let mut fut = fut.chunks_exact_mut(SimdT::lanes()); let mut fut = fut.chunks_exact_mut(SimdT::lanes());
@ -720,94 +472,51 @@ pub(crate) fn diff_op_col_matrix<const M: usize, const N: usize, const D: usize>
} }
} }
} }
*/
#[inline(always)] #[inline(always)]
pub(crate) fn diff_op_row( fn dotproduct<'a>(
block: &'static [&'static [Float]], u: impl IntoIterator<Item = &'a Float>,
diag: &'static [Float], v: impl IntoIterator<Item = &'a Float>,
symmetry: Symmetry, ) -> Float {
optype: OperatorType, u.into_iter().zip(v.into_iter()).fold(0.0, |acc, (&u, &v)| {
) -> impl Fn(ArrayView2<Float>, ArrayViewMut2<Float>) { #[cfg(feature = "fast-float")]
#[inline(always)]
move |prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>| {
assert_eq!(prev.shape(), fut.shape());
let nx = prev.shape()[1];
assert!(nx >= 2 * block.len());
assert_eq!(prev.strides()[1], 1);
assert_eq!(fut.strides()[1], 1);
let dx = if optype == OperatorType::H2 {
1.0 / (nx - 2) as Float
} else {
1.0 / (nx - 1) as Float
};
let idx = 1.0 / dx;
for (prev, mut fut) in prev
.axis_iter(ndarray::Axis(0))
.zip(fut.axis_iter_mut(ndarray::Axis(0)))
{ {
let prev = prev.as_slice().unwrap(); // We do not care about the order of multiplication nor addition
let fut = fut.as_slice_mut().unwrap(); (FastFloat::from(acc) + FastFloat::from(u) * FastFloat::from(v)).into()
assert_eq!(prev.len(), fut.len());
assert!(prev.len() >= 2 * block.len());
for (bl, f) in block.iter().zip(fut.iter_mut()) {
let diff = dotproduct(bl.iter(), prev[..bl.len()].iter());
*f = diff * idx;
} }
#[cfg(not(feature = "fast-float"))]
// The window needs to be aligned to the diagonal elements,
// based on the block size
let window_elems_to_skip = block.len() - ((diag.len() - 1) / 2);
for (window, f) in prev
.windows(diag.len())
.skip(window_elems_to_skip)
.zip(fut.iter_mut().skip(block.len()))
.take(nx - 2 * block.len())
{ {
let diff = dotproduct(diag.iter(), window.iter()); acc + u * v
*f = diff * idx;
}
for (bl, f) in block.iter().zip(fut.iter_mut().rev()) {
let diff = dotproduct(bl.iter(), prev.iter().rev());
*f = idx
* if symmetry == Symmetry::Symmetric {
diff
} else {
-diff
};
}
}
} }
})
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
pub(crate) fn sparse_from_block( pub(crate) fn sparse_from_block<const M: usize, const N: usize, const D: usize>(
block: &[&[Float]], matrix: &BlockMatrix<M, N, D>,
diag: &[Float],
symmetry: Symmetry,
optype: OperatorType, optype: OperatorType,
n: usize, n: usize,
) -> sprs::CsMat<Float> { ) -> sprs::CsMat<Float> {
assert!(n >= 2 * block.len()); assert!(n >= 2 * M);
let nnz = { let nnz = {
let block_elems = block.iter().fold(0, |acc, x| { let blockstart_elems = matrix
acc + x .start
.iter()
.fold(0, |acc, &x| if x != 0.0 { acc + 1 } else { acc })
});
let diag_elems = diag
.iter() .iter()
.fold(0, |acc, &x| if x != 0.0 { acc + 1 } else { acc }); .fold(0, |acc, &x| if x != 0.0 { acc + 1 } else { acc });
2 * block_elems + (n - 2 * block.len()) * diag_elems let diag_elems = matrix
.diag
.iter()
.fold(0, |acc, &x| if x != 0.0 { acc + 1 } else { acc });
let blockend_elems = matrix
.end
.iter()
.fold(0, |acc, &x| if x != 0.0 { acc + 1 } else { acc });
blockstart_elems + (n - 2 * M) * diag_elems + blockend_elems
}; };
let mut mat = sprs::TriMat::with_capacity((n, n), nnz); let mut mat = sprs::TriMat::with_capacity((n, n), nnz);
@ -819,7 +528,7 @@ pub(crate) fn sparse_from_block(
}; };
let idx = 1.0 / dx; let idx = 1.0 / dx;
for (j, bl) in block.iter().enumerate() { for (j, bl) in matrix.start.iter_rows().enumerate() {
for (i, &b) in bl.iter().enumerate() { for (i, &b) in bl.iter().enumerate() {
if b == 0.0 { if b == 0.0 {
continue; continue;
@ -828,9 +537,9 @@ pub(crate) fn sparse_from_block(
} }
} }
for j in block.len()..n - block.len() { for j in M..n - M {
let half_diag_len = diag.len() / 2; let half_diag_len = D / 2;
for (&d, i) in diag.iter().zip(j - half_diag_len..) { for (&d, i) in matrix.diag.iter().zip(j - half_diag_len..) {
if d == 0.0 { if d == 0.0 {
continue; continue;
} }
@ -838,34 +547,35 @@ pub(crate) fn sparse_from_block(
} }
} }
for (bl, j) in block.iter().zip((0..n).rev()).rev() { for (bl, j) in matrix.end.iter_rows().zip(n - M..) {
for (&b, i) in bl.iter().zip((0..n).rev()).rev() { for (&b, i) in bl.iter().zip(n - N..) {
if b == 0.0 { if b == 0.0 {
continue; continue;
} }
if symmetry == Symmetry::AntiSymmetric {
mat.add_triplet(j, i, -b * idx);
} else {
mat.add_triplet(j, i, b * idx); mat.add_triplet(j, i, b * idx);
} }
} }
}
mat.to_csr() mat.to_csr()
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
pub(crate) fn h_matrix(diag: &[Float], n: usize, is_h2: bool) -> sprs::CsMat<Float> { pub(crate) fn h_matrix<const D: usize>(
hmatrix: &DiagonalMatrix<D>,
n: usize,
is_h2: bool,
) -> sprs::CsMat<Float> {
let h = if is_h2 { let h = if is_h2 {
1.0 / (n - 2) as Float 1.0 / (n - 2) as Float
} else { } else {
1.0 / (n - 1) as Float 1.0 / (n - 1) as Float
}; };
let nmiddle = n - 2 * diag.len(); let nmiddle = n - 2 * D;
let iter = diag let iter = hmatrix
.start
.iter() .iter()
.chain(std::iter::repeat(&1.0).take(nmiddle)) .chain(std::iter::repeat(&hmatrix.diag).take(nmiddle))
.chain(diag.iter().rev()) .chain(hmatrix.end.iter())
.map(|&x| h * x); .map(|&x| h * x);
let mut mat = sprs::TriMat::with_capacity((n, n), n); let mut mat = sprs::TriMat::with_capacity((n, n), n);

4
sbp/src/operators/algos/constmatrix.rs
View File

@ -72,6 +72,10 @@ impl<T, const M: usize, const N: usize> Matrix<T, M, N> {
) -> impl ExactSizeIterator<Item = &[T; N]> + DoubleEndedIterator<Item = &[T; N]> { ) -> impl ExactSizeIterator<Item = &[T; N]> + DoubleEndedIterator<Item = &[T; N]> {
self.data.iter() self.data.iter()
} }
#[inline(always)]
pub const fn row(&self, idx: usize) -> &[T; N] {
&self.data[idx]
}
} }
impl<T, const N: usize> ColVector<T, N> { impl<T, const N: usize> ColVector<T, N> {

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

@ -1,4 +1,6 @@
use super::{SbpOperator1d, SbpOperator2d}; use super::{
BlockMatrix, DiagonalMatrix, Matrix, OperatorType, RowVector, SbpOperator1d, SbpOperator2d,
};
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
@ -7,76 +9,61 @@ pub struct SBP4;
impl SBP4 { impl SBP4 {
#[rustfmt::skip] #[rustfmt::skip]
const HBLOCK: &'static [Float] = &[ const H: DiagonalMatrix<4> = DiagonalMatrix::new([
17.0 / 48.0, 59.0 / 48.0, 43.0 / 48.0, 49.0 / 48.0, 17.0 / 48.0, 59.0 / 48.0, 43.0 / 48.0, 49.0 / 48.0,
]; ]);
#[rustfmt::skip] #[rustfmt::skip]
const DIAG: &'static [Float] = &[ const DIFF_DIAG: RowVector<Float, 5> = RowVector::new([[
1.0 / 12.0, -2.0 / 3.0, 0.0, 2.0 / 3.0, -1.0 / 12.0, 1.0 / 12.0, -2.0 / 3.0, 0.0, 2.0 / 3.0, -1.0 / 12.0
]; ]]);
const DIAG_MATRIX: super::RowVector<Float, 5> =
super::RowVector::new([[1.0 / 12.0, -2.0 / 3.0, 0.0, 2.0 / 3.0, -1.0 / 12.0]]);
#[rustfmt::skip] #[rustfmt::skip]
const BLOCK: &'static [&'static [Float]] = &[ const DIFF_BLOCK: Matrix<Float, 4, 6> = Matrix::new([
&[-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0],
&[-1.0/2.0, 0.0, 1.0/2.0],
&[4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0],
&[3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0]
];
#[rustfmt::skip]
const BLOCK_MATRIX: super::Matrix<Float, 4, 6> = super::Matrix::new([
[-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0, 0.0, 0.0], [-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0, 0.0, 0.0],
[-1.0/2.0, 0.0, 1.0/2.0, 0.0, 0.0, 0.0], [-1.0/2.0, 0.0, 1.0/2.0, 0.0, 0.0, 0.0],
[4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0, 0.0], [4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0, 0.0],
[3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0] [3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0]
]); ]);
const BLOCKEND_MATRIX: super::Matrix<Float, 4, 6> = const DIFF_BLOCKEND: super::Matrix<Float, 4, 6> =
super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX))); super::flip_sign(super::flip_ud(super::flip_lr(Self::DIFF_BLOCK)));
const DIFF: BlockMatrix<4, 6, 5> =
BlockMatrix::new(Self::DIFF_BLOCK, Self::DIFF_DIAG, Self::DIFF_BLOCKEND);
#[rustfmt::skip] #[rustfmt::skip]
const D2DIAG: &'static [Float] = &[ const D2DIAG: RowVector<Float, 5> = RowVector::new([[
-1.0 / 12.0, 4.0 / 3.0, -5.0 / 2.0, 4.0 / 3.0, -1.0 / 12.0 -1.0 / 12.0, 4.0 / 3.0, -5.0 / 2.0, 4.0 / 3.0, -1.0 / 12.0
]; ]]);
#[rustfmt::skip] #[rustfmt::skip]
const D2BLOCK: &'static [&'static [Float]] = &[ const D2BLOCK: Matrix<Float, 4, 6> = Matrix::new([
&[2.0, -5.0, 4.0, -1.0], [2.0, -5.0, 4.0, -1.0, 0.0, 0.0],
&[1.0, -2.0, 1.0], [1.0, -2.0, 1.0, 0.0, 0.0, 0.0],
&[-4.0/43.0, 59.0/43.0, -110.0/43.0, 59.0/43.0, -4.0/43.0], [-4.0/43.0, 59.0/43.0, -110.0/43.0, 59.0/43.0, -4.0/43.0, 0.0],
&[-1.0/49.0, 0.0, 59.0/49.0, -118.0/49.0, 64.0/49.0, -4.0/49.0] [-1.0/49.0, 0.0, 59.0/49.0, -118.0/49.0, 64.0/49.0, -4.0/49.0]
]; ]);
const D2: BlockMatrix<4, 6, 5> = BlockMatrix::new(
Self::D2BLOCK,
Self::D2DIAG,
super::flip_ud(super::flip_lr(Self::D2BLOCK)),
);
} }
impl SbpOperator1d for SBP4 { impl SbpOperator1d for SBP4 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d_matrix( super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut)
&Self::BLOCK_MATRIX,
&Self::BLOCKEND_MATRIX,
&Self::DIAG_MATRIX,
super::OperatorType::Normal,
prev,
fut,
)
} }
fn h(&self) -> &'static [Float] { fn h(&self) -> &'static [Float] {
Self::HBLOCK &Self::H.start
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DIFF, super::OperatorType::Normal, n)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::Normal,
n,
)
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::h_matrix(Self::HBLOCK, n, self.is_h2()) super::h_matrix(&Self::H, n, self.is_h2())
} }
fn d2(&self) -> Option<&dyn super::SbpOperator1d2> { fn d2(&self) -> Option<&dyn super::SbpOperator1d2> {
@ -84,76 +71,19 @@ impl SbpOperator1d for SBP4 {
} }
} }
fn diff_op_row_local(prev: ndarray::ArrayView2<Float>, mut fut: ndarray::ArrayViewMut2<Float>) {
for (p, mut f) in prev
.axis_iter(ndarray::Axis(0))
.zip(fut.axis_iter_mut(ndarray::Axis(0)))
{
super::diff_op_1d_slice_matrix(
&SBP4::BLOCK_MATRIX,
&SBP4::BLOCKEND_MATRIX,
&SBP4::DIAG_MATRIX,
super::OperatorType::Normal,
p.as_slice().unwrap(),
f.as_slice_mut().unwrap(),
)
}
}
fn diff_op_col_local(prev: ndarray::ArrayView2<Float>, fut: ndarray::ArrayViewMut2<Float>) {
let optype = super::OperatorType::Normal;
super::diff_op_col_naive_matrix(
&SBP4::BLOCK_MATRIX,
&SBP4::BLOCKEND_MATRIX,
&SBP4::DIAG_MATRIX,
optype,
prev,
fut,
)
}
impl SbpOperator2d for SBP4 { impl SbpOperator2d for SBP4 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * SBP4::BLOCK.len()); super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut)
let symmetry = super::Symmetry::AntiSymmetric;
let optype = super::OperatorType::Normal;
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
//diff_op_row(SBP4::BLOCK, SBP4::DIAG, symmetry, optype)(prev, fut);
diff_op_row_local(prev, fut)
}
([1, _], [1, _]) => {
//diff_op_col(SBP4::BLOCK, SBP4::DIAG, symmetry, optype)(prev, fut);
diff_op_col_local(prev, fut)
}
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
SBP4.diff(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn SbpOperator1d { fn op_xi(&self) -> &dyn SbpOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
} }
impl super::SbpOperator1d2 for SBP4 { impl super::SbpOperator1d2 for SBP4 {
fn diff2(&self, prev: ArrayView1<Float>, mut fut: ArrayViewMut1<Float>) { fn diff2(&self, prev: ArrayView1<Float>, mut fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(&Self::D2, OperatorType::Normal, prev, fut.view_mut());
Self::D2BLOCK,
Self::D2DIAG,
super::Symmetry::Symmetric,
super::OperatorType::Normal,
prev,
fut.view_mut(),
);
let hi = (prev.len() - 1) as Float; let hi = (prev.len() - 1) as Float;
fut.map_inplace(|x| *x *= hi) fut.map_inplace(|x| *x *= hi)
} }
@ -162,13 +92,7 @@ impl super::SbpOperator1d2 for SBP4 {
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diff2_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diff2_matrix(&self, n: usize) -> sprs::CsMat<Float> {
let mut m = super::sparse_from_block( let mut m = super::sparse_from_block(&Self::D2, OperatorType::Normal, n);
Self::D2BLOCK,
Self::D2DIAG,
super::Symmetry::Symmetric,
super::OperatorType::Normal,
n,
);
let hi = (n - 1) as Float; let hi = (n - 1) as Float;
m.map_inplace(|v| v * hi); m.map_inplace(|v| v * hi);
m m

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

@ -1,4 +1,6 @@
use super::{diff_op_col, diff_op_row, Matrix, RowVector, SbpOperator1d, SbpOperator2d}; use super::{
BlockMatrix, DiagonalMatrix, Matrix, OperatorType, RowVector, SbpOperator1d, SbpOperator2d,
};
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
@ -7,30 +9,15 @@ pub struct SBP8;
impl SBP8 { impl SBP8 {
#[rustfmt::skip] #[rustfmt::skip]
const HBLOCK: &'static [Float] = &[ const H: DiagonalMatrix<8> = DiagonalMatrix::new([
2.94890676177879e-01, 1.52572062389771e+00, 2.57452876984127e-01, 1.79811370149912e+00, 4.12708057760141e-01, 1.27848462301587e+00, 9.23295579805997e-01, 1.00933386085916e+00 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] #[rustfmt::skip]
const DIAG: &'static [Float] = &[ const DIFF_DIAG: RowVector<Float, 9> = Matrix::new([[
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 DIAG_MATRIX: RowVector<Float, 9> = Matrix::new([[
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 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] #[rustfmt::skip]
const BLOCK: &'static [&'static [Float]] = &[ const DIFF_BLOCK: Matrix<Float, 8, 12> = Matrix::new([
&[-1.69554360443190e+00, 2.24741246341404e+00, -3.38931922601500e-02, -7.81028168126749e-01, 2.54881486107905e-02, 3.43865227388873e-01, -8.62858162633335e-02, -2.00150583315761e-02],
&[-4.34378988266985e-01, 0.00000000000000e+00, 9.18511925072956e-02, 4.94008626807984e-01, -2.46151762937235e-02, -1.86759403432935e-01, 5.27267838475813e-02, 7.16696483080115e-03],
&[3.88218088704253e-02, -5.44329744454984e-01, 0.00000000000000e+00, 3.89516189693211e-01, 1.36433486528546e-01, 1.03290582800845e-01, -1.79720579323281e-01, 5.59882558852296e-02],
&[1.28088632226564e-01, -4.19172130036008e-01, -5.57707021445779e-02, 0.00000000000000e+00, 1.24714160903055e-01, 2.81285212519100e-01, -3.94470423942641e-02, -1.96981310738430e-02],
&[-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],
&[-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],
&[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],
&[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],
];
#[rustfmt::skip]
const BLOCK_MATRIX: Matrix<Float, 8, 12> = Matrix::new([
[-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.0, 0.0, 0.0, 0.0], [-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.0, 0.0, 0.0, 0.0],
[-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.0, 0.0, 0.0, 0.0], [-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.0, 0.0, 0.0, 0.0],
[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.0, 0.0, 0.0, 0.0], [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.0, 0.0, 0.0, 0.0],
@ -40,107 +27,40 @@ impl SBP8 {
[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.0], [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.0],
[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], [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],
]); ]);
const BLOCKEND_MATRIX: super::Matrix<Float, 8, 12> = const DIFF: BlockMatrix<8, 12, 9> = BlockMatrix::new(
super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX))); Self::DIFF_BLOCK,
Self::DIFF_DIAG,
super::flip_sign(super::flip_ud(super::flip_lr(Self::DIFF_BLOCK))),
);
} }
impl SbpOperator1d for SBP8 { impl SbpOperator1d for SBP8 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d_matrix( super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut)
&Self::BLOCK_MATRIX,
&Self::BLOCKEND_MATRIX,
&Self::DIAG_MATRIX,
super::OperatorType::Normal,
prev,
fut,
)
} }
fn h(&self) -> &'static [Float] { fn h(&self) -> &'static [Float] {
Self::HBLOCK &Self::H.start
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DIFF, super::OperatorType::Normal, n)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::Normal,
n,
)
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::h_matrix(Self::HBLOCK, n, self.is_h2()) super::h_matrix(&Self::H, n, self.is_h2())
} }
} }
fn diff_op_row_local(prev: ndarray::ArrayView2<Float>, mut fut: ndarray::ArrayViewMut2<Float>) {
/*
let mut flipmatrix = SBP8::BLOCK_MATRIX;
flipmatrix *= &-1.0;
*/
for (p, mut f) in prev
.axis_iter(ndarray::Axis(0))
.zip(fut.axis_iter_mut(ndarray::Axis(0)))
{
super::diff_op_1d_slice_matrix(
&SBP8::BLOCK_MATRIX,
&SBP8::BLOCKEND_MATRIX,
//&flipmatrix,
&SBP8::DIAG_MATRIX,
super::OperatorType::Normal,
p.as_slice().unwrap(),
f.as_slice_mut().unwrap(),
)
}
}
fn diff_op_col_local(prev: ndarray::ArrayView2<Float>, fut: ndarray::ArrayViewMut2<Float>) {
let optype = super::OperatorType::Normal;
super::diff_op_col_matrix(
&SBP8::BLOCK_MATRIX,
&SBP8::BLOCKEND_MATRIX,
&SBP8::DIAG_MATRIX,
optype,
prev,
fut,
)
}
impl SbpOperator2d for SBP8 { impl SbpOperator2d for SBP8 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * SBP8::BLOCK.len()); super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut)
let symmetry = super::Symmetry::AntiSymmetric;
let optype = super::OperatorType::Normal;
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
//diff_op_row(SBP8::BLOCK, SBP8::DIAG, symmetry, optype)(prev, fut);
diff_op_row_local(prev, fut);
}
([1, _], [1, _]) => {
//diff_op_col(SBP8::BLOCK, SBP8::DIAG, symmetry, optype)(prev, fut);
diff_op_col_local(prev, fut)
}
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
SBP8.diff(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn SbpOperator1d { fn op_xi(&self) -> &dyn SbpOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
} }
#[test] #[test]

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

@ -1,241 +1,66 @@
use super::{ use super::{
diff_op_col, diff_op_row, Matrix, RowVector, SbpOperator1d, SbpOperator2d, UpwindOperator1d, BlockMatrix, DiagonalMatrix, Matrix, OperatorType, RowVector, SbpOperator1d, SbpOperator2d,
UpwindOperator2d, UpwindOperator1d, UpwindOperator2d,
}; };
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Axis}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
#[derive(Debug, Copy, Clone)] #[derive(Debug, Copy, Clone)]
pub struct Upwind4; pub struct Upwind4;
/// Simdtype used in diff_simd_col
#[cfg(feature = "f32")]
type SimdT = packed_simd::f32x8;
#[cfg(not(feature = "f32"))]
type SimdT = packed_simd::f64x8;
macro_rules! diff_simd_col_7_47 {
($name: ident, $BLOCK: expr, $DIAG: expr, $symmetric: expr) => {
#[inline(never)]
fn $name(prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
use std::slice;
assert_eq!(prev.shape(), fut.shape());
assert_eq!(prev.stride_of(Axis(0)), 1);
assert_eq!(fut.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 Float;
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 Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, 1)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, 2)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, 3)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, 4)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, 5)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, 6)) as *const Float,
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 Float,
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 Float,
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 Float,
SimdT::lanes(),
));
}
}
let a = unsafe {
[
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, nx - 1)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, nx - 2)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, nx - 3)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, nx - 4)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, nx - 5)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, nx - 6)) as *const Float,
SimdT::lanes(),
)),
SimdT::from_slice_unaligned(slice::from_raw_parts(
prev.uget((j, nx - 7)) as *const Float,
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 Float,
SimdT::lanes(),
));
}
}
}
}
};
}
diff_simd_col_7_47!(diff_simd_col, Upwind4::BLOCK, Upwind4::DIAG, false);
diff_simd_col_7_47!(diss_simd_col, Upwind4::DISS_BLOCK, Upwind4::DISS_DIAG, true);
impl Upwind4 { impl Upwind4 {
#[rustfmt::skip] #[rustfmt::skip]
pub const HBLOCK: &'static [Float] = &[ const H: DiagonalMatrix<4> = DiagonalMatrix::new([
49.0 / 144.0, 61.0 / 48.0, 41.0 / 48.0, 149.0 / 144.0 49.0 / 144.0, 61.0 / 48.0, 41.0 / 48.0, 149.0 / 144.0
];
#[rustfmt::skip]
const DIAG: &'static [Float] = &[
-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 DIAG_MATRIX : RowVector<Float, 7> = RowVector::new([
[-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] #[rustfmt::skip]
const BLOCK: &'static [&'static [Float]] = &[ const DIFF_BLOCK: Matrix<Float, 4, 7> = Matrix::new([
&[ -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 BLOCK_MATRIX: Matrix<Float, 4, 7> = Matrix::new([
[ -72.0 / 49.0, 187.0 / 98.0, -20.0 / 49.0, -3.0 / 98.0, 0.0, 0.0, 0.0], [ -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], [-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], [ 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], [ 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],
]); ]);
const BLOCKEND_MATRIX: Matrix<Float, 4, 7> = #[rustfmt::skip]
super::flip_sign(super::flip_ud(super::flip_lr(Self::BLOCK_MATRIX))); const DIFF_DIAG: RowVector<Float, 7> = RowVector::new([[
-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
]]);
const DIFF_BLOCKEND: Matrix<Float, 4, 7> =
super::flip_sign(super::flip_ud(super::flip_lr(Self::DIFF_BLOCK)));
const DIFF: BlockMatrix<4, 7, 7> =
BlockMatrix::new(Self::DIFF_BLOCK, Self::DIFF_DIAG, Self::DIFF_BLOCKEND);
#[rustfmt::skip] #[rustfmt::skip]
const DISS_BLOCK: &'static [&'static [Float]] = &[ const DISS_BLOCK: Matrix<Float, 4, 7> = Matrix::new([
&[-3.0 / 49.0, 9.0 / 49.0, -9.0 / 49.0, 3.0 / 49.0, 0.0, 0.0, 0.0], [-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 / 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 / 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], [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] #[rustfmt::skip]
const DISS_DIAG: &'static [Float; 7] = &[ const DISS_DIAG: RowVector<Float, 7> = Matrix::new([[
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 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
]; ]]);
const DISS_BLOCKEND: Matrix<Float, 4, 7> = super::flip_ud(super::flip_lr(Self::DISS_BLOCK));
const DISS: BlockMatrix<4, 7, 7> =
BlockMatrix::new(Self::DISS_BLOCK, Self::DISS_DIAG, Self::DISS_BLOCKEND);
} }
impl SbpOperator1d for Upwind4 { impl SbpOperator1d for Upwind4 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d_matrix( super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut)
&Self::BLOCK_MATRIX,
&Self::BLOCKEND_MATRIX,
&Self::DIAG_MATRIX,
super::OperatorType::Normal,
prev,
fut,
)
} }
fn h(&self) -> &'static [Float] { fn h(&self) -> &'static [Float] {
Self::HBLOCK &Self::H.start
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DIFF, super::OperatorType::Normal, n)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::Normal,
n,
)
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::h_matrix(Self::HBLOCK, n, self.is_h2()) super::h_matrix(&Self::H, n, self.is_h2())
} }
fn upwind(&self) -> Option<&dyn UpwindOperator1d> { fn upwind(&self) -> Option<&dyn UpwindOperator1d> {
@ -243,53 +68,14 @@ impl SbpOperator1d for Upwind4 {
} }
} }
fn diff_op_row_local(prev: ndarray::ArrayView2<Float>, mut fut: ndarray::ArrayViewMut2<Float>) {
for (p, mut f) in prev
.axis_iter(ndarray::Axis(0))
.zip(fut.axis_iter_mut(ndarray::Axis(0)))
{
super::diff_op_1d_slice_matrix(
&Upwind4::BLOCK_MATRIX,
&Upwind4::BLOCKEND_MATRIX,
&Upwind4::DIAG_MATRIX,
super::OperatorType::Normal,
p.as_slice().unwrap(),
f.as_slice_mut().unwrap(),
)
}
}
impl SbpOperator2d for Upwind4 { impl SbpOperator2d for Upwind4 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind4::BLOCK.len()); super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut);
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => diff_op_row_local(prev, fut),
([1, _], [1, _]) if prev.len_of(Axis(0)) % SimdT::lanes() == 0 => {
diff_simd_col(prev, fut)
}
([1, _], [1, _]) => diff_op_col(
Upwind4::BLOCK,
Upwind4::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::Normal,
)(prev, fut),
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
Upwind4.diff(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn SbpOperator1d { fn op_xi(&self) -> &dyn SbpOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> { fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
Some(Box::new(Self)) Some(Box::new(Self))
} }
@ -385,14 +171,7 @@ fn upwind4_test() {
impl UpwindOperator1d for Upwind4 { impl UpwindOperator1d for Upwind4 {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(&Self::DISS, super::OperatorType::Normal, prev, fut)
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::Normal,
prev,
fut,
)
} }
fn as_sbp(&self) -> &dyn SbpOperator1d { fn as_sbp(&self) -> &dyn SbpOperator1d {
@ -401,53 +180,20 @@ impl UpwindOperator1d for Upwind4 {
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DISS, super::OperatorType::Normal, n)
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::Normal,
n,
)
} }
} }
impl UpwindOperator2d for Upwind4 { impl UpwindOperator2d for Upwind4 {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind4::BLOCK.len());
match (prev.strides(), fut.strides()) { super::diff_op_2d(&Self::DISS, OperatorType::Normal, prev, fut)
([_, 1], [_, 1]) => diff_op_row(
Upwind4::DISS_BLOCK,
Upwind4::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::Normal,
)(prev, fut),
([1, _], [1, _]) if prev.len_of(Axis(0)) % SimdT::lanes() == 0 => {
diss_simd_col(prev, fut);
}
([1, _], [1, _]) => diff_op_col(
Upwind4::DISS_BLOCK,
Upwind4::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::Normal,
)(prev, fut),
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
Upwind4.diss(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn UpwindOperator1d { fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn UpwindOperator1d {
&Self
}
} }
#[test] #[test]

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

@ -1,5 +1,6 @@
use super::{ use super::{
diff_op_col, diff_op_row, SbpOperator1d, SbpOperator2d, UpwindOperator1d, UpwindOperator2d, BlockMatrix, DiagonalMatrix, Matrix, OperatorType, RowVector, SbpOperator1d, SbpOperator2d,
UpwindOperator1d, UpwindOperator2d,
}; };
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
@ -9,49 +10,52 @@ pub struct Upwind4h2;
impl Upwind4h2 { impl Upwind4h2 {
#[rustfmt::skip] #[rustfmt::skip]
const HBLOCK: &'static [Float] = &[ const H: DiagonalMatrix<4> = DiagonalMatrix::new([
0.91e2/0.720e3, 0.325e3/0.384e3, 0.595e3/0.576e3, 0.1909e4/0.1920e4, 0.91e2/0.720e3, 0.325e3/0.384e3, 0.595e3/0.576e3, 0.1909e4/0.1920e4,
]; ]);
#[rustfmt::skip] #[rustfmt::skip]
const DIAG: &'static [Float] = &[ const DIFF_BLOCK: Matrix<Float, 4, 7> = Matrix::new([
[-3.95604395604396e+00, 5.41758241758242e+00, -1.94505494505495e+00, 4.83516483516484e-01, 0.0, 0.0, 0.0],
[-8.09025641025641e-01, 0.00000000000000e+00, 1.03948717948718e+00, -2.47384615384615e-01, 1.69230769230769e-02, 0.0, 0.0],
[2.37983193277311e-01, -8.51680672268908e-01, 0.00000000000000e+00, 7.35966386554622e-01, -1.36134453781513e-01, 1.38655462184874e-02, 0.0],
[-6.14632442814737e-02, 2.10581456259822e-01, -7.64623712240265e-01, 0.00000000000000e+00, 7.42535358826611e-01, -1.41435306443164e-01, 1.44054478784704e-02],
]);
#[rustfmt::skip]
const DIFF_DIAG: RowVector<Float, 7> = RowVector::new([[
-1.43229166666667e-02, 1.40625000000000e-01, -7.38281250000000e-01, 0.00000000000000e+00, 7.38281250000000e-01, -1.40625000000000e-01, 1.43229166666667e-02 -1.43229166666667e-02, 1.40625000000000e-01, -7.38281250000000e-01, 0.00000000000000e+00, 7.38281250000000e-01, -1.40625000000000e-01, 1.43229166666667e-02
]; ]]);
#[rustfmt::skip] const DIFF: BlockMatrix<4, 7, 7> = BlockMatrix::new(
const BLOCK: &'static [&'static [Float]] = &[ Self::DIFF_BLOCK,
&[-3.95604395604396e+00, 5.41758241758242e+00, -1.94505494505495e+00, 4.83516483516484e-01], Self::DIFF_DIAG,
&[-8.09025641025641e-01, 0.00000000000000e+00, 1.03948717948718e+00, -2.47384615384615e-01, 1.69230769230769e-02], super::flip_sign(super::flip_ud(super::flip_lr(Self::DIFF_BLOCK))),
&[2.37983193277311e-01, -8.51680672268908e-01, 0.00000000000000e+00, 7.35966386554622e-01, -1.36134453781513e-01, 1.38655462184874e-02], );
&[-6.14632442814737e-02, 2.10581456259822e-01, -7.64623712240265e-01, 0.00000000000000e+00, 7.42535358826611e-01, -1.41435306443164e-01, 1.44054478784704e-02],
];
#[rustfmt::skip] #[rustfmt::skip]
const DISS_BLOCK: &'static [&'static [Float]] = &[ const DISS_BLOCK: Matrix<Float, 4, 7> = Matrix::new([
&[-2.76989342315976e-01, 5.19355016842454e-01, -3.46236677894969e-01, 1.03871003368491e-01], [-2.76989342315976e-01, 5.19355016842454e-01, -3.46236677894969e-01, 1.03871003368491e-01, 0.0, 0.0, 0.0],
&[7.75570158484731e-02, -1.62342481638964e-01, 1.47715500579822e-01, -7.98531117124082e-02, 1.69230769230769e-02], [7.75570158484731e-02, -1.62342481638964e-01, 1.47715500579822e-01, -7.98531117124082e-02, 1.69230769230769e-02, 0.0, 0.0],
&[-4.23630758836198e-02, 1.21027405937249e-01, -1.91609307039399e-01, 1.82272708078206e-01, -8.31932773109244e-02, 1.38655462184874e-02], [-4.23630758836198e-02, 1.21027405937249e-01, -1.91609307039399e-01, 1.82272708078206e-01, -8.31932773109244e-02, 1.38655462184874e-02, 0.0],
&[1.32037874021759e-02, -6.79734450144910e-02, 1.89370108794365e-01, -2.78654929966754e-01, 2.16081718177056e-01, -8.64326872708224e-02, 1.44054478784704e-02], [1.32037874021759e-02, -6.79734450144910e-02, 1.89370108794365e-01, -2.78654929966754e-01, 2.16081718177056e-01, -8.64326872708224e-02, 1.44054478784704e-02],
]; ]);
#[rustfmt::skip] #[rustfmt::skip]
const DISS_DIAG: &'static [Float; 7] = &[ const DISS_DIAG: RowVector<Float, 7> = RowVector::new([[
1.43229166666667e-02, -8.59375000000000e-02, 2.14843750000000e-01, -2.86458333333333e-01, 2.14843750000000e-01, -8.59375000000000e-02, 1.43229166666667e-02, 1.43229166666667e-02, -8.59375000000000e-02, 2.14843750000000e-01, -2.86458333333333e-01, 2.14843750000000e-01, -8.59375000000000e-02, 1.43229166666667e-02,
]; ]]);
const DISS: BlockMatrix<4, 7, 7> = BlockMatrix::new(
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::flip_ud(super::flip_lr(Self::DISS_BLOCK)),
);
} }
impl SbpOperator1d for Upwind4h2 { impl SbpOperator1d for Upwind4h2 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(&Self::DIFF, OperatorType::H2, prev, fut)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::H2,
prev,
fut,
)
} }
fn h(&self) -> &'static [Float] { fn h(&self) -> &'static [Float] {
Self::HBLOCK &Self::H.start
} }
fn is_h2(&self) -> bool { fn is_h2(&self) -> bool {
true true
@ -59,17 +63,11 @@ impl SbpOperator1d for Upwind4h2 {
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DIFF, OperatorType::H2, n)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::H2,
n,
)
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::h_matrix(Self::HBLOCK, n, self.is_h2()) super::h_matrix(&Self::H, n, self.is_h2())
} }
fn upwind(&self) -> Option<&dyn UpwindOperator1d> { fn upwind(&self) -> Option<&dyn UpwindOperator1d> {
@ -78,80 +76,28 @@ impl SbpOperator1d for Upwind4h2 {
} }
impl SbpOperator2d for Upwind4h2 { impl SbpOperator2d for Upwind4h2 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind4h2::BLOCK.len());
let symmetry = super::Symmetry::AntiSymmetric; super::diff_op_2d(&Self::DIFF, OperatorType::H2, prev, fut)
let optype = super::OperatorType::H2;
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
diff_op_row(Upwind4h2::BLOCK, Upwind4h2::DIAG, symmetry, optype)(prev, fut);
}
([1, _], [1, _]) => {
diff_op_col(Upwind4h2::BLOCK, Upwind4h2::DIAG, symmetry, optype)(prev, fut);
}
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
Upwind4h2.diff(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn SbpOperator1d { fn op_xi(&self) -> &dyn SbpOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> { fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
Some(Box::new(Self)) Some(Box::new(Self))
} }
} }
impl UpwindOperator2d for Upwind4h2 { impl UpwindOperator2d for Upwind4h2 {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind4h2::BLOCK.len());
let symmetry = super::Symmetry::Symmetric; super::diff_op_2d(&Self::DISS, OperatorType::H2, prev, fut)
let optype = super::OperatorType::H2;
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
diff_op_row(
Upwind4h2::DISS_BLOCK,
Upwind4h2::DISS_DIAG,
symmetry,
optype,
)(prev, fut);
}
([1, _], [1, _]) => {
diff_op_col(
Upwind4h2::DISS_BLOCK,
Upwind4h2::DISS_DIAG,
symmetry,
optype,
)(prev, fut);
}
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
Upwind4h2.diss(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn UpwindOperator1d { fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn UpwindOperator1d {
&Self
}
} }
#[test] #[test]
@ -181,14 +127,7 @@ fn upwind4h2_test() {
impl UpwindOperator1d for Upwind4h2 { impl UpwindOperator1d for Upwind4h2 {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(&Self::DISS, super::OperatorType::H2, prev, fut)
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::H2,
prev,
fut,
)
} }
fn as_sbp(&self) -> &dyn SbpOperator1d { fn as_sbp(&self) -> &dyn SbpOperator1d {
@ -197,12 +136,6 @@ impl UpwindOperator1d for Upwind4h2 {
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DISS, super::OperatorType::H2, n)
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::H2,
n,
)
} }
} }

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

@ -1,5 +1,6 @@
use super::{ use super::{
diff_op_col, diff_op_row, SbpOperator1d, SbpOperator2d, UpwindOperator1d, UpwindOperator2d, BlockMatrix, DiagonalMatrix, Matrix, OperatorType, RowVector, SbpOperator1d, SbpOperator2d,
UpwindOperator1d, UpwindOperator2d,
}; };
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
@ -9,72 +10,69 @@ pub struct Upwind9;
impl Upwind9 { impl Upwind9 {
#[rustfmt::skip] #[rustfmt::skip]
const HBLOCK: &'static [Float] = &[ const H: DiagonalMatrix<8> = DiagonalMatrix::new([
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 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] #[rustfmt::skip]
const DIAG: &'static [Float] = &[ const DIFF_DIAG: RowVector<Float, 11> = RowVector::new([[
-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, -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] #[rustfmt::skip]
const BLOCK: &'static [&'static [Float]] = &[ const DIFF_BLOCK: Matrix<Float, 8, 13> = Matrix::new([
&[-1.69567399396458e+00, 2.29023358159400e+00, -2.16473500425698e-01, -5.05879766354449e-01, -1.01161106778154e-01, 2.59147072064383e-01, 1.93922119400659e-02, -4.95844980755642e-02], [-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.0, 0.0, 0.0, 0.0, 0.0],
&[-4.42575354959737e-01, 0.00000000000000e+00, 1.91582959381899e-01, 2.82222626681305e-01, 1.12083989713257e-01, -1.51334868892111e-01, -2.23600502721044e-02, 3.03806983474913e-02], [-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.0, 0.0, 0.0, 0.0,0.0],
&[2.48392603571843e-01, -1.13758367065272e+00, 0.00000000000000e+00, 1.95334726810969e+00, -1.58879011773212e+00, 3.93797129320378e-01, 2.52140821030291e-01, -1.21304033647356e-01], [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.0, 0.0, 0.0, 0.0, 0.0],
&[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], [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.0, 0.0, 0.0, 0.0],
&[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], [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.0, 0.0, 0.0],
&[-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], [-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.0, 0.0],
&[-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], [-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.0],
&[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], [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],
]; ]);
const DIFF: BlockMatrix<8, 13, 11> = BlockMatrix::new(
Self::DIFF_BLOCK,
Self::DIFF_DIAG,
super::flip_sign(super::flip_ud(super::flip_lr(Self::DIFF_BLOCK))),
);
#[rustfmt::skip] #[rustfmt::skip]
const DISS_BLOCK: &'static [&'static [Float]] = &[ const DISS_BLOCK: Matrix<Float, 8, 13> = Matrix::new([
&[-3.99020778658945e-04, 2.05394169917502e-03, -4.24493243399805e-03, 4.38126393542801e-03, -2.18883813216888e-03, 2.98565988131608e-04, 1.38484104084115e-04, -3.94643819928825e-05], [-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.0, 0.0, 0.0, 0.0, 0.0],
&[3.96913216138553e-04, -2.28230530115522e-03, 5.43069719436758e-03, -6.81086901935894e-03, 4.69064759201504e-03, -1.61429862514855e-03, 1.62083873811316e-04, 2.71310693302277e-05], [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.0, 0.0, 0.0, 0.0, 0.0 ],
&[-4.87084939816571e-03, 3.22464611075207e-02, -9.06094757860846e-02, 1.39830191253413e-01, -1.27675500367419e-01, 6.87310321912961e-02, -2.00917702215270e-02, 2.43991122096699e-03], [-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.0, 0.0, 0.0, 0.0, 0.0],
&[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], [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.0, 0.0, 0.0, 0.0],
&[-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], [-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.0, 0.0, 0.0],
&[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], [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.0, 0.0],
&[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], [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.0],
&[-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] [-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] #[rustfmt::skip]
const DISS_DIAG: &'static [Float] = &[ const DISS_DIAG: RowVector<Float, 11> = RowVector::new([[
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, 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,
]; ]]);
const DISS: BlockMatrix<8, 13, 11> = BlockMatrix::new(
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::flip_ud(super::flip_lr(Self::DISS_BLOCK)),
);
} }
impl SbpOperator1d for Upwind9 { impl SbpOperator1d for Upwind9 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::Normal,
prev,
fut,
)
} }
fn h(&self) -> &'static [Float] { fn h(&self) -> &'static [Float] {
Self::HBLOCK &Self::H.start
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DIFF, super::OperatorType::Normal, n)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::Normal,
n,
)
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::h_matrix(Self::HBLOCK, n, self.is_h2()) super::h_matrix(&Self::H, n, self.is_h2())
} }
fn upwind(&self) -> Option<&dyn UpwindOperator1d> { fn upwind(&self) -> Option<&dyn UpwindOperator1d> {
@ -83,35 +81,14 @@ impl SbpOperator1d for Upwind9 {
} }
impl SbpOperator2d for Upwind9 { impl SbpOperator2d for Upwind9 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind9::BLOCK.len());
let symmetry = super::Symmetry::AntiSymmetric; super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut)
let optype = super::OperatorType::Normal;
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
diff_op_row(Upwind9::BLOCK, Upwind9::DIAG, symmetry, optype)(prev, fut);
}
([1, _], [1, _]) => {
diff_op_col(Upwind9::BLOCK, Upwind9::DIAG, symmetry, optype)(prev, fut);
}
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
Upwind9.diff(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn SbpOperator1d { fn op_xi(&self) -> &dyn SbpOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> { fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
Some(Box::new(Self)) Some(Box::new(Self))
} }
@ -119,14 +96,7 @@ impl SbpOperator2d for Upwind9 {
impl UpwindOperator1d for Upwind9 { impl UpwindOperator1d for Upwind9 {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(&Self::DISS, super::OperatorType::Normal, prev, fut)
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::Normal,
prev,
fut,
)
} }
fn as_sbp(&self) -> &dyn SbpOperator1d { fn as_sbp(&self) -> &dyn SbpOperator1d {
@ -135,46 +105,18 @@ impl UpwindOperator1d for Upwind9 {
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DISS, super::OperatorType::Normal, n)
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::Normal,
n,
)
} }
} }
impl UpwindOperator2d for Upwind9 { impl UpwindOperator2d for Upwind9 {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind9::BLOCK.len()); super::diff_op_2d(&Self::DISS, OperatorType::Normal, prev, fut)
let symmetry = super::Symmetry::Symmetric;
let optype = super::OperatorType::Normal;
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
diff_op_row(Upwind9::DISS_BLOCK, Upwind9::DISS_DIAG, symmetry, optype)(prev, fut);
}
([1, _], [1, _]) => {
diff_op_col(Upwind9::DISS_BLOCK, Upwind9::DISS_DIAG, symmetry, optype)(prev, fut);
}
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
Upwind9.diss(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn UpwindOperator1d { fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn UpwindOperator1d {
&Self
}
} }
#[test] #[test]

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

@ -1,5 +1,6 @@
use super::{ use super::{
diff_op_col, diff_op_row, SbpOperator1d, SbpOperator2d, UpwindOperator1d, UpwindOperator2d, BlockMatrix, DiagonalMatrix, Matrix, OperatorType, RowVector, SbpOperator1d, SbpOperator2d,
UpwindOperator1d, UpwindOperator2d,
}; };
use crate::Float; use crate::Float;
use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2};
@ -9,57 +10,60 @@ pub struct Upwind9h2;
impl Upwind9h2 { impl Upwind9h2 {
#[rustfmt::skip] #[rustfmt::skip]
const HBLOCK: &'static [Float] = &[ const H: DiagonalMatrix<8> = DiagonalMatrix::new([
0.13528301e8/0.97297200e8, 0.189103813e9/0.232243200e9, 0.125336387e9/0.116121600e9, 0.24412231e8/0.25804800e8, 0.11976149e8/0.11612160e8, 0.27510113e8/0.27869184e8, 0.142384289e9/0.141926400e9, 0.3018054133e10/0.3019161600e10 0.13528301e8/0.97297200e8, 0.189103813e9/0.232243200e9, 0.125336387e9/0.116121600e9, 0.24412231e8/0.25804800e8, 0.11976149e8/0.11612160e8, 0.27510113e8/0.27869184e8, 0.142384289e9/0.141926400e9, 0.3018054133e10/0.3019161600e10
]; ]);
#[rustfmt::skip] #[rustfmt::skip]
const DIAG: &'static [Float] = &[ const DIFF_DIAG: RowVector<Float, 11> = RowVector::new([[
-7.93650793650794e-04, 9.92063492063492e-03, -5.95238095238095e-02, 2.38095238095238e-01, -8.33333333333333e-01, 0.00000000000000e+00, 8.33333333333333e-01, -2.38095238095238e-01, 5.95238095238095e-02, -9.92063492063492e-03, 7.93650793650794e-04 -7.93650793650794e-04, 9.92063492063492e-03, -5.95238095238095e-02, 2.38095238095238e-01, -8.33333333333333e-01, 0.00000000000000e+00, 8.33333333333333e-01, -2.38095238095238e-01, 5.95238095238095e-02, -9.92063492063492e-03, 7.93650793650794e-04
]; ]]);
#[rustfmt::skip] #[rustfmt::skip]
const BLOCK: &'static [&'static [Float]] = &[ const DIFF_BLOCK: Matrix<Float, 8, 13> = Matrix::new([
&[-3.59606132359119e+00, 4.89833536312447e+00, -1.88949383693789e+00, 5.47742030095865e-01, 1.98323891961440e-01, -1.81919774755160e-01, 9.75536252790282e-03, 1.33182875745630e-02], [-3.59606132359119e+00, 4.89833536312447e+00, -1.88949383693789e+00, 5.47742030095865e-01, 1.98323891961440e-01, -1.81919774755160e-01, 9.75536252790282e-03, 1.33182875745630e-02, 0.0, 0.0, 0.0, 0.0, 0.0],
&[-8.36438764514407e-01, 0.00000000000000e+00, 1.13444534067338e+00, -2.98589126517177e-01, -5.64158604558766e-02, 6.76718494319307e-02, -7.11434545077044e-03, -3.55909316707803e-03], [-8.36438764514407e-01, 0.00000000000000e+00, 1.13444534067338e+00, -2.98589126517177e-01, -5.64158604558766e-02, 6.76718494319307e-02, -7.11434545077044e-03, -3.55909316707803e-03, 0.0, 0.0, 0.0, 0.0, 0.0],
&[2.43402051702672e-01, -8.55808694890095e-01, 0.00000000000000e+00, 7.10023434316597e-01, -8.34706840757189e-02, -2.34098045062645e-02, 9.87675724649922e-03, -6.13059793689851e-04], [2.43402051702672e-01, -8.55808694890095e-01, 0.00000000000000e+00, 7.10023434316597e-01, -8.34706840757189e-02, -2.34098045062645e-02, 9.87675724649922e-03, -6.13059793689851e-04, 0.0, 0.0, 0.0, 0.0, 0.0],
&[-8.05029896098980e-02, 2.56994775122436e-01, -8.10083655220697e-01, 0.00000000000000e+00, 7.86933049074895e-01, -1.76732206311904e-01, 2.52232469669152e-02, -2.67114375632805e-03, 8.38923734582063e-04], [-8.05029896098980e-02, 2.56994775122436e-01, -8.10083655220697e-01, 0.00000000000000e+00, 7.86933049074895e-01, -1.76732206311904e-01, 2.52232469669152e-02, -2.67114375632805e-03, 8.38923734582063e-04, 0.0, 0.0, 0.0, 0.0],
&[-2.67370674924706e-02, 4.45404208225958e-02, 8.73562441688813e-02, -7.21839392529084e-01, 0.00000000000000e+00, 7.74603212862825e-01, -2.00161722480411e-01, 5.10878939438559e-02, -9.61911880020865e-03, 7.69529504016692e-04], [-2.67370674924706e-02, 4.45404208225958e-02, 8.73562441688813e-02, -7.21839392529084e-01, 0.00000000000000e+00, 7.74603212862825e-01, -2.00161722480411e-01, 5.10878939438559e-02, -9.61911880020865e-03, 7.69529504016692e-04, 0.0, 0.0, 0.0],
&[2.56244589039153e-02, -5.58209474181657e-02, 2.55972803937942e-02, 1.69377044771715e-01, -8.09310829929974e-01, 0.00000000000000e+00, 8.33417128898551e-01, -2.39938756878570e-01, 6.03007337701593e-02, -1.00501222950266e-02, 8.04009783602125e-04], [2.56244589039153e-02, -5.58209474181657e-02, 2.55972803937942e-02, 1.69377044771715e-01, -8.09310829929974e-01, 0.00000000000000e+00, 8.33417128898551e-01, -2.39938756878570e-01, 6.03007337701593e-02, -1.00501222950266e-02, 8.04009783602125e-04, 0.0, 0.0],
&[-1.35203347497451e-03, 5.77422023528126e-03, -1.06262408660325e-02, -2.37853245409338e-02, 2.05772021953463e-01, -8.20033622612654e-01, 0.00000000000000e+00, 8.31345778788959e-01, -2.37329555369694e-01, 5.93323888424235e-02, -9.88873147373725e-03, 7.91098517898980e-04], [-1.35203347497451e-03, 5.77422023528126e-03, -1.06262408660325e-02, -2.37853245409338e-02, 2.05772021953463e-01, -8.20033622612654e-01, 0.00000000000000e+00, 8.31345778788959e-01, -2.37329555369694e-01, 5.93323888424235e-02, -9.88873147373725e-03, 7.91098517898980e-04, 0.0],
&[-1.85246767034256e-03, 2.89905176240824e-03, 6.61951741227212e-04, 2.52792141498726e-03, -5.27086038822990e-02, 2.36934258275492e-01, -8.34333946306806e-01, 0.00000000000000e+00, 8.33639122800983e-01, -2.38182606514566e-01, 5.95456516286416e-02, -9.92427527144027e-03, 7.93942021715222e-04] [-1.85246767034256e-03, 2.89905176240824e-03, 6.61951741227212e-04, 2.52792141498726e-03, -5.27086038822990e-02, 2.36934258275492e-01, -8.34333946306806e-01, 0.00000000000000e+00, 8.33639122800983e-01, -2.38182606514566e-01, 5.95456516286416e-02, -9.92427527144027e-03, 7.93942021715222e-04]
]; ]);
const DIFF: BlockMatrix<8, 13, 11> = BlockMatrix::new(
Self::DIFF_BLOCK,
Self::DIFF_DIAG,
super::flip_sign(super::flip_ud(super::flip_lr(Self::DIFF_BLOCK))),
);
#[rustfmt::skip] #[rustfmt::skip]
const DISS_BLOCK: &'static [&'static [Float]] = &[ const DISS_BLOCK: Matrix<Float, 8, 13> = Matrix::new([
&[-9.67684432055993e-03, 2.46709290489728e-02, -3.59750445192099e-02, 3.68391266633410e-02, -2.15642540957128e-02, 6.31810399460400e-03, -5.50815969583136e-04, -6.12008018523449e-05], [-9.67684432055993e-03, 2.46709290489728e-02, -3.59750445192099e-02, 3.68391266633410e-02, -2.15642540957128e-02, 6.31810399460400e-03, -5.50815969583136e-04, -6.12008018523449e-05, 0.0, 0.0, 0.0, 0.0, 0.0],
&[4.21280289799980e-03, -1.11651391003720e-02, 1.78139108670387e-02, -2.04287733194061e-02, 1.39229889120575e-02, -5.16220078628933e-03, 8.08589544070426e-04, -2.17901509897605e-06], [4.21280289799980e-03, -1.11651391003720e-02, 1.78139108670387e-02, -2.04287733194061e-02, 1.39229889120575e-02, -5.16220078628933e-03, 8.08589544070426e-04, -2.17901509897605e-06, 0.0, 0.0, 0.0, 0.0, 0.0],
&[-4.63425679136442e-03, 1.34385494509234e-02, -2.60541352669482e-02, 3.74164435120410e-02, -3.36394001786921e-02, 1.82199905150995e-02, -5.42550576921788e-03, 6.78314528158575e-04], [-4.63425679136442e-03, 1.34385494509234e-02, -2.60541352669482e-02, 3.74164435120410e-02, -3.36394001786921e-02, 1.82199905150995e-02, -5.42550576921788e-03, 6.78314528158575e-04, 0.0, 0.0, 0.0, 0.0, 0.0],
&[5.41433680102586e-03, -1.75829845731036e-02, 4.26893646894443e-02, -7.75663958780684e-02, 9.06986463369770e-02, -6.74666360933304e-02, 3.07997352073168e-02, -7.82499022484376e-03, 8.38923734582063e-04], [5.41433680102586e-03, -1.75829845731036e-02, 4.26893646894443e-02, -7.75663958780684e-02, 9.06986463369770e-02, -6.74666360933304e-02, 3.07997352073168e-02, -7.82499022484376e-03, 8.38923734582063e-04, 0.0, 0.0, 0.0, 0.0],
&[-2.90718839510249e-03, 1.09922241766815e-02, -3.52053141560315e-02, 8.31962208882433e-02, -1.27244413706703e-01, 1.27078393791315e-01, -8.23947157593937e-02, 3.34105586971409e-02, -7.69529504016692e-03, 7.69529504016692e-04], [-2.90718839510249e-03, 1.09922241766815e-02, -3.52053141560315e-02, 8.31962208882433e-02, -1.27244413706703e-01, 1.27078393791315e-01, -8.23947157593937e-02, 3.34105586971409e-02, -7.69529504016692e-03, 7.69529504016692e-04, 0.0, 0.0, 0.0],
&[8.89941714023592e-04, -4.25818033750240e-03, 1.99225160493104e-02, -6.46588399513875e-02, 1.32772390609268e-01, -1.78508501143865e-01, 1.59999344843642e-01, -9.51030239931668e-02, 3.61804402620956e-02, -8.04009783602125e-03, 8.04009783602125e-04], [8.89941714023592e-04, -4.25818033750240e-03, 1.99225160493104e-02, -6.46588399513875e-02, 1.32772390609268e-01, -1.78508501143865e-01, 1.59999344843642e-01, -9.51030239931668e-02, 3.61804402620956e-02, -8.04009783602125e-03, 8.04009783602125e-04, 0.0, 0.0],
&[-7.63397185185940e-05, 6.56275990492236e-04, -5.83721252683346e-03, 2.90439093207059e-02, -8.47041434795334e-02, 1.57429980520300e-01, -1.95480070025582e-01, 1.65419875422483e-01, -9.49318221478776e-02, 3.55994333054541e-02, -7.91098517898980e-03, 7.91098517898980e-04], [-7.63397185185940e-05, 6.56275990492236e-04, -5.83721252683346e-03, 2.90439093207059e-02, -8.47041434795334e-02, 1.57429980520300e-01, -1.95480070025582e-01, 1.65419875422483e-01, -9.49318221478776e-02, 3.55994333054541e-02, -7.91098517898980e-03, 7.91098517898980e-04, 0.0],
&[-8.51254383837181e-06, -1.77491211003807e-06, 7.32410586432031e-04, -7.40542710012764e-03, 3.44704736857857e-02, -9.39121496781828e-02, 1.66014456295034e-01, -1.99926171069111e-01, 1.66727824560197e-01, -9.52730426058266e-02, 3.57273909771850e-02, -7.93942021715222e-03, 7.93942021715222e-04] [-8.51254383837181e-06, -1.77491211003807e-06, 7.32410586432031e-04, -7.40542710012764e-03, 3.44704736857857e-02, -9.39121496781828e-02, 1.66014456295034e-01, -1.99926171069111e-01, 1.66727824560197e-01, -9.52730426058266e-02, 3.57273909771850e-02, -7.93942021715222e-03, 7.93942021715222e-04]
]; ]);
#[rustfmt::skip] #[rustfmt::skip]
const DISS_DIAG: &'static [Float; 11] = &[ const DISS_DIAG: RowVector<Float, 11> = RowVector::new([[
7.93650793650794e-04, -7.93650793650794e-03, 3.57142857142857e-02, -9.52380952380952e-02, 1.66666666666667e-01, -2.00000000000000e-01, 1.66666666666667e-01, -9.52380952380952e-02, 3.57142857142857e-02, -7.93650793650794e-03, 7.93650793650794e-04 7.93650793650794e-04, -7.93650793650794e-03, 3.57142857142857e-02, -9.52380952380952e-02, 1.66666666666667e-01, -2.00000000000000e-01, 1.66666666666667e-01, -9.52380952380952e-02, 3.57142857142857e-02, -7.93650793650794e-03, 7.93650793650794e-04
]; ]]);
const DISS: BlockMatrix<8, 13, 11> = BlockMatrix::new(
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::flip_ud(super::flip_lr(Self::DISS_BLOCK)),
);
} }
impl SbpOperator1d for Upwind9h2 { impl SbpOperator1d for Upwind9h2 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(&Self::DIFF, super::OperatorType::H2, prev, fut)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::H2,
prev,
fut,
)
} }
fn h(&self) -> &'static [Float] { fn h(&self) -> &'static [Float] {
Self::HBLOCK &Self::H.start
} }
fn is_h2(&self) -> bool { fn is_h2(&self) -> bool {
true true
@ -67,17 +71,11 @@ impl SbpOperator1d for Upwind9h2 {
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diff_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DIFF, super::OperatorType::H2, n)
Self::BLOCK,
Self::DIAG,
super::Symmetry::AntiSymmetric,
super::OperatorType::H2,
n,
)
} }
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn h_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::h_matrix(Self::HBLOCK, n, self.is_h2()) super::h_matrix(&Self::H, n, self.is_h2())
} }
fn upwind(&self) -> Option<&dyn UpwindOperator1d> { fn upwind(&self) -> Option<&dyn UpwindOperator1d> {
@ -86,35 +84,14 @@ impl SbpOperator1d for Upwind9h2 {
} }
impl SbpOperator2d for Upwind9h2 { impl SbpOperator2d for Upwind9h2 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind9h2::BLOCK.len());
let symmetry = super::Symmetry::AntiSymmetric; super::diff_op_2d(&Self::DIFF, OperatorType::H2, prev, fut)
let optype = super::OperatorType::H2;
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
diff_op_row(Upwind9h2::BLOCK, Upwind9h2::DIAG, symmetry, optype)(prev, fut);
}
([1, _], [1, _]) => {
diff_op_col(Upwind9h2::BLOCK, Upwind9h2::DIAG, symmetry, optype)(prev, fut);
}
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
Upwind9h2.diff(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn SbpOperator1d { fn op_xi(&self) -> &dyn SbpOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn SbpOperator1d {
&Self
}
fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> { fn upwind(&self) -> Option<Box<dyn UpwindOperator2d>> {
Some(Box::new(Self)) Some(Box::new(Self))
} }
@ -147,14 +124,7 @@ fn upwind9h2_test() {
impl UpwindOperator1d for Upwind9h2 { impl UpwindOperator1d for Upwind9h2 {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) { fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d( super::diff_op_1d(&Self::DISS, super::OperatorType::H2, prev, fut)
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::H2,
prev,
fut,
)
} }
fn as_sbp(&self) -> &dyn SbpOperator1d { fn as_sbp(&self) -> &dyn SbpOperator1d {
self self
@ -162,54 +132,17 @@ impl UpwindOperator1d for Upwind9h2 {
#[cfg(feature = "sparse")] #[cfg(feature = "sparse")]
fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> { fn diss_matrix(&self, n: usize) -> sprs::CsMat<Float> {
super::sparse_from_block( super::sparse_from_block(&Self::DISS, OperatorType::H2, n)
Self::DISS_BLOCK,
Self::DISS_DIAG,
super::Symmetry::Symmetric,
super::OperatorType::H2,
n,
)
} }
} }
impl UpwindOperator2d for Upwind9h2 { impl UpwindOperator2d for Upwind9h2 {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) { fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape()); assert_eq!(prev.shape(), fut.shape());
assert!(prev.shape()[1] >= 2 * Upwind9h2::BLOCK.len());
let symmetry = super::Symmetry::Symmetric; super::diff_op_2d(&Self::DISS, OperatorType::H2, prev, fut)
let optype = super::OperatorType::H2;
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
diff_op_row(
Upwind9h2::DISS_BLOCK,
Upwind9h2::DISS_DIAG,
symmetry,
optype,
)(prev, fut);
}
([1, _], [1, _]) => {
diff_op_col(
Upwind9h2::DISS_BLOCK,
Upwind9h2::DISS_DIAG,
symmetry,
optype,
)(prev, fut);
}
([_, _], [_, _]) => {
// Fallback, work row by row
for (r0, r1) in prev.outer_iter().zip(fut.outer_iter_mut()) {
Upwind9h2.diss(r0, r1);
}
}
_ => unreachable!("Should only be two elements in the strides vectors"),
}
} }
fn op_xi(&self) -> &dyn UpwindOperator1d { fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self &Self
} }
fn op_eta(&self) -> &dyn UpwindOperator1d {
&Self
}
} }