OGS
ThermoHydroMechanicsProcess.cpp
Go to the documentation of this file.
1 
12 
13 #include <cassert>
14 
15 #include "MeshLib/Elements/Utils.h"
17 #include "ProcessLib/Process.h"
22 
23 namespace ProcessLib
24 {
25 namespace ThermoHydroMechanics
26 {
27 template <int DisplacementDim>
29  std::string name,
30  MeshLib::Mesh& mesh,
31  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
32  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
33  unsigned const integration_order,
34  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
35  process_variables,
37  SecondaryVariableCollection&& secondary_variables,
38  bool const use_monolithic_scheme)
39  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
40  integration_order, std::move(process_variables),
41  std::move(secondary_variables), use_monolithic_scheme),
42  _process_data(std::move(process_data))
43 {
44  _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
45  mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
46 
47  _hydraulic_flow = MeshLib::getOrCreateMeshProperty<double>(
48  mesh, "MassFlowRate", MeshLib::MeshItemType::Node, 1);
49 
50  _heat_flux = MeshLib::getOrCreateMeshProperty<double>(
51  mesh, "HeatFlowRate", MeshLib::MeshItemType::Node, 1);
52 
53  _integration_point_writer.emplace_back(
54  std::make_unique<IntegrationPointWriter>(
55  "sigma_ip",
56  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
57  integration_order, _local_assemblers,
59 
60  _integration_point_writer.emplace_back(
61  std::make_unique<IntegrationPointWriter>(
62  "epsilon_ip",
63  static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
64  integration_order, _local_assemblers,
66 }
67 
68 template <int DisplacementDim>
70 {
71  return false;
72 }
73 
74 template <int DisplacementDim>
77  const int process_id) const
78 {
79  // For the monolithic scheme or the M process (deformation) in the staggered
80  // scheme.
81  if (_use_monolithic_scheme || process_id == 2)
82  {
83  auto const& l = *_local_to_global_index_map;
84  return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
85  &l.getGhostIndices(), &this->_sparsity_pattern};
86  }
87 
88  // For staggered scheme and T or H process (pressure).
89  auto const& l = *_local_to_global_index_map_with_base_nodes;
90  return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
91  &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
92 }
93 
94 template <int DisplacementDim>
96 {
97  // Create single component dof in every of the mesh's nodes.
98  _mesh_subset_all_nodes =
99  std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
100  // Create single component dof in the mesh's base nodes.
101  _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
102  _mesh_subset_base_nodes =
103  std::make_unique<MeshLib::MeshSubset>(_mesh, _base_nodes);
104 
105  // TODO move the two data members somewhere else.
106  // for extrapolation of secondary variables of stress or strain
107  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
108  *_mesh_subset_all_nodes};
109  _local_to_global_index_map_single_component =
110  std::make_unique<NumLib::LocalToGlobalIndexMap>(
111  std::move(all_mesh_subsets_single_component),
112  // by location order is needed for output
114 
115  if (_use_monolithic_scheme)
116  {
117  // For temperature, which is the first
118  std::vector<MeshLib::MeshSubset> all_mesh_subsets{
119  *_mesh_subset_base_nodes};
120 
121  // For pressure, which is the second
122  all_mesh_subsets.push_back(*_mesh_subset_base_nodes);
123 
124  // For displacement.
125  const int monolithic_process_id = 0;
126  std::generate_n(std::back_inserter(all_mesh_subsets),
127  getProcessVariables(monolithic_process_id)[2]
128  .get()
129  .getNumberOfGlobalComponents(),
130  [&]() { return *_mesh_subset_all_nodes; });
131 
132  std::vector<int> const vec_n_components{1, 1, DisplacementDim};
133  _local_to_global_index_map =
134  std::make_unique<NumLib::LocalToGlobalIndexMap>(
135  std::move(all_mesh_subsets), vec_n_components,
137  assert(_local_to_global_index_map);
138  }
139  else
140  {
141  // For displacement equation.
142  const int process_id = 2;
143  std::vector<MeshLib::MeshSubset> all_mesh_subsets;
144  std::generate_n(std::back_inserter(all_mesh_subsets),
145  getProcessVariables(process_id)[0]
146  .get()
147  .getNumberOfGlobalComponents(),
148  [&]() { return *_mesh_subset_all_nodes; });
149 
150  std::vector<int> const vec_n_components{DisplacementDim};
151  _local_to_global_index_map =
152  std::make_unique<NumLib::LocalToGlobalIndexMap>(
153  std::move(all_mesh_subsets), vec_n_components,
155 
156  // For pressure equation or temperature equation.
157  // Collect the mesh subsets with base nodes in a vector.
158  std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
159  *_mesh_subset_base_nodes};
160  _local_to_global_index_map_with_base_nodes =
161  std::make_unique<NumLib::LocalToGlobalIndexMap>(
162  std::move(all_mesh_subsets_base_nodes),
163  // by location order is needed for output
165 
166  _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
167  *_local_to_global_index_map_with_base_nodes, _mesh);
168 
169  assert(_local_to_global_index_map);
170  assert(_local_to_global_index_map_with_base_nodes);
171  }
172 }
173 
174 template <int DisplacementDim>
176  NumLib::LocalToGlobalIndexMap const& dof_table,
177  MeshLib::Mesh const& mesh,
178  unsigned const integration_order)
179 {
180  ProcessLib::createLocalAssemblersHM<DisplacementDim,
182  mesh.getElements(), dof_table, _local_assemblers,
183  NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
184  _process_data);
185 
186  _secondary_variables.addSecondaryVariable(
187  "sigma",
189  DisplacementDim>::RowsAtCompileTime,
190  getExtrapolator(), _local_assemblers,
192 
193  _secondary_variables.addSecondaryVariable(
194  "epsilon",
196  DisplacementDim>::RowsAtCompileTime,
197  getExtrapolator(), _local_assemblers,
199 
200  _secondary_variables.addSecondaryVariable(
201  "velocity",
202  makeExtrapolator(mesh.getDimension(), getExtrapolator(),
203  _local_assemblers,
205 
206  _process_data.pressure_interpolated =
207  MeshLib::getOrCreateMeshProperty<double>(
208  const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
210 
211  _process_data.temperature_interpolated =
212  MeshLib::getOrCreateMeshProperty<double>(
213  const_cast<MeshLib::Mesh&>(mesh), "temperature_interpolated",
215 
216  setIPDataInitialConditions(_integration_point_writer, mesh.getProperties(),
217  _local_assemblers);
218 
219  // Initialize local assemblers after all variables have been set.
221  &LocalAssemblerInterface::initialize, _local_assemblers,
222  *_local_to_global_index_map);
223 }
224 
225 template <int DisplacementDim>
227  DisplacementDim>::initializeBoundaryConditions()
228 {
229  if (_use_monolithic_scheme)
230  {
231  const int process_id_of_thermohydromechancs = 0;
232  initializeProcessBoundaryConditionsAndSourceTerms(
233  *_local_to_global_index_map, process_id_of_thermohydromechancs);
234  return;
235  }
236 
237  // Staggered scheme:
238  // for the equations of heat transport
239  const int thermal_process_id = 0;
240  initializeProcessBoundaryConditionsAndSourceTerms(
241  *_local_to_global_index_map_with_base_nodes, thermal_process_id);
242 
243  // for the equations of mass balance
244  const int hydraulic_process_id = 1;
245  initializeProcessBoundaryConditionsAndSourceTerms(
246  *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
247 
248  // for the equations of deformation.
249  const int mechanical_process_id = 2;
250  initializeProcessBoundaryConditionsAndSourceTerms(
251  *_local_to_global_index_map, mechanical_process_id);
252 }
253 
254 template <int DisplacementDim>
256  const double t, double const dt, std::vector<GlobalVector*> const& x,
257  std::vector<GlobalVector*> const& xdot, int const process_id,
259 {
260  DBUG("Assemble the equations for ThermoHydroMechanics");
261 
262  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
263  dof_table = {std::ref(*_local_to_global_index_map)};
264 
265  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
266 
268  _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
269  pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
270  b);
271 }
272 
273 template <int DisplacementDim>
275  assembleWithJacobianConcreteProcess(const double t, double const dt,
276  std::vector<GlobalVector*> const& x,
277  std::vector<GlobalVector*> const& xdot,
278  int const process_id, GlobalMatrix& M,
279  GlobalMatrix& K, GlobalVector& b,
280  GlobalMatrix& Jac)
281 {
282  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
283  dof_tables;
284  // For the monolithic scheme
285  if (_use_monolithic_scheme)
286  {
287  DBUG(
288  "Assemble the Jacobian of ThermoHydroMechanics for the monolithic "
289  "scheme.");
290  dof_tables.emplace_back(*_local_to_global_index_map);
291  }
292  else
293  {
294  // For the staggered scheme
295  if (process_id == 0)
296  {
297  DBUG(
298  "Assemble the Jacobian equations of heat transport process in "
299  "ThermoHydroMechanics for the staggered scheme.");
300  }
301  else if (process_id == 1)
302  {
303  DBUG(
304  "Assemble the Jacobian equations of liquid fluid process in "
305  "ThermoHydroMechanics for the staggered scheme.");
306  }
307  else
308  {
309  DBUG(
310  "Assemble the Jacobian equations of mechanical process in "
311  "ThermoHydroMechanics for the staggered scheme.");
312  }
313  dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
314  dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
315  dof_tables.emplace_back(*_local_to_global_index_map);
316  }
317 
318  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
319 
322  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
323  process_id, M, K, b, Jac);
324 
325  auto copyRhs = [&](int const variable_id, auto& output_vector)
326  {
327  if (_use_monolithic_scheme)
328  {
329  transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
330  output_vector,
331  std::negate<double>());
332  }
333  else
334  {
335  transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
336  output_vector,
337  std::negate<double>());
338  }
339  };
340  if (_use_monolithic_scheme || process_id == 0)
341  {
342  copyRhs(0, *_heat_flux);
343  }
344  if (_use_monolithic_scheme || process_id == 1)
345  {
346  copyRhs(1, *_hydraulic_flow);
347  }
348  if (_use_monolithic_scheme || process_id == 2)
349  {
350  copyRhs(2, *_nodal_forces);
351  }
352 }
353 
354 template <int DisplacementDim>
356  std::vector<GlobalVector*> const& x, double const t, double const dt,
357  const int process_id)
358 {
359  DBUG("PreTimestep ThermoHydroMechanicsProcess.");
360 
361  if (hasMechanicalProcess(process_id))
362  {
363  ProcessLib::ProcessVariable const& pv =
364  getProcessVariables(process_id)[0];
365 
367  &LocalAssemblerInterface::preTimestep, _local_assemblers,
368  pv.getActiveElementIDs(), *_local_to_global_index_map,
369  *x[process_id], t, dt);
370  }
371 }
372 
373 template <int DisplacementDim>
375  std::vector<GlobalVector*> const& x, double const t, double const dt,
376  const int process_id)
377 {
378  if (process_id != 0)
379  {
380  return;
381  }
382 
383  DBUG("PostTimestep ThermoHydroMechanicsProcess.");
384  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
385  auto const n_processes = x.size();
386  dof_tables.reserve(n_processes);
387  for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
388  {
389  dof_tables.push_back(&getDOFTable(process_id));
390  }
391 
392  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
393 
395  &LocalAssemblerInterface::postTimestep, _local_assemblers,
396  pv.getActiveElementIDs(), dof_tables, x, t, dt);
397 }
398 
399 template <int DisplacementDim>
402  GlobalVector const& xdot, const double t,
403  double const dt, const int process_id)
404 {
405  if (!hasMechanicalProcess(process_id))
406  {
407  return;
408  }
409 
410  DBUG("PostNonLinearSolver ThermoHydroMechanicsProcess.");
411  // Calculate strain, stress or other internal variables of mechanics.
412 
413  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
414 
417  pv.getActiveElementIDs(), getDOFTable(process_id), x, xdot, t, dt,
418  _use_monolithic_scheme, process_id);
419 }
420 
421 template <int DisplacementDim>
423  computeSecondaryVariableConcrete(double const t, double const dt,
424  std::vector<GlobalVector*> const& x,
425  GlobalVector const& x_dot,
426  const int process_id)
427 {
428  if (process_id != 0)
429  {
430  return;
431  }
432 
433  DBUG("Compute the secondary variables for ThermoHydroMechanicsProcess.");
434  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
435  auto const n_processes = x.size();
436  dof_tables.reserve(n_processes);
437  for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
438  {
439  dof_tables.push_back(&getDOFTable(process_id));
440  }
441 
442  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
445  pv.getActiveElementIDs(), dof_tables, t, dt, x, x_dot, process_id);
446 }
447 
448 template <int DisplacementDim>
449 std::tuple<NumLib::LocalToGlobalIndexMap*, bool> ThermoHydroMechanicsProcess<
450  DisplacementDim>::getDOFTableForExtrapolatorData() const
451 {
452  const bool manage_storage = false;
453  return std::make_tuple(_local_to_global_index_map_single_component.get(),
454  manage_storage);
455 }
456 
457 template <int DisplacementDim>
460  const int process_id) const
461 {
462  if (hasMechanicalProcess(process_id))
463  {
464  return *_local_to_global_index_map;
465  }
466 
467  // For the equation of pressure
468  return *_local_to_global_index_map_with_base_nodes;
469 }
470 
471 template class ThermoHydroMechanicsProcess<2>;
472 template class ThermoHydroMechanicsProcess<3>;
473 
474 } // namespace ThermoHydroMechanics
475 } // namespace ProcessLib
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:29
Global vector based on Eigen vector.
Definition: EigenVector.h:28
bool isAxiallySymmetric() const
Definition: Mesh.h:131
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:76
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
Properties & getProperties()
Definition: Mesh.h:128
virtual void postTimestep(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
virtual void preTimestep(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void initialize(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void computeSecondaryVariable(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id)
void postNonLinearSolver(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, GlobalVector const &xdot, double const t, double const dt, bool const use_monolithic_scheme, int const process_id)
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer
Definition: Process.h:354
Handles configuration of several secondary variables from the project file.
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
void postNonLinearSolverConcreteProcess(GlobalVector const &x, GlobalVector const &xdot, const double t, double const dt, int const process_id) override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) override
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, const int process_id) override
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
ThermoHydroMechanicsProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable >>> &&process_variables, ThermoHydroMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
std::vector< std::unique_ptr< LocalAssemblerInterface > > _local_assemblers
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
static const double t
std::vector< Node * > getBaseNodes(std::vector< Element * > const &elements)
Definition: Utils.h:26
@ BY_LOCATION
Ordering data by spatial location.
GlobalSparsityPattern computeSparsityPattern(LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh)
Computes a sparsity pattern for the given inputs.
void transformVariableFromGlobalVector(GlobalVector const &input_vector, int const variable_id, NumLib::LocalToGlobalIndexMap const &local_to_global_index_map, MeshLib::PropertyVector< double > &output_vector, Functor mapFunction)
Definition: DOFTableUtil.h:59
void setIPDataInitialConditions(std::vector< std::unique_ptr< IntegrationPointWriter >> const &_integration_point_writer, MeshLib::Properties const &mesh_properties, LocalAssemblersVector &local_assemblers)
void createLocalAssemblersHM(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface >> &local_assemblers, NumLib::IntegrationOrder const integration_order, ExtraCtorArgs &&... extra_ctor_args)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeMemberOnDereferenced(Method method, Container const &container, Args &&... args)
virtual std::vector< double > getSigma() const =0
virtual std::vector< double > getEpsilon() const =0
virtual std::vector< double > const & getIntPtSigma(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtEpsilon(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0