Add assembly benchmark

bem/Cargo.toml

@@ -31,3 +31,11 @@ itertools = "0.10"
 mpi = { version = "0.6.*", optional = true }
 num = "0.4"
 rayon = "1.7"
+
+[dev-dependencies]
+criterion = "0.3"
+
+[[bench]]
+name = "assembly_benchmark"
+harness = false
+

bem/benches/assembly_benchmark.rs

@@ -0,0 +1,52 @@
+use bempp_bem::assembly::{assemble_batched, BoundaryOperator, PDEType};
+use bempp_bem::function_space::SerialFunctionSpace;
+use bempp_element::element::create_element;
+use bempp_grid::shapes::regular_sphere;
+use bempp_tools::arrays::Array2D;
+use bempp_traits::bem::DofMap;
+use bempp_traits::bem::FunctionSpace;
+use bempp_traits::cell::ReferenceCellType;
+use bempp_traits::element::{Continuity, ElementFamily};
+use criterion::{criterion_group, criterion_main, Criterion};
+
+pub fn assembly_benchmark(c: &mut Criterion) {
+    let mut group = c.benchmark_group("assembly");
+    group.sample_size(10);
+
+    for i in 3..5 {
+        let grid = regular_sphere(i);
+        let element = create_element(
+            ElementFamily::Lagrange,
+            ReferenceCellType::Triangle,
+            0,
+            Continuity::Discontinuous,
+        );
+
+        let space = SerialFunctionSpace::new(&grid, &element);
+        let mut matrix =
+            Array2D::<f64>::new((space.dofmap().global_size(), space.dofmap().global_size()));
+
+        group.bench_function(
+            &format!(
+                "Assembly of {}x{} matrix",
+                space.dofmap().global_size(),
+                space.dofmap().global_size()
+            ),
+            |b| {
+                b.iter(|| {
+                    assemble_batched(
+                        &mut matrix,
+                        BoundaryOperator::SingleLayer,
+                        PDEType::Laplace,
+                        &space,
+                        &space,
+                    )
+                })
+            },
+        );
+    }
+    group.finish();
+}
+
+criterion_group!(benches, assembly_benchmark);
+criterion_main!(benches);

bem/examples/assembly_timing.rs

@@ -1,5 +1,3 @@
-// TODO: use Criterion to benchmark properly
-
 use bempp_bem::assembly::{assemble_batched, BoundaryOperator, PDEType};
 use bempp_bem::function_space::SerialFunctionSpace;
 use bempp_element::element::create_element;

bem/src/assembly/batched.rs

@@ -13,7 +13,7 @@ use bempp_traits::element::FiniteElement;
 use bempp_traits::grid::{Geometry, Grid, Topology};
 use bempp_traits::types::Scalar;
 use rayon::prelude::*;
-use std::time::Instant;
+//use std::time::Instant;
 
 fn get_quadrature_rule(
     test_celltype: ReferenceCellType,
@@ -334,7 +334,7 @@ pub fn assemble<'a>(
     trial_space: &SerialFunctionSpace<'a>,
     test_space: &SerialFunctionSpace<'a>,
 ) {
-    let now = Instant::now();
+    // let now = Instant::now();
     // Note: currently assumes that the two grids are the same
     // TODO: implement == and != for grids, then add:
     // if *trial_space.grid() != *test_space.grid() {
@@ -440,8 +440,8 @@ pub fn assemble<'a>(
         }
     }
 
-    println!("Non-adjacent terms: {}ms", now.elapsed().as_millis());
-    let now = Instant::now();
+    // println!("Non-adjacent terms: {}ms", now.elapsed().as_millis());
+    // let now = Instant::now();
 
     if test_space.grid() != trial_space.grid() {
         // If the test and trial grids are different, there are no neighbouring triangles
@@ -572,7 +572,7 @@ pub fn assemble<'a>(
             }
         }
     }
-    println!("Singular terms: {}ms", now.elapsed().as_millis());
+    // println!("Singular terms: {}ms", now.elapsed().as_millis());
 }
 #[cfg(test)]
 mod test {