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