adamantine
adamantine.hh
Go to the documentation of this file.
1 /* SPDX-FileCopyrightText: Copyright (c) 2016 - 2025, the adamantine authors.
2  * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3  */
4 
5 #ifndef ADAMANTINE_HH
6 #define ADAMANTINE_HH
7 
8 #include <Boundary.hh>
9 #include <DataAssimilator.hh>
10 #include <ExperimentalData.hh>
11 #include <Geometry.hh>
12 #include <MaterialProperty.hh>
13 #include <MechanicalPhysics.hh>
14 #include <Microstructure.hh>
15 #include <PointCloud.hh>
16 #include <PostProcessor.hh>
17 #include <RayTracing.hh>
18 #include <ThermalPhysics.hh>
19 #include <ThermalPhysicsInterface.hh>
20 #include <Timer.hh>
21 #include <ensemble_management.hh>
22 #include <experimental_data_utils.hh>
23 #include <material_deposition.hh>
24 #include <types.hh>
25 #include <utils.hh>
26 
27 #include <deal.II/arborx/bvh.h>
28 #include <deal.II/base/index_set.h>
29 #include <deal.II/base/mpi.h>
30 #include <deal.II/base/symmetric_tensor.h>
31 #include <deal.II/base/types.h>
32 #include <deal.II/distributed/cell_data_transfer.templates.h>
33 #include <deal.II/distributed/solution_transfer.h>
34 #include <deal.II/grid/filtered_iterator.h>
35 #include <deal.II/grid/grid_refinement.h>
36 #include <deal.II/lac/trilinos_sparse_matrix.h>
37 #include <deal.II/lac/vector_operation.h>
38 #include <deal.II/numerics/error_estimator.h>
39 
40 #include <boost/algorithm/string.hpp>
41 #include <boost/archive/text_iarchive.hpp>
42 #include <boost/archive/text_oarchive.hpp>
43 #include <boost/property_tree/ptree.hpp>
44 
45 #include <fstream>
46 #include <limits>
47 #include <memory>
48 #include <string>
49 #include <type_traits>
50 #include <unordered_map>
51 
52 #ifdef ADAMANTINE_WITH_CALIPER
53 #include <caliper/cali.h>
54 #endif
55 
56 #include <cmath>
57 #include <iostream>
58 
59 template <int dim, int n_materials, int p_order, typename MaterialStates,
60  typename MemorySpaceType,
61  std::enable_if_t<
62  std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value,
63  int> = 0>
64 void output_pvtu(
65  adamantine::PostProcessor<dim> &post_processor, unsigned int n_time_step,
66  double time,
67  std::unique_ptr<
68  adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>> const
69  &thermal_physics,
70  dealii::LinearAlgebra::distributed::Vector<double, MemorySpaceType>
71  &temperature,
72  std::unique_ptr<adamantine::MechanicalPhysics<
73  dim, n_materials, p_order, MaterialStates, MemorySpaceType>> const
74  &mechanical_physics,
75  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
76  &displacement,
77  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
78  MemorySpaceType> const &material_properties,
79  std::vector<adamantine::Timer> &timers)
80 {
81 #ifdef ADAMANTINE_WITH_CALIPER
82  CALI_CXX_MARK_FUNCTION;
83 #endif
84  timers[adamantine::output].start();
85  if (thermal_physics)
86  {
87  thermal_physics->get_affine_constraints().distribute(temperature);
88  if (mechanical_physics)
89  {
90  mechanical_physics->get_affine_constraints().distribute(displacement);
91  post_processor.template write_output<typename Kokkos::View<
92  double **, typename MemorySpaceType::kokkos_space>::array_layout>(
93  n_time_step, time, temperature, displacement,
94  mechanical_physics->get_stress_tensor(),
95  material_properties.get_state(), material_properties.get_dofs_map(),
96  material_properties.get_dof_handler());
97  }
98  else
99  {
100  post_processor.template write_thermal_output<typename Kokkos::View<
101  double **, typename MemorySpaceType::kokkos_space>::array_layout>(
102  n_time_step, time, temperature, material_properties.get_state(),
103  material_properties.get_dofs_map(),
104  material_properties.get_dof_handler());
105  }
106  }
107  else
108  {
109  mechanical_physics->get_affine_constraints().distribute(displacement);
110  post_processor.template write_mechanical_output<typename Kokkos::View<
111  double **, typename MemorySpaceType::kokkos_space>::array_layout>(
112  n_time_step, time, displacement,
113  mechanical_physics->get_stress_tensor(),
114  material_properties.get_state(), material_properties.get_dofs_map(),
115  material_properties.get_dof_handler());
116  }
117  timers[adamantine::output].stop();
118 }
119 
120 template <int dim, int n_materials, int p_order, typename MaterialStates,
121  typename MemorySpaceType,
122  std::enable_if_t<std::is_same<MemorySpaceType,
123  dealii::MemorySpace::Default>::value,
124  int> = 0>
125 void output_pvtu(
126  adamantine::PostProcessor<dim> &post_processor, unsigned int n_time_step,
127  double time,
128  std::unique_ptr<
129  adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>> const
130  &thermal_physics,
131  dealii::LinearAlgebra::distributed::Vector<double, MemorySpaceType>
132  &temperature,
133  std::unique_ptr<adamantine::MechanicalPhysics<
134  dim, n_materials, p_order, MaterialStates, MemorySpaceType>> const
135  &mechanical_physics,
136  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
137  &displacement,
138  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
139  MemorySpaceType> const &material_properties,
140  std::vector<adamantine::Timer> &timers)
141 {
142 #ifdef ADAMANTINE_WITH_CALIPER
143  CALI_CXX_MARK_FUNCTION;
144 #endif
145  timers[adamantine::output].start();
146  auto state = material_properties.get_state();
147  auto state_host = Kokkos::create_mirror_view_and_copy(
148  dealii::MemorySpace::Host::kokkos_space{}, state);
149  if (thermal_physics)
150  {
151  dealii::LinearAlgebra::distributed::Vector<double,
152  dealii::MemorySpace::Host>
153  temperature_host(temperature.get_partitioner());
154  temperature_host.import(temperature, dealii::VectorOperation::insert);
155  thermal_physics->get_affine_constraints().distribute(temperature_host);
156  if (mechanical_physics)
157  {
158  mechanical_physics->get_affine_constraints().distribute(displacement);
159  post_processor.template write_output<typename Kokkos::View<
160  double **, typename MemorySpaceType::kokkos_space>::array_layout>(
161  n_time_step, time, temperature_host, displacement,
162  mechanical_physics->get_stress_tensor(), state_host,
163  material_properties.get_dofs_map(),
164  material_properties.get_dof_handler());
165  }
166  else
167  {
168  post_processor.template write_thermal_output<typename Kokkos::View<
169  double **, typename MemorySpaceType::kokkos_space>::array_layout>(
170  n_time_step, time, temperature_host, state_host,
171  material_properties.get_dofs_map(),
172  material_properties.get_dof_handler());
173  }
174  }
175  else
176  {
177  mechanical_physics->get_affine_constraints().distribute(displacement);
178  post_processor.template write_mechanical_output<typename Kokkos::View<
179  double **, typename MemorySpaceType::kokkos_space>::array_layout>(
180  n_time_step, time, displacement,
181  mechanical_physics->get_stress_tensor(), state_host,
182  material_properties.get_dofs_map(),
183  material_properties.get_dof_handler());
184  }
185  timers[adamantine::output].stop();
186 }
187 
188 // inlining this function so we can have in the header
189 inline void initialize_timers(MPI_Comm const &communicator,
190  std::vector<adamantine::Timer> &timers)
191 {
192  timers.push_back(adamantine::Timer(communicator, "Main"));
193  timers.push_back(adamantine::Timer(communicator, "Refinement"));
194  timers.push_back(adamantine::Timer(communicator, "Add Material, Search"));
195  timers.push_back(adamantine::Timer(communicator, "Add Material, Activate"));
196  timers.push_back(
197  adamantine::Timer(communicator, "Data Assimilation, Exp. Data"));
198  timers.push_back(
199  adamantine::Timer(communicator, "Data Assimilation, DOF Mapping"));
200  timers.push_back(
201  adamantine::Timer(communicator, "Data Assimilation, Cov. Sparsity"));
202  timers.push_back(
203  adamantine::Timer(communicator, "Data Assimilation, Exp. Cov."));
204  timers.push_back(
205  adamantine::Timer(communicator, "Data Assimilation, Update Ensemble"));
206  timers.push_back(adamantine::Timer(communicator, "Evolve One Time Step"));
207  timers.push_back(adamantine::Timer(
208  communicator, "Evolve One Time Step: evaluate_thermal_physics"));
209  timers.push_back(adamantine::Timer(
210  communicator, "Evolve One Time Step: evaluate_material_properties"));
211  timers.push_back(adamantine::Timer(communicator, "Output"));
212 }
213 
214 template <int dim, int n_materials, int p_order, int fe_degree,
215  typename MaterialStates, typename MemorySpaceType,
216  typename QuadratureType>
217 std::unique_ptr<adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>
218 initialize(
219  MPI_Comm const &communicator, boost::property_tree::ptree const &database,
220  adamantine::Geometry<dim> &geometry, adamantine::Boundary const &boundary,
221  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
222  MemorySpaceType> &material_properties)
223 {
224  return std::make_unique<adamantine::ThermalPhysics<
225  dim, n_materials, p_order, fe_degree, MaterialStates, MemorySpaceType,
226  QuadratureType>>(communicator, database, geometry, boundary,
227  material_properties);
228 }
229 
230 template <int dim, int n_materials, int p_order, int fe_degree,
231  typename MaterialStates, typename MemorySpaceType>
232 std::unique_ptr<adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>
233 initialize_quadrature(
234  std::string const &quadrature_type, MPI_Comm const &communicator,
235  boost::property_tree::ptree const &database,
236  adamantine::Geometry<dim> &geometry, adamantine::Boundary const &boundary,
237  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
238  MemorySpaceType> &material_properties)
239 {
240  if (quadrature_type.compare("gauss") == 0)
241  return initialize<dim, n_materials, p_order, fe_degree, MaterialStates,
242  MemorySpaceType, dealii::QGauss<1>>(
243  communicator, database, geometry, boundary, material_properties);
244  else
245  {
246  adamantine::ASSERT_THROW(quadrature_type.compare("lobatto") == 0,
247  "quadrature should be Gauss or Lobatto.");
248  return initialize<dim, n_materials, p_order, fe_degree, MaterialStates,
249  MemorySpaceType, dealii::QGaussLobatto<1>>(
250  communicator, database, geometry, boundary, material_properties);
251  }
252 }
253 
254 template <int dim, int n_materials, int p_order, typename MaterialStates,
255  typename MemorySpaceType>
256 std::unique_ptr<adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>
257 initialize_thermal_physics(
258  unsigned int fe_degree, std::string const &quadrature_type,
259  MPI_Comm const &communicator, boost::property_tree::ptree const &database,
260  adamantine::Geometry<dim> &geometry, adamantine::Boundary const &boundary,
261  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
262  MemorySpaceType> &material_properties)
263 {
264  switch (fe_degree)
265  {
266  case 1:
267  {
268  return initialize_quadrature<dim, n_materials, p_order, 1, MaterialStates,
269  MemorySpaceType>(quadrature_type, communicator,
270  database, geometry, boundary,
271  material_properties);
272  }
273  case 2:
274  {
275  return initialize_quadrature<dim, n_materials, p_order, 2, MaterialStates,
276  MemorySpaceType>(quadrature_type, communicator,
277  database, geometry, boundary,
278  material_properties);
279  }
280  case 3:
281  {
282  return initialize_quadrature<dim, n_materials, p_order, 3, MaterialStates,
283  MemorySpaceType>(quadrature_type, communicator,
284  database, geometry, boundary,
285  material_properties);
286  }
287  case 4:
288  {
289  return initialize_quadrature<dim, n_materials, p_order, 4, MaterialStates,
290  MemorySpaceType>(quadrature_type, communicator,
291  database, geometry, boundary,
292  material_properties);
293  }
294  default:
295  {
296  adamantine::ASSERT_THROW(fe_degree == 5,
297  "fe_degree should be between 1 and 5.");
298  return initialize_quadrature<dim, n_materials, p_order, 5, MaterialStates,
299  MemorySpaceType>(quadrature_type, communicator,
300  database, geometry, boundary,
301  material_properties);
302  }
303  }
304 }
305 
306 template <int dim, int n_materials, int p_order, typename MaterialStates,
307  typename MemorySpaceType>
308 void refine_and_transfer(
309  std::unique_ptr<adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>
310  &thermal_physics,
311  std::unique_ptr<adamantine::MechanicalPhysics<
312  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>
313  &mechanical_physics,
314  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
315  MemorySpaceType> &material_properties,
316  dealii::DoFHandler<dim> &dof_handler,
317  dealii::LA::distributed::Vector<double, MemorySpaceType> &solution)
318 {
319 #ifdef ADAMANTINE_WITH_CALIPER
320  CALI_CXX_MARK_FUNCTION;
321 #endif
322 
323  dealii::parallel::distributed::Triangulation<dim> &triangulation =
324  dynamic_cast<dealii::parallel::distributed::Triangulation<dim> &>(
325  const_cast<dealii::Triangulation<dim> &>(
326  dof_handler.get_triangulation()));
327 
328  // Update the material state from the ThermalOperator to MaterialProperty
329  // because, for now, we need to use state from MaterialProperty to perform the
330  // transfer to the refined mesh.
331  thermal_physics->set_state_to_material_properties();
332 
333  // Transfer of the solution
334  dealii::parallel::distributed::SolutionTransfer<
335  dim, dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>>
336  solution_transfer(dof_handler);
337 
338  // Transfer material state
339  unsigned int const direction_data_size = 2;
340  unsigned int const phase_history_data_size = 1;
341  unsigned int constexpr n_material_states = MaterialStates::n_material_states;
342  std::vector<std::vector<double>> data_to_transfer;
343  std::vector<double> dummy_cell_data(n_material_states + direction_data_size +
344  phase_history_data_size,
345  std::numeric_limits<double>::infinity());
346  auto state_host = Kokkos::create_mirror_view_and_copy(
347  Kokkos::HostSpace{}, material_properties.get_state());
348  unsigned int cell_id = 0;
349  unsigned int activated_cell_id = 0;
350  for (auto const &cell : dof_handler.active_cell_iterators())
351  {
352  if (cell->is_locally_owned())
353  {
354  std::vector<double> cell_data(n_material_states + direction_data_size +
355  phase_history_data_size);
356  for (unsigned int i = 0; i < n_material_states; ++i)
357  cell_data[i] = state_host(i, cell_id);
358  if (cell->active_fe_index() == 0)
359  {
360  cell_data[n_material_states] =
361  thermal_physics->get_deposition_cos(activated_cell_id);
362  cell_data[n_material_states + 1] =
363  thermal_physics->get_deposition_sin(activated_cell_id);
364 
365  if (thermal_physics->get_has_melted(activated_cell_id))
366  cell_data[n_material_states + direction_data_size] = 1.0;
367  else
368  cell_data[n_material_states + direction_data_size] = 0.0;
369 
370  ++activated_cell_id;
371  }
372  else
373  {
374  cell_data[n_material_states] = std::numeric_limits<double>::infinity();
375  cell_data[n_material_states + 1] =
376  std::numeric_limits<double>::infinity();
377  cell_data[n_material_states + direction_data_size] =
378  std::numeric_limits<double>::infinity();
379  }
380  data_to_transfer.push_back(cell_data);
381  ++cell_id;
382  }
383  else
384  {
385  data_to_transfer.push_back(dummy_cell_data);
386  }
387  }
388 
389  // Prepare the Triangulation and the diffent data transfer objects for
390  // refinement
391  triangulation.prepare_coarsening_and_refinement();
392  // Prepare for refinement of the solution
393  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
394  solution_host;
395  if constexpr (std::is_same_v<MemorySpaceType, dealii::MemorySpace::Host>)
396  {
397  // We need to apply the constraints before the mesh transfer
398  thermal_physics->get_affine_constraints().distribute(solution);
399  // We need to update the ghost values before we can do the interpolation on
400  // the new mesh.
401  solution.update_ghost_values();
402  solution_transfer.prepare_for_coarsening_and_refinement(solution);
403  }
404  else
405  {
406  solution_host.reinit(solution.get_partitioner());
407  solution_host.import(solution, dealii::VectorOperation::insert);
408  // We need to apply the constraints before the mesh transfer
409  thermal_physics->get_affine_constraints().distribute(solution_host);
410  // We need to update the ghost values before we can do the interpolation on
411  // the new mesh.
412  solution_host.update_ghost_values();
413  solution_transfer.prepare_for_coarsening_and_refinement(solution_host);
414  }
415 
416  dealii::parallel::distributed::CellDataTransfer<
417  dim, dim, std::vector<std::vector<double>>>
418  cell_data_trans(triangulation);
419  cell_data_trans.prepare_for_coarsening_and_refinement(data_to_transfer);
420 
421  if (mechanical_physics)
422  {
423  mechanical_physics->prepare_transfer_mpi();
424  }
425 
426 #ifdef ADAMANTINE_WITH_CALIPER
427  CALI_MARK_BEGIN("refine triangulation");
428 #endif
429  // Execute the refinement
430  triangulation.execute_coarsening_and_refinement();
431 #ifdef ADAMANTINE_WITH_CALIPER
432  CALI_MARK_END("refine triangulation");
433 #endif
434 
435  // Update the AffineConstraints and resize the solution
436  thermal_physics->setup_dofs();
437  thermal_physics->initialize_dof_vector(0., solution);
438 
439  // Update MaterialProperty DoFHandler and resize the state vectors
440  material_properties.reinit_dofs();
441 
442  // Interpolate the solution
443  if constexpr (std::is_same_v<MemorySpaceType, dealii::MemorySpace::Host>)
444  {
445  solution_transfer.interpolate(solution);
446  }
447  else
448  {
449  solution_host.reinit(solution.get_partitioner());
450  solution_transfer.interpolate(solution_host);
451  solution.import(solution_host, dealii::VectorOperation::insert);
452  }
453 
454  // Unpack the material state and repopulate the material state
455  std::vector<std::vector<double>> transferred_data(
456  triangulation.n_active_cells(),
457  std::vector<double>(n_material_states + direction_data_size +
458  phase_history_data_size));
459  cell_data_trans.unpack(transferred_data);
460  auto state = material_properties.get_state();
461  state_host = Kokkos::create_mirror_view(state);
462  unsigned int total_cell_id = 0;
463  cell_id = 0;
464  std::vector<double> transferred_cos;
465  std::vector<double> transferred_sin;
466  std::vector<bool> has_melted;
467  for (auto const &cell : dof_handler.active_cell_iterators())
468  {
469  if (cell->is_locally_owned())
470  {
471  for (unsigned int i = 0; i < n_material_states; ++i)
472  {
473  state_host(i, cell_id) = transferred_data[total_cell_id][i];
474  }
475  if (cell->active_fe_index() == 0)
476  {
477  transferred_cos.push_back(
478  transferred_data[total_cell_id][n_material_states]);
479  transferred_sin.push_back(
480  transferred_data[total_cell_id][n_material_states + 1]);
481 
482  // Convert from double back to bool
483  if (transferred_data[total_cell_id]
484  [n_material_states + direction_data_size] > 0.5)
485  has_melted.push_back(true);
486  else
487  has_melted.push_back(false);
488  }
489  ++cell_id;
490  }
491  ++total_cell_id;
492  }
493 
494  // Update the deposition cos and sin
495  thermal_physics->set_material_deposition_orientation(transferred_cos,
496  transferred_sin);
497 
498  // Update the melted indicator
499  thermal_physics->set_has_melted_vector(has_melted);
500 
501  // Copy the data back to material_property
502  Kokkos::deep_copy(state, state_host);
503 
504  // Update the material states in the ThermalOperator
505  thermal_physics->get_state_from_material_properties();
506 
507 #if ADAMANTINE_DEBUG
508  // Check that we are not losing material
509  cell_id = 0;
510  for (auto const &cell : dof_handler.active_cell_iterators())
511  {
512  if (cell->is_locally_owned())
513  {
514  double material_ratio = 0.;
515  for (unsigned int i = 0; i < n_material_states; ++i)
516  {
517  material_ratio += state_host(i, cell_id);
518  }
519  ASSERT(std::abs(material_ratio - 1.) < 1e-14, "Material is lost.");
520  ++cell_id;
521  }
522  }
523 #endif
524 
525  if (mechanical_physics)
526  {
527  mechanical_physics->complete_transfer_mpi();
528  }
529 }
530 
531 template <int dim>
532 std::vector<typename dealii::parallel::distributed::Triangulation<
533  dim>::active_cell_iterator>
534 compute_cells_to_refine(
535  dealii::parallel::distributed::Triangulation<dim> &triangulation,
536  double const time, double const next_refinement_time,
537  unsigned int const n_time_steps,
538  std::vector<std::shared_ptr<adamantine::HeatSource<dim>>> const
539  &heat_sources)
540 {
541  // Compute the position of the beams between time and next_refinement_time and
542  // create a list of cells that will intersect the beam.
543 
544  // Build the bounding boxes associated with the locally owned cells
545  std::vector<dealii::BoundingBox<dim>> cell_bounding_boxes;
546  for (auto const &cell : triangulation.active_cell_iterators() |
547  dealii::IteratorFilters::LocallyOwnedCell())
548  {
549  cell_bounding_boxes.push_back(cell->bounding_box());
550  }
551  dealii::ArborXWrappers::BVH bvh(cell_bounding_boxes);
552 
553  double const bounding_box_scaling = 2.0;
554  std::vector<dealii::BoundingBox<dim>> heat_source_bounding_boxes;
555  for (unsigned int i = 0; i < n_time_steps; ++i)
556  {
557  double const current_time = time + static_cast<double>(i) /
558  static_cast<double>(n_time_steps) *
559  (next_refinement_time - time);
560 
561  // Build the bounding boxes associated with the heat sources
562  for (auto &beam : heat_sources)
563  {
564  heat_source_bounding_boxes.push_back(
565  beam->get_bounding_box(current_time, bounding_box_scaling));
566  }
567  }
568 
569  // Perform the search with ArborX. Since we are only interested in locally
570  // owned cells, we use BVH.
571  dealii::ArborXWrappers::BoundingBoxIntersectPredicate bb_intersect(
572  heat_source_bounding_boxes);
573  auto [indices, offset] = bvh.query(bb_intersect);
574 
575  // Put the indices into a set to get rid of the duplicates and to make it
576  // easier to check if the indices are found.
577  std::unordered_set<int> indices_to_refine;
578  indices_to_refine.insert(indices.begin(), indices.end());
579 
580  std::vector<typename dealii::parallel::distributed::Triangulation<
581  dim>::active_cell_iterator>
582  cells_to_refine;
583  int cell_index = 0;
584  for (auto const &cell : triangulation.active_cell_iterators() |
585  dealii::IteratorFilters::LocallyOwnedCell())
586  {
587  if (indices_to_refine.count(cell_index))
588  {
589  cells_to_refine.push_back(cell);
590  }
591  ++cell_index;
592  }
593 
594  return cells_to_refine;
595 }
596 
597 template <int dim, int n_materials, int p_order, int fe_degree,
598  typename MaterialStates, typename MemorySpaceType>
599 void refine_mesh(
600  std::unique_ptr<adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>
601  &thermal_physics,
602  std::unique_ptr<adamantine::MechanicalPhysics<
603  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>
604  &mechanical_physics,
605  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
606  MemorySpaceType> &material_properties,
607  dealii::LA::distributed::Vector<double, MemorySpaceType> &solution,
608  std::vector<std::shared_ptr<adamantine::HeatSource<dim>>> const
609  &heat_sources,
610  double const time, double const next_refinement_time,
611  unsigned int const time_steps_refinement,
612  boost::property_tree::ptree const &refinement_database)
613 {
614 #ifdef ADAMANTINE_WITH_CALIPER
615  CALI_CXX_MARK_FUNCTION;
616 #endif
617  dealii::DoFHandler<dim> &dof_handler = thermal_physics->get_dof_handler();
618  dealii::parallel::distributed::Triangulation<dim> &triangulation =
619  dynamic_cast<dealii::parallel::distributed::Triangulation<dim> &>(
620  const_cast<dealii::Triangulation<dim> &>(
621  dof_handler.get_triangulation()));
622 
623  // PropertyTreeInput refinement.n_refinements
624  unsigned int const n_refinements =
625  refinement_database.get("n_refinements", 2);
626 
627  for (unsigned int i = 0; i < n_refinements; ++i)
628  {
629  // Compute the cells to be refined.
630  auto cells_to_refine =
631  compute_cells_to_refine(triangulation, time, next_refinement_time,
632  time_steps_refinement, heat_sources);
633 
634  // PropertyTreeInput refinement.coarsen_after_beam
635  const bool coarsen_after_beam =
636  refinement_database.get<bool>("coarsen_after_beam", false);
637 
638  // If coarsening is allowed, set the coarsening flag everywhere
639  if (coarsen_after_beam)
640  {
641  for (auto cell : dealii::filter_iterators(
642  triangulation.active_cell_iterators(),
643  dealii::IteratorFilters::LocallyOwnedCell()))
644  {
645  if (cell->level() > 0)
646  cell->set_coarsen_flag();
647  }
648  }
649 
650  // Flag the cells for refinement.
651  for (auto &cell : cells_to_refine)
652  {
653  if (coarsen_after_beam)
654  cell->clear_coarsen_flag();
655 
656  if (cell->level() < static_cast<int>(n_refinements))
657  cell->set_refine_flag();
658  }
659 
660  // Execute the refinement and transfer the solution onto the new mesh.
661  refine_and_transfer(thermal_physics, mechanical_physics,
662  material_properties, dof_handler, solution);
663  }
664 
665  // Recompute the inverse of the mass matrix
666  thermal_physics->compute_inverse_mass_matrix();
667 }
668 
669 template <int dim, int n_materials, int p_order, typename MaterialStates,
670  typename MemorySpaceType>
671 void refine_mesh(
672  std::unique_ptr<adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>
673  &thermal_physics,
674  std::unique_ptr<adamantine::MechanicalPhysics<
675  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>
676  &mechanical_physics,
677  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
678  MemorySpaceType> &material_properties,
679  dealii::LA::distributed::Vector<double, MemorySpaceType> &solution,
680  std::vector<std::shared_ptr<adamantine::HeatSource<dim>>> const
681  &heat_sources,
682  double const time, double const next_refinement_time,
683  unsigned int const time_steps_refinement,
684  boost::property_tree::ptree const &refinement_database)
685 {
686  if (!thermal_physics)
687  return;
688 
689  switch (thermal_physics->get_fe_degree())
690  {
691  case 1:
692  {
693  refine_mesh<dim, n_materials, p_order, 1, MaterialStates>(
694  thermal_physics, mechanical_physics, material_properties, solution,
695  heat_sources, time, next_refinement_time, time_steps_refinement,
696  refinement_database);
697  break;
698  }
699  case 2:
700  {
701  refine_mesh<dim, n_materials, p_order, 2, MaterialStates>(
702  thermal_physics, mechanical_physics, material_properties, solution,
703  heat_sources, time, next_refinement_time, time_steps_refinement,
704  refinement_database);
705  break;
706  }
707  case 3:
708  {
709  refine_mesh<dim, n_materials, p_order, 3, MaterialStates>(
710  thermal_physics, mechanical_physics, material_properties, solution,
711  heat_sources, time, next_refinement_time, time_steps_refinement,
712  refinement_database);
713  break;
714  }
715  case 4:
716  {
717  refine_mesh<dim, n_materials, p_order, 4, MaterialStates>(
718  thermal_physics, mechanical_physics, material_properties, solution,
719  heat_sources, time, next_refinement_time, time_steps_refinement,
720  refinement_database);
721  break;
722  }
723  case 5:
724  {
725  refine_mesh<dim, n_materials, p_order, 5, MaterialStates>(
726  thermal_physics, mechanical_physics, material_properties, solution,
727  heat_sources, time, next_refinement_time, time_steps_refinement,
728  refinement_database);
729  break;
730  }
731  default:
732  {
733  adamantine::ASSERT_THROW(false, "fe_degree should be between 1 and 5.");
734  }
735  }
736 }
737 
738 template <int dim, int n_materials, int p_order, typename MaterialStates,
739  typename MemorySpaceType>
740 std::pair<dealii::LinearAlgebra::distributed::Vector<double,
741  dealii::MemorySpace::Host>,
742  dealii::LinearAlgebra::distributed::Vector<double,
743  dealii::MemorySpace::Host>>
744 run(MPI_Comm const &communicator, boost::property_tree::ptree const &database,
745  std::vector<adamantine::Timer> &timers)
746 {
747 #ifdef ADAMANTINE_WITH_CALIPER
748  CALI_CXX_MARK_FUNCTION;
749 #endif
750  // Get optional units property tree
751  boost::optional<boost::property_tree::ptree const &> units_optional_database =
752  database.get_child_optional("units");
753 
754  // Create the Geometry
755  boost::property_tree::ptree geometry_database =
756  database.get_child("geometry");
757  adamantine::Geometry<dim> geometry(communicator, geometry_database,
758  units_optional_database);
759 
760  // Create the MaterialProperty
761  boost::property_tree::ptree material_database =
762  database.get_child("materials");
763  adamantine::MaterialProperty<dim, n_materials, p_order, MaterialStates,
764  MemorySpaceType>
765  material_properties(communicator, geometry.get_triangulation(),
766  material_database);
767 
768  // Extract the physics property tree
769  boost::property_tree::ptree physics_database = database.get_child("physics");
770  bool const use_thermal_physics = physics_database.get<bool>("thermal");
771  bool const use_mechanical_physics = physics_database.get<bool>("mechanical");
772 
773  // Extract the discretization property tree
774  boost::property_tree::ptree discretization_database =
775  database.get_child("discretization");
776 
777  // Extract the post-processor property tree
778  boost::property_tree::ptree post_processor_database =
779  database.get_child("post_processor");
780 
781  // Extract the verbosity
782  // PropertyTreeInput verbose_output
783  bool const verbose_output = database.get("verbose_output", false);
784 
785  // Create the Boundary
786  boost::property_tree::ptree boundary_database =
787  database.get_child("boundary");
788  adamantine::Boundary boundary(
789  boundary_database, geometry.get_triangulation().get_boundary_ids());
790 
791  // Create ThermalPhysics if necessary
792  std::unique_ptr<adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>
793  thermal_physics;
794  std::vector<std::shared_ptr<adamantine::HeatSource<dim>>> heat_sources;
795  if (use_thermal_physics)
796  {
797  // PropertyTreeInput discretization.thermal.fe_degree
798  unsigned int const fe_degree =
799  discretization_database.get<unsigned int>("thermal.fe_degree");
800  // PropertyTreeInput discretization.thermal.quadrature
801  std::string quadrature_type =
802  discretization_database.get("thermal.quadrature", "gauss");
803  thermal_physics = initialize_thermal_physics<dim>(
804  fe_degree, quadrature_type, communicator, database, geometry, boundary,
805  material_properties);
806  heat_sources = thermal_physics->get_heat_sources();
807  post_processor_database.put("thermal_output", true);
808  }
809 
810  // PropertyTreeInput materials.initial_temperature
811  double const initial_temperature =
812  material_database.get("initial_temperature", 300.);
813 
814  // Get the reference temperature(s) for the substrate and deposited
815  // material(s)
816  std::vector<double> material_reference_temps;
817  // PropertyTreeInput materials.n_materials
818  unsigned int const n_materials_runtime =
819  material_database.get<unsigned int>("n_materials");
820  for (unsigned int i = 0; i < n_materials_runtime; ++i)
821  {
822  // Use the solidus as the reference temperature, in the future we may want
823  // to use something else
824  // PropertyTreeInput materials.material_n.solidus
825  double const reference_temperature = material_database.get<double>(
826  "material_" + std::to_string(i) + ".solidus");
827  material_reference_temps.push_back(reference_temperature);
828  }
829  material_reference_temps.push_back(initial_temperature);
830 
831  // Create MechanicalPhysics
832  std::unique_ptr<adamantine::MechanicalPhysics<
833  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>
834  mechanical_physics;
835  if (use_mechanical_physics)
836  {
837  // PropertyTreeInput discretization.mechanical.fe_degree
838  unsigned int const fe_degree =
839  discretization_database.get<unsigned int>("mechanical.fe_degree");
840  mechanical_physics = std::make_unique<adamantine::MechanicalPhysics<
841  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>(
842  communicator, fe_degree, geometry, boundary, material_properties,
843  material_reference_temps);
844  post_processor_database.put("mechanical_output", true);
845  }
846 
847  // Create PostProcessor
848  std::unique_ptr<adamantine::PostProcessor<dim>> post_processor;
849  bool thermal_output = post_processor_database.get("thermal_output", false);
850  bool mechanical_output =
851  post_processor_database.get("mechanical_output", false);
852  if ((thermal_output) && (mechanical_output))
853  {
854  post_processor = std::make_unique<adamantine::PostProcessor<dim>>(
855  communicator, post_processor_database,
856  thermal_physics->get_dof_handler(),
857  mechanical_physics->get_dof_handler());
858  }
859  else
860  {
861  post_processor = std::make_unique<adamantine::PostProcessor<dim>>(
862  communicator, post_processor_database,
863  thermal_output ? thermal_physics->get_dof_handler()
864  : mechanical_physics->get_dof_handler());
865  }
866 
867  // Get the checkpoint and restart subtrees
868  boost::optional<boost::property_tree::ptree const &>
869  checkpoint_optional_database = database.get_child_optional("checkpoint");
870  boost::optional<boost::property_tree::ptree const &>
871  restart_optional_database = database.get_child_optional("restart");
872 
873  unsigned int time_steps_checkpoint = std::numeric_limits<unsigned int>::max();
874  std::string checkpoint_filename;
875  bool checkpoint_overwrite = true;
876  if (checkpoint_optional_database)
877  {
878  auto checkpoint_database = checkpoint_optional_database.get();
879  // PropertyTreeInput checkpoint.time_steps_between_checkpoint
880  time_steps_checkpoint =
881  checkpoint_database.get<unsigned int>("time_steps_between_checkpoint");
882  // PropertyTreeInput checkpoint.filename_prefix
883  checkpoint_filename =
884  checkpoint_database.get<std::string>("filename_prefix");
885  // PropertyTreeInput checkpoint.overwrite_files
886  checkpoint_overwrite = checkpoint_database.get<bool>("overwrite_files");
887  }
888 
889  bool restart = false;
890  std::string restart_filename;
891  if (restart_optional_database)
892  {
893  restart = true;
894  auto restart_database = restart_optional_database.get();
895  // PropertyTreeInput restart.filename_prefix
896  restart_filename = restart_database.get<std::string>("filename_prefix");
897  }
898 
899  dealii::LA::distributed::Vector<double, MemorySpaceType> temperature;
900  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
901  displacement;
902  if (use_thermal_physics)
903  {
904  if (restart == false)
905  {
906  thermal_physics->setup();
907  thermal_physics->initialize_dof_vector(initial_temperature, temperature);
908  }
909  else
910  {
911 #ifdef ADAMANTINE_WITH_CALIPER
912  CALI_MARK_BEGIN("restart from file");
913 #endif
914  thermal_physics->load_checkpoint(restart_filename, temperature);
915 #ifdef ADAMANTINE_WITH_CALIPER
916  CALI_MARK_END("restart from file");
917 #endif
918  }
919  }
920 
921  if (use_mechanical_physics)
922  {
923  // We currently do not support restarting the mechanical simulation.
924  adamantine::ASSERT_THROW(
925  restart == false,
926  "Mechanical simulation cannot be restarted from a file");
927  if (use_thermal_physics)
928  {
929  // Thermo-mechanical simulation
930  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
931  temperature_host(temperature.get_partitioner());
932  temperature_host.import(temperature, dealii::VectorOperation::insert);
933  mechanical_physics->setup_dofs(thermal_physics->get_dof_handler(),
934  temperature_host,
935  thermal_physics->get_has_melted_vector());
936  }
937  else
938  {
939  // Mechanical only simulation
940  mechanical_physics->setup_dofs();
941  }
942  displacement = mechanical_physics->solve();
943  }
944 
945  unsigned int progress = 0;
946  unsigned int n_time_step = 0;
947  double time = 0.;
948  double activation_time_end = -1.;
949  unsigned int const rank =
950  dealii::Utilities::MPI::this_mpi_process(communicator);
951  if (restart == true)
952  {
953  if (rank == 0)
954  {
955  std::cout << "Restarting from file" << std::endl;
956  }
957  std::ifstream file{restart_filename + "_time.txt"};
958  boost::archive::text_iarchive ia{file};
959  ia >> time;
960  ia >> n_time_step;
961  }
962  // PropertyTreeInput geometry.deposition_time
963  double const activation_time =
964  geometry_database.get<double>("deposition_time", 0.);
965 
966  // Output the initial solution
967  output_pvtu(*post_processor, n_time_step, time, thermal_physics, temperature,
968  mechanical_physics, displacement, material_properties, timers);
969  ++n_time_step;
970 
971  // Create the bounding boxes used for material deposition
972  auto [material_deposition_boxes, deposition_times, deposition_cos,
973  deposition_sin] =
974  adamantine::create_material_deposition_boxes<dim>(geometry_database,
975  heat_sources);
976  // Extract the time-stepping database
977  boost::property_tree::ptree time_stepping_database =
978  database.get_child("time_stepping");
979  // PropertyTreeInput time_stepping.time_step
980  double time_step = time_stepping_database.get<double>("time_step");
981  // PropertyTreeInput time_stepping.scan_path_for_duration
982  bool const scan_path_for_duration =
983  time_stepping_database.get("scan_path_for_duration", false);
984  // PropertyTreeInput time_stepping.duration
985  double const duration = scan_path_for_duration
986  ? std::numeric_limits<double>::max()
987  : time_stepping_database.get<double>("duration");
988 
989  // Extract the refinement database
990  boost::property_tree::ptree refinement_database =
991  database.get_child("refinement");
992  // PropertyTreeInput refinement.time_steps_between_refinement
993  unsigned int const time_steps_refinement =
994  refinement_database.get("time_steps_between_refinement", 10);
995  // PropertyTreeInput post_processor.time_steps_between_output
996  unsigned int const time_steps_output =
997  post_processor_database.get("time_steps_between_output", 1);
998  // PropertyTreeInput materials.new_material_temperature
999  double const new_material_temperature =
1000  database.get("materials.new_material_temperature", 300.);
1001 
1002  // Extract the microstructure database if present
1003  boost::optional<boost::property_tree::ptree const &>
1004  microstructure_optional_database =
1005  database.get_child_optional("microstructure");
1006  bool const compute_microstructure =
1007  microstructure_optional_database ? true : false;
1008  std::unique_ptr<adamantine::Microstructure<dim>> microstructure;
1009  if (compute_microstructure)
1010  {
1011  auto microstructure_database = microstructure_optional_database.get();
1012  // PropertyTreeInput microstructure.filename_prefix
1013  std::string microstructure_filename =
1014  microstructure_database.get<std::string>("filename_prefix");
1015  microstructure = std::make_unique<adamantine::Microstructure<dim>>(
1016  communicator, microstructure_filename);
1017  }
1018 
1019 #ifdef ADAMANTINE_WITH_CALIPER
1020  CALI_CXX_MARK_LOOP_BEGIN(main_loop_id, "main_loop");
1021 #endif
1022  while (time < duration)
1023  {
1024 #ifdef ADAMANTINE_WITH_CALIPER
1025  CALI_CXX_MARK_LOOP_ITERATION(main_loop_id, n_time_step - 1);
1026 #endif
1027  if ((time + time_step) > duration)
1028  time_step = duration - time;
1029 
1030  // Refine the mesh the first time we get in the loop and after
1031  // time_steps_refinement time steps.
1032  if (((n_time_step == 1) || ((n_time_step % time_steps_refinement) == 0)) &&
1033  use_thermal_physics)
1034  {
1035  timers[adamantine::refine].start();
1036  double next_refinement_time = time + time_steps_refinement * time_step;
1037  refine_mesh(thermal_physics, mechanical_physics, material_properties,
1038  temperature, heat_sources, time, next_refinement_time,
1039  time_steps_refinement, refinement_database);
1040  timers[adamantine::refine].stop();
1041  if ((rank == 0) && (verbose_output == true))
1042  {
1043  std::cout << "n_time_step: " << n_time_step << " time: " << time
1044  << " n_dofs after mesh refinement: "
1045  << thermal_physics->get_dof_handler().n_dofs() << std::endl;
1046  }
1047  }
1048 
1049  // Add material if necessary.
1050  // We use an epsilon to get the "expected" behavior when the deposition
1051  // time and the time match should match exactly but don't because of
1052  // floating point accuracy.
1053  timers[adamantine::add_material_activate].start();
1054  if (time > activation_time_end)
1055  {
1056  // If we use scan_path_for_duration, we may need to read the scan path
1057  // file once again.
1058  if (scan_path_for_duration)
1059  {
1060  // Check if we have reached the end of current scan path.
1061  bool need_updated_scan_path = false;
1062  for (auto &source : heat_sources)
1063  {
1064  if (time > source->get_scan_path().get_segment_list().back().end_time)
1065  {
1066  need_updated_scan_path = true;
1067  break;
1068  }
1069  }
1070 
1071  if (need_updated_scan_path)
1072  {
1073  // Check if we have reached the end of the file. If not, read the
1074  // updated scan path file
1075  bool scan_path_end = true;
1076  for (auto &source : heat_sources)
1077  {
1078  if (!source->get_scan_path().is_finished())
1079  {
1080  scan_path_end = false;
1081  // This functions waits for the scan path file to be updated
1082  // before reading the file.
1083  source->get_scan_path().read_file();
1084  }
1085  }
1086 
1087  // If we have reached the end of scan path file for all the heat
1088  // sources, we just exit.
1089  if (scan_path_end)
1090  {
1091  break;
1092  }
1093 
1094  std::tie(material_deposition_boxes, deposition_times, deposition_cos,
1095  deposition_sin) =
1096  adamantine::create_material_deposition_boxes<dim>(
1097  geometry_database, heat_sources);
1098  }
1099  }
1100 
1101  double const eps = time_step / 1e10;
1102 
1103  auto activation_start =
1104  std::lower_bound(deposition_times.begin(), deposition_times.end(),
1105  time - eps) -
1106  deposition_times.begin();
1107  activation_time_end =
1108  std::min(time + std::max(activation_time, time_step), duration) - eps;
1109  auto activation_end =
1110  std::lower_bound(deposition_times.begin(), deposition_times.end(),
1111  activation_time_end) -
1112  deposition_times.begin();
1113  if (activation_start < activation_end)
1114  {
1115  if (use_thermal_physics)
1116  {
1117 #ifdef ADAMANTINE_WITH_CALIPER
1118  CALI_MARK_BEGIN("add material");
1119 #endif
1120  // Compute the elements to activate.
1121  // TODO Right now, we compute the list of cells that get activated for
1122  // the entire material deposition. We should restrict the list to the
1123  // cells that are activated between activation_start and
1124  // activation_end.
1125  timers[adamantine::add_material_search].start();
1126  auto elements_to_activate = adamantine::get_elements_to_activate(
1127  thermal_physics->get_dof_handler(), material_deposition_boxes);
1128  timers[adamantine::add_material_search].stop();
1129 
1130  // For now assume that all deposited material has never been melted
1131  // (may or may not be reasonable)
1132  std::vector<bool> has_melted(deposition_cos.size(), false);
1133 
1134  thermal_physics->add_material_start(
1135  elements_to_activate, deposition_cos, deposition_sin, has_melted,
1136  activation_start, activation_end, temperature);
1137 
1138  if (use_mechanical_physics)
1139  {
1140  mechanical_physics->prepare_transfer_mpi();
1141  }
1142 
1143 #ifdef ADAMANTINE_WITH_CALIPER
1144  CALI_MARK_BEGIN("refine triangulation");
1145 #endif
1146  dealii::parallel::distributed::Triangulation<dim> &triangulation =
1147  dynamic_cast<dealii::parallel::distributed::Triangulation<dim> &>(
1148  const_cast<dealii::Triangulation<dim> &>(
1149  thermal_physics->get_dof_handler().get_triangulation()));
1150  triangulation.execute_coarsening_and_refinement();
1151 #ifdef ADAMANTINE_WITH_CALIPER
1152  CALI_MARK_END("refine triangulation");
1153 #endif
1154 
1155  thermal_physics->add_material_end(new_material_temperature,
1156  temperature);
1157 
1158  if (use_mechanical_physics)
1159  {
1160  mechanical_physics->complete_transfer_mpi();
1161  }
1162 
1163 #ifdef ADAMANTINE_WITH_CALIPER
1164  CALI_MARK_END("add material");
1165 #endif
1166  }
1167  }
1168 
1169  if ((rank == 0) && (verbose_output == true) &&
1170  (activation_end - activation_start > 0) && use_thermal_physics)
1171  {
1172  std::cout << "n_time_step: " << n_time_step << " time: " << time
1173  << " n_dofs after cell activation: "
1174  << thermal_physics->get_dof_handler().n_dofs() << std::endl;
1175  }
1176  }
1177  timers[adamantine::add_material_activate].stop();
1178 
1179  // If thermomechanics are being solved, mark cells that are above the
1180  // solidus as cells that should have their reference temperature reset.
1181  // This cannot be in the thermomechanics solve because some cells may go
1182  // above the solidus and then back below the solidus in the time between
1183  // thermomechanical solves.
1184  if (use_thermal_physics && use_mechanical_physics)
1185  {
1186  thermal_physics->mark_has_melted(material_reference_temps[0],
1187  temperature);
1188  }
1189 
1190  timers[adamantine::evol_time].start();
1191 
1192  // Solve the thermal problem
1193  if (use_thermal_physics)
1194  {
1195  if (compute_microstructure)
1196  {
1197 #ifdef ADAMANTINE_WITH_CALIPER
1198  CALI_MARK_BEGIN("compute microstructure");
1199 #endif
1200  if constexpr (std::is_same_v<MemorySpaceType,
1201  dealii::MemorySpace::Host>)
1202  {
1203  microstructure->set_old_temperature(temperature);
1204  }
1205  else
1206  {
1207  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
1208  temperature_host(temperature.get_partitioner());
1209  temperature_host.import_elements(temperature,
1210  dealii::VectorOperation::insert);
1211  microstructure->set_old_temperature(temperature_host);
1212  }
1213 #ifdef ADAMANTINE_WITH_CALIPER
1214  CALI_MARK_END("compute microstructure");
1215 #endif
1216  }
1217 
1218  time = thermal_physics->evolve_one_time_step(time, time_step, temperature,
1219  timers);
1220 
1221  if (compute_microstructure)
1222  {
1223 #ifdef ADAMANTINE_WITH_CALIPER
1224  CALI_MARK_BEGIN("compute microstructure");
1225 #endif
1226  if constexpr (std::is_same_v<MemorySpaceType,
1227  dealii::MemorySpace::Host>)
1228  {
1229  microstructure->compute_G_and_R(material_properties,
1230  thermal_physics->get_dof_handler(),
1231  temperature, time_step);
1232  }
1233  else
1234  {
1235  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
1236  temperature_host(temperature.get_partitioner());
1237  temperature_host.import_elements(temperature,
1238  dealii::VectorOperation::insert);
1239  microstructure->compute_G_and_R(material_properties,
1240  thermal_physics->get_dof_handler(),
1241  temperature_host, time_step);
1242  }
1243 #ifdef ADAMANTINE_WITH_CALIPER
1244  CALI_MARK_END("compute microstructure");
1245 #endif
1246  }
1247  }
1248 
1249  // Solve the (thermo-)mechanical problem
1250  if (use_mechanical_physics)
1251  {
1252  // Since there is no history dependence in the model, only calculate
1253  // mechanics when outputting
1254  if (n_time_step % time_steps_output == 0)
1255  {
1256  if (use_thermal_physics)
1257  {
1258  // Update the material state
1259  thermal_physics->set_state_to_material_properties();
1260  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
1261  temperature_host(temperature.get_partitioner());
1262  temperature_host.import(temperature, dealii::VectorOperation::insert);
1263  mechanical_physics->setup_dofs(
1264  thermal_physics->get_dof_handler(), temperature_host,
1265  thermal_physics->get_has_melted_vector());
1266  }
1267  else
1268  {
1269  mechanical_physics->setup_dofs();
1270  }
1271  displacement = mechanical_physics->solve();
1272  }
1273  }
1274 
1275  timers[adamantine::evol_time].stop();
1276 
1277  if (n_time_step % time_steps_checkpoint == 0)
1278  {
1279 #ifdef ADAMANTINE_WITH_CALIPER
1280  CALI_MARK_BEGIN("save checkpoint");
1281 #endif
1282  if (rank == 0)
1283  {
1284  std::cout << "Checkpoint reached" << std::endl;
1285  }
1286 
1287  std::string output_dir =
1288  post_processor_database.get<std::string>("output_dir", "");
1289  std::string filename_prefix =
1290  checkpoint_overwrite
1291  ? checkpoint_filename
1292  : checkpoint_filename + '_' + std::to_string(n_time_step);
1293  thermal_physics->save_checkpoint(filename_prefix, temperature);
1294  std::ofstream file{output_dir + filename_prefix + "_time.txt"};
1295  boost::archive::text_oarchive oa{file};
1296  oa << time;
1297  oa << n_time_step;
1298 #ifdef ADAMANTINE_WITH_CALIPER
1299  CALI_MARK_END("save checkpoint");
1300 #endif
1301  }
1302 
1303  // Output progress on screen
1304  if (rank == 0)
1305  {
1306  double adim_time = time / (duration / 10.);
1307  double int_part = 0;
1308  std::modf(adim_time, &int_part);
1309  if (int_part > progress)
1310  {
1311  std::cout << int_part * 10 << '%' << " completed" << std::endl;
1312  progress = static_cast<unsigned int>(int_part);
1313  }
1314  }
1315 
1316  // Output the solution
1317  if (n_time_step % time_steps_output == 0)
1318  {
1319  if (use_thermal_physics)
1320  {
1321  thermal_physics->set_state_to_material_properties();
1322  }
1323  output_pvtu(*post_processor, n_time_step, time, thermal_physics,
1324  temperature, mechanical_physics, displacement,
1325  material_properties, timers);
1326  }
1327  ++n_time_step;
1328  }
1329 
1330 #ifdef ADAMANTINE_WITH_CALIPER
1331  CALI_CXX_MARK_LOOP_END(main_loop_id);
1332 #endif
1333 
1334  post_processor->write_pvd();
1335 
1336  // This is only used for integration test
1337  if constexpr (std::is_same_v<MemorySpaceType, dealii::MemorySpace::Host>)
1338  {
1339  if (use_thermal_physics)
1340  {
1341  thermal_physics->get_affine_constraints().distribute(temperature);
1342  }
1343  if (use_mechanical_physics)
1344  {
1345  mechanical_physics->get_affine_constraints().distribute(displacement);
1346  }
1347  return std::make_pair(temperature, displacement);
1348  }
1349  else
1350  {
1351  dealii::LinearAlgebra::distributed::Vector<double,
1352  dealii::MemorySpace::Host>
1353  temperature_host(temperature.get_partitioner());
1354  temperature_host.import(temperature, dealii::VectorOperation::insert);
1355  if (use_thermal_physics)
1356  {
1357  thermal_physics->get_affine_constraints().distribute(temperature_host);
1358  }
1359  if (use_mechanical_physics)
1360  {
1361  mechanical_physics->get_affine_constraints().distribute(displacement);
1362  }
1363  return std::make_pair(temperature_host, displacement);
1364  }
1365 }
1366 
1367 void split_global_communicator(MPI_Comm global_communicator,
1368  unsigned int global_ensemble_size,
1369  MPI_Comm &local_communicator,
1370  unsigned int &local_ensemble_size,
1371  unsigned int &first_local_member, int &my_color)
1372 {
1373  unsigned int const global_rank =
1374  dealii::Utilities::MPI::this_mpi_process(global_communicator);
1375  unsigned int const global_n_procs =
1376  dealii::Utilities::MPI::n_mpi_processes(global_communicator);
1377 
1378  double const avg_n_procs_per_member =
1379  static_cast<double>(global_n_procs) /
1380  static_cast<double>(global_ensemble_size);
1381  if (avg_n_procs_per_member > 1)
1382  {
1383  local_ensemble_size = 1;
1384  // We need all the members to use the same partitioning otherwise the dofs
1385  // may be number differently which is problematic for the data assimulation.
1386  adamantine::ASSERT_THROW(
1387  avg_n_procs_per_member == std::floor(avg_n_procs_per_member),
1388  "Number of MPI ranks should be less than the number of ensemble "
1389  "members or a multiple of this number.");
1390  // Assign color
1391  my_color =
1392  global_rank / static_cast<int>(std::floor(avg_n_procs_per_member));
1393 
1394  first_local_member = my_color;
1395  }
1396  else
1397  {
1398  my_color = global_rank;
1399  // Compute the local ensemble size
1400  unsigned int const members_per_proc = global_ensemble_size / global_n_procs;
1401  unsigned int const leftover_members =
1402  global_ensemble_size - global_n_procs * members_per_proc;
1403  // We assign all the leftover members to processor 0. We could do a
1404  // better job distributing the leftover
1405  local_ensemble_size = members_per_proc;
1406  if (global_rank == 0)
1407  {
1408  local_ensemble_size += leftover_members;
1409  first_local_member = 0;
1410  }
1411  else
1412  {
1413  first_local_member = global_rank * local_ensemble_size + leftover_members;
1414  }
1415  }
1416 
1417  MPI_Comm_split(global_communicator, my_color, 0, &local_communicator);
1418 }
1419 
1420 template <int dim, int n_materials, int p_order, typename MaterialStates,
1421  typename MemorySpaceType>
1422 std::vector<dealii::LA::distributed::BlockVector<double>>
1423 run_ensemble(MPI_Comm const &global_communicator,
1424  boost::property_tree::ptree const &database,
1425  std::vector<adamantine::Timer> &timers)
1426 {
1427 #ifdef ADAMANTINE_WITH_CALIPER
1428  CALI_CXX_MARK_FUNCTION;
1429 #endif
1430 
1431  // ------ Extract child property trees -----
1432  // Mandatory subtrees
1433  boost::property_tree::ptree geometry_database =
1434  database.get_child("geometry");
1435  boost::property_tree::ptree boundary_database =
1436  database.get_child("boundary");
1437  boost::property_tree::ptree discretization_database =
1438  database.get_child("discretization");
1439  boost::property_tree::ptree time_stepping_database =
1440  database.get_child("time_stepping");
1441  boost::property_tree::ptree post_processor_database =
1442  database.get_child("post_processor");
1443  boost::property_tree::ptree refinement_database =
1444  database.get_child("refinement");
1445  boost::property_tree::ptree material_database =
1446  database.get_child("materials");
1447 
1448  // Optional subtrees
1449  boost::optional<boost::property_tree::ptree const &> units_optional_database =
1450  database.get_child_optional("units");
1451  boost::optional<const boost::property_tree::ptree &>
1452  experiment_optional_database = database.get_child_optional("experiment");
1453  boost::optional<const boost::property_tree::ptree &>
1454  data_assimilation_optional_database =
1455  database.get_child_optional("data_assimilation");
1456 
1457  // Verbosity
1458  // PropertyTreeInput verbose_output
1459  bool const verbose_output = database.get("verbose_output", false);
1460 
1461  // ------ Get finite element implementation parameters -----
1462  // PropertyTreeInput discretization.fe_degree
1463  unsigned int const fe_degree =
1464  discretization_database.get<unsigned int>("thermal.fe_degree");
1465  // PropertyTreeInput discretization.quadrature
1466  std::string quadrature_type =
1467  discretization_database.get("thermal.quadrature", "gauss");
1468  std::transform(quadrature_type.begin(), quadrature_type.end(),
1469  quadrature_type.begin(),
1470  [](unsigned char c) { return std::tolower(c); });
1471 
1472  // ------ Split MPI communicator -----
1473  // PropertyTreeInput ensemble.ensemble_size
1474  unsigned int const global_ensemble_size =
1475  database.get("ensemble.ensemble_size", 5);
1476  // Distribute the processors among the ensemble members
1477  MPI_Comm local_communicator;
1478  unsigned int local_ensemble_size = std::numeric_limits<unsigned int>::max();
1479  unsigned int first_local_member = std::numeric_limits<unsigned int>::max();
1480  int my_color = -1;
1481  split_global_communicator(global_communicator, global_ensemble_size,
1482  local_communicator, local_ensemble_size,
1483  first_local_member, my_color);
1484 
1485  // ------ Set up the ensemble members -----
1486  // There might be a more efficient way to share some of these objects
1487  // between ensemble members. For now, we'll do the simpler approach of
1488  // duplicating everything.
1489 
1490  // Create a new property tree database for each ensemble member
1491  std::vector<boost::property_tree::ptree> database_ensemble =
1492  adamantine::create_database_ensemble(
1493  database, local_communicator, first_local_member, local_ensemble_size,
1494  global_ensemble_size);
1495 
1496  std::vector<std::unique_ptr<
1497  adamantine::ThermalPhysicsInterface<dim, MemorySpaceType>>>
1498  thermal_physics_ensemble(local_ensemble_size);
1499 
1500  std::vector<std::vector<std::shared_ptr<adamantine::HeatSource<dim>>>>
1501  heat_sources_ensemble(local_ensemble_size);
1502 
1503  std::vector<std::unique_ptr<adamantine::Geometry<dim>>> geometry_ensemble;
1504 
1505  std::vector<std::unique_ptr<adamantine::Boundary>> boundary_ensemble;
1506 
1507  std::vector<std::unique_ptr<adamantine::MaterialProperty<
1508  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>>
1509  material_properties_ensemble;
1510 
1511  std::vector<std::unique_ptr<adamantine::PostProcessor<dim>>>
1512  post_processor_ensemble;
1513 
1514  // Create the vector of augmented state vectors
1515  std::vector<dealii::LA::distributed::BlockVector<double>>
1516  solution_augmented_ensemble(local_ensemble_size);
1517 
1518  // Give names to the blocks in the augmented state vector
1519  int constexpr base_state = 0;
1520  int constexpr augmented_state = 1;
1521 
1522  // ----- Initialize the data assimilation object -----
1523  boost::property_tree::ptree data_assimilation_database;
1524  bool assimilate_data = false;
1525  std::vector<adamantine::AugmentedStateParameters> augmented_state_parameters;
1526 
1527  if (data_assimilation_optional_database)
1528  {
1529  // PropertyTreeInput data_assimilation
1530  data_assimilation_database = data_assimilation_optional_database.get();
1531 
1532  // PropertyTreeInput data_assimilation.assimilate_data
1533  if (data_assimilation_database.get("assimilate_data", false))
1534  {
1535  assimilate_data = true;
1536 
1537  // PropertyTreeInput data_assimilation.augment_with_beam_0_absorption
1538  if (data_assimilation_database.get("augment_with_beam_0_absorption",
1539  false))
1540  {
1541  augmented_state_parameters.push_back(
1542  adamantine::AugmentedStateParameters::beam_0_absorption);
1543  }
1544  // PropertyTreeInput data_assimilation.augment_with_beam_0_max_power
1545  if (data_assimilation_database.get("augment_with_beam_0_max_power",
1546  false))
1547  {
1548  augmented_state_parameters.push_back(
1549  adamantine::AugmentedStateParameters::beam_0_max_power);
1550  }
1551  }
1552  }
1553  adamantine::DataAssimilator data_assimilator(global_communicator,
1554  local_communicator, my_color,
1555  data_assimilation_database);
1556 
1557  // Get the checkpoint and restart subtrees
1558  boost::optional<boost::property_tree::ptree const &>
1559  checkpoint_optional_database = database.get_child_optional("checkpoint");
1560  boost::optional<boost::property_tree::ptree const &>
1561  restart_optional_database = database.get_child_optional("restart");
1562 
1563  unsigned int const global_rank =
1564  dealii::Utilities::MPI::this_mpi_process(global_communicator);
1565  unsigned int time_steps_checkpoint = std::numeric_limits<unsigned int>::max();
1566  std::string checkpoint_filename;
1567  bool checkpoint_overwrite = true;
1568  if (checkpoint_optional_database)
1569  {
1570  auto checkpoint_database = checkpoint_optional_database.get();
1571  // PropertyTreeInput checkpoint.time_steps_between_checkpoint
1572  time_steps_checkpoint =
1573  checkpoint_database.get<unsigned int>("time_steps_between_checkpoint");
1574  // PropertyTreeInput checkpoint.filename_prefix
1575  checkpoint_filename =
1576  checkpoint_database.get<std::string>("filename_prefix");
1577  // PropertyTreeInput checkpoint.overwrite_files
1578  checkpoint_overwrite = checkpoint_database.get<bool>("overwrite_files");
1579  }
1580 
1581  bool restart = false;
1582  std::string restart_filename;
1583  if (restart_optional_database)
1584  {
1585  restart = true;
1586  auto restart_database = restart_optional_database.get();
1587  // PropertyTreeInput restart.filename_prefix
1588  restart_filename = restart_database.get<std::string>("filename_prefix");
1589  }
1590  for (unsigned int member = 0; member < local_ensemble_size; ++member)
1591  {
1592  // Resize the augmented ensemble block vector to have two blocks
1593  solution_augmented_ensemble[member].reinit(2);
1594 
1595  // Edit the database for the ensemble
1596  if (database.get<unsigned int>("sources.n_beams") > 0)
1597  {
1598  // Populate the parameter augmentation block of the augmented state
1599  // ensemble
1600  if (assimilate_data)
1601  {
1602  solution_augmented_ensemble[member]
1603  .block(augmented_state)
1604  .reinit(augmented_state_parameters.size());
1605 
1606  for (unsigned int index = 0; index < augmented_state_parameters.size();
1607  ++index)
1608  {
1609  // FIXME: Need to consider how we want to generalize this. It could
1610  // get unwieldy if we want to specify every parameter of an
1611  // arbitrary number of beams.
1612  if (augmented_state_parameters.at(index) ==
1613  adamantine::AugmentedStateParameters::beam_0_absorption)
1614  {
1615  solution_augmented_ensemble[member].block(augmented_state)[index] =
1616  database_ensemble[member].get<double>(
1617  "sources.beam_0.absorption_efficiency");
1618  }
1619  else if (augmented_state_parameters.at(index) ==
1620  adamantine::AugmentedStateParameters::beam_0_max_power)
1621  {
1622  solution_augmented_ensemble[member].block(augmented_state)[index] =
1623  database_ensemble[member].get<double>(
1624  "sources.beam_0.max_power");
1625  }
1626  }
1627  }
1628  }
1629 
1630  solution_augmented_ensemble[member].collect_sizes();
1631 
1632  geometry_ensemble.push_back(std::make_unique<adamantine::Geometry<dim>>(
1633  local_communicator, geometry_database, units_optional_database));
1634 
1635  boundary_ensemble.push_back(std::make_unique<adamantine::Boundary>(
1636  boundary_database,
1637  geometry_ensemble.back()->get_triangulation().get_boundary_ids()));
1638 
1639  material_properties_ensemble.push_back(
1640  std::make_unique<adamantine::MaterialProperty<
1641  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>(
1642  local_communicator, geometry_ensemble.back()->get_triangulation(),
1643  material_database));
1644 
1645  thermal_physics_ensemble[member] = initialize_thermal_physics<dim>(
1646  fe_degree, quadrature_type, local_communicator,
1647  database_ensemble[member], *geometry_ensemble[member],
1648  *boundary_ensemble[member], *material_properties_ensemble[member]);
1649  heat_sources_ensemble[member] =
1650  thermal_physics_ensemble[member]->get_heat_sources();
1651 
1652  if (restart == false)
1653  {
1654  thermal_physics_ensemble[member]->setup();
1655  thermal_physics_ensemble[member]->initialize_dof_vector(
1656  database_ensemble[member].get("materials.initial_temperature", 300.),
1657  solution_augmented_ensemble[member].block(base_state));
1658  }
1659  else
1660  {
1661 #ifdef ADAMANTINE_WITH_CALIPER
1662  CALI_MARK_BEGIN("restart from file");
1663 #endif
1664  thermal_physics_ensemble[member]->load_checkpoint(
1665  restart_filename + '_' + std::to_string(member),
1666  solution_augmented_ensemble[member].block(base_state));
1667 #ifdef ADAMANTINE_WITH_CALIPER
1668  CALI_MARK_END("restart from file");
1669 #endif
1670  }
1671  solution_augmented_ensemble[member].collect_sizes();
1672 
1673  // For now we only output temperature
1674  post_processor_database.put("thermal_output", true);
1675  post_processor_ensemble.push_back(
1676  std::make_unique<adamantine::PostProcessor<dim>>(
1677  local_communicator, post_processor_database,
1678  thermal_physics_ensemble[member]->get_dof_handler(),
1679  first_local_member + member));
1680  }
1681 
1682  // PostProcessor for outputting the experimental data
1683  boost::property_tree::ptree post_processor_expt_database;
1684  // PropertyTreeInput post_processor.file_name
1685  std::string expt_file_prefix =
1686  post_processor_database.get<std::string>("filename_prefix") + ".expt";
1687  post_processor_expt_database.put("filename_prefix", expt_file_prefix);
1688  post_processor_expt_database.put("thermal_output", true);
1689  std::unique_ptr<adamantine::PostProcessor<dim>> post_processor_expt;
1690  // PropertyTreeInput experiment.output_experiment_on_mesh
1691  bool const output_experiment_on_mesh =
1692  experiment_optional_database
1693  ? experiment_optional_database->get("output_experiment_on_mesh", true)
1694  : false;
1695  // Optionally output the experimental data projected onto the mesh. Since the
1696  // experimental data is unique we may not need all the processors to write the
1697  // output files.
1698  if (output_experiment_on_mesh && (my_color == 0))
1699  {
1700  post_processor_expt = std::make_unique<adamantine::PostProcessor<dim>>(
1701  local_communicator, post_processor_expt_database,
1702  thermal_physics_ensemble[0]->get_dof_handler());
1703  }
1704 
1705  // ----- Initialize time and time stepping counters -----
1706  unsigned int progress = 0;
1707  unsigned int n_time_step = 0;
1708  double time = 0.;
1709  double activation_time_end = -1.;
1710  double da_time = -1.;
1711  if (restart == true)
1712  {
1713  if (global_rank == 0)
1714  {
1715  std::cout << "Restarting from file" << std::endl;
1716  }
1717  std::ifstream file{restart_filename + "_time.txt"};
1718  boost::archive::text_iarchive ia{file};
1719  ia >> time;
1720  ia >> n_time_step;
1721  }
1722 
1723  // PropertyTreeInput geometry.deposition_time
1724  double const activation_time =
1725  geometry_database.get<double>("deposition_time", 0.);
1726 
1727  // ----- Read the experimental data -----
1728  std::vector<std::vector<double>> frame_time_stamps;
1729  std::unique_ptr<adamantine::ExperimentalData<dim>> experimental_data;
1730  unsigned int experimental_frame_index =
1731  std::numeric_limits<unsigned int>::max();
1732 
1733  if (experiment_optional_database)
1734  {
1735  auto experiment_database = experiment_optional_database.get();
1736 
1737  // PropertyTreeInput experiment.read_in_experimental_data
1738  bool const read_in_experimental_data =
1739  experiment_database.get("read_in_experimental_data", false);
1740 
1741  if (read_in_experimental_data)
1742  {
1743  if (global_rank == 0)
1744  std::cout << "Reading the experimental log file..." << std::endl;
1745 
1746  frame_time_stamps =
1747  adamantine::read_frame_timestamps(experiment_database);
1748 
1749  adamantine::ASSERT_THROW(frame_time_stamps.size() > 0,
1750  "Experimental data parsing is activated, but "
1751  "the log shows zero cameras.");
1752  adamantine::ASSERT_THROW(frame_time_stamps[0].size() > 0,
1753  "Experimental data parsing is activated, but "
1754  "the log shows zero data frames.");
1755 
1756  if (global_rank == 0)
1757  std::cout << "Done. Log entries found for " << frame_time_stamps.size()
1758  << " camera(s), with " << frame_time_stamps[0].size()
1759  << " frame(s)." << std::endl;
1760 
1761  // PropertyTreeInput experiment.format
1762  std::string experiment_format =
1763  experiment_database.get<std::string>("format");
1764 
1765  if (boost::iequals(experiment_format, "point_cloud"))
1766  {
1767  experimental_data = std::make_unique<adamantine::PointCloud<dim>>(
1768  adamantine::PointCloud<dim>(experiment_database));
1769  }
1770  else
1771  {
1772  if constexpr (dim == 3)
1773  {
1774  experimental_data =
1775  std::make_unique<adamantine::RayTracing>(adamantine::RayTracing(
1776  experiment_database,
1777  thermal_physics_ensemble[0]->get_dof_handler()));
1778  }
1779  }
1780 
1781  // Advance the experimental frame counter upon restart, if necessary
1782  if (restart)
1783  {
1784  bool found_frame = false;
1785  while (!found_frame)
1786  {
1787  if (frame_time_stamps[0][experimental_frame_index + 1] < time)
1788  {
1789  experimental_frame_index++;
1790  }
1791  else
1792  {
1793  found_frame = true;
1794  }
1795  }
1796  }
1797  }
1798  }
1799 
1800  // ----- Output the initial solution -----
1801  std::unique_ptr<adamantine::MechanicalPhysics<
1802  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>
1803  mechanical_physics;
1804  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
1805  displacement;
1806  for (unsigned int member = 0; member < local_ensemble_size; ++member)
1807  {
1808  output_pvtu(*post_processor_ensemble[member], n_time_step, time,
1809  thermal_physics_ensemble[member],
1810  solution_augmented_ensemble[member].block(base_state),
1811  mechanical_physics, displacement,
1812  *material_properties_ensemble[member], timers);
1813  }
1814 
1815  // ----- Increment the time step -----
1816  ++n_time_step;
1817 
1818  // ----- Get refinement and time stepping parameters -----
1819  // PropertyTreeInput refinement.time_steps_between_refinement
1820  unsigned int const time_steps_refinement =
1821  refinement_database.get("time_steps_between_refinement", 10);
1822  // PropertyTreeInput time_stepping.time_step
1823  double time_step = time_stepping_database.get<double>("time_step");
1824  double const da_time_half_window = 1.01 * time_step;
1825  // PropertyTreeInput time_stepping.scan_path_for_duration
1826  bool const scan_path_for_duration =
1827  time_stepping_database.get("scan_path_for_duration", false);
1828  // PropertyTreeInput time_stepping.duration
1829  double const duration = scan_path_for_duration
1830  ? std::numeric_limits<double>::max()
1831  : time_stepping_database.get<double>("duration");
1832  // PropertyTreeInput post_processor.time_steps_between_output
1833  unsigned int const time_steps_output =
1834  post_processor_database.get("time_steps_between_output", 1);
1835  // PropertyTreeInput post_processor.output_on_data_assimilation
1836  bool const output_on_da =
1837  post_processor_database.get("output_on_data_assimilation", true);
1838 
1839  // ----- Deposit material -----
1840  // For now assume that all ensemble members share the same geometry (they
1841  // have independent adamantine::Geometry objects, but all are constructed
1842  // from identical parameters), base new additions on the 0th ensemble member
1843  // since all the sources use the same scan path.
1844  auto [material_deposition_boxes, deposition_times, deposition_cos,
1845  deposition_sin] =
1846  adamantine::create_material_deposition_boxes<dim>(
1847  geometry_database, heat_sources_ensemble[0]);
1848 
1849  // ----- Compute bounding heat sources -----
1850  // When using AMR, we refine the cells that the heat sources intersect. Since
1851  // each ensemble members can have slightly different sources, we create a new
1852  // set of heat sources that encompasses the member sources. We use these new
1853  // sources when refining the mesh.
1854  auto bounding_heat_sources = adamantine::get_bounding_heat_sources<dim>(
1855  database_ensemble, global_communicator);
1856 
1857  // ----- Main time stepping loop -----
1858  if (global_rank == 0)
1859  std::cout << "Starting the main time stepping loop..." << std::endl;
1860 
1861 #ifdef ADAMANTINE_WITH_CALIPER
1862  CALI_CXX_MARK_LOOP_BEGIN(main_loop_id, "main_loop");
1863 #endif
1864  while (time < duration)
1865  {
1866 #ifdef ADAMANTINE_WITH_CALIPER
1867  CALI_CXX_MARK_LOOP_ITERATION(main_loop_id, n_time_step - 1);
1868 #endif
1869  if ((time + time_step) > duration)
1870  time_step = duration - time;
1871 
1872  // ----- Refine the mesh if necessary -----
1873  // Refine the mesh the first time we get in the loop and after
1874  // time_steps_refinement time steps.
1875  if ((n_time_step == 1) || ((n_time_step % time_steps_refinement) == 0))
1876  {
1877  timers[adamantine::refine].start();
1878  double const next_refinement_time =
1879  time + time_steps_refinement * time_step;
1880 
1881  for (unsigned int member = 0; member < local_ensemble_size; ++member)
1882  {
1883  // FIXME
1884  std::unique_ptr<adamantine::MechanicalPhysics<
1885  dim, n_materials, p_order, MaterialStates, MemorySpaceType>>
1886  dummy;
1887  refine_mesh(thermal_physics_ensemble[member], dummy,
1888  *material_properties_ensemble[member],
1889  solution_augmented_ensemble[member].block(base_state),
1890  bounding_heat_sources, time, next_refinement_time,
1891  time_steps_refinement, refinement_database);
1892  solution_augmented_ensemble[member].collect_sizes();
1893  }
1894 
1895  timers[adamantine::refine].stop();
1896  if ((global_rank == 0) && (verbose_output == true))
1897  {
1898  std::cout << "n_time_step: " << n_time_step << " time: " << time
1899  << " n_dofs: "
1900  << thermal_physics_ensemble[0]->get_dof_handler().n_dofs()
1901  << std::endl;
1902  }
1903  }
1904 
1905  // We use an epsilon to get the "expected" behavior when the deposition
1906  // time and the time match should match exactly but don't because of
1907  // floating point accuracy.
1908  timers[adamantine::add_material_activate].start();
1909  if (time > activation_time_end)
1910  {
1911  // If we use scan_path_for_duration, we may need to read the scan path
1912  // file once again.
1913  if (scan_path_for_duration)
1914  {
1915  // Check if we have reached the end of current scan path.
1916  bool need_updated_scan_path = false;
1917  {
1918  for (auto &source : heat_sources_ensemble[0])
1919  {
1920  if (time >
1921  source->get_scan_path().get_segment_list().back().end_time)
1922  {
1923  need_updated_scan_path = true;
1924  break;
1925  }
1926  }
1927  }
1928 
1929  if (need_updated_scan_path)
1930  {
1931  // Check if we have reached the end of the file. If not, read the
1932  // updated scan path file. We assume that the scan paths are identical
1933  // for all the heat sources and thus, we can use
1934  // heat_sources_ensemble[0] to get the material deposition boxes and
1935  // times. We still need every ensemble member to read the scan path in
1936  // order to compute the correct heat sources.
1937  bool scan_path_end = true;
1938  for (unsigned int member = 0; member < local_ensemble_size; ++member)
1939  {
1940  for (auto &source : heat_sources_ensemble[member])
1941  {
1942  if (!source->get_scan_path().is_finished())
1943  {
1944  scan_path_end = false;
1945  // This functions waits for the scan path file to be updated
1946  // before reading the file.
1947  source->get_scan_path().read_file();
1948  }
1949  }
1950  }
1951 
1952  // If we have reached the end of scan path file for all the heat
1953  // sources, we just exit.
1954  if (scan_path_end)
1955  {
1956  break;
1957  }
1958 
1959  std::tie(material_deposition_boxes, deposition_times, deposition_cos,
1960  deposition_sin) =
1961  adamantine::create_material_deposition_boxes<dim>(
1962  geometry_database, heat_sources_ensemble[0]);
1963  }
1964  }
1965 
1966  double const eps = time_step / 1e10;
1967  auto activation_start =
1968  std::lower_bound(deposition_times.begin(), deposition_times.end(),
1969  time - eps) -
1970  deposition_times.begin();
1971  activation_time_end =
1972  std::min(time + std::max(activation_time, time_step), duration) - eps;
1973  auto activation_end =
1974  std::lower_bound(deposition_times.begin(), deposition_times.end(),
1975  activation_time_end) -
1976  deposition_times.begin();
1977  if (activation_start < activation_end)
1978  {
1979 #ifdef ADAMANTINE_WITH_CALIPER
1980  CALI_MARK_BEGIN("add material");
1981 #endif
1982  for (unsigned int member = 0; member < local_ensemble_size; ++member)
1983  {
1984  // Compute the elements to activate.
1985  // TODO Right now, we compute the list of cells that get activated
1986  // for the entire material deposition. We should restrict the list
1987  // to the cells that are activated between activation_start and
1988  // activation_end.
1989  timers[adamantine::add_material_search].start();
1990  auto elements_to_activate = adamantine::get_elements_to_activate(
1991  thermal_physics_ensemble[member]->get_dof_handler(),
1992  material_deposition_boxes);
1993  timers[adamantine::add_material_search].stop();
1994  // For now assume that all deposited material has never been
1995  // melted (may or may not be reasonable)
1996  std::vector<bool> has_melted(deposition_cos.size(), false);
1997 
1998  thermal_physics_ensemble[member]->add_material_start(
1999  elements_to_activate, deposition_cos, deposition_sin, has_melted,
2000  activation_start, activation_end,
2001  solution_augmented_ensemble[member].block(base_state));
2002 
2003 #ifdef ADAMANTINE_WITH_CALIPER
2004  CALI_MARK_BEGIN("refine triangulation");
2005 #endif
2006  dealii::DoFHandler<dim> &dof_handler =
2007  thermal_physics_ensemble[member]->get_dof_handler();
2008  dealii::parallel::distributed::Triangulation<dim> &triangulation =
2009  dynamic_cast<dealii::parallel::distributed::Triangulation<dim> &>(
2010  const_cast<dealii::Triangulation<dim> &>(
2011  dof_handler.get_triangulation()));
2012 
2013  triangulation.execute_coarsening_and_refinement();
2014 #ifdef ADAMANTINE_WITH_CALIPER
2015  CALI_MARK_END("refine triangulation");
2016 #endif
2017 
2018  thermal_physics_ensemble[member]->add_material_end(
2019  database_ensemble[member].get(
2020  "materials.new_material_temperature", 300.),
2021  solution_augmented_ensemble[member].block(base_state));
2022 
2023  solution_augmented_ensemble[member].collect_sizes();
2024  }
2025 #ifdef ADAMANTINE_WITH_CALIPER
2026  CALI_MARK_END("add material");
2027 #endif
2028  }
2029 
2030  if ((global_rank == 0) && (verbose_output == true) &&
2031  (activation_end - activation_start > 0))
2032  {
2033  std::cout << "n_time_step: " << n_time_step << " time: " << time
2034  << " n_dofs after cell activation: "
2035  << thermal_physics_ensemble[0]->get_dof_handler().n_dofs()
2036  << std::endl;
2037  }
2038  }
2039  timers[adamantine::add_material_activate].stop();
2040 
2041  // ----- Evolve the solution by one time step -----
2042  double const old_time = time;
2043  timers[adamantine::evol_time].start();
2044 
2045  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2046  {
2047  time = thermal_physics_ensemble[member]->evolve_one_time_step(
2048  old_time, time_step,
2049  solution_augmented_ensemble[member].block(base_state), timers);
2050  }
2051  timers[adamantine::evol_time].stop();
2052 
2053  // ----- Perform data assimilation -----
2054  if (assimilate_data)
2055  {
2056  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2057  {
2058  thermal_physics_ensemble[member]->get_affine_constraints().distribute(
2059  solution_augmented_ensemble[member].block(base_state));
2060  }
2061 
2062  // Currently assume that all frames are synced so that the 0th camera
2063  // frame time is the relevant time
2064  double frame_time = std::numeric_limits<double>::max();
2065  if ((experimental_frame_index + 1) < frame_time_stamps[0].size())
2066  {
2067  if (time > da_time + da_time_half_window)
2068  {
2069  da_time = frame_time_stamps[0][experimental_frame_index + 1];
2070  }
2071  if (frame_time_stamps[0][experimental_frame_index + 1] <= time)
2072  {
2073  experimental_frame_index = experimental_data->read_next_frame();
2074  frame_time = frame_time_stamps[0][experimental_frame_index];
2075  }
2076  }
2077 
2078  if (frame_time <= time)
2079  {
2080  adamantine::ASSERT_THROW(
2081  frame_time > old_time || n_time_step == 1,
2082  "Unexpectedly missed a data assimilation frame.");
2083  if (global_rank == 0)
2084  std::cout << "Performing data assimilation at time " << time << "..."
2085  << std::endl;
2086 
2087  // Print out the augmented parameters
2088  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2089  {
2090  std::cout << "Rank: " << global_rank
2091  << " | Old parameters for member "
2092  << first_local_member + member << ": ";
2093  for (auto param : solution_augmented_ensemble[member].block(1))
2094  std::cout << param << " ";
2095 
2096  std::cout << std::endl;
2097  }
2098 
2099  timers[adamantine::da_experimental_data].start();
2100 #ifdef ADAMANTINE_WITH_CALIPER
2101  CALI_MARK_BEGIN("da_experimental_data");
2102 #endif
2103  auto points_values = experimental_data->get_points_values();
2104  auto const &thermal_dof_handler =
2105  thermal_physics_ensemble[0]->get_dof_handler();
2106  auto expt_to_dof_mapping = adamantine::get_expt_to_dof_mapping(
2107  points_values, thermal_dof_handler);
2108  std::cout << "Rank: " << global_rank
2109  << " | Number expt sites mapped to DOFs: "
2110  << expt_to_dof_mapping.first.size() << std::endl;
2111 #ifdef ADAMANTINE_WITH_CALIPER
2112  CALI_MARK_END("da_experimental_data");
2113 #endif
2114  timers[adamantine::da_experimental_data].stop();
2115 
2116  if (post_processor_expt)
2117  {
2118  dealii::LA::distributed::Vector<double, MemorySpaceType>
2119  temperature_expt;
2120  temperature_expt.reinit(
2121  solution_augmented_ensemble[0].block(base_state));
2122  temperature_expt.add(1.0e10);
2123  adamantine::set_with_experimental_data(
2124  global_communicator, points_values, expt_to_dof_mapping,
2125  temperature_expt, verbose_output);
2126 
2127  thermal_physics_ensemble[0]->get_affine_constraints().distribute(
2128  temperature_expt);
2129  post_processor_expt
2130  ->template write_thermal_output<typename Kokkos::View<
2131  double **,
2132  typename MemorySpaceType::kokkos_space>::array_layout>(
2133  n_time_step, time, temperature_expt,
2134  material_properties_ensemble[0]->get_state(),
2135  material_properties_ensemble[0]->get_dofs_map(),
2136  material_properties_ensemble[0]->get_dof_handler());
2137  }
2138 
2139  // Only continue data assimilation if some of the observations are
2140  // mapped to DOFs
2141  if (expt_to_dof_mapping.first.size() > 0)
2142  {
2143  // NOTE: As is, this updates the dof mapping and covariance sparsity
2144  // pattern for every data assimilation operation. Strictly, this is
2145  // only necessary if the mesh changes (both updates) or the locations
2146  // of observations changes (the dof mapping). In practice, changes to
2147  // the mesh due to deposition likely cause the updates to be required
2148  // for each operation. If this is a bottleneck, it can be fixed in the
2149  // future.
2150  timers[adamantine::da_dof_mapping].start();
2151 #ifdef ADAMANTINE_WITH_CALIPER
2152  CALI_MARK_BEGIN("da_dof_mapping");
2153 #endif
2154  data_assimilator.update_dof_mapping<dim>(expt_to_dof_mapping);
2155 #ifdef ADAMANTINE_WITH_CALIPER
2156  CALI_MARK_END("da_dof_mapping");
2157 #endif
2158  timers[adamantine::da_dof_mapping].stop();
2159 
2160  timers[adamantine::da_covariance_sparsity].start();
2161 #ifdef ADAMANTINE_WITH_CALIPER
2162  CALI_MARK_BEGIN("da_covariance_sparsity");
2163 #endif
2164  data_assimilator.update_covariance_sparsity_pattern<dim>(
2165  thermal_dof_handler,
2166  solution_augmented_ensemble[0].block(augmented_state).size());
2167 #ifdef ADAMANTINE_WITH_CALIPER
2168  CALI_MARK_END("da_covariance_sparsity");
2169 #endif
2170  timers[adamantine::da_covariance_sparsity].start();
2171 
2172  unsigned int experimental_data_size = points_values.values.size();
2173 
2174  // Create the R matrix (the observation covariance matrix)
2175  // PropertyTreeInput experiment.estimated_uncertainty
2176  timers[adamantine::da_obs_covariance].start();
2177 #ifdef ADAMANTINE_WITH_CALIPER
2178  CALI_MARK_BEGIN("da_obs_covariance");
2179 #endif
2180  double variance_entries = experiment_optional_database.get().get(
2181  "estimated_uncertainty", 0.0);
2182  variance_entries = variance_entries * variance_entries;
2183 
2184  dealii::SparsityPattern pattern(experimental_data_size,
2185  experimental_data_size, 1);
2186  for (unsigned int i = 0; i < experimental_data_size; ++i)
2187  {
2188  pattern.add(i, i);
2189  }
2190  pattern.compress();
2191 
2192  dealii::SparseMatrix<double> R(pattern);
2193  for (unsigned int i = 0; i < experimental_data_size; ++i)
2194  {
2195  R.add(i, i, variance_entries);
2196  }
2197 #ifdef ADAMANTINE_WITH_CALIPER
2198  CALI_MARK_END("da_obs_covariance");
2199 #endif
2200  timers[adamantine::da_obs_covariance].stop();
2201 
2202  // Perform data assimilation to update the augmented state ensemble
2203  timers[adamantine::da_update_ensemble].start();
2204 #ifdef ADAMANTINE_WITH_CALIPER
2205  CALI_MARK_BEGIN("da_update_ensemble");
2206 #endif
2207  data_assimilator.update_ensemble(solution_augmented_ensemble,
2208  points_values.values, R);
2209 #ifdef ADAMANTINE_WITH_CALIPER
2210  CALI_MARK_END("da_update_ensemble");
2211 #endif
2212  timers[adamantine::da_update_ensemble].stop();
2213 
2214  // Extract the parameters from the augmented state
2215  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2216  {
2217  for (unsigned int index = 0;
2218  index < augmented_state_parameters.size(); ++index)
2219  {
2220  // FIXME: Need to consider how we want to generalize this. It
2221  // could get unwieldy if we want to specify every parameter of an
2222  // arbitrary number of beams.
2223  if (augmented_state_parameters.at(index) ==
2224  adamantine::AugmentedStateParameters::beam_0_absorption)
2225  {
2226  database_ensemble[member].put(
2227  "sources.beam_0.absorption_efficiency",
2228  solution_augmented_ensemble[member].block(
2229  augmented_state)[index]);
2230  }
2231  else if (augmented_state_parameters.at(index) ==
2232  adamantine::AugmentedStateParameters::beam_0_max_power)
2233  {
2234  database_ensemble[member].put(
2235  "sources.beam_0.max_power",
2236  solution_augmented_ensemble[member].block(
2237  augmented_state)[index]);
2238  }
2239  }
2240  }
2241 
2242  if (global_rank == 0)
2243  std::cout << "Done." << std::endl;
2244 
2245  // Print out the augmented parameters
2246  if (solution_augmented_ensemble[0].block(1).size() > 0)
2247  {
2248  for (unsigned int member = 0; member < local_ensemble_size;
2249  ++member)
2250  {
2251  std::cout << "Rank: " << global_rank
2252  << " | New parameters for member "
2253  << first_local_member + member << ": ";
2254  for (auto param : solution_augmented_ensemble[member].block(1))
2255  std::cout << param << " ";
2256 
2257  std::cout << std::endl;
2258  }
2259  }
2260  }
2261  else
2262  {
2263  if (global_rank == 0)
2264  std::cout
2265  << "WARNING: NO EXPERIMENTAL DATA POINTS MAPPED ONTO THE "
2266  "SIMULATION MESH. SKIPPING DATA ASSIMILATION OPERATION."
2267  << std::endl;
2268  }
2269  }
2270 
2271  // Update the heat source in the ThermalPhysics objects
2272  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2273  {
2274  thermal_physics_ensemble[member]->update_physics_parameters(
2275  database_ensemble[member].get_child("sources"));
2276  }
2277  }
2278 
2279  // ----- Checkpoint the ensemble members -----
2280  if (n_time_step % time_steps_checkpoint == 0)
2281  {
2282 #ifdef ADAMANTINE_WITH_CALIPER
2283  CALI_MARK_BEGIN("save checkpoint");
2284 #endif
2285  if (global_rank == 0)
2286  {
2287  std::cout << "Checkpoint reached" << std::endl;
2288  }
2289 
2290  std::string output_dir =
2291  post_processor_database.get<std::string>("output_dir", "");
2292  std::string filename_prefix =
2293  checkpoint_overwrite
2294  ? checkpoint_filename
2295  : checkpoint_filename + '_' + std::to_string(n_time_step);
2296  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2297  {
2298  thermal_physics_ensemble[member]->save_checkpoint(
2299  filename_prefix + '_' + std::to_string(first_local_member + member),
2300  solution_augmented_ensemble[member].block(base_state));
2301  }
2302  std::ofstream file{output_dir + filename_prefix + "_time.txt"};
2303  boost::archive::text_oarchive oa{file};
2304  oa << time;
2305  oa << n_time_step;
2306 #ifdef ADAMANTINE_WITH_CALIPER
2307  CALI_MARK_END("save checkpoint");
2308 #endif
2309  }
2310 
2311  // ----- Output progress on screen -----
2312  if (global_rank == 0)
2313  {
2314  double adim_time = time / (duration / 10.);
2315  double int_part = 0;
2316  std::modf(adim_time, &int_part);
2317  if (int_part > progress)
2318  {
2319  std::cout << int_part * 10 << '%' << " completed" << std::endl;
2320  ++progress;
2321  }
2322  }
2323 
2324  // ----- Output the solution -----
2325  if ((n_time_step % time_steps_output == 0) ||
2326  (output_on_da && time > (da_time - da_time_half_window) &&
2327  time < (da_time + da_time_half_window)))
2328  {
2329  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2330  {
2331  thermal_physics_ensemble[member]->set_state_to_material_properties();
2332  output_pvtu(*post_processor_ensemble[member], n_time_step, time,
2333  thermal_physics_ensemble[member],
2334  solution_augmented_ensemble[member].block(base_state),
2335  mechanical_physics, displacement,
2336  *material_properties_ensemble[member], timers);
2337  }
2338  }
2339 
2340  ++n_time_step;
2341  }
2342 
2343 #ifdef ADAMANTINE_WITH_CALIPER
2344  CALI_CXX_MARK_LOOP_END(main_loop_id);
2345 #endif
2346 
2347  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2348  {
2349  post_processor_ensemble[member]->write_pvd();
2350  }
2351 
2352  // This is only used for integration test
2353  if constexpr (std::is_same_v<MemorySpaceType, dealii::MemorySpace::Host>)
2354  {
2355  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2356  {
2357  thermal_physics_ensemble[member]->get_affine_constraints().distribute(
2358  solution_augmented_ensemble[member].block(base_state));
2359  }
2360  return solution_augmented_ensemble;
2361  }
2362  else
2363  {
2364  // NOTE: Currently unused. Added for the future case where run_ensemble is
2365  // functional on the device.
2366  std::vector<dealii::LA::distributed::BlockVector<double>>
2367  solution_augmented_ensemble_host(local_ensemble_size);
2368 
2369  for (unsigned int member = 0; member < local_ensemble_size; ++member)
2370  {
2371  solution_augmented_ensemble[member].reinit(2);
2372 
2373  solution_augmented_ensemble_host[member].block(0).reinit(
2374  solution_augmented_ensemble[member].block(0).get_partitioner());
2375 
2376  solution_augmented_ensemble_host[member].block(0).import(
2377  solution_augmented_ensemble[member].block(0),
2378  dealii::VectorOperation::insert);
2379  thermal_physics_ensemble[member]->get_affine_constraints().distribute(
2380  solution_augmented_ensemble_host[member].block(0));
2381 
2382  solution_augmented_ensemble_host[member].block(1).reinit(
2383  solution_augmented_ensemble[member].block(0).get_partitioner());
2384 
2385  solution_augmented_ensemble_host[member].block(1).import(
2386  solution_augmented_ensemble[member].block(1),
2387  dealii::VectorOperation::insert);
2388  }
2389 
2390  return solution_augmented_ensemble_host;
2391  }
2392 }
2393 #endif
PostProcessor.hh
run_ensemble
std::vector< dealii::LA::distributed::BlockVector< double > > run_ensemble(MPI_Comm const &global_communicator, boost::property_tree::ptree const &database, std::vector< adamantine::Timer > &timers)
Definition: adamantine.hh:1423
compute_cells_to_refine
std::vector< typename dealii::parallel::distributed::Triangulation< dim >::active_cell_iterator > 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 >>> const &heat_sources)
Definition: adamantine.hh:534
refine_and_transfer
void refine_and_transfer(std::unique_ptr< adamantine::ThermalPhysicsInterface< dim, MemorySpaceType >> &thermal_physics, std::unique_ptr< adamantine::MechanicalPhysics< dim, n_materials, p_order, MaterialStates, MemorySpaceType >> &mechanical_physics, adamantine::MaterialProperty< dim, n_materials, p_order, MaterialStates, MemorySpaceType > &material_properties, dealii::DoFHandler< dim > &dof_handler, dealii::LA::distributed::Vector< double, MemorySpaceType > &solution)
Definition: adamantine.hh:308
split_global_communicator
void split_global_communicator(MPI_Comm global_communicator, unsigned int global_ensemble_size, MPI_Comm &local_communicator, unsigned int &local_ensemble_size, unsigned int &first_local_member, int &my_color)
Definition: adamantine.hh:1367
initialize_thermal_physics
std::unique_ptr< adamantine::ThermalPhysicsInterface< dim, MemorySpaceType > > initialize_thermal_physics(unsigned int fe_degree, std::string const &quadrature_type, MPI_Comm const &communicator, boost::property_tree::ptree const &database, adamantine::Geometry< dim > &geometry, adamantine::Boundary const &boundary, adamantine::MaterialProperty< dim, n_materials, p_order, MaterialStates, MemorySpaceType > &material_properties)
Definition: adamantine.hh:257
initialize
std::unique_ptr< adamantine::ThermalPhysicsInterface< dim, MemorySpaceType > > initialize(MPI_Comm const &communicator, boost::property_tree::ptree const &database, adamantine::Geometry< dim > &geometry, adamantine::Boundary const &boundary, adamantine::MaterialProperty< dim, n_materials, p_order, MaterialStates, MemorySpaceType > &material_properties)
Definition: adamantine.hh:218
output_pvtu
void output_pvtu(adamantine::PostProcessor< dim > &post_processor, unsigned int n_time_step, double time, std::unique_ptr< adamantine::ThermalPhysicsInterface< dim, MemorySpaceType >> const &thermal_physics, dealii::LinearAlgebra::distributed::Vector< double, MemorySpaceType > &temperature, std::unique_ptr< adamantine::MechanicalPhysics< dim, n_materials, p_order, MaterialStates, MemorySpaceType >> const &mechanical_physics, dealii::LA::distributed::Vector< double, dealii::MemorySpace::Host > &displacement, adamantine::MaterialProperty< dim, n_materials, p_order, MaterialStates, MemorySpaceType > const &material_properties, std::vector< adamantine::Timer > &timers)
Definition: adamantine.hh:64
refine_mesh
void refine_mesh(std::unique_ptr< adamantine::ThermalPhysicsInterface< dim, MemorySpaceType >> &thermal_physics, std::unique_ptr< adamantine::MechanicalPhysics< dim, n_materials, p_order, MaterialStates, MemorySpaceType >> &mechanical_physics, adamantine::MaterialProperty< dim, n_materials, p_order, MaterialStates, MemorySpaceType > &material_properties, dealii::LA::distributed::Vector< double, MemorySpaceType > &solution, std::vector< std::shared_ptr< adamantine::HeatSource< 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)
Definition: adamantine.hh:599
initialize_timers
void initialize_timers(MPI_Comm const &communicator, std::vector< adamantine::Timer > &timers)
Definition: adamantine.hh:189
run
std::pair< dealii::LinearAlgebra::distributed::Vector< double, dealii::MemorySpace::Host >, dealii::LinearAlgebra::distributed::Vector< double, dealii::MemorySpace::Host > > run(MPI_Comm const &communicator, boost::property_tree::ptree const &database, std::vector< adamantine::Timer > &timers)
Definition: adamantine.hh:744
initialize_quadrature
std::unique_ptr< adamantine::ThermalPhysicsInterface< dim, MemorySpaceType > > initialize_quadrature(std::string const &quadrature_type, MPI_Comm const &communicator, boost::property_tree::ptree const &database, adamantine::Geometry< dim > &geometry, adamantine::Boundary const &boundary, adamantine::MaterialProperty< dim, n_materials, p_order, MaterialStates, MemorySpaceType > &material_properties)
Definition: adamantine.hh:233
adamantine::DataAssimilator::update_dof_mapping
void update_dof_mapping(std::pair< std::vector< int >, std::vector< int >> const &expt_to_dof_mapping)
Definition: DataAssimilator.cc:500
adamantine::DataAssimilator::update_covariance_sparsity_pattern
void update_covariance_sparsity_pattern(dealii::DoFHandler< dim > const &dof_handler, const unsigned int parameter_size)
Definition: DataAssimilator.cc:508
adamantine::DataAssimilator::update_ensemble
void update_ensemble(std::vector< dealii::LA::distributed::BlockVector< double >> &augmented_state_ensemble, std::vector< double > const &expt_data, dealii::SparseMatrix< double > const &R)
Definition: DataAssimilator.cc:90
adamantine::Geometry::get_triangulation
dealii::parallel::distributed::Triangulation< dim > & get_triangulation()
Definition: Geometry.hh:48
adamantine::MaterialProperty::get
double get(dealii::types::material_id material_id, Property prop) const
Definition: MaterialProperty.hh:440
adamantine::PostProcessor
Definition: PostProcessor.hh:41
ensemble_management.hh
experimental_data_utils.hh
material_deposition.hh
adamantine::get_expt_to_dof_mapping
std::pair< std::vector< int >, std::vector< int > > get_expt_to_dof_mapping(PointsValues< dim > const &points_values, dealii::DoFHandler< dim > const &dof_handler)
Definition: experimental_data_utils.cc:70
adamantine::output
@ output
Definition: types.hh:146
adamantine::evol_time
@ evol_time
Definition: types.hh:143
adamantine::da_obs_covariance
@ da_obs_covariance
Definition: types.hh:141
adamantine::da_experimental_data
@ da_experimental_data
Definition: types.hh:138
adamantine::da_dof_mapping
@ da_dof_mapping
Definition: types.hh:139
adamantine::add_material_search
@ add_material_search
Definition: types.hh:136
adamantine::da_covariance_sparsity
@ da_covariance_sparsity
Definition: types.hh:140
adamantine::da_update_ensemble
@ da_update_ensemble
Definition: types.hh:142
adamantine::refine
@ refine
Definition: types.hh:135
adamantine::add_material_activate
@ add_material_activate
Definition: types.hh:137
adamantine::create_database_ensemble
std::vector< boost::property_tree::ptree > create_database_ensemble(boost::property_tree::ptree const &database, MPI_Comm local_communicator, unsigned int first_local_member, unsigned int local_ensemble_size, unsigned int global_ensemble_size)
Definition: ensemble_management.cc:213
adamantine::ASSERT_THROW
void ASSERT_THROW(bool cond, std::string const &message)
Definition: utils.hh:70
adamantine::get_elements_to_activate
std::vector< std::vector< typename dealii::DoFHandler< dim >::active_cell_iterator > > get_elements_to_activate(dealii::DoFHandler< dim > const &dof_handler, std::vector< dealii::BoundingBox< dim >> const &material_deposition_boxes)
Definition: material_deposition.cc:314
adamantine::set_with_experimental_data
void set_with_experimental_data(MPI_Comm const &communicator, PointsValues< dim > const &points_values, std::pair< std::vector< int >, std::vector< int >> &expt_to_dof_mapping, dealii::LinearAlgebra::distributed::Vector< double > &temperature, bool verbose_output)
Definition: experimental_data_utils.cc:109
adamantine::read_frame_timestamps
std::vector< std::vector< double > > read_frame_timestamps(boost::property_tree::ptree const &experiment_database)
Definition: experimental_data_utils.cc:137
ASSERT
#define ASSERT(condition, message)
Definition: utils.hh:68
Generated by doxygen 1.9.1