From 8cc7a3de4f296b7661bd2d620c599bd1ad385bd8 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 28 Aug 2024 10:45:12 +0100 Subject: [PATCH] cell pair assemblers with caching --- src/assembly/boundary/assemblers.rs | 8 +- src/assembly/boundary/cell_pair_assemblers.rs | 508 +++++++++++++++++- 2 files changed, 513 insertions(+), 3 deletions(-) diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index 99ada30f..1b753b98 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -4,7 +4,10 @@ pub mod double_layer; pub mod hypersingular; pub mod single_layer; -use super::cell_pair_assemblers::{NonsingularCellPairAssembler, SingularCellPairAssembler}; +use super::cell_pair_assemblers::{ + NonsingularCellPairAssembler, NonsingularCellPairAssemblerWithTestCaching, + SingularCellPairAssembler, +}; use crate::assembly::common::{equal_grids, RawData2D, RlstArray, SparseMatrixData}; use crate::quadrature::duffy::{ quadrilateral_duffy, quadrilateral_triangle_duffy, triangle_duffy, triangle_quadrilateral_duffy, @@ -264,10 +267,11 @@ fn assemble_batch_nonadjacent< let test_evaluator = test_grid.geometry_map(test_cell_type, test_points.data()); let trial_evaluator = trial_grid.geometry_map(trial_cell_type, trial_points.data()); - let mut a = NonsingularCellPairAssembler::new( + let mut a = NonsingularCellPairAssemblerWithTestCaching::new( npts_test, npts_trial, deriv_size, + test_cells, &assembler.integrand, &assembler.kernel, test_evaluator, diff --git a/src/assembly/boundary/cell_pair_assemblers.rs b/src/assembly/boundary/cell_pair_assemblers.rs index 45849b63..d01eeced 100644 --- a/src/assembly/boundary/cell_pair_assemblers.rs +++ b/src/assembly/boundary/cell_pair_assemblers.rs @@ -234,7 +234,6 @@ impl< } } } -// TODO: make version of this with trial caching impl< 'a, T: RlstScalar, @@ -317,3 +316,510 @@ impl< } } } + +/// Assembler for the contributions from pairs of non-neighbouring cells with trial geometry caching +pub struct NonsingularCellPairAssemblerWithTrialCaching< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, +> { + integrand: &'a I, + kernel: &'a K, + test_evaluator: TestG, + trial_evaluator: TrialG, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + k: RlstArray, + test_mapped_pts: RlstArray, + trial_mapped_pts: Vec>, + test_normals: RlstArray, + trial_normals: Vec>, + test_jacobians: RlstArray, + trial_jacobians: Vec>, + test_jdet: Vec, + trial_jdet: Vec>, + test_weights: &'a [T::Real], + trial_weights: &'a [T::Real], + test_cell: usize, + trial_cell: usize, + trial_cells: &'a [usize], + npts_trial: usize, +} + +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, + > NonsingularCellPairAssemblerWithTrialCaching<'a, T, I, TestG, TrialG, K> +{ + // Allow dead code as this currently is only used in testing + #[allow(dead_code)] + #[allow(clippy::too_many_arguments)] + /// Create new + pub fn new( + npts_test: usize, + npts_trial: usize, + deriv_size: usize, + trial_cells: &'a [usize], + integrand: &'a I, + kernel: &'a K, + test_evaluator: TestG, + trial_evaluator: TrialG, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + test_weights: &'a [T::Real], + trial_weights: &'a [T::Real], + ) -> Self { + Self { + integrand, + kernel, + test_evaluator, + trial_evaluator, + test_table, + trial_table, + k: rlst_dynamic_array3!(T, [deriv_size, npts_test, npts_trial]), + test_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts_test]), + trial_mapped_pts: trial_cells + .iter() + .map(|_| rlst_dynamic_array2!(T::Real, [3, npts_trial])) + .collect::>(), + test_normals: rlst_dynamic_array2!(T::Real, [3, npts_test]), + trial_normals: trial_cells + .iter() + .map(|_| rlst_dynamic_array2!(T::Real, [3, npts_trial])) + .collect::>(), + test_jacobians: rlst_dynamic_array2!(T::Real, [6, npts_test]), + trial_jacobians: trial_cells + .iter() + .map(|_| rlst_dynamic_array2!(T::Real, [6, npts_trial])) + .collect::>(), + test_jdet: vec![T::Real::zero(); npts_test], + trial_jdet: vec![vec![]; trial_cells.len()], + test_weights, + trial_weights, + test_cell: 0, + trial_cell: 0, + trial_cells, + npts_trial, + } + } +} +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, + > CellPairAssembler + for NonsingularCellPairAssemblerWithTrialCaching<'a, T, I, TestG, TrialG, K> +{ + type T = T; + fn set_test_cell(&mut self, test_cell: usize) { + self.test_cell = test_cell; + self.test_evaluator + .points(test_cell, self.test_mapped_pts.data_mut()); + self.test_evaluator.jacobians_dets_normals( + test_cell, + self.test_jacobians.data_mut(), + &mut self.test_jdet, + self.test_normals.data_mut(), + ); + } + fn set_trial_cell(&mut self, trial_cell: usize) { + self.trial_cell = self + .trial_cells + .iter() + .position(|r| *r == trial_cell) + .unwrap(); + if self.trial_jdet[self.trial_cell].is_empty() { + self.trial_jdet[self.trial_cell] = vec![T::Real::zero(); self.npts_trial]; + self.trial_evaluator.points( + trial_cell, + self.trial_mapped_pts[self.trial_cell].data_mut(), + ); + self.trial_evaluator.jacobians_dets_normals( + trial_cell, + self.trial_jacobians[self.trial_cell].data_mut(), + &mut self.trial_jdet[self.trial_cell], + self.trial_normals[self.trial_cell].data_mut(), + ); + } + } + fn assemble(&mut self, local_mat: &mut RlstArray) { + self.kernel.assemble_st( + self.test_mapped_pts.data(), + self.trial_mapped_pts[self.trial_cell].data(), + self.k.data_mut(), + ); + + let test_geometry = AssemblerGeometry::new( + &self.test_mapped_pts, + &self.test_normals, + &self.test_jacobians, + &self.test_jdet, + ); + let trial_geometry = AssemblerGeometry::new( + &self.trial_mapped_pts[self.trial_cell], + &self.trial_normals[self.trial_cell], + &self.trial_jacobians[self.trial_cell], + &self.trial_jdet[self.trial_cell], + ); + + for (trial_i, mut col) in local_mat.col_iter_mut().enumerate() { + for (test_i, entry) in col.iter_mut().enumerate() { + *entry = T::zero(); + for (test_index, test_wt) in self.test_weights.iter().enumerate() { + for (trial_index, trial_wt) in self.trial_weights.iter().enumerate() { + *entry += unsafe { + self.integrand.evaluate_nonsingular( + self.test_table, + self.trial_table, + test_index, + trial_index, + test_i, + trial_i, + &self.k, + &test_geometry, + &trial_geometry, + ) * num::cast::( + *test_wt + * self.test_jdet[test_index] + * *trial_wt + * self.trial_jdet[self.trial_cell][trial_index], + ) + .unwrap() + }; + } + } + } + } + } +} + +/// Assembler for the contributions from pairs of non-neighbouring cells with trial geometry caching +pub struct NonsingularCellPairAssemblerWithTestCaching< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, +> { + integrand: &'a I, + kernel: &'a K, + test_evaluator: TestG, + trial_evaluator: TrialG, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + k: RlstArray, + test_mapped_pts: Vec>, + trial_mapped_pts: RlstArray, + test_normals: Vec>, + trial_normals: RlstArray, + test_jacobians: Vec>, + trial_jacobians: RlstArray, + test_jdet: Vec>, + trial_jdet: Vec, + test_weights: &'a [T::Real], + trial_weights: &'a [T::Real], + test_cell: usize, + trial_cell: usize, + test_cells: &'a [usize], + npts_test: usize, +} + +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, + > NonsingularCellPairAssemblerWithTestCaching<'a, T, I, TestG, TrialG, K> +{ + #[allow(clippy::too_many_arguments)] + /// Create new + pub fn new( + npts_test: usize, + npts_trial: usize, + deriv_size: usize, + test_cells: &'a [usize], + integrand: &'a I, + kernel: &'a K, + test_evaluator: TestG, + trial_evaluator: TrialG, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + test_weights: &'a [T::Real], + trial_weights: &'a [T::Real], + ) -> Self { + Self { + integrand, + kernel, + test_evaluator, + trial_evaluator, + test_table, + trial_table, + k: rlst_dynamic_array3!(T, [deriv_size, npts_test, npts_trial]), + test_mapped_pts: test_cells + .iter() + .map(|_| rlst_dynamic_array2!(T::Real, [3, npts_test])) + .collect::>(), + trial_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts_trial]), + test_normals: test_cells + .iter() + .map(|_| rlst_dynamic_array2!(T::Real, [3, npts_test])) + .collect::>(), + trial_normals: rlst_dynamic_array2!(T::Real, [3, npts_trial]), + test_jacobians: test_cells + .iter() + .map(|_| rlst_dynamic_array2!(T::Real, [6, npts_test])) + .collect::>(), + trial_jacobians: rlst_dynamic_array2!(T::Real, [6, npts_trial]), + test_jdet: vec![vec![]; test_cells.len()], + trial_jdet: vec![T::Real::zero(); npts_trial], + test_weights, + trial_weights, + test_cell: 0, + trial_cell: 0, + test_cells, + npts_test, + } + } +} +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, + > CellPairAssembler + for NonsingularCellPairAssemblerWithTestCaching<'a, T, I, TestG, TrialG, K> +{ + type T = T; + fn set_test_cell(&mut self, test_cell: usize) { + self.test_cell = self + .test_cells + .iter() + .position(|r| *r == test_cell) + .unwrap(); + if self.test_jdet[self.test_cell].is_empty() { + self.test_jdet[self.test_cell] = vec![T::Real::zero(); self.npts_test]; + self.test_evaluator + .points(test_cell, self.test_mapped_pts[self.test_cell].data_mut()); + self.test_evaluator.jacobians_dets_normals( + test_cell, + self.test_jacobians[self.test_cell].data_mut(), + &mut self.test_jdet[self.test_cell], + self.test_normals[self.test_cell].data_mut(), + ); + } + } + fn set_trial_cell(&mut self, trial_cell: usize) { + self.trial_cell = trial_cell; + self.trial_evaluator + .points(trial_cell, self.trial_mapped_pts.data_mut()); + self.trial_evaluator.jacobians_dets_normals( + trial_cell, + self.trial_jacobians.data_mut(), + &mut self.trial_jdet, + self.trial_normals.data_mut(), + ); + } + fn assemble(&mut self, local_mat: &mut RlstArray) { + self.kernel.assemble_st( + self.test_mapped_pts[self.test_cell].data(), + self.trial_mapped_pts.data(), + self.k.data_mut(), + ); + + let test_geometry = AssemblerGeometry::new( + &self.test_mapped_pts[self.test_cell], + &self.test_normals[self.test_cell], + &self.test_jacobians[self.test_cell], + &self.test_jdet[self.test_cell], + ); + let trial_geometry = AssemblerGeometry::new( + &self.trial_mapped_pts, + &self.trial_normals, + &self.trial_jacobians, + &self.trial_jdet, + ); + + for (trial_i, mut col) in local_mat.col_iter_mut().enumerate() { + for (test_i, entry) in col.iter_mut().enumerate() { + *entry = T::zero(); + for (test_index, test_wt) in self.test_weights.iter().enumerate() { + for (trial_index, trial_wt) in self.trial_weights.iter().enumerate() { + *entry += unsafe { + self.integrand.evaluate_nonsingular( + self.test_table, + self.trial_table, + test_index, + trial_index, + test_i, + trial_i, + &self.k, + &test_geometry, + &trial_geometry, + ) * num::cast::( + *test_wt + * self.test_jdet[self.test_cell][test_index] + * *trial_wt + * self.trial_jdet[trial_index], + ) + .unwrap() + }; + } + } + } + } + } +} + +#[cfg(test)] +mod test { + use super::*; + use crate::{ + assembly::{ + boundary::integrands::SingleLayerBoundaryIntegrand, common::GreenKernelEvalType, + kernels::KernelEvaluator, + }, + quadrature::simplex_rules::simplex_rule, + }; + use approx::*; + use itertools::izip; + use ndelement::{ + ciarlet::LagrangeElementFamily, + traits::{ElementFamily, FiniteElement}, + types::{Continuity, ReferenceCellType}, + }; + use ndgrid::{shapes::regular_sphere, traits::Grid}; + use rlst::{rlst_dynamic_array4, DefaultIterator, RandomAccessMut}; + + #[test] + fn test_non_singular() { + let grid = regular_sphere::(2); + let family = LagrangeElementFamily::::new(2, Continuity::Standard); + let element = family.element(ReferenceCellType::Triangle); + + let mut test_cells = vec![]; + let mut trial_cells = vec![]; + + for i in 0..grid.entity_count(ReferenceCellType::Triangle) { + if i % 2 == 0 { + test_cells.push(i); + } else { + trial_cells.push(i); + } + } + + let npts_test = 3; + let qrule_test = simplex_rule(ReferenceCellType::Triangle, npts_test).unwrap(); + let mut test_points = rlst_dynamic_array2!(f64, [2, npts_test]); + for i in 0..npts_test { + for j in 0..2 { + *test_points.get_mut([j, i]).unwrap() = qrule_test.points[2 * i + j]; + } + } + let test_weights = qrule_test.weights; + + let npts_trial = 6; + let qrule_trial = simplex_rule(ReferenceCellType::Triangle, npts_trial).unwrap(); + let mut trial_points = rlst_dynamic_array2!(f64, [2, npts_trial]); + for i in 0..npts_trial { + for j in 0..2 { + *trial_points.get_mut([j, i]).unwrap() = qrule_trial.points[2 * i + j]; + } + } + let trial_weights = qrule_trial.weights; + + let mut test_table = rlst_dynamic_array4!(f64, element.tabulate_array_shape(0, npts_test)); + element.tabulate(&test_points, 0, &mut test_table); + let mut trial_table = + rlst_dynamic_array4!(f64, element.tabulate_array_shape(0, npts_trial)); + element.tabulate(&trial_points, 0, &mut trial_table); + + let integrand = SingleLayerBoundaryIntegrand::new(); + let kernel = KernelEvaluator::new_laplace(GreenKernelEvalType::Value); + + let mut a0 = NonsingularCellPairAssembler::new( + npts_test, + npts_trial, + 1, + &integrand, + &kernel, + grid.geometry_map(ReferenceCellType::Triangle, test_points.data()), + grid.geometry_map(ReferenceCellType::Triangle, trial_points.data()), + &test_table, + &trial_table, + &test_weights, + &trial_weights, + ); + + let mut a1 = NonsingularCellPairAssemblerWithTestCaching::new( + npts_test, + npts_trial, + 1, + &test_cells, + &integrand, + &kernel, + grid.geometry_map(ReferenceCellType::Triangle, test_points.data()), + grid.geometry_map(ReferenceCellType::Triangle, trial_points.data()), + &test_table, + &trial_table, + &test_weights, + &trial_weights, + ); + + let mut a2 = NonsingularCellPairAssemblerWithTrialCaching::new( + npts_test, + npts_trial, + 1, + &trial_cells, + &integrand, + &kernel, + grid.geometry_map(ReferenceCellType::Triangle, test_points.data()), + grid.geometry_map(ReferenceCellType::Triangle, trial_points.data()), + &test_table, + &trial_table, + &test_weights, + &trial_weights, + ); + + let mut result0 = rlst_dynamic_array2!(f64, [element.dim(), element.dim()]); + let mut result1 = rlst_dynamic_array2!(f64, [element.dim(), element.dim()]); + let mut result2 = rlst_dynamic_array2!(f64, [element.dim(), element.dim()]); + + for test_cell in &test_cells { + a0.set_test_cell(*test_cell); + a1.set_test_cell(*test_cell); + a2.set_test_cell(*test_cell); + for trial_cell in &trial_cells { + a0.set_trial_cell(*trial_cell); + a1.set_trial_cell(*trial_cell); + a2.set_trial_cell(*trial_cell); + + a0.assemble(&mut result0); + a1.assemble(&mut result1); + a2.assemble(&mut result2); + for (col0, col1, col2) in + izip!(result0.col_iter(), result1.col_iter(), result2.col_iter()) + { + for (value0, value1, value2) in izip!(col0.iter(), col1.iter(), col2.iter()) { + assert_relative_eq!(value0, value1); + assert_relative_eq!(value0, value2); + } + } + } + } + } +}