Introduce HeatSources

application/adamantine.hh

@@ -520,7 +520,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<dealii::MemorySpace::Host, dim> &heat_sources,
     double const current_source_height, double const refinement_beam_cutoff)
 {
 
@@ -537,25 +537,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.max_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;
         }
       }
     }
@@ -572,7 +569,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<MemorySpaceType, dim> const &heat_sources,
     double const time, double const next_refinement_time,
     unsigned int const time_steps_refinement,
     boost::property_tree::ptree const &refinement_database)
@@ -609,14 +606,16 @@ void refine_mesh(
   const double refinement_beam_cutoff =
       refinement_database.get<double>("beam_cutoff", 1.0e-15);
 
+  adamantine::HeatSources<dealii::MemorySpace::Host, dim> host_heat_sources = heat_sources.copy_to_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);
+            host_heat_sources, current_source_height, refinement_beam_cutoff);
 
     // PropertyTreeInput refinement.coarsen_after_beam
     const bool coarsen_after_beam =
@@ -661,7 +660,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<MemorySpaceType, dim> const &heat_sources,
     double const time, double const next_refinement_time,
     unsigned int const time_steps_refinement,
     boost::property_tree::ptree const &refinement_database)
@@ -759,7 +758,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<MemorySpaceType, dim> heat_sources;
   if (use_thermal_physics)
   {
     // PropertyTreeInput discretization.thermal.fe_degree
@@ -936,11 +935,13 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database,
               mechanical_physics, displacement, material_properties, timers);
   ++n_time_step;
 
+  adamantine::HeatSources<dealii::MemorySpace::Host, dim> host_heat_sources = heat_sources.copy_to_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);
+                                                        host_heat_sources);
   // Extract the time-stepping database
   boost::property_tree::ptree time_stepping_database =
       database.get_child("time_stepping");
@@ -1364,7 +1365,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<MemorySpaceType, dim>>
       heat_sources_ensemble(local_ensemble_size);
 
   std::vector<std::unique_ptr<adamantine::Geometry<dim>>> geometry_ensemble;
@@ -1718,7 +1719,7 @@ run_ensemble(MPI_Comm const &global_communicator,
         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,
+                    heat_sources_ensemble[member].copy_to_host(), time, next_refinement_time,
                     time_steps_refinement, refinement_database);
         solution_augmented_ensemble[member].collect_sizes();
       }

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

@@ -9,11 +9,12 @@
 #include <instantiation.hh>
 #include <types.hh>
 
+#include <deal.II/base/memory_space.h>
+
 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 final
 {
 public:
   /**
@@ -38,18 +40,28 @@ 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;
+
+  /**
+   * (Re)sets the BeamHeatSourceProperties member variable, necessary if the
+   * beam parameters vary in time (e.g. due to data assimilation).
+   */
+  void set_beam_properties(boost::property_tree::ptree const &database);
+
+  /**
+   * Return the beam properties.
+   */
+  BeamHeatSourceProperties get_beam_properties() const;
 
 private:
   bool _source_on = false;
@@ -58,7 +70,23 @@ private:
   double _value;
   dealii::Point<dim> _min_point;
   dealii::Point<dim> _max_point;
+  double _alpha;
+  BeamHeatSourceProperties _beam;
 };
+
+template <int dim>
+void CubeHeatSource<dim>::set_beam_properties(
+    boost::property_tree::ptree const &database)
+{
+  _beam.set_from_database(database);
+}
+
+template <int dim>
+BeamHeatSourceProperties CubeHeatSource<dim>::get_beam_properties() const
+{
+  return _beam;
+}
+
 } // namespace adamantine
 
 #endif

source/DataAssimilator.hh

@@ -14,7 +14,7 @@
 #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/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)