From 36739a737a2e432bc332ed484e58202701da8cc0 Mon Sep 17 00:00:00 2001
From: Dominik Kern <dominik.kern.ifgt@gmail.com>
Date: Tue, 16 Nov 2021 09:33:19 +0100
Subject: [PATCH] add parameter to config, use it in HM process and document it
---
.../t_coupling_scheme_parameter.md | 1 +
.../CreateHydroMechanicsProcess.cpp | 54 ++++++++++++++++---
.../HydroMechanics/HydroMechanicsFEM-impl.h | 8 ++-
.../HydroMechanicsProcessData.h | 7 +++
.../docs/userguide/basics/conventions.md | 44 +++++++++++++++
web/data/bib_other.yaml | 24 ++++++++-
6 files changed, 127 insertions(+), 11 deletions(-)
create mode 100644 Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_coupling_scheme_parameter.md
diff --git a/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_coupling_scheme_parameter.md b/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_coupling_scheme_parameter.md
new file mode 100644
index 00000000000..0e7606dac27
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/t_coupling_scheme_parameter.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::HydroMechanics::HydroMechanicsProcessData::fixed_stress_stabilization_parameter
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
index 5bf16cb6b0e..746cc7a898e 100644
--- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -61,6 +61,45 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
const bool use_monolithic_scheme =
!(coupling_scheme && (*coupling_scheme == "staggered"));
+ auto coupling_scheme_parameter_optional =
+ //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__coupling_scheme_parameter}
+ config.getConfigParameterOptional<double>("coupling_scheme_parameter");
+ double coupling_scheme_parameter =
+ 0.5; // default value recommended [Mikelic & Wheeler]
+
+ if (coupling_scheme_parameter_optional)
+ {
+ if (use_monolithic_scheme)
+ {
+ WARN(
+ "Monolithic scheme ignores coupling scheme parameter set in "
+ "project file.");
+ }
+ else
+ {
+ coupling_scheme_parameter =
+ coupling_scheme_parameter_optional.value();
+ // optimum not a-priori known, but within certain interval [Storvik
+ // & Nordbotten]
+ double const csp_min = 1.0 / 6.0;
+ double const csp_max = 1.0;
+ if (coupling_scheme_parameter < csp_min ||
+ coupling_scheme_parameter > csp_max)
+ {
+ WARN(
+ "Value of coupling scheme parameter = {:g} is out of "
+ "reasonable range ({:g}, {:g}).",
+ coupling_scheme_parameter, csp_min, csp_max);
+ }
+ }
+ }
+
+ if (!use_monolithic_scheme)
+ {
+ DBUG("Using value {:g} for coupling parameter of staggered scheme.",
+ coupling_scheme_parameter);
+ }
+
/// \section processvariableshm Process Variables
//! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__process_variables}
@@ -186,14 +225,13 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
variable_u->getShapeFunctionOrder();
HydroMechanicsProcessData<DisplacementDim> process_data{
- materialIDs(mesh),
- std::move(media_map),
- std::move(solid_constitutive_relations),
- initial_stress,
- specific_body_force,
- mass_lumping,
- hydraulic_process_id,
- mechanics_related_process_id,
+ materialIDs(mesh), std::move(media_map),
+ std::move(solid_constitutive_relations), initial_stress,
+ specific_body_force, mass_lumping,
+ coupling_scheme_parameter, // this parameter gets its specific meaning
+ // in the process depending on implemented
+ // coupling scheme
+ hydraulic_process_id, mechanics_related_process_id,
use_taylor_hood_elements};
SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index e016db67841..de158b88e22 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -554,6 +554,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const& identity2 = Invariants::identity2;
+ auto const fixed_stress_stabilization_parameter =
+ _process_data.fixed_stress_stabilization_parameter;
+
int const n_integration_points = _integration_method.getNumberOfPoints();
for (int ip = 0; ip < n_integration_points; ip++)
{
@@ -623,8 +626,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const K_over_mu = K / mu;
- // optimal coupling parameter for fixed-stress split [Wheeler]
- auto const beta_FS = 0.5 * alpha_b * alpha_b / K_S;
+ // artificial compressibility to account for coupling
+ auto const beta_FS =
+ fixed_stress_stabilization_parameter * alpha_b * alpha_b / K_S;
laplace.noalias() +=
rho_fr * dNdx_p.transpose() * K_over_mu * dNdx_p * w;
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
index 63ac7cd7e90..1cbdb7b6a8f 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h
@@ -58,6 +58,13 @@ struct HydroMechanicsProcessData
/// If set mass lumping will be applied to the pressure equation.
bool const mass_lumping;
+ /// An optional input to set an algorithmic parameter of the staggered
+ /// scheme. In the HydroMechanics process the fixed-stress split has been
+ /// implemented as staggered scheme with a stabilization parameter to be
+ /// set. For more details see [user guide -
+ /// conventions](https://www.opengeosys.org/docs/userguide/basics/conventions/#fixed-stress-split-for-hydro-mechanical-processes).
+ double const fixed_stress_stabilization_parameter;
+
/// ID of hydraulic process.
int const hydraulic_process_id;
diff --git a/web/content/docs/userguide/basics/conventions.md b/web/content/docs/userguide/basics/conventions.md
index b5bb818dc89..7b2181e8665 100644
--- a/web/content/docs/userguide/basics/conventions.md
+++ b/web/content/docs/userguide/basics/conventions.md
@@ -95,3 +95,47 @@ For 2D, the Kelvin-Vector of the stress tensor looks like $\sigma=(\sigma_{xx},\
For Kelvin mapping also see the [conversion function documentation](https://doxygen.opengeosys.org/d6/dce/namespacemathlib_1_1kelvinvector#ad78b122c10e91732e95181b6c9a92486).
+## Staggered Scheme
+
+A staggered scheme solves coupled problems by alternating on separate physical domains (e.g. thermal and mechanical) in contrast to monolithic schemes which solve all domains simultaneously (e.g. thermomechanical).
+Thus staggered schemes add another level of iterations, however this may pay off since the subproblems are smaller and they enable a finer tuning of the specific solvers.
+
+### Fixed-stress Split for Hydro-mechanical Processes
+
+For hydro-mechanical processes the fixed-stress split has been implemented, since it turned out advantageous [[1]](#1).
+For sake of brevity, we do not describe the scheme itself, but intend to provide guidance for its stabilization parameter.
+On this parameter depends how many coupling iterations are needed and thus how long it takes to obtain a solution.
+The optimal value of this parameter is not a-priori known, only that it lies within a certain interval [[2]](#2)
+\begin{equation}
+\frac{1}{2}\frac{\alpha^2}{K_\mathrm{1D}} \le \beta_\mathrm{FS} \le \frac{\alpha^2}{K_\mathrm{ph}},
+\end{equation}
+where $\alpha$ denotes Biot's coefficient, and $K_\mathrm{1D}$ and $K_\mathrm{ph}$ are the bulk moduli described next.
+By $K_\mathrm{ph}$ we mean the bulk modulus adjusted to the spatial dimension, so in three dimensions it coincides with the drained bulk modulus, whereas in lower dimensions it becomes a constrained bulk modulus (2D plane strain, 1D uniaxial strain).
+Assuming isotropic, linear elasticity we have
+\begin{eqnarray}
+K_\mathrm{3D} &=& \lambda + \frac{2}{3}\mu, \\
+K_\mathrm{2D} &=& \lambda + \frac{2}{2}\mu, \\
+K_\mathrm{1D} &=& \lambda + \frac{2}{1}\mu. \\
+\end{eqnarray}
+OGS sets the stabilization parameter, which corresponds to a coupling compressibility, via the `coupling_scheme_parameter`
+\begin{equation}
+\beta_\mathrm{FS} = p_\mathrm{FS} \frac{\alpha^2}{K_\mathrm{3D}},
+\end{equation}
+by default to $p_\mathrm{FS}=\frac{1}{2}$.
+For isotropic, linear elasticity we provide the interval [[2]](#2) and the recommended value [[3]](#3) in dependence on Poisson's ratio $\nu$ (note $\frac{\lambda}{\mu}=\frac{2\nu}{1-2\nu}$).
+
+| | 2D | 3D |
+| ------ | ------ | ------ |
+| $p_\mathrm{FS}^\mathrm{min}$ | $\frac{1}{6}\frac{1+\nu}{1-\nu}$ | same as 2D |
+| $p_\mathrm{FS}^\mathrm{MW}$ | $\frac{1+\nu}{3}$ | $\frac{1}{2}$ |
+| $p_\mathrm{FS}^\mathrm{max}$ | $\frac{2(1+\nu)}{3}$ | $1$ |
+
+## References
+<a id="1">[1]</a>
+{{< bib "kimtchjua2009" >}}
+
+<a id="2">[2]</a>
+{{< bib "stonor2019" >}}
+
+<a id="3">[3]</a>
+{{< bib "mikwhe2013" >}}
\ No newline at end of file
diff --git a/web/data/bib_other.yaml b/web/data/bib_other.yaml
index d2f8020b813..4e8a1ad3855 100644
--- a/web/data/bib_other.yaml
+++ b/web/data/bib_other.yaml
@@ -175,8 +175,30 @@ entries:
last: Kim
- first: H.A.
last: Tchelepi
- - first: R
+ - first: R.
last: Juanes
+ stonor2019:
+ type: article
+ title: 'On the optimization of the fixed-stress splitting for Biots equations'
+ journal: International Journal for Numerical Methods in Engineering
+ volume: '120'
+ number: 02
+ pages: 179--194
+ year: '2019'
+ note: ''
+ doi: https://doi.org/10.1002/nme.6130
+ url: https://onlinelibrary.wiley.com/doi/full/10.1002/nme.6130
+ author:
+ - first: E.
+ last: Storvik
+ - first: J.W.
+ last: Both
+ - first: K.
+ last: Kumar
+ - first: J.M.
+ last: Nordbotten
+ - first: F.A.
+ last: Radu
mikwhe2013:
type: article
title: 'Convergence of iterative coupling for coupled flow
--
GitLab