OGS 6.1.0-1699-ge946d4c5f
CreateHydroMechanicsProcess.cpp
Go to the documentation of this file.
1 
11 
12 #include <cassert>
13 
18 
20 #include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
21 
22 #include "HydroMechanicsProcess.h"
24 
25 namespace ProcessLib
26 {
27 namespace LIE
28 {
29 namespace HydroMechanics
30 {
31 template <unsigned GlobalDim>
32 std::unique_ptr<Process> createHydroMechanicsProcess(
33  MeshLib::Mesh& mesh,
34  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
35  std::vector<ProcessVariable> const& variables,
36  std::vector<std::unique_ptr<ParameterBase>> const& parameters,
37  unsigned const integration_order,
38  BaseLib::ConfigTree const& config)
39 {
41  config.checkConfigParameter("type", "HYDRO_MECHANICS_WITH_LIE");
42  DBUG("Create HydroMechanicsProcess with LIE.");
43  auto const staggered_scheme =
45  config.getConfigParameterOptional<std::string>("coupling_scheme");
46  const bool use_monolithic_scheme =
47  !(staggered_scheme && (*staggered_scheme == "staggered"));
48 
49  // Process variables
51  auto const pv_conf = config.getConfigSubtree("process_variables");
52  auto range =
54  pv_conf.getConfigParameterList<std::string>("process_variable");
55  std::vector<std::reference_wrapper<ProcessVariable>> p_u_process_variables;
56  std::vector<std::reference_wrapper<ProcessVariable>> p_process_variables;
57  std::vector<std::reference_wrapper<ProcessVariable>> u_process_variables;
58  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
59  process_variables;
60  for (std::string const& pv_name : range)
61  {
62  if (pv_name != "pressure" && pv_name != "displacement" &&
63  pv_name.find("displacement_jump") != 0)
64  {
65  OGS_FATAL(
66  "Found a process variable name '%s'. It should be "
67  "'displacement' or 'displacement_jumpN' or 'pressure'");
68  }
69  auto variable = std::find_if(variables.cbegin(), variables.cend(),
70  [&pv_name](ProcessVariable const& v) {
71  return v.getName() == pv_name;
72  });
73 
74  if (variable == variables.end())
75  {
76  OGS_FATAL(
77  "Could not find process variable '%s' in the provided "
78  "variables "
79  "list for config tag <%s>.",
80  pv_name.c_str(), "process_variable");
81  }
82  DBUG("Found process variable '%s' for config tag <%s>.",
83  variable->getName().c_str(), "process_variable");
84 
85  if (pv_name.find("displacement") != std::string::npos &&
86  variable->getNumberOfComponents() != GlobalDim)
87  {
88  OGS_FATAL(
89  "Number of components of the process variable '%s' is "
90  "different "
91  "from the displacement dimension: got %d, expected %d",
92  variable->getName().c_str(),
93  variable->getNumberOfComponents(),
94  GlobalDim);
95  }
96 
97  if (!use_monolithic_scheme)
98  {
99  if (pv_name == "pressure")
100  {
101  p_process_variables.emplace_back(
102  const_cast<ProcessVariable&>(*variable));
103  }
104  else
105  {
106  u_process_variables.emplace_back(
107  const_cast<ProcessVariable&>(*variable));
108  }
109  }
110  else
111  {
112  p_u_process_variables.emplace_back(
113  const_cast<ProcessVariable&>(*variable));
114  }
115  }
116 
117  if (p_u_process_variables.size() > 3 || u_process_variables.size() > 2)
118  {
119  OGS_FATAL("Currently only one displacement jump is supported");
120  }
121 
122  if (!use_monolithic_scheme)
123  {
124  process_variables.push_back(std::move(p_process_variables));
125  process_variables.push_back(std::move(u_process_variables));
126  }
127  else
128  {
129  process_variables.push_back(std::move(p_u_process_variables));
130  }
131 
132  auto solid_constitutive_relations =
133  MaterialLib::Solids::createConstitutiveRelations<GlobalDim>(parameters,
134  config);
135 
136  // Intrinsic permeability
137  auto& intrinsic_permeability = findParameter<double>(
138  config,
140  "intrinsic_permeability", parameters, 1);
141 
142  DBUG("Use '%s' as intrinsic permeabiltiy parameter.",
143  intrinsic_permeability.name.c_str());
144 
145  // Storage coefficient
146  auto& specific_storage = findParameter<double>(
147  config,
149  "specific_storage", parameters, 1);
150 
151  DBUG("Use '%s' as specific storage parameter.",
152  specific_storage.name.c_str());
153 
154  // Fluid viscosity
155  auto& fluid_viscosity = findParameter<double>(
156  config,
158  "fluid_viscosity", parameters, 1);
159  DBUG("Use '%s' as fluid viscosity parameter.",
160  fluid_viscosity.name.c_str());
161 
162  // Fluid density
163  auto& fluid_density = findParameter<double>(
164  config,
166  "fluid_density", parameters, 1);
167  DBUG("Use '%s' as fluid density parameter.", fluid_density.name.c_str());
168 
169  // Biot coefficient
170  auto& biot_coefficient = findParameter<double>(
171  config,
173  "biot_coefficient", parameters, 1);
174  DBUG("Use '%s' as Biot coefficient parameter.",
175  biot_coefficient.name.c_str());
176 
177  // Porosity
178  auto& porosity = findParameter<double>(
179  config,
181  "porosity", parameters, 1);
182  DBUG("Use '%s' as porosity parameter.", porosity.name.c_str());
183 
184  // Solid density
185  auto& solid_density = findParameter<double>(
186  config,
188  "solid_density", parameters, 1);
189  DBUG("Use '%s' as solid density parameter.", solid_density.name.c_str());
190 
191  // Specific body force
192  Eigen::Matrix<double, GlobalDim, 1> specific_body_force;
193  {
194  std::vector<double> const b =
196  config.getConfigParameter<std::vector<double>>(
197  "specific_body_force");
198  if (b.size() != GlobalDim)
199  {
200  OGS_FATAL(
201  "The size of the specific body force vector does not match the "
202  "displacement dimension. Vector size is %d, displacement "
203  "dimension is %d",
204  b.size(), GlobalDim);
205  }
206 
207  std::copy_n(b.data(), b.size(), specific_body_force.data());
208  }
209 
210  // Fracture constitutive relation.
211  std::unique_ptr<MaterialLib::Fracture::FractureModelBase<GlobalDim>>
212  fracture_model = nullptr;
213  auto const opt_fracture_model_config =
215  config.getConfigSubtreeOptional("fracture_model");
216  if (opt_fracture_model_config)
217  {
218  auto& fracture_model_config = *opt_fracture_model_config;
219 
220  auto const frac_type =
222  fracture_model_config.peekConfigParameter<std::string>("type");
223 
224  if (frac_type == "LinearElasticIsotropic")
225  {
226  fracture_model =
227  MaterialLib::Fracture::createLinearElasticIsotropic<GlobalDim>(
228  parameters, fracture_model_config);
229  }
230  else if (frac_type == "MohrCoulomb")
231  {
232  fracture_model =
233  MaterialLib::Fracture::createMohrCoulomb<GlobalDim>(
234  parameters, fracture_model_config);
235  }
236  else if (frac_type == "CohesiveZoneModeI")
237  {
238  fracture_model = MaterialLib::Fracture::CohesiveZoneModeI::
239  createCohesiveZoneModeI<GlobalDim>(parameters,
240  fracture_model_config);
241  }
242  else
243  {
244  OGS_FATAL(
245  "Cannot construct fracture constitutive relation of given type "
246  "'%s'.",
247  frac_type.c_str());
248  }
249  }
250 
251  // Fracture properties
252  std::unique_ptr<FracturePropertyHM> frac_prop = nullptr;
253  auto opt_fracture_properties_config =
255  config.getConfigSubtreeOptional("fracture_properties");
256  if (opt_fracture_properties_config)
257  {
258  auto& fracture_properties_config = *opt_fracture_properties_config;
259 
260  frac_prop = std::make_unique<ProcessLib::LIE::FracturePropertyHM>(
261  0 /*fracture_id*/,
263  fracture_properties_config.getConfigParameter<int>("material_id"),
264  ProcessLib::findParameter<double>(
266  fracture_properties_config, "initial_aperture", parameters, 1),
267  ProcessLib::findParameter<double>(
269  fracture_properties_config, "specific_storage", parameters, 1),
270  ProcessLib::findParameter<double>(
272  fracture_properties_config, "biot_coefficient", parameters, 1));
273  if (frac_prop->aperture0.isTimeDependent())
274  {
275  OGS_FATAL(
276  "The initial aperture parameter '%s' must not be "
277  "time-dependent.",
278  frac_prop->aperture0.name.c_str());
279  }
280  }
281 
282  // initial effective stress in matrix
283  auto& initial_effective_stress = findParameter<double>(
284  config,
286  "initial_effective_stress", parameters,
288  DBUG("Use '%s' as initial effective stress parameter.",
289  initial_effective_stress.name.c_str());
290 
291  // initial effective stress in fracture
292  auto& initial_fracture_effective_stress = findParameter<double>(
293  config,
295  "initial_fracture_effective_stress", parameters, GlobalDim);
296  DBUG("Use '%s' as initial fracture effective stress parameter.",
297  initial_fracture_effective_stress.name.c_str());
298 
299  // deactivation of matrix elements in flow
300  auto opt_deactivate_matrix_in_flow =
302  config.getConfigParameterOptional<bool>("deactivate_matrix_in_flow");
303  bool const deactivate_matrix_in_flow =
304  opt_deactivate_matrix_in_flow && *opt_deactivate_matrix_in_flow;
305  ;
306  if (deactivate_matrix_in_flow)
307  INFO("Deactivate matrix elements in flow calculation.");
308 
309  // Reference temperature
310  const auto& reference_temperature =
312  config.getConfigParameter<double>(
313  "reference_temperature", std::numeric_limits<double>::quiet_NaN());
314 
316  materialIDs(mesh),
317  std::move(solid_constitutive_relations),
318  intrinsic_permeability,
319  specific_storage,
323  porosity,
324  solid_density,
325  specific_body_force,
326  std::move(fracture_model),
327  std::move(frac_prop),
328  initial_effective_stress,
329  initial_fracture_effective_stress,
330  deactivate_matrix_in_flow,
332 
333  SecondaryVariableCollection secondary_variables;
334 
335  NumLib::NamedFunctionCaller named_function_caller(
336  {"HydroMechanics_displacement"});
337 
338  ProcessLib::createSecondaryVariables(config, secondary_variables,
339  named_function_caller);
340 
341  return std::make_unique<HydroMechanicsProcess<GlobalDim>>(
342  mesh, std::move(jacobian_assembler), parameters, integration_order,
343  std::move(process_variables), std::move(process_data),
344  std::move(secondary_variables), std::move(named_function_caller),
345  use_monolithic_scheme);
346 }
347 
348 template std::unique_ptr<Process> createHydroMechanicsProcess<2>(
349  MeshLib::Mesh& mesh,
350  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
351  std::vector<ProcessVariable> const& variables,
352  std::vector<std::unique_ptr<ParameterBase>> const& parameters,
353  unsigned const integration_order,
354  BaseLib::ConfigTree const& config);
355 template std::unique_ptr<Process> createHydroMechanicsProcess<3>(
356  MeshLib::Mesh& mesh,
357  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
358  std::vector<ProcessVariable> const& variables,
359  std::vector<std::unique_ptr<ParameterBase>> const& parameters,
360  unsigned const integration_order,
361  BaseLib::ConfigTree const& config);
362 
363 } // namespace HydroMechanics
364 } // namespace LIE
365 } // namespace ProcessLib
void createSecondaryVariables(BaseLib::ConfigTree const &config, SecondaryVariableCollection &secondary_variables, NumLib::NamedFunctionCaller &named_function_caller)
template std::unique_ptr< Process > createHydroMechanicsProcess< 2 >(MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterBase >> const &parameters, unsigned const integration_order, BaseLib::ConfigTree const &config)
Kelvin vector dimensions for given displacement dimension.
Definition: KelvinVector.h:23
double fluid_viscosity(const double p, const double T, const double x)
Builds expression trees of named functions dynamically at runtime.
void checkConfigParameter(std::string const &param, T const &value) const
double fluid_density(const double p, const double T, const double x)
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition: Mesh.cpp:342
template std::unique_ptr< Process > createHydroMechanicsProcess< 3 >(MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterBase >> const &parameters, unsigned const integration_order, BaseLib::ConfigTree const &config)
ConfigTree getConfigSubtree(std::string const &root) const
Definition: ConfigTree.cpp:146
Handles configuration of several secondary variables from the project file.
boost::optional< T > getConfigParameterOptional(std::string const &param) const
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
Range< ValueIterator< T > > getConfigParameterList(std::string const &param) const
std::unique_ptr< Process > createHydroMechanicsProcess(MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterBase >> const &parameters, unsigned const integration_order, BaseLib::ConfigTree const &config)