Introduce HeatSources

application/adamantine.hh

@@ -521,7 +521,7 @@ compute_cells_to_refine(
     dealii::parallel::distributed::Triangulation<dim> &triangulation,
     double const time, double const next_refinement_time,
     unsigned int const n_time_steps,
-    std::vector<std::shared_ptr<adamantine::HeatSource<dim>>> &heat_sources,
+    adamantine::HeatSources<dim, dealii::MemorySpace::Host> &heat_sources,
     double const current_source_height, double const refinement_beam_cutoff)
 {
 
@@ -538,25 +538,22 @@ compute_cells_to_refine(
     double const current_time = time + static_cast<double>(i) /
                                            static_cast<double>(n_time_steps) *
                                            (next_refinement_time - time);
-    for (auto &beam : heat_sources)
+    heat_sources.update_time(current_time);
+    for (auto cell :
+         dealii::filter_iterators(triangulation.active_cell_iterators(),
+                                  dealii::IteratorFilters::LocallyOwnedCell()))
     {
-      beam->update_time(current_time);
-      for (auto cell : dealii::filter_iterators(
-               triangulation.active_cell_iterators(),
-               dealii::IteratorFilters::LocallyOwnedCell()))
+      // Check the value at the center of the cell faces. For most cases this
+      // should be sufficient, but if the beam is small compared to the
+      // coarsest mesh we may need to add other points to check (e.g.
+      // quadrature points, vertices).
+      for (unsigned int f = 0; f < cell->reference_cell().n_faces(); ++f)
       {
-        // Check the value at the center of the cell faces. For most cases this
-        // should be sufficient, but if the beam is small compared to the
-        // coarsest mesh we may need to add other points to check (e.g.
-        // quadrature points, vertices).
-        for (unsigned int f = 0; f < cell->reference_cell().n_faces(); ++f)
+        if (heat_sources.value(cell->face(f)->center(), current_source_height) >
+            refinement_beam_cutoff)
         {
-          if (beam->value(cell->face(f)->center(), current_source_height) >
-              refinement_beam_cutoff)
-          {
-            cells_to_refine.push_back(cell);
-            break;
-          }
+          cells_to_refine.push_back(cell);
+          break;
         }
       }
     }
@@ -573,7 +570,7 @@ void refine_mesh(
     adamantine::MaterialProperty<dim, p_order, MaterialStates, MemorySpaceType>
         &material_properties,
     dealii::LA::distributed::Vector<double, MemorySpaceType> &solution,
-    std::vector<std::shared_ptr<adamantine::HeatSource<dim>>> &heat_sources,
+    adamantine::HeatSources<dim, MemorySpaceType> const &heat_sources,
     double const time, double const next_refinement_time,
     unsigned int const time_steps_refinement,
     boost::property_tree::ptree const &refinement_database)
@@ -610,14 +607,17 @@ void refine_mesh(
   const double refinement_beam_cutoff =
       refinement_database.get<double>("beam_cutoff", 1.0e-15);
 
+  adamantine::HeatSources<dim, dealii::MemorySpace::Host> heat_sources_host =
+      heat_sources.copy_to(dealii::MemorySpace::Host{});
+
   for (unsigned int i = 0; i < n_refinements; ++i)
   {
     // Compute the cells to be refined.
     std::vector<typename dealii::parallel::distributed::Triangulation<
         dim>::active_cell_iterator>
         cells_to_refine = compute_cells_to_refine(
             triangulation, time, next_refinement_time, time_steps_refinement,
-            heat_sources, current_source_height, refinement_beam_cutoff);
+            heat_sources_host, current_source_height, refinement_beam_cutoff);
 
     // PropertyTreeInput refinement.coarsen_after_beam
     const bool coarsen_after_beam =
@@ -662,7 +662,7 @@ void refine_mesh(
     adamantine::MaterialProperty<dim, p_order, MaterialStates, MemorySpaceType>
         &material_properties,
     dealii::LA::distributed::Vector<double, MemorySpaceType> &solution,
-    std::vector<std::shared_ptr<adamantine::HeatSource<dim>>> &heat_sources,
+    adamantine::HeatSources<dim, MemorySpaceType> const &heat_sources,
     double const time, double const next_refinement_time,
     unsigned int const time_steps_refinement,
     boost::property_tree::ptree const &refinement_database)
@@ -760,7 +760,7 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database,
   // Create ThermalPhysics if necessary
   std::unique_ptr<adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>
       thermal_physics;
-  std::vector<std::shared_ptr<adamantine::HeatSource<dim>>> heat_sources;
+  adamantine::HeatSources<dim, MemorySpaceType> heat_sources;
   std::vector<std::pair<double, bool>> scan_path_end;
   if (use_thermal_physics)
   {
@@ -774,12 +774,16 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database,
         fe_degree, quadrature_type, communicator, database, geometry,
         material_properties);
     heat_sources = thermal_physics->get_heat_sources();
+    adamantine::HeatSources<dim, dealii::MemorySpace::Host> heat_sources_host =
+      heat_sources.copy_to(dealii::MemorySpace::Host{});
+
     // Store the current end time of each heat source and set a flag that the
     // scan path has changed
-    for (auto const &source : heat_sources)
+    auto const &scan_paths = heat_sources_host.get_scan_paths();
+    for (auto const &scan_path : scan_paths)
     {
-      scan_path_end.emplace_back(
-          source->get_scan_path().get_segment_list().back().end_time, true);
+      scan_path_end.emplace_back(scan_path.get_segment_list().back().end_time,
+                                 true);
     }
     post_processor_database.put("thermal_output", true);
   }
@@ -945,11 +949,14 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database,
               mechanical_physics, displacement, material_properties, timers);
   ++n_time_step;
 
+  adamantine::HeatSources<dim, dealii::MemorySpace::Host> heat_sources_host =
+      heat_sources.copy_to(dealii::MemorySpace::Host{});
+
   // Create the bounding boxes used for material deposition
   auto [material_deposition_boxes, deposition_times, deposition_cos,
         deposition_sin] =
       adamantine::create_material_deposition_boxes<dim>(geometry_database,
-                                                        heat_sources);
+                                                        heat_sources_host);
   // Extract the time-stepping database
   boost::property_tree::ptree time_stepping_database =
       database.get_child("time_stepping");
