diff --git a/sbp/examples/multigrid/bin.rs b/sbp/examples/multigrid/bin.rs index f0d7f47..3d5aa2d 100644 --- a/sbp/examples/multigrid/bin.rs +++ b/sbp/examples/multigrid/bin.rs @@ -1,3 +1,4 @@ +#![feature(str_strip)] use sbp::utils::json_to_grids; use sbp::*; use structopt::StructOpt; @@ -143,7 +144,7 @@ impl System { s.scope(|s| { let fields = &self.fnext; - let bt = euler::extract_boundaries( + let bt = euler::extract_boundaries::( &fields, &mut self.bt, &mut self.eb, @@ -197,10 +198,17 @@ fn main() { let vortexparams = utils::json_to_vortex(json["vortex"].clone()); let mut bt = Vec::with_capacity(jgrids.len()); - let determine_bc = |dir| match dir { + let determine_bc = |dir: Option<&String>| match dir { Some(dir) => { if dir == "vortex" { euler::BoundaryCharacteristic::Vortex(vortexparams) + } else if let Some(grid) = dir.strip_prefix("interpolate:") { + euler::BoundaryCharacteristic::Interpolate( + jgrids + .iter() + .position(|other| other.name.as_ref().map_or(false, |name| name == grid)) + .unwrap(), + ) } else { euler::BoundaryCharacteristic::Grid( jgrids @@ -279,7 +287,7 @@ fn main() { bar.inc(1); sys.advance(dt, &pool); } - bar.finish(); + bar.finish_and_clear(); if let Some(timer) = timer { let duration = timer.elapsed(); @@ -288,7 +296,7 @@ fn main() { output.add_timestep(ntime, &sys.fnow); if opt.error { - let time = ntime as f64 * dt; + let time = ntime as Float * dt; let mut e = 0.0; for (fmod, grid) in sys.fnow.iter().zip(&sys.grids) { let mut fvort = fmod.clone(); diff --git a/sbp/src/euler.rs b/sbp/src/euler.rs index b5993fd..99b74ab 100644 --- a/sbp/src/euler.rs +++ b/sbp/src/euler.rs @@ -1,6 +1,6 @@ use super::grid::{Grid, Metrics}; use super::integrate; -use super::operators::{SbpOperator, UpwindOperator}; +use super::operators::{InterpolationOperator, SbpOperator, UpwindOperator}; use super::Float; use ndarray::azip; use ndarray::prelude::*; @@ -582,6 +582,7 @@ pub enum BoundaryCharacteristic { Grid(usize), Vortex(VortexParameters), // Vortices(Vec), + Interpolate(usize), } #[derive(Clone, Debug)] @@ -601,27 +602,35 @@ fn boundary_extractor<'a>( north: match bc.north { BoundaryCharacteristic::This => field.south(), BoundaryCharacteristic::Vortex(_params) => todo!(), - BoundaryCharacteristic::Grid(_) => panic!("Only working on self grid"), + BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_) => { + panic!("Only working on self grid") + } }, south: match bc.south { BoundaryCharacteristic::This => field.north(), BoundaryCharacteristic::Vortex(_params) => todo!(), - BoundaryCharacteristic::Grid(_) => panic!("Only working on self grid"), + BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_) => { + panic!("Only working on self grid") + } }, west: match bc.west { BoundaryCharacteristic::This => field.east(), BoundaryCharacteristic::Vortex(_params) => todo!(), - BoundaryCharacteristic::Grid(_) => panic!("Only working on self grid"), + BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_) => { + panic!("Only working on self grid") + } }, east: match bc.east { BoundaryCharacteristic::This => field.west(), BoundaryCharacteristic::Vortex(_params) => todo!(), - BoundaryCharacteristic::Grid(_) => panic!("Only working on self grid"), + BoundaryCharacteristic::Grid(_) | BoundaryCharacteristic::Interpolate(_) => { + panic!("Only working on self grid") + } }, } } -pub fn extract_boundaries<'a>( +pub fn extract_boundaries<'a, IO: InterpolationOperator>( fields: &'a [Field], bt: &[BoundaryCharacteristics], eb: &'a mut [BoundaryStorage], @@ -641,6 +650,19 @@ pub fn extract_boundaries<'a>( vortexify(field.view_mut(), grid.north(), v, time); field.view() } + BoundaryCharacteristic::Interpolate(g) => { + let to = eb.n.as_mut().unwrap(); + let fine2coarse = field.nx() < fields[g].nx(); + + for (mut to, from) in to.outer_iter_mut().zip(fields[g].south().outer_iter()) { + if fine2coarse { + IO::fine2coarse(from.view(), to.view_mut()); + } else { + IO::coarse2fine(from.view(), to.view_mut()); + } + } + to.view() + } }, south: match bt.south { BoundaryCharacteristic::This => field.north(), @@ -650,6 +672,19 @@ pub fn extract_boundaries<'a>( vortexify(field.view_mut(), grid.south(), v, time); field.view() } + BoundaryCharacteristic::Interpolate(g) => { + let to = eb.s.as_mut().unwrap(); + let fine2coarse = field.nx() < fields[g].nx(); + + for (mut to, from) in to.outer_iter_mut().zip(fields[g].north().outer_iter()) { + if fine2coarse { + IO::fine2coarse(from.view(), to.view_mut()); + } else { + IO::coarse2fine(from.view(), to.view_mut()); + } + } + to.view() + } }, west: match bt.west { BoundaryCharacteristic::This => field.east(), @@ -659,6 +694,19 @@ pub fn extract_boundaries<'a>( vortexify(field.view_mut(), grid.west(), v, time); field.view() } + BoundaryCharacteristic::Interpolate(g) => { + let to = eb.w.as_mut().unwrap(); + let fine2coarse = field.ny() < fields[g].ny(); + + for (mut to, from) in to.outer_iter_mut().zip(fields[g].east().outer_iter()) { + if fine2coarse { + IO::fine2coarse(from.view(), to.view_mut()); + } else { + IO::coarse2fine(from.view(), to.view_mut()); + } + } + to.view() + } }, east: match bt.east { BoundaryCharacteristic::This => field.west(), @@ -668,6 +716,19 @@ pub fn extract_boundaries<'a>( vortexify(field.view_mut(), grid.east(), v, time); field.view() } + BoundaryCharacteristic::Interpolate(g) => { + let to = eb.e.as_mut().unwrap(); + let fine2coarse = field.ny() < fields[g].ny(); + + for (mut to, from) in to.outer_iter_mut().zip(fields[g].west().outer_iter()) { + if fine2coarse { + IO::fine2coarse(from.view(), to.view_mut()); + } else { + IO::coarse2fine(from.view(), to.view_mut()); + } + } + to.view() + } }, }) .collect() @@ -686,19 +747,27 @@ impl BoundaryStorage { pub fn new(bt: &BoundaryCharacteristics, grid: &Grid) -> Self { Self { n: match bt.north { - BoundaryCharacteristic::Vortex(_) => Some(ndarray::Array2::zeros((4, grid.nx()))), + BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_) => { + Some(ndarray::Array2::zeros((4, grid.nx()))) + } _ => None, }, s: match bt.south { - BoundaryCharacteristic::Vortex(_) => Some(ndarray::Array2::zeros((4, grid.nx()))), + BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_) => { + Some(ndarray::Array2::zeros((4, grid.nx()))) + } _ => None, }, e: match bt.east { - BoundaryCharacteristic::Vortex(_) => Some(ndarray::Array2::zeros((4, grid.ny()))), + BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_) => { + Some(ndarray::Array2::zeros((4, grid.ny()))) + } _ => None, }, w: match bt.west { - BoundaryCharacteristic::Vortex(_) => Some(ndarray::Array2::zeros((4, grid.ny()))), + BoundaryCharacteristic::Vortex(_) | BoundaryCharacteristic::Interpolate(_) => { + Some(ndarray::Array2::zeros((4, grid.ny()))) + } _ => None, }, } diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs index 8109d01..f311d19 100644 --- a/sbp/src/operators.rs +++ b/sbp/src/operators.rs @@ -32,6 +32,11 @@ pub trait UpwindOperator: SbpOperator { } } +pub trait InterpolationOperator: Send + Sync { + fn fine2coarse(fine: ArrayView1, coarse: ArrayViewMut1); + fn coarse2fine(coarse: ArrayView1, fine: ArrayViewMut1); +} + #[macro_export] macro_rules! diff_op_1d { ($name: ident, $BLOCK: expr, $DIAG: expr) => { @@ -96,6 +101,9 @@ pub use traditional4::SBP4; mod traditional8; pub use traditional8::SBP8; +mod interpolation; +pub use interpolation::Interpolation4; + #[cfg(test)] pub(crate) mod testing { use super::*; diff --git a/sbp/src/operators/interpolation.rs b/sbp/src/operators/interpolation.rs new file mode 100644 index 0000000..690ec9d --- /dev/null +++ b/sbp/src/operators/interpolation.rs @@ -0,0 +1,59 @@ +use super::InterpolationOperator; +use crate::Float; +use ndarray::{ArrayView1, ArrayView2, ArrayViewMut1}; + +mod interpolation4; +pub use interpolation4::Interpolation4; +mod interpolation9; +pub use interpolation9::Interpolation9; + +fn interpolate( + input: ArrayView1, + mut output: ArrayViewMut1, + block: ArrayView2, + diag: ArrayView2, + jump: (usize, usize), +) { + use ndarray::Axis; + + output.fill(0.0); + let noutput = output.len(); + + for (i, out) in output.iter_mut().enumerate() { + if i < block.len_of(Axis(0)) { + for k in 0..block.len_of(Axis(1)) { + *out += input[k] * block[(i, k)]; + } + } else if noutput - i <= block.len_of(Axis(0)) { + let row = noutput - i - 1; + let index = input.len() - block.len_of(Axis(1)); + + for k in 0..block.len_of(Axis(1)) { + let col = block.len_of(Axis(1)) - k - 1; + *out += input[index + k] * block[(row, col)]; + } + } else { + let n = i - block.len_of(Axis(0)); + let index = jump.0 + jump.1 * (n / diag.len_of(Axis(0))); + let row = n % diag.len_of(Axis(0)); + + for k in 0..diag.len_of(Axis(1)) { + *out += input[index + k] * diag[(row, k)]; + } + } + } +} + +#[cfg(test)] +pub(crate) fn test_interpolation_operator() { + let x_c = ndarray::Array1::linspace(0.0, 1.0, 101); + let x_f = ndarray::Array1::linspace(0.0, 1.0, 2 * x_c.len() - 1); + + let mut ix_f = ndarray::Array1::zeros(x_f.raw_dim()); + IO::coarse2fine(x_c.view(), ix_f.view_mut()); + approx::assert_abs_diff_eq!(ix_f, x_f, epsilon = 1e-2); + + let mut ix_c = ndarray::Array1::zeros(x_c.raw_dim()); + IO::fine2coarse(x_f.view(), ix_c.view_mut()); + approx::assert_abs_diff_eq!(ix_c, x_c, epsilon = 1e-2); +} diff --git a/sbp/src/operators/interpolation/interpolation4.rs b/sbp/src/operators/interpolation/interpolation4.rs new file mode 100644 index 0000000..e01b67e --- /dev/null +++ b/sbp/src/operators/interpolation/interpolation4.rs @@ -0,0 +1,99 @@ +use super::*; + +pub struct Interpolation4 {} + +impl Interpolation4 { + const F2C_DIAG: &'static [[Float; 7]] = &[[ + -1.0 / 32.0, + 0.0, + 9.0 / 32.0, + 1.0 / 2.0, + 9.0 / 32.0, + 0.0, + -1.0 / 32.0, + ]]; + const F2C_BLOCK: &'static [[Float; 10]] = &[ + [ + 0.15745752473496450946063e23 / 0.30759034429175640424160e23, + 0.74943057537715447001883e23 / 0.104580717059197177442144e24, + 0.1733236963142549444853e22 / 0.522903585295985887210720e24, + -0.2994697914858864417917e22 / 0.26145179264799294360536e23, + -0.641055020554612856232e21 / 0.16340737040499558975335e23, + -0.31200598857765495755439e23 / 0.522903585295985887210720e24, + -0.422055004077148591739e21 / 0.6536294816199823590134e22, + -0.16670527152025920991407e23 / 0.522903585295985887210720e24, + 0.53414892210822113304e20 / 0.3268147408099911795067e22, + 0.32274131131983102388471e23 / 0.522903585295985887210720e24, + ], + [ + -0.6225999401446722477711e22 / 0.907391515660681392512720e24, + 0.1036706267989606133475e22 / 0.3075903442917564042416e22, + 0.328926383150495585114867e24 / 0.907391515660681392512720e24, + 0.89897263376416194502629e23 / 0.362956606264272557005088e24, + 0.1282110041109225712464e22 / 0.56711969728792587032045e23, + 0.24913624504973179743933e23 / 0.1814783031321362785025440e25, + 0.1036353973336644936129e22 / 0.22684787891517034812818e23, + 0.41992032737257490852109e23 / 0.1814783031321362785025440e25, + -0.106829784421644226608e21 / 0.11342393945758517406409e23, + -0.65509482089877943540197e23 / 0.1814783031321362785025440e25, + ], + [ + 0.6225999401446722477711e22 / 0.1322638480454552538238880e25, + -0.15796094028352692249389e23 / 0.264527696090910507647776e24, + 0.40307836352152312671e20 / 0.30759034429175640424160e23, + 0.41110375639957470563661e23 / 0.132263848045455253823888e24, + 0.22428220801327117461888e23 / 0.41332452514204766819965e23, + 0.366744409069694352232131e24 / 0.1322638480454552538238880e25, + -0.806542934441844097041e21 / 0.16532981005681906727986e23, + -0.80111035662200679366237e23 / 0.1322638480454552538238880e25, + 0.53414892210822113304e20 / 0.8266490502840953363993e22, + 0.34196570783806579914981e23 / 0.1322638480454552538238880e25, + ], + ]; + + const C2F_DIAG: &'static [[Float; 4]] = &[ + [0.0, 1.0, 0.0, 0.0], + [-1.0 / 16.0, 9.0 / 16.0, 9.0 / 16.0, -1.0 / 16.0], + ]; + #[rustfmt::skip] + const C2F_BLOCK: &'static [[Float; 7]] = &[ + [2305422334878123.0/2251799813685248.0, -116682606115696273.0/2449958197289549824.0, 116682606115696271.0/4899916394579099648.0, 0.0, 0.0, 0.0, 0.0], + [13716032226970759.0/33214047251857408.0, 6071594962398535.0/9007199254740992.0, -370046899066617177.0/4251398048237748224.0, 0.0, 0.0, 0.0, 0.0], + [129931657769535863.0/49575624698094419968.0, 192639692900606809.0/193654783976931328.0, 377708307469581.0/144115188075855872.0, 0.0, 0.0, 0.0, 0.0], + [-70155234265236545.0/882705526964617216.0, 263247068288746603.0/441352763482308608.0, 30095982505478077.0/55169095435288576.0, -1.0/16.0, 0.0, 0.0, 0.0], + [-48056522745236605.0/1729382256910270464.0, 64075363660315477.0/1152921504606846976.0, 52541429192657303.0/54043195528445952.0, 0.0, 0.0, 0.0, 0.0], + [-1522750502383427.0/36028797018963968.0, 7295481936885541.0/216172782113783808.0, 71596091562607345.0/144115188075855872.0, 147.0/256.0, -1.0/16.0, 0.0, 0.0], + [-39549054396350587.0/864691128455135232.0, 64741671534788765.0/576460752303423488.0, -25192617138438183.0/288230376151711744.0, 49.0/48.0, 0.0, 0.0, 0.0], + [-19526589507951803.0/864691128455135232.0, 131163408473299267.0/2305843009213693952.0, -15639330559927497.0/144115188075855872.0, 147.0/256.0, 9.0/16.0, -1.0/16.0, 0.0], + [80084669807060717.0/6917529027641081856.0, -320338679228242837.0/13835058055282163712.0, 106779559742747605.0/9223372036854775808.0, 0.0, 1.0, 0.0, 0.0], + [37803466236726359.0/864691128455135232.0, -76732832380295345.0/864691128455135232.0, 53407021400548805.0/1152921504606846976.0, -49.0/768.0, 9.0/16.0, 9.0/16.0, -1.0/16.0], + ]; +} + +impl InterpolationOperator for Interpolation4 { + fn fine2coarse(fine: ArrayView1, coarse: ArrayViewMut1) { + assert_eq!(fine.len(), 2 * coarse.len() - 1); + super::interpolate( + fine, + coarse, + ndarray::arr2(Self::F2C_BLOCK).view(), + ndarray::arr2(Self::F2C_DIAG).view(), + (3, 2), + ) + } + fn coarse2fine(coarse: ArrayView1, fine: ArrayViewMut1) { + assert_eq!(fine.len(), 2 * coarse.len() - 1); + super::interpolate( + coarse, + fine, + ndarray::arr2(Self::C2F_BLOCK).view(), + ndarray::arr2(Self::C2F_DIAG).view(), + (4, 1), + ) + } +} + +#[test] +fn test_inter4() { + test_interpolation_operator::(); +} diff --git a/sbp/src/operators/interpolation/interpolation9.rs b/sbp/src/operators/interpolation/interpolation9.rs new file mode 100644 index 0000000..9f99291 --- /dev/null +++ b/sbp/src/operators/interpolation/interpolation9.rs @@ -0,0 +1,82 @@ +use super::*; + +pub struct Interpolation9 {} + +impl Interpolation9 { + #[rustfmt::skip] + const F2C_DIAG: &'static [[Float; 15]] = &[ + [-5.0/4096.0, 0.0, 49.0/4096.0, 0.0, -245.0/4096.0, 0.0, 1225.0/4096.0, 1.0/2.0, 1225.0/4096.0, 0.0, -245.0/4096.0, 0.0, 49.0/4096.0, 0.0, -5.0/4096.0], + ]; + #[rustfmt::skip] + const F2C_BLOCK: &'static [[Float; 22]] = &[ + [5.2686284188785748e-01, 8.4682941042184923e-01, -2.0727251976345169e-03, -3.5467776280365609e-01, -1.5291230463645894e-02, -9.0976491902894502e-02, -2.9298328682108545e-02, -5.2179807595463586e-03, 3.4062250834683325e-02, 7.3813898085096877e-02, 6.5912558521007691e-02, -2.5455425182224583e-02, -7.8858135268769738e-02, 1.0511170158482339e-02, 1.0297466246544207e-01, 5.6949108895613200e-02, -3.4308921157670390e-02, -5.8997482788315750e-02, -4.2966539044152929e-02, -2.5124065411472413e-02, 1.7139909678724394e-03, 4.3615196424187055e-02], + [-2.3468446435194214e-02, 4.3131420786340957e-01, 8.5980850312858598e-02, 4.7603940875129624e-01, 1.3210073514275335e-02, 6.0881411955347041e-02, 2.5585810411810836e-02, 1.0600733302542303e-02, -2.8158836917230867e-02, -6.7994087660688027e-02, -5.4799898209469740e-02, 5.0238370377932259e-02, 1.0850940737329789e-01, 1.3660999523576111e-05, -1.1424280311019673e-01, -7.0594704953755610e-02, 2.6209331260080626e-02, 5.0182418961172676e-02, 3.2103343148941477e-02, 1.7470268584954918e-02, -7.0745228443207762e-04, -2.8373067246476109e-02], + [2.4916714306945922e-01, -4.7775189410021041e-01, -1.8369140990199855e-02, 4.1977084808920612e-01, 6.6316613209558584e-01, 3.7116423563000828e-01, -2.7151118047839545e-01, -2.9085560694898455e-01, 2.7796013161440170e-01, 8.1515123444641679e-01, 5.4525214595381100e-01, -1.0040515073752119e+00, -1.8219288006779970e+00, -1.6683214390616796e-01, 1.6046462777797821e+00, 1.0846974719520888e+00, -2.2882651660792802e-01, -5.2721786590791042e-01, -2.6947479334214830e-01, -1.2454712416724613e-01, -2.8643203612198978e-03, 1.7325527423286030e-01], + [-3.1374744107225773e-02, -3.8560277245237022e-02, 2.2373357921453515e-03, 1.6624250373884275e-01, 1.7101480843281450e-02, 3.4736941670549754e-01, 2.9074915472830576e-01, 1.8037070362251098e-01, -3.1643970637017987e-02, -1.4439700205357936e-01, -6.2814859725063579e-02, 2.3466168189679312e-01, 3.7917743408128213e-01, 6.1631222196174403e-02, -2.8284935028445379e-01, -2.0614120963189675e-01, 2.0913586845593374e-02, 6.8695794722380743e-02, 2.2528020642632283e-02, 5.4061692073344636e-03, 2.0185034965858510e-03, -1.3215948348859834e-03], + [5.9296747167331777e-02, 1.4378307696751280e-01, -4.0409834787451700e-03, -3.7929003879766199e-01, -3.1546591422357252e-02, -4.2678477551108618e-01, -6.4519096894777528e-02, 6.8916457275965848e-01, 1.2653355304352822e+00, 1.0015090142279532e+00, 1.0424390935875409e-01, -9.7829622252487647e-01, -1.3208748944090689e+00, -2.8300008595522019e-01, 8.1422088300754603e-01, 6.3440130137928941e-01, -1.9955760198721691e-02, -1.5064386074684394e-01, -1.4326889571202749e-02, 1.7711792232011369e-02, -7.9217161324989811e-03, -4.8465911892278418e-02], + [-3.2260941454839102e-03, -9.3712947419508016e-03, 2.0590133951067365e-04, 2.1085040742079481e-02, 1.6586905667752524e-03, 2.3510098029135212e-02, 3.5062071002191428e-03, -4.1732089994723873e-02, -2.1826444408900507e-03, 2.0649123576151260e-01, 3.8634802477549085e-01, 3.4657285748507582e-01, 1.7693083329005391e-01, 1.2957209274305518e-02, -7.7913009746446549e-02, -5.8609417199037818e-02, -3.7071852226009880e-04, 1.1759795657464859e-02, -1.3229500781068283e-03, -3.4952049532427699e-03, 7.3662735648173173e-04, 6.4609024440376705e-03], + ]; + + #[rustfmt::skip] + const C2F_DIAG: &'static [[Float; 8]] = &[ + [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0], + [-5.0/2048.0, 49.0/2048.0, -245.0/2048.0, 1225.0/2048.0, 1225.0/2048.0, -245.0/2048.0, 49.0/2048.0, -5.0/2048.0], + ]; + #[rustfmt::skip] + const C2F_BLOCK: &'static [[Float; 15]] = &[ + [1.0537256837757150e+00, -2.4288846421921273e-01, 4.3429701911970164e-01, -3.8281710980097783e-01, 1.6566860024112703e-01, -2.7985729116353025e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [3.2729048243799186e-01, 8.6262841572681914e-01, -1.6091848728719577e-01, -9.0919856879246821e-02, 7.7629100522844718e-02, -1.5709654521213201e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-4.7566986935300764e-03, 1.0210763172289250e+00, -3.6738281980399710e-02, 3.1323929502949190e-02, -1.2954788512749336e-02, 2.0495224548048151e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-1.1627405111828265e-01, 8.0757767164055649e-01, 1.1992982462060080e-01, 3.3248500747768550e-01, -1.7369991978798591e-01, 2.9981467167425814e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-2.1892385771606702e-02, 9.7869742354703762e-02, 8.2744488829725205e-01, 1.4937073869608852e-01, -6.3093182844714504e-02, 1.0300199268276950e-02, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-4.1949770425397409e-02, 1.4527057994032164e-01, 1.4915321003780291e-01, 9.7717676895352013e-01, -2.7490881010963303e-01, 4.7020196058270423e-02, -1.7621744548846038e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-1.8716920677698350e-02, 8.4583038536706734e-02, -1.5116294736984345e-01, 1.1331598176662732e+00, -5.7578343981351678e-02, 9.7153558259133310e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-3.0487057786247208e-03, 3.2050993871316637e-02, -1.4810051238794766e-01, 6.4292482016111874e-01, 5.6249066806053438e-01, -1.0575788798746898e-01, 2.1882030311071686e-02, -2.4414062500000000e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [2.0087735588280062e-02, -8.5933975772489596e-02, 1.4285854720715777e-01, -1.1384914286933631e-01, 1.0424198692657576e+00, -5.5830334193693381e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [4.3530713066204312e-02, -2.0750154911869487e-01, 4.1894972646178857e-01, -5.1951365728643595e-01, 8.2507198333460330e-01, 5.2818839773705129e-01, -1.1043385164237317e-01, 2.4149643697856386e-02, -2.4414062500000000e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [3.8871008670068943e-02, -1.6723606656444842e-01, 2.8023417955710428e-01, -2.2599622598531163e-01, 8.5879136206759521e-02, 9.8824796811582738e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-1.5011980647711863e-02, 1.5331538428225466e-01, -5.1603565889716185e-01, 8.4426924973105621e-01, -8.0594861667778428e-01, 8.8650620748676312e-01, 5.5216925821186580e-01, -1.2074821848928192e-01, 2.3925781250000000e-02, -2.4414062500000000e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-4.6505481330402991e-02, 3.3114452886362927e-01, -9.3638645250322083e-01, 1.3642101479846502e+00, -1.0881747976148743e+00, 4.5257520495295200e-01, 9.2313684964726628e-01, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [6.1988154538518501e-03, 4.1690074257599231e-05, -8.5743943088044547e-02, 2.2173771747897486e-01, -2.3314362515542061e-01, 3.3143525828105390e-02, 5.5216925821186580e-01, 6.0374109244640961e-01, -1.1962890625000000e-01, 2.3925781250000000e-02, -2.4414062500000000e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [6.0727865634526521e-02, -3.4864146923076744e-01, 8.2471336696228448e-01, -1.0176395516358734e+00, 6.7077862432757074e-01, -1.9929537265386069e-01, 0.0000000000000000e+00, 1.0093565365961199e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [3.3584939733883018e-02, -2.1543800558928278e-01, 5.5748392441156636e-01, -7.4165787523468796e-01, 5.2263807167283938e-01, -1.4991829579835428e-01, -1.1043385164237317e-01, 6.0374109244640961e-01, 5.9814453125000000e-01, -1.1962890625000000e-01, 2.3925781250000000e-02, -2.4414062500000000e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-2.0233205958094688e-02, 7.9984554906766037e-02, -1.1760616004611729e-01, 7.5243210278702510e-02, -1.6440130255643871e-02, -9.4826892561270318e-04, 0.0000000000000000e+00, 0.0000000000000000e+00, 1.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-3.4792939561676173e-02, 1.5314463405892736e-01, -2.7096540049760082e-01, 2.4715474039537103e-01, -1.2410475312535753e-01, 3.0080635641143383e-02, 2.2086770328474635e-02, -1.2074821848928192e-01, 5.9814453125000000e-01, 5.9814453125000000e-01, -1.1962890625000000e-01, 2.3925781250000000e-02, -2.4414062500000000e-03, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-2.5338914907631935e-02, 9.7971657014317953e-02, -1.3849747898095244e-01, 8.1051643933268969e-02, -1.1802904442792750e-02, -3.3840026162097964e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 1.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [-1.4816565861655355e-02, 5.3315044287206965e-02, -6.4011414565182903e-02, 1.9450395069625474e-02, 1.4591484787125306e-02, -8.9404603406417257e-03, -2.2537520743341464e-03, 2.4149643697856386e-02, -1.1962890625000000e-01, 5.9814453125000000e-01, 5.9814453125000000e-01, -1.1962890625000000e-01, 2.3925781250000000e-02, -2.4414062500000000e-03, 0.0000000000000000e+00], + [1.0108021789406768e-03, -2.1589736695899393e-03, -1.4721271110470103e-03, 7.2622015612738991e-03, -6.5261380057503940e-03, 1.8842350461727676e-03, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 1.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00, 0.0000000000000000e+00], + [2.5721451516875750e-02, -8.6587755045305664e-02, 8.9045132584768621e-02, -4.7548533304570878e-03, -3.9927614710396069e-02, 1.6526482091438560e-02, 0.0000000000000000e+00, -2.4642493569241209e-03, 2.3925781250000000e-02, -1.1962890625000000e-01, 5.9814453125000000e-01, 5.9814453125000000e-01, -1.1962890625000000e-01, 2.3925781250000000e-02, -2.4414062500000000e-03], + ]; +} + +impl InterpolationOperator for Interpolation9 { + fn fine2coarse(fine: ArrayView1, coarse: ArrayViewMut1) { + assert_eq!(fine.len(), 2 * coarse.len() - 1); + use ndarray::prelude::*; + use std::iter::FromIterator; + let block = Array::from_iter(Self::F2C_BLOCK.iter().flatten().copied()) + .into_shape((Self::F2C_BLOCK.len(), Self::F2C_BLOCK[0].len())) + .unwrap(); + let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied()) + .into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len())) + .unwrap(); + super::interpolate(fine, coarse, block.view(), diag.view(), (5, 2)) + } + fn coarse2fine(coarse: ArrayView1, fine: ArrayViewMut1) { + assert_eq!(fine.len(), 2 * coarse.len() - 1); + use ndarray::prelude::*; + use std::iter::FromIterator; + let block = Array::from_iter(Self::C2F_BLOCK.iter().flatten().copied()) + .into_shape((Self::C2F_BLOCK.len(), Self::C2F_BLOCK[0].len())) + .unwrap(); + let diag = Array::from_iter(Self::C2F_DIAG.iter().flatten().copied()) + .into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len())) + .unwrap(); + super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1)) + } +} + +#[test] +fn test_inter9() { + test_interpolation_operator::(); +}