wip

Author: jeremylt

examples/deal.II/CMakeLists.txt

@@ -11,13 +11,21 @@ IF(NOT ${deal.II_FOUND})
     )
 ENDIF()
 
-DEAL_II_INITIALIZE_CACHED_VARIABLES()
-PROJECT("bps")
+FILE(GLOB SOURCE_FILES "*.cc")
 
-DEAL_II_INITIALIZE_CACHED_VARIABLES()
+FOREACH ( source_file ${SOURCE_FILES} )
+  GET_FILENAME_COMPONENT(file_name ${source_file} NAME)
+  STRING( REPLACE ".cc" "" exec ${file_name} )
 
-ADD_EXECUTABLE(bps bps.cc)
-DEAL_II_SETUP_TARGET(bps)
+  DEAL_II_INITIALIZE_CACHED_VARIABLES()
+  PROJECT(${exec})
 
-TARGET_INCLUDE_DIRECTORIES(bps PUBLIC ${CEED_DIR}/include)
-TARGET_LINK_LIBRARIES(bps ${CEED_DIR}/lib/libceed.so)
+  DEAL_II_INITIALIZE_CACHED_VARIABLES()
+
+  ADD_EXECUTABLE(${exec} ${source_file})
+  DEAL_II_SETUP_TARGET(${exec})
+
+  TARGET_INCLUDE_DIRECTORIES(${exec} PUBLIC ${CEED_DIR}/include)
+  TARGET_LINK_LIBRARIES(${exec} ${CEED_DIR}/lib/libceed.so)
+
+ENDFOREACH ( source_file ${SOURCE_FILES} )

examples/deal.II/README.md

@@ -15,7 +15,11 @@ make
 To run the executable, write:
 
 ```
-./bps
+./bps_cpu
+```
+
+```
+./bps_kokkos
 ```
 
 Optional command-line arguments are shown by adding the command-line argument "--help".

examples/deal.II/bps.cc examples/deal.II/bps-cpu.ccexamples/deal.II/bps.cc renamed to examples/deal.II/bps-cpu.cc

@@ -47,6 +47,7 @@
 
 // include operators
 #include "bps.h"
+#include "bps-cpu.h"
 
 // Test cases
 //TESTARGS(name="BP1") --resource {ceed_resource} --bp BP1 --fe_degree 2 --print_timings 0
@@ -61,7 +62,7 @@ struct Parameters
   unsigned int n_global_refinements = 1;
   unsigned int fe_degree            = 2;
   bool         print_timings        = true;
-  std::string  libCEED_resource      = "/cpu/self";
+  std::string  libCEED_resource     = "/cpu/self";
 
   bool
   parse(int argc, char *argv[])

examples/deal.II/bps-cpu.h