@@ -1025,25 +1032,26 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database,
       }
       if (scan_path_updated)
       {
+        adamantine::HeatSources<dim, dealii::MemorySpace::Host>
+            heat_sources_host =
+                heat_sources.copy_to(dealii::MemorySpace::Host{});
+        heat_sources_host.update_scan_paths();
+        heat_sources = heat_sources_host.copy_to(MemorySpaceType{});
+
+        auto const &scan_paths = heat_sources.get_scan_paths();
         for (unsigned int s = 0; s < scan_path_end.size(); ++s)
         {
-          // Read the scan path file
-          heat_sources[s]->get_scan_path().read_file();
-
           // Update scan_path_end
-          double new_end_time = heat_sources[s]
-                                    ->get_scan_path()
-                                    .get_segment_list()
-                                    .back()
-                                    .end_time;
+          double new_end_time =
+              scan_paths[s].get_segment_list().back().end_time;
           scan_path_end[s].second = scan_path_end[s].first != new_end_time;
           scan_path_end[s].first = new_end_time;
         }
 
         std::tie(material_deposition_boxes, deposition_times, deposition_cos,
                  deposition_sin) =
-            adamantine::create_material_deposition_boxes<dim>(geometry_database,
-                                                              heat_sources);
+            adamantine::create_material_deposition_boxes<dim>(
+                geometry_database, heat_sources_host);
       }
 
       auto activation_start =
