Add Portable MatrixFree initialization for GPU builds

CMakeLists.txt

@@ -118,6 +118,10 @@ if(PRISMS_PF_UNIT_TESTS OR PRISMS_PF_REGRESSION_TESTS OR PRISMS_PF_PERFORMANCE_T
   include(deps/catch2)
 endif()
 
+if(PRISMS_PF_GPU AND TARGET libassert-lib)
+  target_compile_options(libassert-lib PRIVATE --host-only)
+endif()
+
 # =========================================================
 # Configure targets
 # =========================================================
@@ -135,6 +139,14 @@ if(PRISMS_PF_ADDITIONAL_DEGREES)
 else()
   set(_additional_degrees "")
 endif()
+
+# deal.II Portable::MatrixFree is instantiated only for 2D and 3D.
+if(PRISMS_PF_GPU)
+  set(_space_dimensions "2; 3")
+else()
+  set(_space_dimensions "1; 2; 3")
+endif()
+
 configure_file(cmake/templates.in cmake/templates @ONLY)
 
 # TODO: Move this

cmake/functions/prisms_pf_configure_target.cmake

@@ -45,6 +45,29 @@ function(prisms_pf_configure_targets TARGETS)
           ON
     )
 
+    # Keep CUDA frontend forwarding options as pairs. Otherwise,
+    # target_compile_options() may de-duplicate repeated -Xcudafe flags.
+    if(PRISMS_PF_GPU)
+      set(_grouped_cxx_flags)
+      while(PRISMS_PF_CXX_FLAGS_LIST)
+        list(POP_FRONT PRISMS_PF_CXX_FLAGS_LIST _flag)
+
+        if(_flag STREQUAL "-Xcudafe")
+          if(NOT PRISMS_PF_CXX_FLAGS_LIST)
+            message(FATAL_ERROR "-Xcudafe is missing its argument")
+          endif()
+
+          list(POP_FRONT PRISMS_PF_CXX_FLAGS_LIST _cudafe_flag)
+          list(APPEND _grouped_cxx_flags "SHELL:-Xcudafe ${_cudafe_flag}")
+        else()
+          list(APPEND _grouped_cxx_flags "${_flag}")
+        endif()
+      endwhile()
+
+      set(PRISMS_PF_CXX_FLAGS_LIST ${_grouped_cxx_flags})
+    endif()
+
+
     # Add the compile flags, which we inherit from deal.II
     # TODO: We really shouldn't do this and have suggested
     # flags that get inherited from deal.II

cmake/scripts/CMakeLists.txt

@@ -2,6 +2,7 @@
 add_executable(expand_template_instantiations_exe expand_template_instantiations.cc)
 
 target_compile_features(expand_template_instantiations_exe PRIVATE cxx_std_20)
+target_compile_options(expand_template_instantiations_exe PRIVATE --host-only)
 
 set_target_properties(
   expand_template_instantiations_exe

cmake/templates.in

@@ -5,7 +5,7 @@
 REAL_SCALARS := { float; double }
 
 // Space dimensions
-SPACE_DIMENSIONS := { 1; 2; 3 }
+SPACE_DIMENSIONS := { @_space_dimensions@ }
 
 // Element degrees
 ELEMENT_DEGREE := { 1; 2@_additional_degrees@ }

include/prismspf/core/matrix_free_manager.h

@@ -18,6 +18,10 @@
 
 #include <prismspf/config.h>
 
+#ifdef PRISMS_PF_GPU
+#  include <deal.II/matrix_free/portable_matrix_free.h>
+#endif
+
 PRISMS_PF_BEGIN_NAMESPACE
 
 /**
@@ -115,6 +119,18 @@ class MatrixFreeManager
    * @brief Generic Matrix-free object with a scalar and vector entry on each level.
    */
   std::vector<MatrixFree<dim, number>> generic_matrix_free_levels;
+
+#ifdef PRISMS_PF_GPU
+  /**
+   * @brief Portable MatrixFree object for each field on each level.
+   *
+   * deal.II 9.6 Portable::MatrixFree is initialized with one DoFHandler,
+   * so the outer vector is indexed by relative level and the inner vector
+   * is indexed by field index.
+   */
+  std::vector<std::vector<dealii::Portable::MatrixFree<dim, number>>>
+    portable_shared_matrix_free_levels;
+#endif
 };
 
 template <unsigned int dim, typename number>
