Skip to content

Commit

Permalink
WIP: Testing of Laplace evaluator
Browse files Browse the repository at this point in the history
  • Loading branch information
@tbetcke
tbetcke committed Jan 5, 2025
1 parent c3ac7ef commit b2c2029
Show file tree
Hide file tree
Showing 4 changed files with 320 additions and 55 deletions.
5 changes: 5 additions & 0 deletions Cargo.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,13 +30,18 @@ itertools = "0.13.*"
mpi = { version = "0.8.*"}
num = "0.4"
ndelement = { git="https://github.com/bempp/ndelement.git", features = ["mpi"]}
# ndelement = { path = "../ndelement", features = ["mpi"] }
ndgrid = { git="https://github.com/bempp/ndgrid.git", features = ["serde"] }
# ndgrid = { path = "../ndgrid", features = ["serde"]}
# ndgrid = { path = "../ndgrid" }
rayon = "1.9"
rlst = { git = "https://github.com/linalg-rs/rlst.git", features = ["mpi"] }
# rlst = { path = "../rlst", features = ["mpi"] }
green-kernels = { git = "https://github.com/bempp/green-kernels.git", features = ["mpi"] }
# green-kernels = { path = "../green-kernels", features = ["mpi"] }
# c-api-tools = { version = "0.1.0" }
kifmm = { git = "https://github.com/bempp/kifmm.git", features = ["mpi"] }
# kifmm = { path = "../kifmm/kifmm", features = ["mpi"] }
bempp-distributed-tools = { git = "https://github.com/bempp/distributed_tools.git"}
rand = "0.8"
rand_chacha = "0.3"
Expand Down
140 changes: 125 additions & 15 deletions examples/laplace_evaluator.rs
View file
Original file line number Diff line number Diff line change
@@ -1,28 +1,50 @@
//! This file implements an example Laplace evaluator test the different involved operators.
use bempp::{boundary_assemblers::BoundaryAssemblerOptions, function::LocalFunctionSpaceTrait};
use bempp_quadrature::types::NumericalQuadratureGenerator;
use bempp::{
boundary_assemblers::BoundaryAssemblerOptions, evaluator_tools::NeighbourEvaluator,
function::LocalFunctionSpaceTrait,
};
use green_kernels::laplace_3d::Laplace3dKernel;
use mpi::traits::Communicator;
use ndelement::{ciarlet::LagrangeElementFamily, types::ReferenceCellType};
use ndgrid::traits::Grid;
use rlst::operator::interface::DistributedArrayVectorSpace;
use ndgrid::{
traits::{Entity, GeometryMap, Grid},
types::Ownership,
};
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;
use rlst::{
assert_array_relative_eq,
operator::interface::{
distributed_sparse_operator::DistributedCsrMatrixOperator, DistributedArrayVectorSpace,
},
rlst_dynamic_array1, AsApply, Element, IndexLayout, LinearSpace, MultInto, OperatorBase,
};

use rlst::prelude::*;

