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