@@ -132,9 +148,19 @@ MatrixFreeManager<dim, number>::reinit(
     0,
     dealii::update_values | dealii::update_gradients | dealii::update_hessians |
       dealii::update_JxW_values | dealii::update_quadrature_points);
-
+#ifdef PRISMS_PF_GPU
+  using PortableAdditionalData =
+    typename dealii::Portable::MatrixFree<dim, number>::AdditionalData;
+  static const PortableAdditionalData portable_additional_data(
+    dealii::update_values | dealii::update_gradients | dealii::update_hessians |
+    dealii::update_JxW_values | dealii::update_quadrature_points);
+#endif
   const std::array<dealii::DoFHandler<dim>, 2> &generic_dof_handlers =
     dof_manager.get_dof_handlers();
+  const std::vector<std::array<dealii::DoFHandler<dim>, 2>> &generic_dof_handlers_levels =
+    dof_manager.get_dof_handlers_levels();
+  const std::vector<std::array<dealii::AffineConstraints<number>, 2>>
+    &generic_constraints_levels = constraint_manager.get_generic_constraints_levels();
 
   {
     const std::array<dealii::AffineConstraints<number>, 2> &generic_constraints =
@@ -158,6 +184,9 @@ MatrixFreeManager<dim, number>::reinit(
   const unsigned int num_levels = dof_manager.has_mg() ? dof_manager.num_levels() : 0;
   shared_matrix_free_levels.resize(num_levels);
   generic_matrix_free_levels.resize(num_levels);
+#ifdef PRISMS_PF_GPU
+  portable_shared_matrix_free_levels.resize(num_levels);
+#endif
   for (unsigned int relative_level = 0; relative_level < num_levels; ++relative_level)
     {
       const unsigned int       fine_level = num_levels - 1;
@@ -184,6 +213,29 @@ MatrixFreeManager<dim, number>::reinit(
       AdditionalData generic_additional_data;
       generic_additional_data.mg_level = level;
 
+#ifdef PRISMS_PF_GPU
+      // Reinit Portable MatrixFree. deal.II 9.6 initializes one DoFHandler
+      // per Portable::MatrixFree object, so create one object for each field.
+      const auto &field_dof_handlers = dof_manager.get_field_dof_handlers(relative_level);
+      const auto  field_constraints =
+        constraint_manager.get_field_constraints(relative_level);
+
+      auto &portable_shared_matrix_free =
+        portable_shared_matrix_free_levels[relative_level];
+      portable_shared_matrix_free.resize(field_dof_handlers.size());
+
+      for (unsigned int field_index = 0; field_index < field_dof_handlers.size();
+           ++field_index)
+        {
+          portable_shared_matrix_free[field_index].reinit(
+            SystemWide<dim, degree>::mapping,
+            *field_dof_handlers[field_index],
+            *field_constraints[field_index],
+            dealii::QGaussLobatto<1>(degree + 1),
+            portable_additional_data);
+        }
+#endif
+
       // Reinit generic MatrixFree
       generic_mg_matrix_free.reinit(
         SystemWide<dim, degree>::mapping,

include/prismspf/core/system_wide.h

@@ -3,6 +3,7 @@
 
 #pragma once
 
+#include <deal.II/base/mg_level_object.h>
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_system.h>
@@ -14,6 +15,8 @@
 
 #include <prismspf/config.h>
 
+#include <array>
+
 PRISMS_PF_BEGIN_NAMESPACE
 
 /**
@@ -37,9 +40,16 @@ class SystemWide
   /**
    * @brief Scalar and vector FE systems
    */
-  inline static const std::array<const dealii::FESystem<dim>, 2> fe_systems = {
-    dealii::FESystem<dim>(dealii::FE_Q<dim>(dealii::QGaussLobatto<1>(degree + 1)), 1),
-    dealii::FESystem<dim>(dealii::FE_Q<dim>(dealii::QGaussLobatto<1>(degree + 1)), dim)};
+  static const std::array<dealii::FESystem<dim>, 2> &
+  fe_systems()
+  {
+    static const std::array<dealii::FESystem<dim>, 2> value = {
+      dealii::FESystem<dim>(dealii::FE_Q<dim>(dealii::QGaussLobatto<1>(degree + 1)), 1),
+      dealii::FESystem<dim>(dealii::FE_Q<dim>(dealii::QGaussLobatto<1>(degree + 1)),
+                            dim)};
+
+    return value;
+  }
 
   /**
    * @brief Mappings to and from reference cell

include/prismspf/nucleation/nucleation_manager.h

@@ -101,7 +101,7 @@ NucleationManager<dim, degree, number>::attempt_nucleation(
   const double delta_t = nuc_params.nucleation_period * time_info.get_timestep();
   auto        &rng     = user_inputs.misc_parameters.rng;
 
-  static dealii::FEValues<dim> fe_values(SystemWide<dim, degree>::fe_systems[0],
+  static dealii::FEValues<dim> fe_values(SystemWide<dim, degree>::fe_systems()[0],
                                          SystemWide<dim, degree>::quadrature,
                                          dealii::UpdateFlags::update_values |
                                            dealii::UpdateFlags::update_JxW_values);

include/prismspf/nucleation/nucleus.h

@@ -4,6 +4,7 @@
 #pragma once
 
 #include <deal.II/base/point.h>
+#include <deal.II/base/vectorization.h>
 
 #include <prismspf/core/types.h>
 

src/core/dof_manager.cc

@@ -17,7 +17,7 @@ DoFManager<dim, degree>::reinit(const TriangulationManager<dim> &triangulation_m
     {
       dealii::DoFHandler<dim> &dof_handler = level_dof_handlers.at(rank);
       dof_handler.reinit(triangulation_manager.get_triangulation());
-      dof_handler.distribute_dofs(SystemWide<dim, degree>::fe_systems.at(rank));
+      dof_handler.distribute_dofs(SystemWide<dim, degree>::fe_systems().at(rank));
       if (init_mg)
         {
           dof_handler.distribute_mg_dofs();

src/core/refinement_manager.cc

@@ -241,12 +241,12 @@ RefinementManager<dim, degree, number>::mark_cells_for_refinement_and_coarsening
   // TODO: Better off using lazy evaluation
   static dealii::FEValues<dim, dim> scalar_fe_values(
     SystemWide<dim, degree>::mapping,
-    SystemWide<dim, degree>::fe_systems[0],
+    SystemWide<dim, degree>::fe_systems()[0],
     SystemWide<dim, degree>::quadrature,
     dealii::UpdateFlags::update_values | dealii::UpdateFlags::update_gradients);
   static dealii::FEValues<dim, dim> vector_fe_values(
     SystemWide<dim, degree>::mapping,
-    SystemWide<dim, degree>::fe_systems[1],
+    SystemWide<dim, degree>::fe_systems()[1],
     SystemWide<dim, degree>::quadrature,
     dealii::UpdateFlags::update_values | dealii::UpdateFlags::update_gradients);
 

src/grains/flood_filler.cc

@@ -133,7 +133,7 @@ FloodFiller<dim, degree, number>::recursive_flood_fill(
 
       cell->set_user_flag();
 
-      dealii::FEValues<dim> fe_values(SystemWide<dim, degree>::fe_systems[0],
+      dealii::FEValues<dim> fe_values(SystemWide<dim, degree>::fe_systems()[0],
                                       SystemWide<dim, degree>::quadrature,
                                       dealii::update_values);
 

src/grains/grain_reassignment.cc

@@ -113,7 +113,7 @@ GrainReassignmentManager<dim, degree, number>::reassign_grains(
     simplified_grain_representations,
     solve_context.get_solution_indexer(),
     solve_context.get_dof_manager().get_field_dof_handler(last_remapped_field_index),
-    SystemWide<dim, degree>::fe_systems[0].dofs_per_cell);
+    SystemWide<dim, degree>::fe_systems()[0].dofs_per_cell);
 
   Logger::instance() << "Reassigning grains completed.\n\n";
 }