OGS
HydroMechanicsProcess.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14
15#include "HydroMechanicsFEM.h"
20#include "ProcessLib/Process.h"
22
23namespace ProcessLib
24{
25namespace HydroMechanics
26{
27template <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 _integration_point_writer.emplace_back(
51 std::make_unique<IntegrationPointWriter>(
52 "sigma_ip",
53 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
54 integration_order, _local_assemblers, &LocalAssemblerIF::getSigma));
55
56 _integration_point_writer.emplace_back(
57 std::make_unique<IntegrationPointWriter>(
58 "epsilon_ip",
59 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
60 integration_order, _local_assemblers,
62}
63
64template <int DisplacementDim>
66{
67 return false;
68}
69
70template <int DisplacementDim>
73 const int process_id) const
74{
75 // For the monolithic scheme or the M process (deformation) in the staggered
76 // scheme.
77 if (_use_monolithic_scheme || process_id == 1)
78 {
79 auto const& l = *_local_to_global_index_map;
80 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
81 &l.getGhostIndices(), &this->_sparsity_pattern};
82 }
83
84 // For staggered scheme and H process (pressure).
85 auto const& l = *_local_to_global_index_map_with_base_nodes;
86 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
87 &l.getGhostIndices(), &_sparsity_pattern_with_linear_element};
88}
89
90template <int DisplacementDim>
92{
93 // Create single component dof in every of the mesh's nodes.
94 _mesh_subset_all_nodes = std::make_unique<MeshLib::MeshSubset>(
95 _mesh, _mesh.getNodes(), _process_data.use_taylor_hood_elements);
96
97 // Create single component dof in the mesh's base nodes.
98 _base_nodes = MeshLib::getBaseNodes(_mesh.getElements());
99 _mesh_subset_base_nodes = std::make_unique<MeshLib::MeshSubset>(
100 _mesh, _base_nodes, _process_data.use_taylor_hood_elements);
101
102 // TODO move the two data members somewhere else.
103 // for extrapolation of secondary variables of stress or strain
104 std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
105 *_mesh_subset_all_nodes};
106 _local_to_global_index_map_single_component =
107 std::make_unique<NumLib::LocalToGlobalIndexMap>(
108 std::move(all_mesh_subsets_single_component),
109 // by location order is needed for output
111
112 if (_use_monolithic_scheme)
113 {
114 // For pressure, which is the first
115 std::vector<MeshLib::MeshSubset> all_mesh_subsets{
116 *_mesh_subset_base_nodes};
117
118 // For displacement.
119 const int monolithic_process_id = 0;
120 std::generate_n(std::back_inserter(all_mesh_subsets),
121 getProcessVariables(monolithic_process_id)[1]
122 .get()
123 .getNumberOfGlobalComponents(),
124 [&]() { return *_mesh_subset_all_nodes; });
125
126 std::vector<int> const vec_n_components{1, DisplacementDim};
127 _local_to_global_index_map =
128 std::make_unique<NumLib::LocalToGlobalIndexMap>(
129 std::move(all_mesh_subsets), vec_n_components,
131 assert(_local_to_global_index_map);
132 }
133 else
134 {
135 // For displacement equation.
136 const int process_id = 1;
137 std::vector<MeshLib::MeshSubset> all_mesh_subsets;
138 std::generate_n(std::back_inserter(all_mesh_subsets),
139 getProcessVariables(process_id)[0]
140 .get()
141 .getNumberOfGlobalComponents(),
142 [&]() { return *_mesh_subset_all_nodes; });
143
144 std::vector<int> const vec_n_components{DisplacementDim};
145 _local_to_global_index_map =
146 std::make_unique<NumLib::LocalToGlobalIndexMap>(
147 std::move(all_mesh_subsets), vec_n_components,
149
150 // For pressure equation.
151 // Collect the mesh subsets with base nodes in a vector.
152 std::vector<MeshLib::MeshSubset> all_mesh_subsets_base_nodes{
153 *_mesh_subset_base_nodes};
154 _local_to_global_index_map_with_base_nodes =
155 std::make_unique<NumLib::LocalToGlobalIndexMap>(
156 std::move(all_mesh_subsets_base_nodes),
157 // by location order is needed for output
159
160 _sparsity_pattern_with_linear_element = NumLib::computeSparsityPattern(
161 *_local_to_global_index_map_with_base_nodes, _mesh);
162
163 assert(_local_to_global_index_map);
164 assert(_local_to_global_index_map_with_base_nodes);
165 }
166}
167
168template <int DisplacementDim>
170 NumLib::LocalToGlobalIndexMap const& dof_table,
171 MeshLib::Mesh const& mesh,
172 unsigned const integration_order)
173{
176 mesh.getElements(), dof_table, _local_assemblers,
177 mesh.isAxiallySymmetric(), integration_order, _process_data);
178
179 auto add_secondary_variable = [&](std::string const& name,
180 int const num_components,
181 auto get_ip_values_function)
182 {
183 _secondary_variables.addSecondaryVariable(
184 name,
185 makeExtrapolator(num_components, getExtrapolator(),
186 _local_assemblers,
187 std::move(get_ip_values_function)));
188 };
189
190 add_secondary_variable("sigma",
192 DisplacementDim>::RowsAtCompileTime,
193 &LocalAssemblerIF::getIntPtSigma);
194
195 add_secondary_variable("epsilon",
197 DisplacementDim>::RowsAtCompileTime,
198 &LocalAssemblerIF::getIntPtEpsilon);
199
200 add_secondary_variable("velocity", DisplacementDim,
201 &LocalAssemblerIF::getIntPtDarcyVelocity);
202
203 //
204 // enable output of internal variables defined by material models
205 //
207 LocalAssemblerIF>(_process_data.solid_materials,
208 add_secondary_variable);
209
210 _process_data.pressure_interpolated =
211 MeshLib::getOrCreateMeshProperty<double>(
212 const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
214
215 _process_data.principal_stress_vector[0] =
216 MeshLib::getOrCreateMeshProperty<double>(
217 const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_1",
219
220 _process_data.principal_stress_vector[1] =
221 MeshLib::getOrCreateMeshProperty<double>(
222 const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_2",
224
225 _process_data.principal_stress_vector[2] =
226 MeshLib::getOrCreateMeshProperty<double>(
227 const_cast<MeshLib::Mesh&>(mesh), "principal_stress_vector_3",
229
230 _process_data.principal_stress_values =
231 MeshLib::getOrCreateMeshProperty<double>(
232 const_cast<MeshLib::Mesh&>(mesh), "principal_stress_values",
234
235 _process_data.permeability = MeshLib::getOrCreateMeshProperty<double>(
236 const_cast<MeshLib::Mesh&>(mesh), "permeability",
239
240 // Set initial conditions for integration point data.
241 for (auto const& ip_writer : _integration_point_writer)
242 {
243 // Find the mesh property with integration point writer's name.
244 auto const& name = ip_writer->name();
245 if (!mesh.getProperties().existsPropertyVector<double>(name))
246 {
247 continue;
248 }
249 auto const& mesh_property =
250 *mesh.getProperties().template getPropertyVector<double>(name);
251
252 // The mesh property must be defined on integration points.
253 if (mesh_property.getMeshItemType() !=
255 {
256 continue;
257 }
258
259 auto const ip_meta_data =
261
262 // Check the number of components.
263 if (ip_meta_data.n_components !=
264 mesh_property.getNumberOfGlobalComponents())
265 {
266 OGS_FATAL(
267 "Different number of components in meta data ({:d}) than in "
268 "the integration point field data for '{:s}': {:d}.",
269 ip_meta_data.n_components, name,
270 mesh_property.getNumberOfGlobalComponents());
271 }
272
273 // Now we have a properly named vtk's field data array and the
274 // corresponding meta data.
275 std::size_t position = 0;
276 for (auto& local_asm : _local_assemblers)
277 {
278 std::size_t const integration_points_read =
279 local_asm->setIPDataInitialConditions(
280 name, &mesh_property[position],
281 ip_meta_data.integration_order);
282 if (integration_points_read == 0)
283 {
284 OGS_FATAL(
285 "No integration points read in the integration point "
286 "initial conditions set function.");
287 }
288 position += integration_points_read * ip_meta_data.n_components;
289 }
290 }
291
292 // Initialize local assemblers after all variables have been set.
293 GlobalExecutor::executeMemberOnDereferenced(&LocalAssemblerIF::initialize,
294 _local_assemblers,
295 *_local_to_global_index_map);
296}
297
298template <int DisplacementDim>
300{
301 if (_use_monolithic_scheme)
302 {
303 const int process_id_of_hydromechanics = 0;
304 initializeProcessBoundaryConditionsAndSourceTerms(
305 *_local_to_global_index_map, process_id_of_hydromechanics);
306 return;
307 }
308
309 // Staggered scheme:
310 // for the equations of pressure
311 const int hydraulic_process_id = 0;
312 initializeProcessBoundaryConditionsAndSourceTerms(
313 *_local_to_global_index_map_with_base_nodes, hydraulic_process_id);
314
315 // for the equations of deformation.
316 const int mechanical_process_id = 1;
317 initializeProcessBoundaryConditionsAndSourceTerms(
318 *_local_to_global_index_map, mechanical_process_id);
319}
320
321template <int DisplacementDim>
323 const double t, double const dt, std::vector<GlobalVector*> const& x,
324 std::vector<GlobalVector*> const& xdot, int const process_id,
326{
327 DBUG("Assemble the equations for HydroMechanics");
328
329 // Note: This assembly function is for the Picard nonlinear solver. Since
330 // only the Newton-Raphson method is employed to simulate coupled HM
331 // processes in this class, this function is actually not used so far.
332
333 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
334 dof_table = {std::ref(*_local_to_global_index_map)};
335
336 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
337
339 _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
340 pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
341 b);
342}
343
344template <int DisplacementDim>
346 assembleWithJacobianConcreteProcess(const double t, double const dt,
347 std::vector<GlobalVector*> const& x,
348 std::vector<GlobalVector*> const& xdot,
349 int const process_id, GlobalMatrix& M,
351 GlobalMatrix& Jac)
352{
353 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
354 dof_tables;
355 // For the monolithic scheme
356 if (_use_monolithic_scheme)
357 {
358 DBUG(
359 "Assemble the Jacobian of HydroMechanics for the monolithic "
360 "scheme.");
361 dof_tables.emplace_back(*_local_to_global_index_map);
362 }
363 else
364 {
365 // For the staggered scheme
366 if (process_id == 0)
367 {
368 DBUG(
369 "Assemble the Jacobian equations of liquid fluid process in "
370 "HydroMechanics for the staggered scheme.");
371 }
372 else
373 {
374 DBUG(
375 "Assemble the Jacobian equations of mechanical process in "
376 "HydroMechanics for the staggered scheme.");
377 }
378 dof_tables.emplace_back(*_local_to_global_index_map_with_base_nodes);
379 dof_tables.emplace_back(*_local_to_global_index_map);
380 }
381
382 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
383
386 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
387 process_id, M, K, b, Jac);
388
389 auto copyRhs = [&](int const variable_id, auto& output_vector)
390 {
391 if (_use_monolithic_scheme)
392 {
393 transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
394 output_vector,
395 std::negate<double>());
396 }
397 else
398 {
399 transformVariableFromGlobalVector(b, 0, dof_tables[process_id],
400 output_vector,
401 std::negate<double>());
402 }
403 };
404 if (_use_monolithic_scheme || process_id == 0)
405 {
406 copyRhs(0, *_hydraulic_flow);
407 }
408 if (_use_monolithic_scheme || process_id == 1)
409 {
410 copyRhs(1, *_nodal_forces);
411 }
412}
413
414template <int DisplacementDim>
416 std::vector<GlobalVector*> const& x, double const t, double const dt,
417 const int process_id)
418{
419 DBUG("PreTimestep HydroMechanicsProcess.");
420
421 if (hasMechanicalProcess(process_id))
422 {
424 getProcessVariables(process_id)[0];
426 &LocalAssemblerIF::preTimestep, _local_assemblers,
427 pv.getActiveElementIDs(), *_local_to_global_index_map,
428 *x[process_id], t, dt);
429 }
430}
431
432template <int DisplacementDim>
434 std::vector<GlobalVector*> const& x, double const t, double const dt,
435 const int process_id)
436{
437 if (process_id != 0)
438 {
439 return;
440 }
441
442 DBUG("PostTimestep HydroMechanicsProcess.");
443 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
444 auto const n_processes = x.size();
445 dof_tables.reserve(n_processes);
446 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
447 {
448 dof_tables.push_back(&getDOFTable(process_id));
449 }
450
451 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
453 &LocalAssemblerIF::postTimestep, _local_assemblers,
454 pv.getActiveElementIDs(), dof_tables, x, t, dt);
455}
456
457template <int DisplacementDim>
459 GlobalVector const& x, GlobalVector const& xdot, const double t,
460 double const dt, const int process_id)
461{
462 DBUG("PostNonLinearSolver HydroMechanicsProcess.");
463 // Calculate strain, stress or other internal variables of mechanics.
464 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
466 &LocalAssemblerIF::postNonLinearSolver, _local_assemblers,
467 pv.getActiveElementIDs(), getDOFTable(process_id), x, xdot, t, dt,
468 _use_monolithic_scheme, process_id);
469}
470
471template <int DisplacementDim>
473 setInitialConditionsConcreteProcess(std::vector<GlobalVector*>& x,
474 double const t,
475 int const process_id)
476{
477 if (process_id != _process_data.hydraulic_process_id)
478 {
479 return;
480 }
481
482 DBUG("Set initial conditions of HydroMechanicsProcess.");
483
484 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
486 &LocalAssemblerIF::setInitialConditions, _local_assemblers,
487 pv.getActiveElementIDs(), getDOFTable(process_id), *x[process_id], t,
488 _use_monolithic_scheme, process_id);
489}
490
491template <int DisplacementDim>
493 double const t, double const dt, std::vector<GlobalVector*> const& x,
494 GlobalVector const& x_dot, const int process_id)
495{
496 if (process_id != 0)
497 {
498 return;
499 }
500
501 DBUG("Compute the secondary variables for HydroMechanicsProcess.");
502 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
503 auto const n_processes = x.size();
504 dof_tables.reserve(n_processes);
505 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
506 {
507 dof_tables.push_back(&getDOFTable(process_id));
508 }
509
510 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
512 &LocalAssemblerIF::computeSecondaryVariable, _local_assemblers,
513 pv.getActiveElementIDs(), dof_tables, t, dt, x, x_dot, process_id);
514}
515
516template <int DisplacementDim>
517std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
519{
520 const bool manage_storage = false;
521 return std::make_tuple(_local_to_global_index_map_single_component.get(),
522 manage_storage);
523}
524
525template <int DisplacementDim>
528{
529 if (hasMechanicalProcess(process_id))
530 {
531 return *_local_to_global_index_map;
532 }
533
534 // For the equation of pressure
535 return *_local_to_global_index_map_with_base_nodes;
536}
537
538template class HydroMechanicsProcess<2>;
539template class HydroMechanicsProcess<3>;
540
541} // namespace HydroMechanics
542} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:29
Global vector based on Eigen vector.
Definition: EigenVector.h:28
bool isAxiallySymmetric() const
Definition: Mesh.h:131
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:76
Properties & getProperties()
Definition: Mesh.h:128
bool existsPropertyVector(std::string const &name) const
void postNonLinearSolverConcreteProcess(GlobalVector const &x, GlobalVector const &xdot, const double t, double const dt, int const process_id) override
std::vector< std::unique_ptr< LocalAssemblerIF > > _local_assemblers
HydroMechanicsProcess(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, HydroMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) override
void assembleWithJacobianConcreteProcess(const double t, double const, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
MeshLib::PropertyVector< double > * _nodal_forces
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) override
MeshLib::PropertyVector< double > * _hydraulic_flow
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
void assembleConcreteProcess(const double t, double const, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, const int process_id) override
std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const override
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer
Definition: Process.h:352
Handles configuration of several secondary variables from the project file.
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)
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)
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
Definition: KelvinVector.h:48
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
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 solidMaterialInternalToSecondaryVariables(std::map< int, std::unique_ptr< MaterialLib::Solids::MechanicsBase< DisplacementDim > > > const &solid_materials, AddSecondaryVariableCallback const &add_secondary_variable)
IntegrationPointMetaData getIntegrationPointMetaData(MeshLib::Properties const &properties, std::string const &name)
void createLocalAssemblersHM(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface > > &local_assemblers, 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)
std::pair< std::vector< GlobalVector * >, std::vector< GlobalVector * > > setInitialConditions(double const t0, std::vector< std::unique_ptr< ProcessData > > const &per_process_data)
Definition: TimeLoop.cpp:171
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