From 78cc4091ba81af26381e6fa6ae191a34b84c6059 Mon Sep 17 00:00:00 2001 From: RussellBentley Date: Sat, 28 Dec 2024 00:05:26 -0500 Subject: [PATCH] Create Domain trait, add owned and slice implementations, update all the things, some misc cleanup as well --- examples/gen_1d.rs | 7 +- examples/gen_2d.rs | 9 +-- examples/heat_1d_ap_direct.rs | 8 +-- examples/heat_1d_p_direct.rs | 9 +-- examples/heat_1d_p_fft.rs | 10 +-- examples/heat_2d_p_direct.rs | 8 +-- examples/heat_2d_p_fft.rs | 8 +-- src/domain/bc/periodic.rs | 39 ++++++----- src/domain/bc_check.rs | 0 src/domain/domain_view.rs | 123 ---------------------------------- src/domain/gather_args.rs | 24 +++---- src/domain/mod.rs | 4 +- src/domain/view/chunk.rs | 34 ++++++++++ src/domain/view/mod.rs | 69 +++++++++++++++++++ src/domain/view/owned.rs | 43 ++++++++++++ src/domain/view/slice.rs | 42 ++++++++++++ src/image.rs | 7 +- src/par_stencil.rs | 42 ++++++------ src/solver/direct.rs | 33 +++++---- src/solver/periodic_direct.rs | 34 +++++----- src/solver/periodic_plan.rs | 32 ++++----- src/solver/trapezoid.rs | 16 ++--- tests/base_solver_compare.rs | 67 ++++++------------ 23 files changed, 337 insertions(+), 331 deletions(-) delete mode 100644 src/domain/bc_check.rs delete mode 100644 src/domain/domain_view.rs create mode 100644 src/domain/view/chunk.rs create mode 100644 src/domain/view/mod.rs create mode 100644 src/domain/view/owned.rs create mode 100644 src/domain/view/slice.rs diff --git a/examples/gen_1d.rs b/examples/gen_1d.rs index d490e78..726b030 100644 --- a/examples/gen_1d.rs +++ b/examples/gen_1d.rs @@ -12,11 +12,8 @@ fn main() { // Create domains let grid_bound = args.grid_bounds(); let buffer_size = grid_bound.buffer_size(); - let mut grid_input = vec![0.0; buffer_size]; - let mut input_domain = Domain::new(grid_bound, &mut grid_input); - - let mut grid_output = vec![0.0; buffer_size]; - let mut output_domain = Domain::new(grid_bound, &mut grid_output); + let mut input_domain = OwnedDomain::new(grid_bound); + let mut output_domain = OwnedDomain::new(grid_bound); // Fill in with IC values (use normal dist for spike in the middle) let n_f = buffer_size as f64; diff --git a/examples/gen_2d.rs b/examples/gen_2d.rs index a14509b..9c1237a 100644 --- a/examples/gen_2d.rs +++ b/examples/gen_2d.rs @@ -19,13 +19,8 @@ fn main() { let chunk_size = 1000; // Create domains - let buffer_size = grid_bound.buffer_size(); - let mut grid_input = vec![0.0; buffer_size]; - let mut input_domain: Domain<2> = Domain::new(grid_bound, &mut grid_input); - - let mut grid_output = vec![0.0; buffer_size]; - let mut output_domain: Domain<2> = - Domain::new(grid_bound, &mut grid_output); + let mut input_domain = OwnedDomain::new(grid_bound); + let mut output_domain = OwnedDomain::new(grid_bound); // Fill in with IC values (use normal dist for spike in the middle) let width_f = grid_bound.bounds[(0, 1)] as f64 + 1.0; diff --git a/examples/heat_1d_ap_direct.rs b/examples/heat_1d_ap_direct.rs index 3d4ee96..9daf299 100644 --- a/examples/heat_1d_ap_direct.rs +++ b/examples/heat_1d_ap_direct.rs @@ -9,12 +9,8 @@ fn main() { // Create domains let grid_bound = args.grid_bounds(); - let buffer_size = grid_bound.buffer_size(); - let mut grid_input = vec![0.0; buffer_size]; - let mut input_domain = Domain::new(grid_bound, &mut grid_input); - - let mut grid_output = vec![0.0; buffer_size]; - let mut output_domain = Domain::new(grid_bound, &mut grid_output); + let mut input_domain = OwnedDomain::new(grid_bound); + let mut output_domain = OwnedDomain::new(grid_bound); // Create BC let bc = ConstantCheck::new(1.0, grid_bound); diff --git a/examples/heat_1d_p_direct.rs b/examples/heat_1d_p_direct.rs index b906b68..cb9ac8f 100644 --- a/examples/heat_1d_p_direct.rs +++ b/examples/heat_1d_p_direct.rs @@ -10,15 +10,12 @@ fn main() { // Create domains let grid_bound = args.grid_bounds(); - let buffer_size = grid_bound.buffer_size(); - let mut grid_input = vec![0.0; buffer_size]; - let mut input_domain = Domain::new(grid_bound, &mut grid_input); + let mut input_domain = OwnedDomain::new(grid_bound); - let mut grid_output = vec![0.0; buffer_size]; - let mut output_domain = Domain::new(grid_bound, &mut grid_output); + let mut output_domain = OwnedDomain::new(grid_bound); // Fill in with IC values (use normal dist for spike in the middle) - let n_f = buffer_size as f64; + let n_f = grid_bound.buffer_size() as f64; let sigma_sq: f64 = (n_f / 25.0) * (n_f / 25.0); let ic_gen = |world_coord: Coord<1>| { let x = (world_coord[0] as f64) - (n_f / 2.0); diff --git a/examples/heat_1d_p_fft.rs b/examples/heat_1d_p_fft.rs index 529a353..00b9f2f 100644 --- a/examples/heat_1d_p_fft.rs +++ b/examples/heat_1d_p_fft.rs @@ -1,4 +1,3 @@ -use fftw::array::AlignedVec; use nhls::domain::*; use nhls::image_1d_example::*; use nhls::util::*; @@ -10,15 +9,12 @@ fn main() { // Create domains let grid_bound = args.grid_bounds(); - let buffer_size = grid_bound.buffer_size(); - let mut grid_input = AlignedVec::new(buffer_size); - let mut input_domain = Domain::new(grid_bound, &mut grid_input); + let mut input_domain = OwnedDomain::new(grid_bound); - let mut grid_output = AlignedVec::new(buffer_size); - let mut output_domain = Domain::new(grid_bound, &mut grid_output); + let mut output_domain = OwnedDomain::new(grid_bound); // Fill in with IC values (use normal dist for spike in the middle) - let n_f = buffer_size as f64; + let n_f = grid_bound.buffer_size() as f64; let sigma_sq: f64 = (n_f / 25.0) * (n_f / 25.0); let ic_gen = |world_coord: Coord<1>| { let x = (world_coord[0] as f64) - (n_f / 2.0); diff --git a/examples/heat_2d_p_direct.rs b/examples/heat_2d_p_direct.rs index 3da132f..ce05029 100644 --- a/examples/heat_2d_p_direct.rs +++ b/examples/heat_2d_p_direct.rs @@ -1,4 +1,3 @@ -use fftw::array::AlignedVec; use nhls::domain::*; use nhls::image::*; use nhls::util::*; @@ -20,11 +19,8 @@ fn main() { let stencil = nhls::standard_stencils::heat_2d(1.0, 1.0, 1.0, 0.2, 0.2); // Create domains - let buffer_size = grid_bound.buffer_size(); - let mut input_buffer = AlignedVec::new(buffer_size); - let mut output_buffer = AlignedVec::new(buffer_size); - let mut input_domain = Domain::new(grid_bound, &mut input_buffer); - let mut output_domain = Domain::new(grid_bound, &mut output_buffer); + let mut input_domain = OwnedDomain::new(grid_bound); + let mut output_domain = OwnedDomain::new(grid_bound); // Fill in with IC values (use normal dist for spike in the middle) // Write out initial frame diff --git a/examples/heat_2d_p_fft.rs b/examples/heat_2d_p_fft.rs index 85907ce..4d11f0a 100644 --- a/examples/heat_2d_p_fft.rs +++ b/examples/heat_2d_p_fft.rs @@ -1,4 +1,3 @@ -use fftw::array::AlignedVec; use nhls::domain::*; use nhls::image::*; use nhls::util::*; @@ -20,11 +19,8 @@ fn main() { let stencil = nhls::standard_stencils::heat_2d(1.0, 1.0, 1.0, 0.2, 0.2); // Create domains - let buffer_size = grid_bound.buffer_size(); - let mut input_buffer = AlignedVec::new(buffer_size); - let mut output_buffer = AlignedVec::new(buffer_size); - let mut input_domain = Domain::new(grid_bound, &mut input_buffer); - let mut output_domain = Domain::new(grid_bound, &mut output_buffer); + let mut input_domain = OwnedDomain::new(grid_bound); + let mut output_domain = OwnedDomain::new(grid_bound); // Fill in with IC values (use normal dist for spike in the middle) // Write out initial frame diff --git a/src/domain/bc/periodic.rs b/src/domain/bc/periodic.rs index 35306f8..a7cfdf9 100644 --- a/src/domain/bc/periodic.rs +++ b/src/domain/bc/periodic.rs @@ -1,19 +1,30 @@ -use crate::domain::bc::BCCheck; -use crate::domain::Domain; +use crate::domain::*; use crate::util::*; -pub struct PeriodicCheck<'a, const GRID_DIMENSION: usize> { - domain: &'a Domain<'a, GRID_DIMENSION>, +pub struct PeriodicCheck< + 'a, + const GRID_DIMENSION: usize, + DomainType: DomainView, +> { + domain: &'a DomainType, } -impl<'a, const GRID_DIMENSION: usize> PeriodicCheck<'a, GRID_DIMENSION> { - pub fn new(domain: &'a Domain<'a, GRID_DIMENSION>) -> Self { +impl< + 'a, + const GRID_DIMENSION: usize, + DomainType: DomainView, + > PeriodicCheck<'a, GRID_DIMENSION, DomainType> +{ + pub fn new(domain: &'a DomainType) -> Self { PeriodicCheck { domain } } } -impl BCCheck - for PeriodicCheck<'_, GRID_DIMENSION> +impl< + const GRID_DIMENSION: usize, + DomainType: DomainView + Sync, + > BCCheck + for PeriodicCheck<'_, GRID_DIMENSION, DomainType> { fn check(&self, world_coord: &Coord) -> Option { let p_coord = &self.domain.aabb().periodic_coord(world_coord); @@ -34,14 +45,11 @@ mod unit_tests { fn periodic_check_test() { { let aabb = AABB::new(matrix![0, 10]); - let n_r = aabb.buffer_size(); - let mut buffer = fftw::array::AlignedVec::new(n_r); - for i in 0..n_r { - buffer.as_slice_mut()[i] = i as f64; - } - let domain = Domain::new(aabb, buffer.as_slice_mut()); + + let mut domain = OwnedDomain::new(aabb); + domain.par_set_values(|coord| coord[0] as f64, 1); let bc = PeriodicCheck::new(&domain); - for i in 0..n_r { + for (i, _) in domain.buffer().iter().enumerate() { let v = bc.check(&vector![i as i32]); assert_eq!(v, None); } @@ -70,5 +78,6 @@ mod unit_tests { buffer_a[i] = i as f64; buffer_b[i] = (n_r - i) as f64; } + // TODO: Not sure what I wanted to do here, but clearly didn't finish } } diff --git a/src/domain/bc_check.rs b/src/domain/bc_check.rs deleted file mode 100644 index e69de29..0000000 diff --git a/src/domain/domain_view.rs b/src/domain/domain_view.rs deleted file mode 100644 index 8bb2647..0000000 --- a/src/domain/domain_view.rs +++ /dev/null @@ -1,123 +0,0 @@ -use crate::util::*; -use rayon::prelude::*; - -pub struct Domain<'a, const GRID_DIMENSION: usize> { - aabb: AABB, - buffer: &'a mut [f64], -} - -impl<'a, const GRID_DIMENSION: usize> Domain<'a, GRID_DIMENSION> { - pub fn aabb(&self) -> &AABB { - &self.aabb - } - - pub fn set_aabb(&mut self, aabb: AABB) { - debug_assert!(aabb.buffer_size() <= self.buffer.len()); - // TODO: should we re-slice here? - self.aabb = aabb; - } - - pub fn buffer(&self) -> &[f64] { - self.buffer - } - - pub fn buffer_mut(&mut self) -> &mut [f64] { - self.buffer - } - - pub fn new(aabb: AABB, buffer: &'a mut [f64]) -> Self { - debug_assert_eq!(buffer.len(), aabb.buffer_size()); - Domain { aabb, buffer } - } - - pub fn view(&self, world_coord: &Coord) -> f64 { - debug_assert!(self.aabb.contains(world_coord)); - let index = self.aabb.coord_to_linear(world_coord); - self.buffer[index] - } - - pub fn modify(&mut self, world_coord: &Coord, value: f64) { - debug_assert!(self.aabb.contains(world_coord)); - let index = self.aabb.coord_to_linear(world_coord); - self.buffer[index] = value; - } - - pub fn par_modify_access( - &mut self, - chunk_size: usize, - ) -> impl ParallelIterator> { - par_modify_access_impl(self.buffer, &self.aabb, chunk_size) - } - - pub fn par_set_values< - F: FnOnce(Coord) -> f64 + Send + Sync + Copy, - >( - &mut self, - f: F, - chunk_size: usize, - ) { - self.par_modify_access(chunk_size).for_each( - |mut d: DomainChunk<'_, GRID_DIMENSION>| { - d.coord_iter_mut().for_each( - |(world_coord, value_mut): ( - Coord, - &mut f64, - )| { - *value_mut = f(world_coord); - }, - ) - }, - ); - } -} - -/// Why not just put this into Domain::par_modify_access? -/// Rust compiler can't figure out how to borrow aabb and buffer -/// at the same time in this way. -/// By putting their borrows into one function call first we work around it. -fn par_modify_access_impl<'a, const GRID_DIMENSION: usize>( - buffer: &'a mut [f64], - aabb: &'a AABB, - chunk_size: usize, -) -> impl ParallelIterator> + 'a { - buffer[0..aabb.buffer_size()] - .par_chunks_mut(chunk_size) - .enumerate() - .map(move |(i, buffer_chunk): (usize, &mut [f64])| { - let offset = i * chunk_size; - DomainChunk::new(offset, aabb, buffer_chunk) - }) -} - -pub struct DomainChunk<'a, const GRID_DIMENSION: usize> { - offset: usize, - aabb: &'a AABB, - buffer: &'a mut [f64], -} - -impl<'a, const GRID_DIMENSION: usize> DomainChunk<'a, GRID_DIMENSION> { - pub fn new( - offset: usize, - aabb: &'a AABB, - buffer: &'a mut [f64], - ) -> Self { - DomainChunk { - offset, - aabb, - buffer, - } - } - - pub fn coord_iter_mut( - &mut self, - ) -> impl Iterator, &mut f64)> { - self.buffer - .iter_mut() - .enumerate() - .map(|(i, v): (usize, &mut f64)| { - let linear_index = self.offset + i; - let coord = self.aabb.linear_to_coord(linear_index); - (coord, v) - }) - } -} diff --git a/src/domain/gather_args.rs b/src/domain/gather_args.rs index 4bdeb9b..dbfed8f 100644 --- a/src/domain/gather_args.rs +++ b/src/domain/gather_args.rs @@ -7,10 +7,11 @@ pub fn gather_args< Operation, const GRID_DIMENSION: usize, const NEIGHBORHOOD_SIZE: usize, + DomainType: DomainView, >( stencil: &StencilF64, bc: &BC, - input: &Domain, + input: &DomainType, world_coord: &Coord, ) -> [f64; NEIGHBORHOOD_SIZE] where @@ -23,12 +24,6 @@ where result[i] = bc .check(&n_world_coord) .unwrap_or_else(|| input.view(&n_world_coord)); - /* - println!( - "wc: {:?}, ni: {:?}, nwc: {:?}, r: {}\n", - world_coord, n_i, n_world_coord, result[i] - ); - */ } result } @@ -42,13 +37,11 @@ mod unit_tests { #[test] fn gather_args_test_const() { let bound = AABB::new(matrix![0, 9; 0, 9]); - let n_r = bound.buffer_size(); - let mut buffer = fftw::array::AlignedVec::new(n_r); - for i in 0..n_r { - let coord = bound.linear_to_coord(i); - buffer.as_slice_mut()[i] = (coord[0] + 3 * coord[1]) as f64; - } - let domain = Domain::new(bound, buffer.as_slice_mut()); + let mut domain = OwnedDomain::new(bound); + domain.par_set_values( + |coord: Coord<2>| (coord[0] + 3 * coord[1]) as f64, + 1, + ); let bc = ConstantCheck::new(-4.0, bound); let stencil = Stencil::new( [[0, -1], [0, 1], [1, 0], [-1, 0], [0, 0]], @@ -82,7 +75,8 @@ mod unit_tests { buffer.as_slice_mut()[i] ); } - let domain = Domain::new(bound, buffer.as_slice_mut()); + let mut domain = OwnedDomain::new(bound); + domain.par_set_values(|coord| (coord[0] + 3 * coord[1]) as f64, 1); let bc = PeriodicCheck::new(&domain); let stencil = Stencil::new( [[0, -1], [0, 1], [1, 0], [-1, 0], [0, 0]], diff --git a/src/domain/mod.rs b/src/domain/mod.rs index 415ebaf..0970440 100644 --- a/src/domain/mod.rs +++ b/src/domain/mod.rs @@ -6,9 +6,9 @@ //! and translate from world coordinates into view coordinates. mod bc; -mod domain_view; mod gather_args; +mod view; pub use bc::*; -pub use domain_view::*; pub use gather_args::*; +pub use view::*; diff --git a/src/domain/view/chunk.rs b/src/domain/view/chunk.rs new file mode 100644 index 0000000..38659f4 --- /dev/null +++ b/src/domain/view/chunk.rs @@ -0,0 +1,34 @@ +use crate::util::*; + +pub struct DomainChunk<'a, const GRID_DIMENSION: usize> { + offset: usize, + aabb: &'a AABB, + buffer: &'a mut [f64], +} + +impl<'a, const GRID_DIMENSION: usize> DomainChunk<'a, GRID_DIMENSION> { + pub fn new( + offset: usize, + aabb: &'a AABB, + buffer: &'a mut [f64], + ) -> Self { + DomainChunk { + offset, + aabb, + buffer, + } + } + + pub fn coord_iter_mut( + &mut self, + ) -> impl Iterator, &mut f64)> { + self.buffer + .iter_mut() + .enumerate() + .map(|(i, v): (usize, &mut f64)| { + let linear_index = self.offset + i; + let coord = self.aabb.linear_to_coord(linear_index); + (coord, v) + }) + } +} diff --git a/src/domain/view/mod.rs b/src/domain/view/mod.rs new file mode 100644 index 0000000..dff08fe --- /dev/null +++ b/src/domain/view/mod.rs @@ -0,0 +1,69 @@ +mod chunk; +mod owned; +mod slice; + +pub use chunk::*; +pub use owned::*; +pub use slice::*; + +use crate::util::*; +use rayon::prelude::*; + +pub trait DomainView { + fn aabb(&self) -> &AABB; + + fn set_aabb(&mut self, aabb: AABB); + + fn buffer(&self) -> &[f64]; + + fn buffer_mut(&mut self) -> (&AABB, &mut [f64]); + + fn view(&self, world_coord: &Coord) -> f64; + + fn par_modify_access<'a>( + &'a mut self, + chunk_size: usize, + ) -> impl ParallelIterator> { + let (aabb, buffer) = self.buffer_mut(); + par_modify_access_impl(buffer, aabb, chunk_size) + } + + fn par_set_values< + F: FnOnce(Coord) -> f64 + Send + Sync + Copy, + >( + &mut self, + f: F, + chunk_size: usize, + ) { + self.par_modify_access(chunk_size).for_each( + |mut d: DomainChunk<'_, GRID_DIMENSION>| { + d.coord_iter_mut().for_each( + |(world_coord, value_mut): ( + Coord, + &mut f64, + )| { + *value_mut = f(world_coord); + }, + ) + }, + ); + } +} + +/// Why not just put this into Domain::par_modify_access? +/// Rust compiler can't figure out how to borrow aabb and buffer +/// at the same time in this way. +/// By putting their borrows into one function call first we work around it. +fn par_modify_access_impl<'a, const GRID_DIMENSION: usize>( + buffer: &'a mut [f64], + aabb: &'a AABB, + chunk_size: usize, +) -> impl ParallelIterator> + 'a { + buffer[0..aabb.buffer_size()] + .par_chunks_mut(chunk_size) + .enumerate() + .map(move |(i, buffer_chunk): (usize, &mut [f64])| { + let offset = i * chunk_size; + DomainChunk::new(offset, aabb, buffer_chunk) + }) +} diff --git a/src/domain/view/owned.rs b/src/domain/view/owned.rs new file mode 100644 index 0000000..54136db --- /dev/null +++ b/src/domain/view/owned.rs @@ -0,0 +1,43 @@ +use super::*; +use crate::util::*; +use fftw::array::*; + +pub struct OwnedDomain { + aabb: AABB, + buffer: AlignedVec, +} + +impl OwnedDomain { + pub fn new(aabb: AABB) -> Self { + let buffer = AlignedVec::new(aabb.buffer_size()); + OwnedDomain { aabb, buffer } + } +} + +impl DomainView + for OwnedDomain +{ + fn aabb(&self) -> &AABB { + &self.aabb + } + + fn set_aabb(&mut self, aabb: AABB) { + debug_assert!(aabb.buffer_size() <= self.buffer.len()); + // TODO: should we re-slice here? + self.aabb = aabb; + } + + fn buffer(&self) -> &[f64] { + &self.buffer + } + + fn buffer_mut(&mut self) -> (&AABB, &mut [f64]) { + (&self.aabb, &mut self.buffer) + } + + fn view(&self, world_coord: &Coord) -> f64 { + debug_assert!(self.aabb.contains(world_coord)); + let index = self.aabb.coord_to_linear(world_coord); + self.buffer[index] + } +} diff --git a/src/domain/view/slice.rs b/src/domain/view/slice.rs new file mode 100644 index 0000000..a171c28 --- /dev/null +++ b/src/domain/view/slice.rs @@ -0,0 +1,42 @@ +use super::*; +use crate::util::*; + +pub struct SliceDomain<'a, const GRID_DIMENSION: usize> { + aabb: AABB, + buffer: &'a mut [f64], +} + +impl<'a, const GRID_DIMENSION: usize> SliceDomain<'a, GRID_DIMENSION> { + pub fn new(aabb: AABB, buffer: &'a mut [f64]) -> Self { + debug_assert_eq!(buffer.len(), aabb.buffer_size()); + SliceDomain { aabb, buffer } + } +} + +impl<'a, const GRID_DIMENSION: usize> DomainView + for SliceDomain<'a, GRID_DIMENSION> +{ + fn aabb(&self) -> &AABB { + &self.aabb + } + + fn set_aabb(&mut self, aabb: AABB) { + debug_assert!(aabb.buffer_size() <= self.buffer.len()); + // TODO: should we re-slice here? + self.aabb = aabb; + } + + fn buffer(&self) -> &[f64] { + self.buffer + } + + fn buffer_mut(&mut self) -> (&AABB, &mut [f64]) { + (&self.aabb, self.buffer) + } + + fn view(&self, world_coord: &Coord) -> f64 { + debug_assert!(self.aabb.contains(world_coord)); + let index = self.aabb.coord_to_linear(world_coord); + self.buffer[index] + } +} diff --git a/src/image.rs b/src/image.rs index 62235a5..f72b13b 100644 --- a/src/image.rs +++ b/src/image.rs @@ -1,4 +1,4 @@ -use crate::domain::Domain; +use crate::domain::*; use crate::util::*; pub struct Image1D { @@ -30,7 +30,10 @@ impl Image1D { } } -pub fn image2d>(domain: &Domain<2>, s: &F) { +pub fn image2d, DomainType: DomainView<2>>( + domain: &DomainType, + s: &F, +) { let aabb = domain.aabb(); let exclusive_bounds = aabb.exclusive_bounds(); let gradient = colorous::TURBO; diff --git a/src/par_stencil.rs b/src/par_stencil.rs index 4d1388a..ddd660f 100644 --- a/src/par_stencil.rs +++ b/src/par_stencil.rs @@ -8,11 +8,12 @@ pub fn apply< Operation, const GRID_DIMENSION: usize, const NEIGHBORHOOD_SIZE: usize, + DomainType: DomainView + Sync, >( bc: &BC, stencil: &StencilF64, - input: &Domain, - output: &mut Domain, + input: &DomainType, + output: &mut DomainType, chunk_size: usize, ) where Operation: StencilOperation, @@ -45,31 +46,27 @@ mod unit_test { fn par_stencil_test_1d_simple() { let stencil = Stencil::new([[0]], |args: &[f64; 1]| args[0]); let bound = AABB::new(matrix![0, 99]); - let n_r = bound.buffer_size(); + let chunk_size = 1; { - let mut input_buffer = vec![1.0; n_r]; - let input_domain = Domain::new(bound, &mut input_buffer); - - let mut output_buffer = vec![2.0; n_r]; - let mut output_domain = Domain::new(bound, &mut output_buffer); + let mut input_domain = OwnedDomain::new(bound); + let mut output_domain = OwnedDomain::new(bound); + input_domain.par_set_values(|_| 1.0, chunk_size); let bc = PeriodicCheck::new(&input_domain); - apply(&bc, &stencil, &input_domain, &mut output_domain, 1); - for x in &output_buffer { + apply(&bc, &stencil, &input_domain, &mut output_domain, chunk_size); + for x in output_domain.buffer() { assert_approx_eq!(f64, *x, 1.0); } } { - let mut input_buffer = vec![2.0; n_r]; - let input_domain = Domain::new(bound, &mut input_buffer); - - let mut output_buffer = vec![1.0; n_r]; - let mut output_domain = Domain::new(bound, &mut output_buffer); + let mut input_domain = OwnedDomain::new(bound); + let mut output_domain = OwnedDomain::new(bound); + input_domain.par_set_values(|_| 2.0, chunk_size); let bc = PeriodicCheck::new(&input_domain); apply(&bc, &stencil, &input_domain, &mut output_domain, 1); - for x in &output_buffer { + for x in output_domain.buffer() { assert_approx_eq!(f64, *x, 2.0); } } @@ -88,6 +85,7 @@ mod unit_test { #[test] fn par_stencil_trapezoid_test_1d_simple() { + let chunk_size = 2; let stencil = Stencil::new([[-1], [0], [1]], |args| { let mut r = 0.0; for a in args { @@ -99,11 +97,9 @@ mod unit_test { let input_bound = AABB::new(matrix![0, 10]); let output_bound = AABB::new(matrix![1, 9]); - let mut input_buffer = vec![1.0; input_bound.buffer_size()]; - let mut output_buffer = vec![0.0; output_bound.buffer_size()]; - - let input_domain = Domain::new(input_bound, &mut input_buffer); - let mut output_domain = Domain::new(output_bound, &mut output_buffer); + let mut input_domain = OwnedDomain::new(input_bound); + let mut output_domain = OwnedDomain::new(output_bound); + input_domain.par_set_values(|_| 1.0, chunk_size); let bc = ErrorCheck { bounds: input_bound, @@ -111,8 +107,8 @@ mod unit_test { apply(&bc, &stencil, &input_domain, &mut output_domain, 2); - for i in output_buffer { - assert_approx_eq!(f64, i, 1.0); + for i in output_domain.buffer() { + assert_approx_eq!(f64, *i, 1.0); } } } diff --git a/src/solver/direct.rs b/src/solver/direct.rs index 96364d1..0c14b4e 100644 --- a/src/solver/direct.rs +++ b/src/solver/direct.rs @@ -3,16 +3,16 @@ use crate::par_stencil; use crate::stencil::*; pub fn box_apply< - 'a, BC, Operation, const GRID_DIMENSION: usize, const NEIGHBORHOOD_SIZE: usize, + DomainType: DomainView + Sync, >( bc: &BC, stencil: &StencilF64, - input: &mut Domain<'a, GRID_DIMENSION>, - output: &mut Domain<'a, GRID_DIMENSION>, + input: &mut DomainType, + output: &mut DomainType, steps: usize, chunk_size: usize, ) where @@ -51,12 +51,11 @@ mod unit_tests { { let chunk_size = 3; assert_approx_eq!(f64, stencil.apply(&[1.0; NEIGHBORHOOD_SIZE]), 1.0); - let n_r = bound.buffer_size(); - let mut input_buffer = vec![1.0; n_r]; - let mut output_buffer = vec![2.0; n_r]; - let mut input_domain = Domain::new(*bound, &mut input_buffer); - let mut output_domain = Domain::new(*bound, &mut output_buffer); + let mut input_domain = OwnedDomain::new(*bound); + let mut output_domain = OwnedDomain::new(*bound); + + input_domain.par_set_values(|_| 1.0, chunk_size); box_apply( bc_lookup, @@ -67,7 +66,7 @@ mod unit_tests { chunk_size, ); - for x in &output_buffer[0..n_r] { + for x in output_domain.buffer() { assert_approx_eq!(f64, *x, 1.0); } } @@ -155,21 +154,21 @@ mod unit_tests { #[test] fn shifter() { + let chunk_size = 1; let stencil = Stencil::new([[-1]], |args: &[f64; 1]| args[0]); let max_size = AABB::new(matrix![0, 9]); let mut input_buffer = AlignedVec::new(10); for i in 0..10 { input_buffer[i] = i as f64; } - let mut output_buffer = AlignedVec::new(10); - let mut input_domain = - Domain::new(max_size, input_buffer.as_slice_mut()); - let mut output_domain = - Domain::new(max_size, output_buffer.as_slice_mut()); + let mut input_domain = OwnedDomain::new(max_size); + let mut output_domain = OwnedDomain::new(max_size); + + input_domain + .par_set_values(|coord: Coord<1>| coord[0] as f64, chunk_size); let bc_lookup = ConstantCheck::new(-1.0, max_size); - let chunk_size = 1; let steps = 3; box_apply( @@ -181,10 +180,10 @@ mod unit_tests { chunk_size, ); for i in 0..3 { - assert_approx_eq!(f64, output_buffer[i], -1.0); + assert_approx_eq!(f64, output_domain.buffer()[i], -1.0); } for i in 3..10 { - assert_approx_eq!(f64, output_buffer[i], (i - 3) as f64); + assert_approx_eq!(f64, output_domain.buffer()[i], (i - 3) as f64); } } } diff --git a/src/solver/periodic_direct.rs b/src/solver/periodic_direct.rs index 3404895..5d5b5b4 100644 --- a/src/solver/periodic_direct.rs +++ b/src/solver/periodic_direct.rs @@ -3,14 +3,14 @@ use crate::par_stencil; use crate::stencil::*; pub fn direct_periodic_apply< - 'a, Operation, const GRID_DIMENSION: usize, const NEIGHBORHOOD_SIZE: usize, + DomainType: DomainView + Sync, >( stencil: &StencilF64, - input: &mut Domain<'a, GRID_DIMENSION>, - output: &mut Domain<'a, GRID_DIMENSION>, + input: &mut DomainType, + output: &mut DomainType, steps: usize, chunk_size: usize, ) where @@ -24,7 +24,7 @@ pub fn direct_periodic_apply< } std::mem::swap(input, output); } - //println!("final t"); + let bc = PeriodicCheck::new(input); par_stencil::apply(&bc, stencil, input, output, chunk_size); } @@ -50,12 +50,10 @@ mod unit_tests { { let chunk_size = 3; assert_approx_eq!(f64, stencil.apply(&[1.0; NEIGHBORHOOD_SIZE]), 1.0); - let n_r = bound.buffer_size(); - let mut input_buffer = vec![1.0; n_r]; - let mut output_buffer = vec![2.0; n_r]; - let mut input_domain = Domain::new(*bound, &mut input_buffer); - let mut output_domain = Domain::new(*bound, &mut output_buffer); + let mut input_domain = OwnedDomain::new(*bound); + let mut output_domain = OwnedDomain::new(*bound); + input_domain.par_set_values(|_| 1.0, chunk_size); direct_periodic_apply( stencil, &mut input_domain, @@ -64,7 +62,7 @@ mod unit_tests { chunk_size, ); - for x in &output_buffer[0..n_r] { + for x in output_domain.buffer() { assert_approx_eq!(f64, *x, 1.0); } } @@ -140,18 +138,18 @@ mod unit_tests { #[test] fn shifter() { + let chunk_size = 1; let stencil = Stencil::new([[-1]], |args: &[f64; 1]| args[0]); let bound = AABB::new(matrix![0, 9]); let mut input_buffer = AlignedVec::new(10); for i in 0..10 { input_buffer[i] = i as f64; } - let mut input_domain = Domain::new(bound, input_buffer.as_slice_mut()); + let mut input_domain = OwnedDomain::new(bound); + let mut output_domain = OwnedDomain::new(bound); + input_domain + .par_set_values(|coord: Coord<1>| coord[0] as f64, chunk_size); - let mut output_buffer = AlignedVec::new(10); - let mut output_domain = - Domain::new(bound, output_buffer.as_slice_mut()); - let chunk_size = 1; let n = 1; direct_periodic_apply( &stencil, @@ -161,7 +159,11 @@ mod unit_tests { chunk_size, ); for i in 0..10 { - assert_approx_eq!(f64, output_buffer[(i + n) % 10], i as f64); + assert_approx_eq!( + f64, + output_domain.buffer()[(i + n) % 10], + i as f64 + ); } } } diff --git a/src/solver/periodic_plan.rs b/src/solver/periodic_plan.rs index e8a7c82..6ad07e6 100644 --- a/src/solver/periodic_plan.rs +++ b/src/solver/periodic_plan.rs @@ -1,4 +1,4 @@ -use crate::domain::Domain; +use crate::domain::*; use crate::par_slice; use crate::solver::fft_plan::*; use crate::stencil::*; @@ -74,10 +74,10 @@ where } } - pub fn apply( + pub fn apply>( &mut self, - input: &mut Domain, - output: &mut Domain, + input: &mut DomainType, + output: &mut DomainType, steps: usize, chunk_size: usize, ) { @@ -99,7 +99,7 @@ where let n_c = input.aabb().complex_buffer_size(); fft_plan .forward_plan - .r2c(input.buffer_mut(), &mut self.complex_buffer[0..n_c]) + .r2c(input.buffer_mut().1, &mut self.complex_buffer[0..n_c]) .unwrap(); par_slice::multiply_by( &mut self.complex_buffer[0..n_c], @@ -108,9 +108,9 @@ where ); fft_plan .backward_plan - .c2r(&mut self.complex_buffer[0..n_c], output.buffer_mut()) + .c2r(&mut self.complex_buffer[0..n_c], output.buffer_mut().1) .unwrap(); - par_slice::div(output.buffer_mut(), n_r as f64, chunk_size); + par_slice::div(output.buffer_mut().1, n_r as f64, chunk_size); } fn new_convolution( @@ -188,20 +188,13 @@ mod unit_tests { ) where Operation: StencilOperation, { - let chunk_size = 1; + let chunk_size = 3; assert_approx_eq!(f64, stencil.apply(&[1.0; NEIGHBORHOOD_SIZE]), 1.0); - let rbs = bound.buffer_size(); - let mut input_x = fftw::array::AlignedVec::new(rbs); - for x in input_x.as_slice_mut() { - *x = 1.0f64; - } - let input_copy = input_x.clone(); - let mut input_domain = Domain::new(bound, input_x.as_slice_mut()); + let mut input_domain = OwnedDomain::new(bound); + let mut output_domain = OwnedDomain::new(bound); - let mut result_buffer = fftw::array::AlignedVec::new(rbs); - let mut output_domain = - Domain::new(bound, result_buffer.as_slice_mut()); + input_domain.par_set_values(|_| 1.0, chunk_size); plan_library.apply( &mut input_domain, @@ -209,10 +202,9 @@ mod unit_tests { steps, chunk_size, ); - for x in &result_buffer[0..rbs] { + for x in output_domain.buffer() { assert_approx_eq!(f64, *x, 1.0); } - assert_eq!(input_x.as_slice(), input_copy.as_slice()); } #[test] diff --git a/src/solver/trapezoid.rs b/src/solver/trapezoid.rs index f4b8a5e..0da3603 100644 --- a/src/solver/trapezoid.rs +++ b/src/solver/trapezoid.rs @@ -18,16 +18,16 @@ pub fn trapezoid_input_region( } pub fn trapezoid_apply< - 'a, BC, Operation, const GRID_DIMENSION: usize, const NEIGHBORHOOD_SIZE: usize, + DomainType: DomainView + Sync, >( bc: &BC, stencil: &StencilF64, - input: &mut Domain<'a, GRID_DIMENSION>, - output: &mut Domain<'a, GRID_DIMENSION>, + input: &mut DomainType, + output: &mut DomainType, sloped_sides: &Bounds, stencil_slopes: &Bounds, steps: usize, @@ -168,11 +168,11 @@ mod unit_tests { { let sloped_sides = matrix![1, 1]; let input_bound = AABB::new(matrix![10, 40]); - let mut input_buffer = vec![1.0; input_bound.buffer_size()]; - let mut output_buffer = vec![1.0; input_bound.buffer_size()]; - let mut input_domain = Domain::new(input_bound, &mut input_buffer); - let mut output_domain = - Domain::new(input_bound, &mut output_buffer); + let mut input_domain = OwnedDomain::new(input_bound); + let mut output_domain = OwnedDomain::new(input_bound); + + input_domain.par_set_values(|_| 1.0, chunk_size); + let bc = ConstantCheck::new(1.0, input_bound); trapezoid_apply( &bc, diff --git a/tests/base_solver_compare.rs b/tests/base_solver_compare.rs index 9caff5c..01c72d2 100644 --- a/tests/base_solver_compare.rs +++ b/tests/base_solver_compare.rs @@ -8,42 +8,23 @@ use float_cmp::assert_approx_eq; use nalgebra::matrix; #[test] -fn thermal_1d_compare() { +fn heat_1d_compare() { // Grid size let grid_bound = AABB::new(matrix![0, 999]); let n_steps = 16; - // Step size t - let dt: f64 = 1.0; - - // Step size x - let dx: f64 = 1.0; - - // Heat transfer coefficient - let k: f64 = 0.5; - let chunk_size = 100; - let stencil = Stencil::new([[-1], [0], [1]], |args: &[f64; 3]| { - let left = args[0]; - let middle = args[1]; - let right = args[2]; - middle + (k * dt / (dx * dx)) * (left - 2.0 * middle + right) - }); + let stencil = nhls::standard_stencils::heat_1d(1.0, 1.0, 0.5); // Create domains let buffer_size = grid_bound.buffer_size(); - let mut grid_input = vec![0.0; buffer_size]; - let mut direct_input_domain = Domain::new(grid_bound, &mut grid_input); + let mut direct_input_domain = OwnedDomain::new(grid_bound); + let mut direct_output_domain = OwnedDomain::new(grid_bound); - let mut grid_output = vec![0.0; buffer_size]; - let mut direct_output_domain = Domain::new(grid_bound, &mut grid_output); - - let mut fft_input = AlignedVec::new(buffer_size); - let mut fft_output = AlignedVec::new(buffer_size); - let mut fft_input_domain = Domain::new(grid_bound, &mut fft_input); - let mut fft_output_domain = Domain::new(grid_bound, &mut fft_output); + let mut fft_input_domain = OwnedDomain::new(grid_bound); + let mut fft_output_domain = OwnedDomain::new(grid_bound); // Fill in with IC values (use normal dist for spike in the middle) let n_f = buffer_size as f64; @@ -53,9 +34,7 @@ fn thermal_1d_compare() { let exp = -x * x / (2.0 * sigma_sq); exp.exp() }; - direct_input_domain.par_set_values(ic_gen, chunk_size); - fft_input_domain.par_set_values(ic_gen, chunk_size); let mut periodic_library = @@ -80,8 +59,12 @@ fn thermal_1d_compare() { ); for i in 0..buffer_size { - // TODO THIS IS PRETTY BAD - assert_approx_eq!(f64, fft_output[i], grid_output[i], epsilon = 0.001); + assert_approx_eq!( + f64, + fft_output_domain.buffer()[i], + direct_output_domain.buffer()[i], + epsilon = 0.0000000000001 + ); } } @@ -104,15 +87,14 @@ fn periodic_compare() { for i in 0..n_r { input_a[i] = i as f64; } - let mut input_b = input_a.clone(); - let mut domain_a_input = Domain::new(bound, input_a.as_slice_mut()); - let mut domain_b_input = Domain::new(bound, input_b.as_slice_mut()); + let mut domain_a_input = OwnedDomain::new(bound); + let mut domain_b_input = OwnedDomain::new(bound); + let mut domain_a_output = OwnedDomain::new(bound); + let mut domain_b_output = OwnedDomain::new(bound); - let mut output_a = AlignedVec::new(n_r); - let mut output_b = AlignedVec::new(n_r); - let mut domain_a_output = Domain::new(bound, output_a.as_slice_mut()); - let mut domain_b_output = Domain::new(bound, output_b.as_slice_mut()); + domain_a_input.par_set_values(|coord| coord[0] as f64, chunk_size); + domain_b_input.par_set_values(|coord| coord[0] as f64, chunk_size); direct_periodic_apply( &stencil, @@ -131,19 +113,10 @@ fn periodic_compare() { ); for i in 0..n_r { - println!( - "n: {}, p: {}, d: {}", - output_a[i], - output_b[i], - (output_a[i] - output_b[i]).abs() - ); - } - for i in 0..n_r { - // TODO THIS IS BROKE assert_approx_eq!( f64, - output_a[i], - output_b[i], + domain_a_output.buffer()[i], + domain_b_output.buffer()[i], epsilon = 0.0000000000001 ); }