OGS
ThermoRichardsMechanicsProcess.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14
15#include "BaseLib/Error.h"
19#include "ProcessLib/Process.h"
22
23namespace ProcessLib
24{
25namespace ThermoRichardsMechanics
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 heat_flux_ = MeshLib::getOrCreateMeshProperty<double>(
51 mesh, "HeatFlowRate", MeshLib::MeshItemType::Node, 1);
52
53 // TODO (naumov) remove ip suffix. Probably needs modification of the mesh
54 // properties, s.t. there is no "overlapping" with cell/point data.
55 // See getOrCreateMeshProperty.
56 _integration_point_writer.emplace_back(
57 std::make_unique<IntegrationPointWriter>(
58 "sigma_ip",
59 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
60 integration_order, local_assemblers_, &LocalAssemblerIF::getSigma));
61
62 _integration_point_writer.emplace_back(
63 std::make_unique<IntegrationPointWriter>(
64 "saturation_ip", 1 /*n components*/, integration_order,
66
67 _integration_point_writer.emplace_back(
68 std::make_unique<IntegrationPointWriter>(
69 "porosity_ip", 1 /*n components*/, integration_order,
71
72 _integration_point_writer.emplace_back(
73 std::make_unique<IntegrationPointWriter>(
74 "transport_porosity_ip", 1 /*n components*/, integration_order,
76
77 _integration_point_writer.emplace_back(
78 std::make_unique<IntegrationPointWriter>(
79 "swelling_stress_ip",
80 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
81 integration_order, local_assemblers_,
83
84 _integration_point_writer.emplace_back(
85 std::make_unique<IntegrationPointWriter>(
86 "epsilon_ip",
87 static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
88 integration_order, local_assemblers_,
90}
91
92template <int DisplacementDim>
94{
95 return false;
96}
97
98template <int DisplacementDim>
101 const int /*process_id*/) const
102{
103 auto const& l = *_local_to_global_index_map;
104 return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
105 &l.getGhostIndices(), &this->_sparsity_pattern};
106}
107
108template <int DisplacementDim>
110{
111 // Create single component dof in every of the mesh's nodes.
112 _mesh_subset_all_nodes = std::make_unique<MeshLib::MeshSubset>(
113 _mesh, _mesh.getNodes(), process_data_.use_TaylorHood_elements);
114 // Create single component dof in the mesh's base nodes.
115 base_nodes_ = MeshLib::getBaseNodes(_mesh.getElements());
116 mesh_subset_base_nodes_ = std::make_unique<MeshLib::MeshSubset>(
117 _mesh, base_nodes_, process_data_.use_TaylorHood_elements);
118
119 // TODO move the two data members somewhere else.
120 // for extrapolation of secondary variables of stress or strain
121 std::vector<MeshLib::MeshSubset> all__meshsubsets_single_component{
122 *_mesh_subset_all_nodes};
123 local_to_global_index_map_single_component_ =
124 std::make_unique<NumLib::LocalToGlobalIndexMap>(
125 std::move(all__meshsubsets_single_component),
126 // by location order is needed for output
128
129 // For temperature, which is the first variable.
130 std::vector<MeshLib::MeshSubset> all__meshsubsets{*mesh_subset_base_nodes_};
131
132 // For pressure, which is the second variable
133 all__meshsubsets.push_back(*mesh_subset_base_nodes_);
134
135 // For displacement.
136 const int monolithic_process_id = 0;
137 std::generate_n(std::back_inserter(all__meshsubsets),
138 getProcessVariables(monolithic_process_id)[2]
139 .get()
140 .getNumberOfGlobalComponents(),
141 [&]() { return *_mesh_subset_all_nodes; });
142
143 std::vector<int> const vec_n_components{1, 1, DisplacementDim};
144 _local_to_global_index_map =
145 std::make_unique<NumLib::LocalToGlobalIndexMap>(
146 std::move(all__meshsubsets), vec_n_components,
148 assert(_local_to_global_index_map);
149}
150
151template <int DisplacementDim>
153 NumLib::LocalToGlobalIndexMap const& dof_table,
154 MeshLib::Mesh const& mesh,
155 unsigned const integration_order)
156{
157 createLocalAssemblersHM<DisplacementDim,
159 mesh.getElements(), dof_table, local_assemblers_,
160 mesh.isAxiallySymmetric(), integration_order, process_data_);
161
162 auto add_secondary_variable = [&](std::string const& name,
163 int const num_components,
164 auto get_ip_values_function)
165 {
166 _secondary_variables.addSecondaryVariable(
167 name,
168 makeExtrapolator(num_components, getExtrapolator(),
169 local_assemblers_,
170 std::move(get_ip_values_function)));
171 };
172
173 add_secondary_variable("sigma",
175 DisplacementDim>::RowsAtCompileTime,
176 &LocalAssemblerIF::getIntPtSigma);
177
178 add_secondary_variable("swelling_stress",
180 DisplacementDim>::RowsAtCompileTime,
181 &LocalAssemblerIF::getIntPtSwellingStress);
182
183 add_secondary_variable("epsilon",
185 DisplacementDim>::RowsAtCompileTime,
186 &LocalAssemblerIF::getIntPtEpsilon);
187
188 add_secondary_variable("velocity", DisplacementDim,
189 &LocalAssemblerIF::getIntPtDarcyVelocity);
190
191 add_secondary_variable("saturation", 1,
192 &LocalAssemblerIF::getIntPtSaturation);
193
194 add_secondary_variable("porosity", 1, &LocalAssemblerIF::getIntPtPorosity);
195
196 add_secondary_variable("transport_porosity", 1,
197 &LocalAssemblerIF::getIntPtTransportPorosity);
198
199 add_secondary_variable("dry_density_solid", 1,
200 &LocalAssemblerIF::getIntPtDryDensitySolid);
201
202 add_secondary_variable("liquid_density", 1,
203 &LocalAssemblerIF::getIntPtLiquidDensity);
204
205 add_secondary_variable("viscosity", 1,
206 &LocalAssemblerIF::getIntPtViscosity);
207
208 //
209 // enable output of internal variables defined by material models
210 //
212 LocalAssemblerIF>(process_data_.solid_materials,
213 add_secondary_variable);
214
215 process_data_.element_saturation = MeshLib::getOrCreateMeshProperty<double>(
216 const_cast<MeshLib::Mesh&>(mesh), "saturation_avg",
218
219 process_data_.element_porosity = MeshLib::getOrCreateMeshProperty<double>(
220 const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
222
223 process_data_.element_liquid_density =
224 MeshLib::getOrCreateMeshProperty<double>(
225 const_cast<MeshLib::Mesh&>(mesh), "liquid_density_avg",
227
228 process_data_.element_viscosity = MeshLib::getOrCreateMeshProperty<double>(
229 const_cast<MeshLib::Mesh&>(mesh), "viscosity_avg",
231
232 process_data_.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
233 const_cast<MeshLib::Mesh&>(mesh), "stress_avg",
236 DisplacementDim>::RowsAtCompileTime);
237
238 process_data_.pressure_interpolated =
239 MeshLib::getOrCreateMeshProperty<double>(
240 const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
242 process_data_.temperature_interpolated =
243 MeshLib::getOrCreateMeshProperty<double>(
244 const_cast<MeshLib::Mesh&>(mesh), "temperature_interpolated",
246
247 // Set initial conditions for integration point data.
248 for (auto const& ip_writer : _integration_point_writer)
249 {
250 // Find the mesh property with integration point writer's name.
251 auto const& name = ip_writer->name();
252 if (!mesh.getProperties().existsPropertyVector<double>(name))
253 {
254 continue;
255 }
256 auto const& _meshproperty =
257 *mesh.getProperties().template getPropertyVector<double>(name);
258
259 // The mesh property must be defined on integration points.
260 if (_meshproperty.getMeshItemType() !=
262 {
263 continue;
264 }
265
266 auto const ip_meta_data =
268
269 // Check the number of components.
270 if (ip_meta_data.n_components !=
271 _meshproperty.getNumberOfGlobalComponents())
272 {
273 OGS_FATAL(
274 "Different number of components in meta data ({:d}) than in "
275 "the integration point field data for '{:s}': {:d}.",
276 ip_meta_data.n_components, name,
277 _meshproperty.getNumberOfGlobalComponents());
278 }
279
280 // Now we have a properly named vtk's field data array and the
281 // corresponding meta data.
282 std::size_t position = 0;
283 for (auto& local_asm : local_assemblers_)
284 {
285 std::size_t const integration_points_read =
286 local_asm->setIPDataInitialConditions(
287 name, &_meshproperty[position],
288 ip_meta_data.integration_order);
289 if (integration_points_read == 0)
290 {
291 OGS_FATAL(
292 "No integration points read in the integration point "
293 "initial conditions set function.");
294 }
295 position += integration_points_read * ip_meta_data.n_components;
296 }
297 }
298
299 // Initialize local assemblers after all variables have been set.
300 GlobalExecutor::executeMemberOnDereferenced(&LocalAssemblerIF::initialize,
301 local_assemblers_,
302 *_local_to_global_index_map);
303}
304
305template <int DisplacementDim>
307 DisplacementDim>::initializeBoundaryConditions()
308{
309 const int process_id = 0;
310 initializeProcessBoundaryConditionsAndSourceTerms(
311 *_local_to_global_index_map, process_id);
312}
313
314template <int DisplacementDim>
316 setInitialConditionsConcreteProcess(std::vector<GlobalVector*>& x,
317 double const t,
318 int const process_id)
319{
320 DBUG("SetInitialConditions ThermoRichardsMechanicsProcess.");
321
323 &LocalAssemblerIF::setInitialConditions, local_assemblers_,
324 *_local_to_global_index_map, *x[process_id], t, _use_monolithic_scheme,
325 process_id);
326}
327
328template <int DisplacementDim>
330 const double /*t*/, double const /*dt*/,
331 std::vector<GlobalVector*> const& /*x*/,
332 std::vector<GlobalVector*> const& /*xdot*/, int const /*process_id*/,
333 GlobalMatrix& /*M*/, GlobalMatrix& /*K*/, GlobalVector& /*b*/)
334{
335 OGS_FATAL(
336 "The Picard method or the Newton-Raphson method with numerical "
337 "Jacobian is not implemented for ThermoRichardsMechanics with the full "
338 "monolithic coupling scheme");
339}
340
341template <int DisplacementDim>
343 assembleWithJacobianConcreteProcess(const double t, double const dt,
344 std::vector<GlobalVector*> const& x,
345 std::vector<GlobalVector*> const& xdot,
346 int const process_id, GlobalMatrix& M,
348 GlobalMatrix& Jac)
349{
350 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
351 dof_tables;
352
353 DBUG(
354 "Assemble the Jacobian of ThermoRichardsMechanics for the monolithic "
355 "scheme.");
356 dof_tables.emplace_back(*_local_to_global_index_map);
357
358 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
359
362 local_assemblers_, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
363 process_id, M, K, b, Jac);
364
365 auto copyRhs = [&](int const variable_id, auto& output_vector)
366 {
367 transformVariableFromGlobalVector(b, variable_id, dof_tables[0],
368 output_vector, std::negate<double>());
369 };
370
371 copyRhs(0, *heat_flux_);
372 copyRhs(1, *hydraulic_flow_);
373 copyRhs(2, *nodal_forces_);
374}
375
376template <int DisplacementDim>
378 postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
379 double const t, double const dt,
380 const int process_id)
381{
382 DBUG("PostTimestep ThermoRichardsMechanicsProcess.");
383
384 auto const dof_tables = getDOFTables(x.size());
385
386 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
388 &LocalAssemblerIF::postTimestep, local_assemblers_,
389 pv.getActiveElementIDs(), dof_tables, x, t, dt);
390}
391
392template <int DisplacementDim>
394 computeSecondaryVariableConcrete(const double t, const double dt,
395 std::vector<GlobalVector*> const& x,
396 GlobalVector const& x_dot,
397 int const process_id)
398{
399 DBUG("Compute the secondary variables for ThermoRichardsMechanicsProcess.");
400
401 auto const dof_tables = getDOFTables(x.size());
402
403 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
404
406 &LocalAssemblerIF::computeSecondaryVariable, local_assemblers_,
407 pv.getActiveElementIDs(), dof_tables, t, dt, x, x_dot, process_id);
408}
409
410template <int DisplacementDim>
411std::tuple<NumLib::LocalToGlobalIndexMap*, bool> ThermoRichardsMechanicsProcess<
412 DisplacementDim>::getDOFTableForExtrapolatorData() const
413{
414 const bool manage_storage = false;
415 return std::make_tuple(local_to_global_index_map_single_component_.get(),
416 manage_storage);
417}
418
419template <int DisplacementDim>
422 const int /*process_id*/) const
423{
424 return *_local_to_global_index_map;
425}
426
427template <int DisplacementDim>
428std::vector<NumLib::LocalToGlobalIndexMap const*>
430 int const number_of_processes) const
431{
432 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
433 dof_tables.reserve(number_of_processes);
434 std::generate_n(std::back_inserter(dof_tables), number_of_processes,
435 [&]() { return &getDOFTable(dof_tables.size()); });
436 return dof_tables;
437}
438
441
442} // namespace ThermoRichardsMechanics
443} // 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
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.
Global assembler for the monolithic scheme of the non-isothermal Richards flow coupled with mechanics...
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 computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id) override
ThermoRichardsMechanicsProcess(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, ThermoRichardsMechanicsProcessData< DisplacementDim > &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme)
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
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 postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, double const t, double const dt, const int process_id) override
std::vector< NumLib::LocalToGlobalIndexMap const * > getDOFTables(const int number_of_processes) const
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) override
NumLib::LocalToGlobalIndexMap const & getDOFTable(const int process_id) const 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)
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.
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 > getEpsilon() const =0
virtual std::vector< double > getTransportPorosity() const =0
virtual std::vector< double > getSwellingStress() const =0
virtual std::vector< double > getPorosity() const =0
virtual std::vector< double > getSaturation() const =0
virtual std::vector< double > getSigma() const =0