@@ -1410,7 +1418,7 @@ run_ensemble(MPI_Comm const &global_communicator,
       adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>>
       thermal_physics_ensemble(local_ensemble_size);
 
-  std::vector<std::vector<std::shared_ptr<adamantine::HeatSource<dim>>>>
+  std::vector<adamantine::HeatSources<dim, MemorySpaceType>>
       heat_sources_ensemble(local_ensemble_size);
 
   std::vector<std::unique_ptr<adamantine::Geometry<dim>>> geometry_ensemble;
@@ -1761,11 +1769,13 @@ run_ensemble(MPI_Comm const &global_communicator,
 
       for (unsigned int member = 0; member < local_ensemble_size; ++member)
       {
-        refine_mesh(thermal_physics_ensemble[member],
-                    *material_properties_ensemble[member],
-                    solution_augmented_ensemble[member].block(base_state),
-                    heat_sources_ensemble[member], time, next_refinement_time,
-                    time_steps_refinement, refinement_database);
+        refine_mesh(
+            thermal_physics_ensemble[member],
+            *material_properties_ensemble[member],
+            solution_augmented_ensemble[member].block(base_state),
+            heat_sources_ensemble[member].copy_to(dealii::MemorySpace::Host{}),
+            time, next_refinement_time, time_steps_refinement,
+            refinement_database);
         solution_augmented_ensemble[member].collect_sizes();
       }
 

indent

@@ -2,3 +2,4 @@ clang-format -style=file -i source/*.cc
 clang-format -style=file -i source/*.hh
 clang-format -style=file -i tests/*.cc
 clang-format -style=file -i application/*.cc
+clang-format -style=file -i application/*.hh

source/CMakeLists.txt

@@ -7,7 +7,7 @@ set(Adamantine_HEADERS
     ${CMAKE_CURRENT_SOURCE_DIR}/ExperimentalData.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/Geometry.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/GoldakHeatSource.hh
-    ${CMAKE_CURRENT_SOURCE_DIR}/HeatSource.hh
+    ${CMAKE_CURRENT_SOURCE_DIR}/HeatSources.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/ImplicitOperator.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/MaterialProperty.hh
     ${CMAKE_CURRENT_SOURCE_DIR}/MaterialProperty.templates.hh

source/CubeHeatSource.cc

@@ -13,7 +13,6 @@ namespace adamantine
 {
 template <int dim>
 CubeHeatSource<dim>::CubeHeatSource(boost::property_tree::ptree const &database)
-    : HeatSource<dim>()
 {
   _start_time = database.get<double>("start_time");
   _end_time = database.get<double>("end_time");

source/CubeHeatSource.hh

@@ -8,7 +8,9 @@
 #ifndef CUBE_HEAT_SOURCE_HH
 #define CUBE_HEAT_SOURCE_HH
 
-#include <HeatSource.hh>
+#include <BeamHeatSourceProperties.hh>
+
+#include <deal.II/base/point.h>
 
 namespace adamantine
 {
@@ -17,7 +19,7 @@ namespace adamantine
  * used for verification purpose.
  */
 template <int dim>
-class CubeHeatSource final : public HeatSource<dim>
+class CubeHeatSource
 {
 public:
   /**
@@ -38,18 +40,17 @@ public:
   /**
    * Set the time variable.
    */
-  void update_time(double time) final;
+  void update_time(double time);
 
   /**
    * Return the value of the source for a given point and time.
    */
-  double value(dealii::Point<dim> const &point,
-               double const /*height*/) const final;
+  double value(dealii::Point<dim> const &point, double const /*height*/) const;
   /**
    * Compute the current height of the where the heat source meets the material
    * (i.e. the current scan path height).
    */
-  double get_current_height(double const time) const final;
+  double get_current_height(double const time) const;
 
 private:
   bool _source_on = false;
@@ -58,7 +59,9 @@ private:
   double _value;
   dealii::Point<dim> _min_point;
   dealii::Point<dim> _max_point;
+  double _alpha;
 };
+
 } // namespace adamantine
 
 #endif

source/DataAssimilator.hh

@@ -14,7 +14,6 @@
 #include <deal.II/fe/mapping_q1_eulerian.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/la_parallel_block_vector.h>
-#include <deal.II/lac/la_vector.h>
 #include <deal.II/lac/solver_gmres.h>
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/trilinos_sparse_matrix.h>

source/ElectronBeamHeatSource.cc

@@ -9,30 +9,34 @@
 #include <instantiation.hh>
 #include <types.hh>
 
+#include <deal.II/base/memory_space.h>
+
 namespace adamantine
 {
 
-template <int dim>
-ElectronBeamHeatSource<dim>::ElectronBeamHeatSource(
-    boost::property_tree::ptree const &database)
-    : HeatSource<dim>(database)
+template <int dim, typename MemorySpaceType>
+ElectronBeamHeatSource<dim, MemorySpaceType>::ElectronBeamHeatSource(
+    BeamHeatSourceProperties const &beam,
+    ScanPath<MemorySpaceType> const &scan_path)
+    : _beam(beam), _scan_path(scan_path)
 {
 }
 
-template <int dim>
-void ElectronBeamHeatSource<dim>::update_time(double time)
+template <int dim, typename MemorySpaceType>
+void ElectronBeamHeatSource<dim, MemorySpaceType>::update_time(double time)
 {
+  static const double log_01 = std::log(0.1);
   _beam_center = this->_scan_path.value(time);
   double segment_power_modifier = this->_scan_path.get_power_modifier(time);
   _alpha =
       -this->_beam.absorption_efficiency * this->_beam.max_power *
-      segment_power_modifier * _log_01 /
+      segment_power_modifier * log_01 /
       (dealii::numbers::PI * this->_beam.radius_squared * this->_beam.depth);
 }
 
-template <int dim>
-double ElectronBeamHeatSource<dim>::value(dealii::Point<dim> const &point,
-                                          double const height) const
+template <int dim, typename MemorySpaceType>
+double ElectronBeamHeatSource<dim, MemorySpaceType>::value(
+    dealii::Point<dim> const &point, double const height) const
 {
   double const z = point[axis<dim>::z] - height;
   if ((z + this->_beam.depth) < 0.)
@@ -52,14 +56,17 @@ double ElectronBeamHeatSource<dim>::value(dealii::Point<dim> const &point,
           std::pow(point[axis<dim>::y] - _beam_center[axis<dim>::y], 2);
     }
 
+    static const double log_01 = std::log(0.1);
+
     // Electron beam heat source equation
     double heat_source =
-        _alpha * std::exp(_log_01 * xpy_squared / this->_beam.radius_squared) *
+        _alpha * std::exp(log_01 * xpy_squared / this->_beam.radius_squared) *
         distribution_z;
 
     return heat_source;
   }
 }
 } // namespace adamantine
 
-INSTANTIATE_DIM(ElectronBeamHeatSource)
+INSTANTIATE_DIM_DEVICE(ElectronBeamHeatSource)
+INSTANTIATE_DIM_HOST(ElectronBeamHeatSource)