From 66842576053db6423399248800ec0fb66a04497c Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 13 Aug 2018 11:01:32 +0200
Subject: [PATCH] [BGRaCreep] Changed the type of parameters of the creep
 model.

---
 MaterialLib/SolidModels/CreateCreepBGRa.cpp   | 27 +++++++-------
 MaterialLib/SolidModels/CreepBGRa.cpp         | 37 +++++++++++++------
 MaterialLib/SolidModels/CreepBGRa.h           | 25 ++++++++-----
 .../CreepAfterExcavation.prj                  | 29 +++++++++++++--
 .../SimpleAxisymmetricCreep.prj               | 29 +++++++++++++--
 ...eAxisymmetricCreepWithAnalyticSolution.prj | 29 +++++++++++++--
 6 files changed, 129 insertions(+), 47 deletions(-)

diff --git a/MaterialLib/SolidModels/CreateCreepBGRa.cpp b/MaterialLib/SolidModels/CreateCreepBGRa.cpp
index 9c63f0ea152..3ba7d751ad2 100644
--- a/MaterialLib/SolidModels/CreateCreepBGRa.cpp
+++ b/MaterialLib/SolidModels/CreateCreepBGRa.cpp
@@ -22,6 +22,7 @@
 #include "BaseLib/Error.h"
 
 #include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/ProcessUtils.h"  // required for findParameter
 
 namespace MaterialLib
 {
@@ -39,27 +40,27 @@ createCreepBGRa(
     config.checkConfigParameter("type", "CreepBGRa");
     DBUG("Create CreepBGRa material");
 
-    // Read elastic data frist.
+    // Read elastic data first.
     const bool skip_type_checking = true;
     auto elastic_data =
         MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
             parameters, config, skip_type_checking);
 
-    //! \ogs_file_param_special{material__solid__constitutive_relation__CreepBGRa__a}
-    const auto A = config.getConfigParameter<double>("a");
-    DBUG("CreepBGRa parameter A=%g", A);
+    auto& A = ProcessLib::findParameter<double>(
+        //! \ogs_file_param_special{material__solid__constitutive_relation__CreepBGRa__a}
+        config, "a", parameters, 1);
 
-    //! \ogs_file_param_special{material__solid__constitutive_relation__CreepBGRa__n}
-    const auto n = config.getConfigParameter<double>("n");
-    DBUG("CreepBGRa parameter n=%g", n);
+    auto& n = ProcessLib::findParameter<double>(
+        //! \ogs_file_param_special{material__solid__constitutive_relation__CreepBGRa__n}
+        config, "n", parameters, 1);
 
-    //! \ogs_file_param_special{material__solid__constitutive_relation__CreepBGRa__sigma0}
-    const auto sigma0 = config.getConfigParameter<double>("sigma0");
-    DBUG("CreepBGRa parameter sigma0=%g", sigma0);
+    auto& sigma0 = ProcessLib::findParameter<double>(
+        //! \ogs_file_param_special{material__solid__constitutive_relation__CreepBGRa__sigma0}
+        config, "sigma0", parameters, 1);
 
-    //! \ogs_file_param_special{material__solid__constitutive_relation__CreepBGRa__q}
-    const auto Q = config.getConfigParameter<double>("q");
-    DBUG("CreepBGRa parameter Q=%g", Q);
+    auto& Q = ProcessLib::findParameter<double>(
+        //! \ogs_file_param_special{material__solid__constitutive_relation__CreepBGRa__q}
+        config, "q", parameters, 1);
 
     auto const nonlinear_solver_parameters =
         createNewtonRaphsonSolverParameters(config);
diff --git a/MaterialLib/SolidModels/CreepBGRa.cpp b/MaterialLib/SolidModels/CreepBGRa.cpp
index 5448fcf4355..d09b7aff663 100644
--- a/MaterialLib/SolidModels/CreepBGRa.cpp
+++ b/MaterialLib/SolidModels/CreepBGRa.cpp
@@ -43,9 +43,16 @@ CreepBGRa<DisplacementDim>::integrateStress(
     KelvinVector sigma_try = sigma_prev + C * (eps - eps_prev);
     ResidualVectorType solution = sigma_try;
 
+    const double A = _a(t, x)[0];
+    const double n = _n(t, x)[0];
+    const double sigma0 = _sigma_f(t, x)[0];
+    const double Q = _q(t, x)[0];
+
+    const double coef = A * std::pow(1.5, 0.5 * (1 + n)) / std::pow(sigma0, n);
+
     const double b =
-        dt * _coef *
-        std::exp(-_q / (MaterialLib::PhysicalConstant::IdealGasConstant * T));
+        dt * coef *
+        std::exp(-Q / (MaterialLib::PhysicalConstant::IdealGasConstant * T));
     auto const& deviatoric_matrix = Invariants::deviatoric_projection;
 
     // In newton_solver.solve(), the Jacobian is calculated first, and then
@@ -65,12 +72,12 @@ CreepBGRa<DisplacementDim>::integrateStress(
         double const norm_s_n1 = Invariants::FrobeniusNorm(s_n1);
         // side effect
         pow_norm_s_n1_n_minus_one_2b_G =
-            2.0 * b * this->_mp.mu(t, x) * std::pow(norm_s_n1, _n - 1);
-        jacobian = KelvinMatrix::Identity() +
-                   pow_norm_s_n1_n_minus_one_2b_G *
-                       (deviatoric_matrix +
-                        ((_n - 1) / (norm_s_n1 * norm_s_n1)) * s_n1 *
-                            s_n1.transpose());
+            2.0 * b * this->_mp.mu(t, x) * std::pow(norm_s_n1, n - 1);
+        jacobian =
+            KelvinMatrix::Identity() +
+            pow_norm_s_n1_n_minus_one_2b_G *
+                (deviatoric_matrix +
+                 ((n - 1) / (norm_s_n1 * norm_s_n1)) * s_n1 * s_n1.transpose());
     };
 
     auto const update_residual = [&](ResidualVectorType& r) {
@@ -110,11 +117,17 @@ double CreepBGRa<DisplacementDim>::getTemperatureRelatedCoefficient(
     double const t, double const dt, ProcessLib::SpatialPosition const& x,
     double const T, double const deviatoric_stress_norm) const
 {
-    return 2.0 * _coef *
-           std::exp(-_q /
+    const double A = _a(t, x)[0];
+    const double n = _n(t, x)[0];
+    const double sigma0 = _sigma_f(t, x)[0];
+    const double Q = _q(t, x)[0];
+
+    const double coef = A * std::pow(1.5, 0.5 * (1 + n)) / std::pow(sigma0, n);
+    return 2.0 * coef *
+           std::exp(-Q /
                     (MaterialLib::PhysicalConstant::IdealGasConstant * T)) *
-           this->_mp.mu(t, x) * std::pow(deviatoric_stress_norm, _n - 1) * dt *
-           _q / (MaterialLib::PhysicalConstant::IdealGasConstant * T * T);
+           this->_mp.mu(t, x) * std::pow(deviatoric_stress_norm, n - 1) * dt *
+           Q / (MaterialLib::PhysicalConstant::IdealGasConstant * T * T);
 }
 
 template class CreepBGRa<2>;
diff --git a/MaterialLib/SolidModels/CreepBGRa.h b/MaterialLib/SolidModels/CreepBGRa.h
index a243d308cc7..ee29d0e23c4 100644
--- a/MaterialLib/SolidModels/CreepBGRa.h
+++ b/MaterialLib/SolidModels/CreepBGRa.h
@@ -18,6 +18,7 @@
 #include "LinearElasticIsotropic.h"
 #include "MathLib/KelvinVector.h"
 #include "NumLib/NewtonRaphson.h"
+#include "ProcessLib/Parameter/Parameter.h"
 
 namespace MaterialLib
 {
@@ -50,6 +51,8 @@ public:
     using KelvinMatrix =
         MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
 
+    using Parameter = ProcessLib::Parameter<double>;
+
     std::unique_ptr<
         typename MechanicsBase<DisplacementDim>::MaterialStateVariables>
     createMaterialStateVariables() override
@@ -61,19 +64,21 @@ public:
     CreepBGRa(
         typename LinearElasticIsotropic<DisplacementDim>::MaterialProperties mp,
         NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters,
-        const double A, const double n, const double sigma0, const double Q)
+        Parameter const& A, Parameter const& n, Parameter const& sigma_f,
+        Parameter const& Q)
         : LinearElasticIsotropic<DisplacementDim>(std::move(mp)),
           _nonlinear_solver_parameters(std::move(nonlinear_solver_parameters)),
+          _a(A),
           _n(n),
-          _coef(A * std::pow(1.5, 0.5 * (1 + _n)) / std::pow(sigma0, _n)),
+          _sigma_f(sigma_f),
           _q(Q)
     {
     }
 
-    boost::optional<std::tuple<KelvinVector,
-                               std::unique_ptr<typename MechanicsBase<
-                                   DisplacementDim>::MaterialStateVariables>,
-                               KelvinMatrix>>
+    boost::optional<
+        std::tuple<KelvinVector, std::unique_ptr<typename MechanicsBase<
+                                     DisplacementDim>::MaterialStateVariables>,
+                   KelvinMatrix>>
     integrateStress(
         double const t, ProcessLib::SpatialPosition const& x, double const dt,
         KelvinVector const& eps_prev, KelvinVector const& eps,
@@ -94,10 +99,10 @@ public:
 private:
     NumLib::NewtonRaphsonSolverParameters const _nonlinear_solver_parameters;
 
-    const double _n;  /// Creep rate exponent n.
-    /// \f$A\left(\frac{3}{2}\right)^{n/2+1}/\sigma_{eff}^n \f$
-    const double _coef;
-    const double _q;  /// Activation energy
+    Parameter const& _a;        /// A parameter determined by experiment.
+    Parameter const& _n;        /// Creep rate exponent n.
+    Parameter const& _sigma_f;  /// A stress scaling factor.
+    Parameter const& _q;        /// Activation energy
 };
 
 extern template class CreepBGRa<2>;
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj
index de8f4d01129..afd7a74da54 100644
--- a/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/CreepAfterExcavation/CreepAfterExcavation.prj
@@ -15,10 +15,10 @@
                 <type>CreepBGRa</type>
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
-                <a>2.0833333333333334e-06</a>
-                <n>4.9</n>
-                <sigma0>1e6</sigma0>
-                <q>54000</q>
+                <a>A</a>
+                <n>n</n>
+                <sigma0>sigma_f</sigma0>
+                <q>Q</q>
                 <nonlinear_solver>
                     <maximum_iterations>1000</maximum_iterations>
                     <error_tolerance>2e-8</error_tolerance>
@@ -109,6 +109,27 @@
     </time_loop>
 
     <parameters>
+        <parameter>
+            <name>A</name>
+            <type>Constant</type>
+            <value>2.0833333333333334e-06</value>
+        </parameter>
+        <parameter>
+            <name>n</name>
+            <type>Constant</type>
+            <value>4.9</value>
+        </parameter>
+        <parameter>
+            <name>sigma_f</name>
+            <type>Constant</type>
+            <value>1.0e+6</value>
+        </parameter>
+        <parameter>
+            <name>Q</name>
+            <type>Constant</type>
+            <value>54000</value>
+        </parameter>
+
         <parameter>
             <name>E</name>
             <type>MeshElement</type>
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/SimpleAxisymmetricCreep/SimpleAxisymmetricCreep.prj b/Tests/Data/ThermoMechanics/CreepBGRa/SimpleAxisymmetricCreep/SimpleAxisymmetricCreep.prj
index 86f983a7484..e7752f3824a 100644
--- a/Tests/Data/ThermoMechanics/CreepBGRa/SimpleAxisymmetricCreep/SimpleAxisymmetricCreep.prj
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/SimpleAxisymmetricCreep/SimpleAxisymmetricCreep.prj
@@ -11,10 +11,10 @@
                 <type>CreepBGRa</type>
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
-                <a>0.18</a>
-                <n>5.0</n>
-                <sigma0>1</sigma0>
-                <q>54000</q>
+                <a>A</a>
+                <n>n</n>
+                <sigma0>sigma_f</sigma0>
+                <q>Q</q>
                 <nonlinear_solver>
                     <maximum_iterations>1000</maximum_iterations>
                     <error_tolerance>1e-8</error_tolerance>
@@ -105,6 +105,27 @@
     </time_loop>
 
     <parameters>
+        <parameter>
+            <name>A</name>
+            <type>Constant</type>
+            <value>0.18</value>
+        </parameter>
+        <parameter>
+            <name>n</name>
+            <type>Constant</type>
+            <value>5.0</value>
+        </parameter>
+        <parameter>
+            <name>sigma_f</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>Q</name>
+            <type>Constant</type>
+            <value>54000</value>
+        </parameter>
+
         <parameter>
             <name>E</name>
             <type>Constant</type>
diff --git a/Tests/Data/ThermoMechanics/CreepBGRa/SimpleAxisymmetricCreep/SimpleAxisymmetricCreepWithAnalyticSolution.prj b/Tests/Data/ThermoMechanics/CreepBGRa/SimpleAxisymmetricCreep/SimpleAxisymmetricCreepWithAnalyticSolution.prj
index cef8ba5bf9b..7c66a11b8bb 100644
--- a/Tests/Data/ThermoMechanics/CreepBGRa/SimpleAxisymmetricCreep/SimpleAxisymmetricCreepWithAnalyticSolution.prj
+++ b/Tests/Data/ThermoMechanics/CreepBGRa/SimpleAxisymmetricCreep/SimpleAxisymmetricCreepWithAnalyticSolution.prj
@@ -11,10 +11,10 @@
                 <type>CreepBGRa</type>
                 <youngs_modulus>E</youngs_modulus>
                 <poissons_ratio>nu</poissons_ratio>
-                <a>0.18</a>
-                <n>5.0</n>
-                <sigma0>1</sigma0>
-                <q>54000</q>
+                <a>A</a>
+                <n>n</n>
+                <sigma0>sigma_f</sigma0>
+                <q>Q</q>
                 <nonlinear_solver>
                     <maximum_iterations>1000</maximum_iterations>
                     <error_tolerance>1e-8</error_tolerance>
@@ -93,6 +93,27 @@
     </time_loop>
 
     <parameters>
+        <parameter>
+            <name>A</name>
+            <type>Constant</type>
+            <value>0.18</value>
+        </parameter>
+        <parameter>
+            <name>n</name>
+            <type>Constant</type>
+            <value>5.0</value>
+        </parameter>
+        <parameter>
+            <name>sigma_f</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>Q</name>
+            <type>Constant</type>
+            <value>54000</value>
+        </parameter>
+
         <parameter>
             <name>E</name>
             <type>Constant</type>
-- 
GitLab