OGS
CreateHydroMechanicsProcess.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14#include <range/v3/numeric/accumulate.hpp>
15#include <range/v3/view/transform.hpp>
16
30#include "ParameterLib/Utils.h"
32
33namespace ProcessLib
34{
35namespace LIE
36{
37namespace HydroMechanics
38{
39template <int DisplacementDim>
40std::unique_ptr<Process> createHydroMechanicsProcess(
41 std::string const& name,
42 MeshLib::Mesh& mesh,
43 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
44 std::vector<ProcessVariable> const& variables,
45 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
46 std::optional<ParameterLib::CoordinateSystem> const&
47 local_coordinate_system,
48 unsigned const integration_order,
49 BaseLib::ConfigTree const& config,
50 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
51{
53 config.checkConfigParameter("type", "HYDRO_MECHANICS_WITH_LIE");
54 DBUG("Create HydroMechanicsProcess with LIE.");
55 auto const coupling_scheme =
57 config.getConfigParameterOptional<std::string>("coupling_scheme");
58 const bool use_monolithic_scheme =
59 !(coupling_scheme && (*coupling_scheme == "staggered"));
60
63 auto const pv_conf = config.getConfigSubtree("process_variables");
65 auto range =
67 pv_conf.getConfigParameterList<std::string>("process_variable");
68 std::vector<std::reference_wrapper<ProcessVariable>> p_u_process_variables;
69 std::vector<std::reference_wrapper<ProcessVariable>> p_process_variables;
70 std::vector<std::reference_wrapper<ProcessVariable>> u_process_variables;
71 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
72 process_variables;
73 for (std::string const& pv_name : range)
74 {
75 if (pv_name != "pressure" && pv_name != "displacement" &&
76 !pv_name.starts_with("displacement_jump"))
77 {
79 "Found a process variable name '{}'. It should be "
80 "'displacement' or 'displacement_jumpN' or 'pressure'",
81 pv_name);
82 }
83 auto variable = std::find_if(variables.cbegin(), variables.cend(),
84 [&pv_name](ProcessVariable const& v)
85 { return v.getName() == pv_name; });
86
87 if (variable == variables.end())
88 {
90 "Could not find process variable '{:s}' in the provided "
91 "variables list for config tag <{:s}>.",
92 pv_name, "process_variable");
93 }
94 DBUG("Found process variable '{:s}' for config tag <{:s}>.",
95 variable->getName(), "process_variable");
96
97 if (pv_name.find("displacement") != std::string::npos &&
98 variable->getNumberOfGlobalComponents() != DisplacementDim)
99 {
100 OGS_FATAL(
101 "Number of components of the process variable '{:s}' is "
102 "different from the displacement dimension: got {:d}, expected "
103 "{:d}",
104 variable->getName(),
105 variable->getNumberOfGlobalComponents(),
106 DisplacementDim);
107 }
108
109 if (!use_monolithic_scheme)
110 {
111 if (pv_name == "pressure")
112 {
113 p_process_variables.emplace_back(
114 const_cast<ProcessVariable&>(*variable));
115 }
116 else
117 {
118 u_process_variables.emplace_back(
119 const_cast<ProcessVariable&>(*variable));
120 }
121 }
122 else
123 {
124 p_u_process_variables.emplace_back(
125 const_cast<ProcessVariable&>(*variable));
126 }
127 }
128
129 if (p_u_process_variables.size() > 3 || u_process_variables.size() > 2)
130 {
131 OGS_FATAL("Currently only one displacement jump is supported");
132 }
133
134 if (!use_monolithic_scheme)
135 {
136 process_variables.push_back(std::move(p_process_variables));
137 process_variables.push_back(std::move(u_process_variables));
138 }
139 else
140 {
141 process_variables.push_back(std::move(p_u_process_variables));
142 }
143
145 auto solid_constitutive_relations =
147 parameters, local_coordinate_system, materialIDs(mesh), config);
148
149 // Specific body force
150 Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
151 {
152 std::vector<double> const b =
154 config.getConfigParameter<std::vector<double>>(
155 "specific_body_force");
156 if (b.size() != DisplacementDim)
157 {
158 OGS_FATAL(
159 "The size of the specific body force vector does not match the "
160 "displacement dimension. Vector size is {:d}, displacement "
161 "dimension is {:d}",
162 b.size(), DisplacementDim);
163 }
164
165 std::copy_n(b.data(), b.size(), specific_body_force.data());
166 }
167
168 // Fracture constitutive relation.
169 std::unique_ptr<MaterialLib::Fracture::FractureModelBase<DisplacementDim>>
170 fracture_model = nullptr;
171 auto const opt_fracture_model_config =
173 config.getConfigSubtreeOptional("fracture_model");
174 if (opt_fracture_model_config)
175 {
176 auto& fracture_model_config = *opt_fracture_model_config;
177
178 auto const frac_type =
180 fracture_model_config.peekConfigParameter<std::string>("type");
181
182 if (frac_type == "LinearElasticIsotropic")
183 {
184 fracture_model =
186 DisplacementDim>(parameters, fracture_model_config);
187 }
188 else if (frac_type == "Coulomb")
189 {
190 fracture_model =
192 parameters, fracture_model_config);
193 }
194 else if (frac_type == "CohesiveZoneModeI")
195 {
198 fracture_model_config);
199 }
200 else
201 {
202 OGS_FATAL(
203 "Cannot construct fracture constitutive relation of given type "
204 "'{:s}'.",
205 frac_type);
206 }
207 }
208
209 // Fracture properties
210 std::vector<FractureProperty> fracture_properties;
211 for (
212 auto fracture_properties_config :
214 config.getConfigSubtreeList("fracture_properties"))
215 {
216 fracture_properties.emplace_back(
217 fracture_properties.size(),
219 fracture_properties_config.getConfigParameter<int>("material_id"),
222 fracture_properties_config, "initial_aperture", parameters, 1,
223 &mesh));
224 }
225
226 std::size_t const n_var_du =
227 ranges::accumulate(
228 process_variables | ranges::views::transform(ranges::size),
229 std::size_t{0}) -
230 2 /* for pressure and displacement */;
231 if (n_var_du != fracture_properties.size())
232 {
233 OGS_FATAL(
234 "The number of displacement jumps {} and the number of "
235 "<fracture_properties> {} are not consistent.",
236 n_var_du,
237 fracture_properties.size());
238 }
239
240 // initial effective stress in matrix
241 auto& initial_effective_stress = ParameterLib::findParameter<double>(
242 config,
244 "initial_effective_stress", parameters,
246 &mesh);
247 DBUG("Use '{:s}' as initial effective stress parameter.",
248 initial_effective_stress.name);
249
250 // initial effective stress in fracture
251 auto& initial_fracture_effective_stress = ParameterLib::findParameter<
252 double>(
253 config,
255 "initial_fracture_effective_stress", parameters, DisplacementDim,
256 &mesh);
257 DBUG("Use '{:s}' as initial fracture effective stress parameter.",
258 initial_fracture_effective_stress.name);
259
260 // deactivation of matrix elements in flow
261 auto opt_deactivate_matrix_in_flow =
263 config.getConfigParameterOptional<bool>("deactivate_matrix_in_flow");
264 bool const deactivate_matrix_in_flow =
265 opt_deactivate_matrix_in_flow && *opt_deactivate_matrix_in_flow;
266
267 if (deactivate_matrix_in_flow)
268 INFO("Deactivate matrix elements in flow calculation.");
269
271 auto const use_b_bar = config.getConfigParameter<bool>("use_b_bar", false);
272
273 auto media_map =
275
276 std::array const requiredMediumProperties = {
280 std::array const requiredFluidProperties = {MaterialPropertyLib::viscosity,
282 std::array const requiredSolidProperties = {MaterialPropertyLib::density};
283
285
286 for (auto const& medium : media_map.media())
287 {
288 checkRequiredProperties(*medium, requiredMediumProperties);
289 checkRequiredProperties(fluidPhase(*medium), requiredFluidProperties);
290 checkRequiredProperties(medium->phase("Solid"),
291 requiredSolidProperties);
292 }
293
294 // Check whether fracture permeability is given as a scalar value.
295 for (auto const& element_id : mesh.getElements() | MeshLib::views::ids)
296 {
297 media_map.checkElementHasMedium(element_id);
298 auto const& medium = *media_map.getMedium(element_id);
299
300 // For fracture element
301 if (mesh.getElement(element_id)->getDimension() != DisplacementDim)
302 {
305 auto const permeability =
307 .value(variables, x_position, 0.0 /*t*/, 0.0 /*dt*/);
308 if (!std::holds_alternative<double>(permeability))
309 {
310 OGS_FATAL(
311 "The permeability model for the fracture must be "
312 "isotropic, and it must return a scalar value.");
313 }
314 }
315 }
316
318 materialIDs(mesh), std::move(solid_constitutive_relations),
319 std::move(media_map), specific_body_force,
320 std::move(fracture_model), std::move(fracture_properties),
321 initial_effective_stress, initial_fracture_effective_stress,
322 deactivate_matrix_in_flow, use_b_bar};
323
324 SecondaryVariableCollection secondary_variables;
325
326 ProcessLib::createSecondaryVariables(config, secondary_variables);
327
328 return std::make_unique<HydroMechanicsProcess<DisplacementDim>>(
329 std::move(name), mesh, std::move(jacobian_assembler), parameters,
330 integration_order, std::move(process_variables),
331 std::move(process_data), std::move(secondary_variables),
332 use_monolithic_scheme);
333}
334
335template std::unique_ptr<Process> createHydroMechanicsProcess<2>(
336 std::string const& name,
337 MeshLib::Mesh& mesh,
338 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
339 std::vector<ProcessVariable> const& variables,
340 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
341 std::optional<ParameterLib::CoordinateSystem> const&
342 local_coordinate_system,
343 unsigned const integration_order,
344 BaseLib::ConfigTree const& config,
345 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
346
347template std::unique_ptr<Process> createHydroMechanicsProcess<3>(
348 std::string const& name,
349 MeshLib::Mesh& mesh,
350 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
351 std::vector<ProcessVariable> const& variables,
352 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
353 std::optional<ParameterLib::CoordinateSystem> const&
354 local_coordinate_system,
355 unsigned const integration_order,
356 BaseLib::ConfigTree const& config,
357 std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media);
358
359} // namespace HydroMechanics
360} // namespace LIE
361} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
std::optional< ConfigTree > getConfigSubtreeOptional(std::string const &root) const
std::optional< T > getConfigParameterOptional(std::string const &param) const
T getConfigParameter(std::string const &param) const
Range< SubtreeIterator > getConfigSubtreeList(std::string const &root) const
ConfigTree getConfigSubtree(std::string const &root) const
Range< ValueIterator< T > > getConfigParameterList(std::string const &param) const
void checkConfigParameter(std::string const &param, std::string_view const value) const
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
const Element * getElement(std::size_t idx) const
Get the element with the given index.
Definition Mesh.h:96
Handles configuration of several secondary variables from the project file.
std::unique_ptr< FractureModelBase< DisplacementDim > > createCohesiveZoneModeI(std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, BaseLib::ConfigTree const &config)
std::unique_ptr< FractureModelBase< DisplacementDim > > createLinearElasticIsotropic(std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, BaseLib::ConfigTree const &config)
std::unique_ptr< FractureModelBase< DisplacementDim > > createCoulomb(std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, BaseLib::ConfigTree const &config)
std::map< int, std::shared_ptr< MaterialLib::Solids::MechanicsBase< DisplacementDim > > > createConstitutiveRelations(std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, MeshLib::PropertyVector< int > const *const material_ids, BaseLib::ConfigTree const &config)
void checkMPLPhasesForSinglePhaseFlow(MeshLib::Mesh const &mesh, MaterialPropertyLib::MaterialSpatialDistributionMap const &media_map)
MaterialSpatialDistributionMap createMaterialSpatialDistributionMap(std::map< int, std::shared_ptr< Medium > > const &media, MeshLib::Mesh const &mesh)
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
Definition Mesh.h:227
OGS_NO_DANGLING Parameter< ParameterDataType > & findParameter(std::string const &parameter_name, std::vector< std::unique_ptr< ParameterBase > > const &parameters, int const num_components, MeshLib::Mesh const *const mesh=nullptr)
Definition Utils.h:102
template std::unique_ptr< Process > createHydroMechanicsProcess< 2 >(std::string const &name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
std::unique_ptr< Process > createHydroMechanicsProcess(std::string const &name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
template std::unique_ptr< Process > createHydroMechanicsProcess< 3 >(std::string const &name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
void createSecondaryVariables(BaseLib::ConfigTree const &config, SecondaryVariableCollection &secondary_variables)