@@ -0,0 +1,213 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+//  Authors: Peter Munch, Martin Kronbichler
+//
+// ---------------------------------------------------------------------
+
+// deal.II includes
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/mapping.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/shape_info.h>
+#include <deal.II/matrix_free/tools.h>
+
+// libCEED includes
+#include <ceed.h>
+#include <ceed/backend.h>
+
+// QFunction source
+#include "bps-qfunctions.h"
+
+using namespace dealii;
+
+
+template <int dim, typename Number>
+class OperatorDealii : public OperatorBase<Number>
+{
+public:
+  using VectorType = typename OperatorBase<Number>::VectorType;
+
+  /**
+   * Constructor.
+   */
+  OperatorDealii(const Mapping<dim>              &mapping,
+                 const DoFHandler<dim>           &dof_handler,
+                 const AffineConstraints<Number> &constraints,
+                 const Quadrature<dim>           &quadrature,
+                 const BPType                    &bp)
+    : mapping(mapping)
+    , dof_handler(dof_handler)
+    , constraints(constraints)
+    , quadrature(quadrature)
+    , bp(bp)
+  {
+    reinit();
+  }
+
+  /**
+   * Destructor.
+   */
+  ~OperatorDealii() = default;
+
+  /**
+   * Initialized internal data structures, particularly, MatrixFree.
+   */
+  void
+  reinit() override
+  {
+    // configure MatrixFree
+    typename MatrixFree<dim, Number>::AdditionalData additional_data;
+    additional_data.tasks_parallel_scheme =
+      MatrixFree<dim, Number>::AdditionalData::TasksParallelScheme::none;
+
+    // create MatrixFree
+    matrix_free.reinit(mapping, dof_handler, constraints, quadrature, additional_data);
+  }
+
+  /**
+   * Matrix-vector product.
+   */
+  void
+  vmult(VectorType &dst, const VectorType &src) const override
+  {
+    if (dof_handler.get_fe().n_components() == 1)
+      {
+        matrix_free.cell_loop(&OperatorDealii::do_cell_integral_range<1>, this, dst, src, true);
+      }
+    else
+      {
+        AssertThrow(dof_handler.get_fe().n_components() == dim, ExcInternalError());
+
+        matrix_free.cell_loop(&OperatorDealii::do_cell_integral_range<dim>, this, dst, src, true);
+      }
+  }
+
+  /**
+   * Initialize vector.
+   */
+  void
+  initialize_dof_vector(VectorType &vec) const override
+  {
+    matrix_free.initialize_dof_vector(vec);
+  }
+
+  /**
+   * Compute inverse of diagonal.
+   */
+  void
+  compute_inverse_diagonal(VectorType &diagonal) const override
+  {
+    this->initialize_dof_vector(diagonal);
+
+    if (dof_handler.get_fe().n_components() == 1)
+      {
+        MatrixFreeTools::compute_diagonal(matrix_free,
+                                          diagonal,
+                                          &OperatorDealii::do_cell_integral_local<1>,
+                                          this);
+      }
+    else
+      {
+        AssertThrow(dof_handler.get_fe().n_components() == dim, ExcInternalError());
+
+        MatrixFreeTools::compute_diagonal(matrix_free,
+                                          diagonal,
+                                          &OperatorDealii::do_cell_integral_local<dim>,
+                                          this);
+      }
+
+    for (auto &i : diagonal)
+      i = (std::abs(i) > 1.0e-10) ? (1.0 / i) : 1.0;
+  }
+
+private:
+  /**
+   * Cell integral without vector access.
+   */
+  template <int n_components>
+  void
+  do_cell_integral_local(FEEvaluation<dim, -1, 0, n_components, Number> &phi) const
+  {
+    if (bp <= BPType::BP2) // mass matrix
+      {
+        phi.evaluate(EvaluationFlags::values);
+        for (const auto q : phi.quadrature_point_indices())
+          phi.submit_value(phi.get_value(q), q);
+        phi.integrate(EvaluationFlags::values);
+      }
+    else // Poisson operator
+      {
+        phi.evaluate(EvaluationFlags::gradients);
+        for (const auto q : phi.quadrature_point_indices())
+          phi.submit_gradient(phi.get_gradient(q), q);
+        phi.integrate(EvaluationFlags::gradients);
+      }
+  }
+
+  /**
+   * Cell integral on a range of cells.
+   */
+  template <int n_components>
+  void
+  do_cell_integral_range(const MatrixFree<dim, Number>               &matrix_free,
+                         VectorType                                  &dst,
+                         const VectorType                            &src,
+                         const std::pair<unsigned int, unsigned int> &range) const
+  {
+    FEEvaluation<dim, -1, 0, n_components, Number> phi(matrix_free, range);
+
+    for (unsigned cell = range.first; cell < range.second; ++cell)
+      {
+        phi.reinit(cell);
+        phi.read_dof_values(src);            // read source vector
+        do_cell_integral_local(phi);         // cell integral
+        phi.distribute_local_to_global(dst); // write to destination vector
+      }
+  }
+
+  /**
+   * Mapping object passed to the constructor.
+   */
+  const Mapping<dim> &mapping;
+
+  /**
+   * DoFHandler object passed to the constructor.
+   */
+  const DoFHandler<dim> &dof_handler;
+
+  /**
+   * Constraints object passed to the constructor.
+   */
+  const AffineConstraints<Number> &constraints;
+
+  /**
+   * Quadrature rule object passed to the constructor.
+   */
+  const Quadrature<dim> &quadrature;
+
+  /**
+   * Selected BP.
+   */
+  const BPType bp;
+
+  /**
+   * MatrixFree object.
+   */
+  MatrixFree<dim, Number> matrix_free;
+};