fn main() {
let universe = mpi::initialize().unwrap();
let world = universe.world();

let grid = bempp::shapes::regular_sphere::<f64, _>(3, 1, &world);
let mut rng = ChaCha8Rng::seed_from_u64(world.rank() as u64);

let grid = bempp::shapes::regular_sphere::<f64, _>(6, 1, &world);

let quad_degree = 6;
// Get the number of cells in the grid.

let n_cells = grid.entity_iter(2).count();

println!("Number of cells: {}", n_cells);

let space = bempp::function::FunctionSpace::new(
&grid,
&LagrangeElementFamily::<f64>::new(1, ndelement::types::Continuity::Standard),
&LagrangeElementFamily::<f64>::new(0, ndelement::types::Continuity::Discontinuous),
);

let mut options = BoundaryAssemblerOptions::default();
options.set_regular_quadrature_degree(ReferenceCellType::Triangle, 6);
options.set_regular_quadrature_degree(ReferenceCellType::Triangle, quad_degree);

let quad_degree = options
.get_regular_quadrature_degree(ReferenceCellType::Triangle)
Expand All @@ -39,26 +61,114 @@ fn main() {
let space_layout =
bempp_distributed_tools::SingleProcessIndexLayout::new(0, space.global_size(), &world);

let point_layout = bempp_distributed_tools::SingleProcessIndexLayout::new(
0,
quad_degree * space.global_size(),
&world,
);
let point_layout =
bempp_distributed_tools::SingleProcessIndexLayout::new(0, quad_degree * n_cells, &world);

// Instantiate function spaces.

let array_function_space = DistributedArrayVectorSpace::<_, f64>::new(&space_layout);
let point_function_space = DistributedArrayVectorSpace::<_, f64>::new(&point_layout);

let qrule = bempp_quadrature::simplex_rules::simplex_rule_triangle(6).unwrap();
let qrule = bempp_quadrature::simplex_rules::simplex_rule_triangle(quad_degree).unwrap();

println!("Here");
let space_to_point = bempp::evaluator_tools::basis_to_point_map(
&space,
&array_function_space,
&point_function_space,
&qrule.points,
&qrule.weights,
false,
);

let point_to_space = bempp::evaluator_tools::point_to_basis_map(
&space,
&point_function_space,
&array_function_space,
&qrule.points,
&qrule.weights,
);

// We now have to get all the points from the grid. We do this by iterating through all cells.

let mut points = vec![0 as f64; 3 * quad_degree * n_cells];

let geometry_map = grid.geometry_map(ReferenceCellType::Triangle, &qrule.points);

for cell in grid
.entity_iter(2)
.filter(|e| matches!(e.ownership(), Ownership::Owned))
{
let start_index = 3 * point_function_space
.index_layout()
.global2local(world.rank() as usize, qrule.npoints * cell.global_index())
.unwrap();
geometry_map.points(
cell.local_index(),
&mut points[start_index..start_index + 3 * qrule.npoints],
);
}

// let kernel_evaluator = bempp::greens_function_evaluators::dense_evaluator::DenseEvaluator::new(
// &points,
// &points,
// green_kernels::types::GreenKernelEvalType::Value,
// true,
// Laplace3dKernel::default(),
// &point_function_space,
// &point_function_space,
// );

let kernel_evaluator = bempp::greens_function_evaluators::kifmm_evaluator::KiFmmEvaluator::new(
&points,
&points,
1,
3,
5,
&point_function_space,
&point_function_space,
);

let correction = NeighbourEvaluator::new(
&qrule.points,
Laplace3dKernel::default(),
green_kernels::types::GreenKernelEvalType::Value,
&point_function_space,
&point_function_space,
&grid,
);

let corrected_evaluator = kernel_evaluator.r().sum(correction.r().scale(-1.0));

let prod1 = corrected_evaluator.r().product(space_to_point.r());

let singular_operator = DistributedCsrMatrixOperator::new(
assembler.assemble_singular(
&space,
array_function_space.index_layout(),
&space,
array_function_space.index_layout(),
),
&array_function_space,
&array_function_space,
);

let laplace_evaluator = (point_to_space.product(prod1)).sum(singular_operator.r());

let mut x = array_function_space.zero();
x.view_mut()
.local_mut()
.fill_from_equally_distributed(&mut rng);

let res = laplace_evaluator.apply(&x);

let res_local = res.view().local();

let mut expected = rlst_dynamic_array1!(f64, [space.global_size()]);

expected
.r_mut()
.simple_mult_into(dense_matrix.r(), x.view().local().r());

let rel_diff = (expected.r() - res_local.r()).norm_2() / expected.r().norm_2();

println!("Relative difference: {}", rel_diff);
}
49 changes: 41 additions & 8 deletions examples/neighbour_evaluator.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,15 +4,16 @@ use bempp::evaluator_tools::NeighbourEvaluator;
use bempp_distributed_tools::IndexLayoutFromLocalCounts;
use green_kernels::laplace_3d::Laplace3dKernel;
use mpi::traits::Communicator;
use ndelement::types::ReferenceCellType;
use ndgrid::{
traits::{Entity, Grid},
traits::{Entity, GeometryMap, Grid},
types::Ownership,
};
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;
use rlst::{
operator::interface::DistributedArrayVectorSpace, rlst_dynamic_array2, AsApply, Element,
LinearSpace, RawAccess,
IndexLayout, LinearSpace, RandomAccessMut, RawAccess,
};

fn main() {
Expand All @@ -39,7 +40,7 @@ fn main() {

let space = DistributedArrayVectorSpace::<_, f64>::new(&index_layout);

let evaluator = NeighbourEvaluator::new(
let neighbour_evaluator = NeighbourEvaluator::new(
points.data(),
Laplace3dKernel::default(),
green_kernels::types::GreenKernelEvalType::Value,
Expand All @@ -48,13 +49,45 @@ fn main() {
&grid,
);

// Create an element in the space.
// We now manually test the evaluator. For that we first create a dense evaluator so that we have comparison.

// For the evaluator we need all the points.

let mut physical_points = vec![0 as f64; 3 * n_points * n_cells];

let geometry_map = grid.geometry_map(ReferenceCellType::Triangle, points.data());

for cell in grid
.entity_iter(2)
.filter(|e| matches!(e.ownership(), Ownership::Owned))
{
let start_index = 3 * index_layout
.global2local(world.rank() as usize, n_points * cell.global_index())
.unwrap();
geometry_map.points(
cell.local_index(),
&mut physical_points[start_index..start_index + 3 * n_points],
);
}

let kernel_evaluator = bempp::greens_function_evaluators::dense_evaluator::DenseEvaluator::new(
&physical_points,
&physical_points,
green_kernels::types::GreenKernelEvalType::Value,
true,
Laplace3dKernel::default(),
&space,
&space,
);

let mut x = space.zero();

x.view_mut()
.local_mut()
.fill_from_equally_distributed(&mut rng);
*x.view_mut().local_mut().get_mut([0]).unwrap() = 1.0;
*x.view_mut().local_mut().get_mut([1]).unwrap() = 2.0;

let actual = neighbour_evaluator.apply(&x);
let expected = kernel_evaluator.apply(&x);

let _res = evaluator.apply(&x);
println!("Actual: {}", actual.view().local()[[0]]);
println!("Expected: {}", expected.view().local()[[0]]);
}
Loading

0 comments on commit b2c2029

Please sign in to comment.