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

working d1 SBP4

This commit is contained in:
Magnus Ulimoen 2021-01-28 21:27:25 +01:00
parent 94e8fb5b7c
commit 30c563c19d
3 changed files with 36 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
sbp/src/lib.rs
View File

@ -1,5 +1,6 @@
#![feature(core_intrinsics)]
#![feature(array_windows)]
#![feature(array_chunks)]
/// Type used for floats, configure with the `f32` feature
#[cfg(feature = "f32")]

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

@ -114,7 +114,7 @@ pub(crate) mod constmatrix {
let mut v = Self::default();
for i in 0..M {
for j in 0..N {
v[(i, j)] = self[(N - 1 - i, M - 1 - j)].clone()
v[(i, j)] = self[(M - 1 - i, N - 1 - j)].clone()
}
}
v
@ -276,12 +276,11 @@ pub(crate) fn diff_op_1d_slice_matrix<const M: usize, const N: usize, const D: u
for (window, f) in prev
.array_windows::<D>()
.skip(window_elems_to_skip)
.zip(fut.iter_mut().skip(M))
.zip(fut.array_chunks_mut::<1>().skip(M))
.take(nx - 2 * M)
{
let fut = ColVector::<Float, 1>::map_to_col_mut(unsafe {
std::mem::transmute::<&mut Float, &mut [Float; 1]>(f)
});
// impl From here?
let fut = ColVector::<Float, 1>::map_to_col_mut(f);
let prev = ColVector::<_, D>::map_to_col(window);
diag.matmul_into(prev, fut);
@ -297,7 +296,8 @@ 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[nx - N..]).try_into().unwrap());
let prev = prev.array_windows::<N>().last().unwrap();
let prev = ColVector::<_, N>::map_to_col(prev);
let fut = ColVector::<_, M>::map_to_col_mut((&mut fut[nx - M..]).try_into().unwrap());
flipped.matmul_into(prev, fut);

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

@ -10,10 +10,14 @@ impl SBP4 {
const HBLOCK: &'static [Float] = &[
17.0 / 48.0, 59.0 / 48.0, 43.0 / 48.0, 49.0 / 48.0,
];
#[rustfmt::skip]
const DIAG: &'static [Float] = &[
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]
const BLOCK: &'static [&'static [Float]] = &[
&[-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0],
@ -21,6 +25,13 @@ impl SBP4 {
&[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],
[-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],
[3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0]
]);
#[rustfmt::skip]
const D2DIAG: &'static [Float] = &[
@ -71,6 +82,22 @@ 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::DIAG_MATRIX,
super::Symmetry::AntiSymmetric,
super::OperatorType::Normal,
p.as_slice().unwrap(),
f.as_slice_mut().unwrap(),
)
}
}
impl SbpOperator2d for SBP4 {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
@ -81,7 +108,8 @@ impl SbpOperator2d for SBP4 {
match (prev.strides(), fut.strides()) {
([_, 1], [_, 1]) => {
diff_op_row(SBP4::BLOCK, SBP4::DIAG, symmetry, optype)(prev, fut);
//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);