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, dealii::MemorySpace::Host> &heat_sources,
     double const time, double const next_refinement_time,
     unsigned int const time_steps_refinement,
     boost::property_tree::ptree const &refinement_database)
@@ -662,7 +659,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, dealii::MemorySpace::Host> &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 +757,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, dealii::MemorySpace::Host> heat_sources;
   if (use_thermal_physics)
   {
     // PropertyTreeInput discretization.thermal.fe_degree
@@ -1011,9 +1008,10 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database,
       {
         // Check if we have reached the end of current scan path.
         bool need_updated_scan_path = false;
-        for (auto &source : heat_sources)
+        auto const &scan_paths = heat_sources.get_scan_paths();
+        for (auto &scan_path : scan_paths)
         {
-          if (time > source->get_scan_path().get_segment_list().back().end_time)
+          if (time > scan_path.get_segment_list().back().end_time)
           {
             need_updated_scan_path = true;
             break;
@@ -1024,17 +1022,7 @@ run(MPI_Comm const &communicator, boost::property_tree::ptree const &database,
         {
           // Check if we have reached the end of the file. If not, read the
           // updated scan path file
-          bool scan_path_end = true;
-          for (auto &source : heat_sources)
-          {
-            if (!source->get_scan_path().is_finished())
-            {
-              scan_path_end = false;
-              // This functions waits for the scan path file to be updated
-              // before reading the file.
-              source->get_scan_path().read_file();
-            }
-          }
+          bool scan_path_end = heat_sources.update_scan_paths();
 
           // If we have reached the end of scan path file for all the heat
           // sources, we just exit.
@@ -1416,7 +1404,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;
@@ -1819,10 +1807,10 @@ run_ensemble(MPI_Comm const &global_communicator,
         // Check if we have reached the end of current scan path.
         bool need_updated_scan_path = false;
         {
-          for (auto &source : heat_sources_ensemble[0])
+          auto const &scan_paths = heat_sources_ensemble[0].get_scan_paths();
+          for (auto &scan_path : scan_paths)
           {
-            if (time >
-                source->get_scan_path().get_segment_list().back().end_time)
+            if (time > scan_path.get_segment_list().back().end_time)
             {
               need_updated_scan_path = true;
               break;
@@ -1841,16 +1829,7 @@ run_ensemble(MPI_Comm const &global_communicator,
           bool scan_path_end = true;
           for (unsigned int member = 0; member < local_ensemble_size; ++member)
           {
-            for (auto &source : heat_sources_ensemble[member])
-            {
-              if (!source->get_scan_path().is_finished())
-              {
-                scan_path_end = false;
-                // This functions waits for the scan path file to be updated
-                // before reading the file.
-                source->get_scan_path().read_file();
-              }
-            }
+            scan_path_end &= heat_sources_ensemble[member].update_scan_paths();
           }
 
           // If we have reached the end of scan path file for all the heat

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,24 +9,29 @@
 #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)
+    : _beam(database),
+      _scan_path(database.get<std::string>("scan_path_file"),
+                 database.get<std::string>("scan_path_file_format"))
 {
 }
 
 template <int dim>
 void ElectronBeamHeatSource<dim>::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);
 }
 
@@ -52,9 +57,11 @@ 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;

source/ElectronBeamHeatSource.hh

@@ -8,19 +8,20 @@
 #ifndef ELECTRON_BEAM_HEAT_SOURCE_HH
 #define ELECTRON_BEAM_HEAT_SOURCE_HH
 
-#include <HeatSource.hh>
+#include <BeamHeatSourceProperties.hh>
+#include <ScanPath.hh>
 
 #include <limits>
 
 namespace adamantine
 {
 /**
- * A derived class from HeatSource for a model of an electron beam heat source.
+ * A model of an electron beam heat source.
  * The form of the heat source model is taken from the following reference:
  * Raghavan et al, Acta Materilia, 112, 2016, pp 303-314.
  */
 template <int dim>
-class ElectronBeamHeatSource final : public HeatSource<dim>
+class ElectronBeamHeatSource
 {
 public:
   /**
@@ -38,20 +39,71 @@ public:
   /**
    * Set the time variable.
    */
-  void update_time(double time) final;
+  void update_time(double time);
 
   /**
    * Returns the value of an electron beam heat source at a specified 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;
+
+  /**
+   * Return the scan path.
+   */
+  ScanPath &get_scan_path();
+
+  void set_scan_path(ScanPath const scan_path) { _scan_path = scan_path; }
+
+  /**
+   * 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;
+
+  /**
+   * (Re)set 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 const &get_beam_properties() const;
 
 private:
   dealii::Point<3> _beam_center;
   double _alpha = std::numeric_limits<double>::signaling_NaN();
-  double const _log_01 = std::log(0.1);
+  BeamHeatSourceProperties _beam;
+  ScanPath _scan_path;
 };
+
+template <int dim>
+ScanPath &ElectronBeamHeatSource<dim>::get_scan_path()
+{
+  return _scan_path;
+}
+
+template <int dim>
+double ElectronBeamHeatSource<dim>::get_current_height(double const time) const
+{
+  return _scan_path.value(time)[2];
+}
+
+template <int dim>
+void ElectronBeamHeatSource<dim>::set_beam_properties(
+    boost::property_tree::ptree const &database)
+{
+  _beam.set_from_database(database);
+}
+
+template <int dim>
+BeamHeatSourceProperties const &
+ElectronBeamHeatSource<dim>::get_beam_properties() const
+{
+  return _beam;
+}
+
 } // namespace adamantine
 
 #endif