From 821cad8dff874c2ab92f91af972598a59f2c7b9a Mon Sep 17 00:00:00 2001 From: Timo Betcke Date: Tue, 5 Nov 2024 11:09:57 +0000 Subject: [PATCH] WiP: Clean up assembler --- benches/assembly_benchmark.rs | 2 +- src/assembly/boundary/assemblers.rs | 259 ++++++++++------------------ src/assembly/boundary/integrands.rs | 2 +- src/function.rs | 2 +- src/helmholtz.rs | 8 +- src/laplace.rs | 8 +- 6 files changed, 106 insertions(+), 175 deletions(-) diff --git a/benches/assembly_benchmark.rs b/benches/assembly_benchmark.rs index f3b06a24..0d60a5a7 100644 --- a/benches/assembly_benchmark.rs +++ b/benches/assembly_benchmark.rs @@ -41,7 +41,7 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) { space.global_size(), space.global_size() ), - |b| b.iter(|| black_box(assembler.assemble_singular_into_csr(&space, &space))), + |b| b.iter(|| black_box(assembler.assemble_singular(&space, &space))), ); } group.finish(); diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index 4691d453..274661e0 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -32,35 +32,6 @@ use rlst::{ }; use std::collections::HashMap; -fn neighbours( - test_grid: &TestGrid, - trial_grid: &TrialGrid, - test_cell: usize, - trial_cell: usize, -) -> bool { - if !equal_grids(test_grid, trial_grid) { - false - } else { - let test_vertices = trial_grid - .entity(2, test_cell) - .unwrap() - .topology() - .sub_entity_iter(0) - .collect::>(); - for v in trial_grid - .entity(2, trial_cell) - .unwrap() - .topology() - .sub_entity_iter(0) - { - if test_vertices.contains(&v) { - return true; - } - } - false - } -} - /// Options for a boundary assembler #[derive(Clone)] pub struct BoundaryAssemblerOptions { @@ -129,15 +100,73 @@ pub struct BoundaryAssembler< pub(crate) table_derivs: usize, } -unsafe impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: Kernel> - Sync for BoundaryAssembler<'o, T, Integrand, K> -{ -} - impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: Kernel> BoundaryAssembler<'o, T, Integrand, K> { - /// Create new + /// Assemble the singular part into a CSR matrix. + pub fn assemble_singular>( + &self, + trial_space: &Space, + test_space: &Space, + ) -> CsrMatrix { + let shape = [test_space.global_size(), trial_space.global_size()]; + let sparse_matrix = self.assemble_singular_part(shape, trial_space, test_space); + + CsrMatrix::::from_aij( + sparse_matrix.shape, + &sparse_matrix.rows, + &sparse_matrix.cols, + &sparse_matrix.data, + ) + .unwrap() + } + + /// Assemble into a dense matrix. + pub fn assemble< + Space: FunctionSpace, + Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut, + >( + &self, + output: &mut Array2, + trial_space: &Space, + test_space: &Space, + ) { + let test_colouring = test_space.cell_colouring(); + let trial_colouring = trial_space.cell_colouring(); + + if !trial_space.is_serial() || !test_space.is_serial() { + panic!("Dense assembly can only be used for function spaces stored in serial"); + } + if output.shape()[0] != test_space.global_size() + || output.shape()[1] != trial_space.global_size() + { + panic!("Matrix has wrong shape"); + } + + let output_raw = RawData2D { + data: output.data_mut().as_mut_ptr(), + shape: output.shape(), + }; + + self.assemble_nonsingular_part( + &output_raw, + trial_space, + test_space, + &trial_colouring, + &test_colouring, + ); + + let sparse_matrix = self.assemble_singular_part(output.shape(), trial_space, test_space); + + let data = sparse_matrix.data; + let rows = sparse_matrix.rows; + let cols = sparse_matrix.cols; + for ((i, j), value) in rows.iter().zip(cols.iter()).zip(data.iter()) { + *output.get_mut([*i, *j]).unwrap() += *value; + } + } + + /// Create new Boundary assembler pub fn new( integrand: Integrand, kernel: KernelEvaluator, @@ -154,52 +183,8 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: } } - // /// Set (non-singular) quadrature degree for a cell type - // pub fn set_quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) { - // *self.options.quadrature_degrees.get_mut(&cell).unwrap() = degree; - // } - - // /// Get (non-singular) quadrature degree for a cell type - // pub fn quadrature_degree(&self, cell: ReferenceCellType) -> Option { - // self.options.quadrature_degrees.get(&cell).copied() - // } - - // /// Set singular quadrature degree for a pair of cell types - // pub fn set_singular_quadrature_degree( - // &mut self, - // cells: (ReferenceCellType, ReferenceCellType), - // degree: usize, - // ) { - // *self - // .options - // .singular_quadrature_degrees - // .get_mut(&cells) - // .unwrap() = degree; - // } - - // /// Get singular quadrature degree for a pair of cell types - // pub fn singular_quadrature_degree( - // &self, - // cells: (ReferenceCellType, ReferenceCellType), - // ) -> Option { - // self.options - // .singular_quadrature_degrees - // .get(&cells) - // .copied() - // } - - // /// Set the maximum size of a batch of cells to send to an assembly function - // pub fn set_batch_size(&mut self, size: usize) { - // self.options.batch_size = size; - // } - - // /// Get the maximum size of a batch of cells to send to an assembly function - // pub fn batch_size(&self) -> usize { - // self.options.batch_size - // } - /// Assemble the singular contributions - pub(crate) fn assemble_singular_part + Sync>( + fn assemble_singular_part>( &self, shape: [usize; 2], trial_space: &Space, @@ -361,7 +346,7 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: } /// Assemble the non-singular contributions into a dense matrix - pub(crate) fn assemble_nonsingular + Sync>( + fn assemble_nonsingular_part>( &self, output: &RawData2D, trial_space: &Space, @@ -487,89 +472,6 @@ impl<'o, T: RlstScalar + MatrixInverse, Integrand: BoundaryIntegrand, K: } } } - - /// Assemble the singular part into a CSR matrix. - pub fn assemble_singular_into_csr + Sync>( - &self, - trial_space: &Space, - test_space: &Space, - ) -> CsrMatrix { - let shape = [test_space.global_size(), trial_space.global_size()]; - let sparse_matrix = self.assemble_singular_part(shape, trial_space, test_space); - - CsrMatrix::::from_aij( - sparse_matrix.shape, - &sparse_matrix.rows, - &sparse_matrix.cols, - &sparse_matrix.data, - ) - .unwrap() - } - - /// Assemble into a dense matrix. - pub fn assemble_into_dense< - Space: FunctionSpace + Sync, - Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut, - >( - &self, - output: &mut Array2, - trial_space: &Space, - test_space: &Space, - ) { - let test_colouring = test_space.cell_colouring(); - let trial_colouring = trial_space.cell_colouring(); - - self.assemble_nonsingular_into_dense( - output, - trial_space, - test_space, - &trial_colouring, - &test_colouring, - ); - - let sparse_matrix = self.assemble_singular_part(output.shape(), trial_space, test_space); - let data = sparse_matrix.data; - let rows = sparse_matrix.rows; - let cols = sparse_matrix.cols; - for ((i, j), value) in rows.iter().zip(cols.iter()).zip(data.iter()) { - *output.get_mut([*i, *j]).unwrap() += *value; - } - } - - /// Assemble the nonsingular part into a dense matrix. - pub fn assemble_nonsingular_into_dense< - Space: FunctionSpace + Sync, - Array2: RandomAccessMut<2, Item = T> + Shape<2> + RawAccessMut, - >( - &self, - output: &mut Array2, - trial_space: &Space, - test_space: &Space, - trial_colouring: &HashMap>>, - test_colouring: &HashMap>>, - ) { - if !trial_space.is_serial() || !test_space.is_serial() { - panic!("Dense assembly can only be used for function spaces stored in serial"); - } - if output.shape()[0] != test_space.global_size() - || output.shape()[1] != trial_space.global_size() - { - panic!("Matrix has wrong shape"); - } - - let output_raw = RawData2D { - data: output.data_mut().as_mut_ptr(), - shape: output.shape(), - }; - - self.assemble_nonsingular( - &output_raw, - trial_space, - test_space, - trial_colouring, - test_colouring, - ) - } } fn get_singular_quadrature_rule( @@ -841,3 +743,32 @@ fn get_pairs_if_smallest( } Some(pairs) } + +fn neighbours( + test_grid: &TestGrid, + trial_grid: &TrialGrid, + test_cell: usize, + trial_cell: usize, +) -> bool { + if !equal_grids(test_grid, trial_grid) { + false + } else { + let test_vertices = trial_grid + .entity(2, test_cell) + .unwrap() + .topology() + .sub_entity_iter(0) + .collect::>(); + for v in trial_grid + .entity(2, trial_cell) + .unwrap() + .topology() + .sub_entity_iter(0) + { + if test_vertices.contains(&v) { + return true; + } + } + false + } +} diff --git a/src/assembly/boundary/integrands.rs b/src/assembly/boundary/integrands.rs index 06670ea3..762e2af3 100644 --- a/src/assembly/boundary/integrands.rs +++ b/src/assembly/boundary/integrands.rs @@ -175,7 +175,7 @@ impl<'a, T: RlstScalar, G: CellGeometry> GeometryAccess for Geometr } } -pub unsafe trait BoundaryIntegrand { +pub unsafe trait BoundaryIntegrand: Sync { //! Integrand //! //! # Safety diff --git a/src/function.rs b/src/function.rs index 61f3688e..2771dd5e 100644 --- a/src/function.rs +++ b/src/function.rs @@ -14,7 +14,7 @@ use rlst::RlstScalar; use std::collections::HashMap; /// A function space -pub trait FunctionSpace { +pub trait FunctionSpace: Sync { /// Scalar type type T: RlstScalar; /// The grid type diff --git a/src/helmholtz.rs b/src/helmholtz.rs index 08a1ec4b..10488d83 100644 --- a/src/helmholtz.rs +++ b/src/helmholtz.rs @@ -44,7 +44,7 @@ pub mod assembler { let assembler = BoundaryAssembler::new(SingleLayerBoundaryIntegrand::new(), kernel, options, 1, 0); - assembler.assemble_into_dense(&mut output, trial_space, test_space); + assembler.assemble(&mut output, trial_space, test_space); output } @@ -72,7 +72,7 @@ pub mod assembler { let assembler = BoundaryAssembler::new(DoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0); - assembler.assemble_into_dense(&mut output, trial_space, test_space); + assembler.assemble(&mut output, trial_space, test_space); output } @@ -105,7 +105,7 @@ pub mod assembler { 0, ); - assembler.assemble_into_dense(&mut output, trial_space, test_space); + assembler.assemble(&mut output, trial_space, test_space); output } @@ -140,7 +140,7 @@ pub mod assembler { let assembler = BoundaryAssembler::new(integrand, kernel, options, 4, 1); - assembler.assemble_into_dense(&mut output, trial_space, test_space); + assembler.assemble(&mut output, trial_space, test_space); output } diff --git a/src/laplace.rs b/src/laplace.rs index 6731dc4d..e97fea6c 100644 --- a/src/laplace.rs +++ b/src/laplace.rs @@ -38,7 +38,7 @@ pub mod assembler { let assembler = BoundaryAssembler::new(SingleLayerBoundaryIntegrand::new(), kernel, options, 1, 0); - assembler.assemble_into_dense(&mut output, trial_space, test_space); + assembler.assemble(&mut output, trial_space, test_space); output } @@ -62,7 +62,7 @@ pub mod assembler { let assembler = BoundaryAssembler::new(DoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0); - assembler.assemble_into_dense(&mut output, trial_space, test_space); + assembler.assemble(&mut output, trial_space, test_space); output } @@ -91,7 +91,7 @@ pub mod assembler { 0, ); - assembler.assemble_into_dense(&mut output, trial_space, test_space); + assembler.assemble(&mut output, trial_space, test_space); output } @@ -120,7 +120,7 @@ pub mod assembler { 1, ); - assembler.assemble_into_dense(&mut output, trial_space, test_space); + assembler.assemble(&mut output, trial_space, test_space); output }