Merge branch 'total_stress0_TH2M' into 'master'

[TH2M] Initial total stress input See merge request ogs/ogs!4930
ufz · Feb 29, 2024 · 7784143 · 7784143
2 parents 9da4d22 + 17de5c8
commit 7784143
Show file tree

Hide file tree

Showing 13 changed files with 567 additions and 19 deletions.
diff --git a/Documentation/ProjectFile/prj/processes/process/TH2M/initial_stress/a_type.md b/Documentation/ProjectFile/prj/processes/process/TH2M/initial_stress/a_type.md
@@ -0,0 +1,2 @@
+An optional attribute to indicate the type of the initial stress, `total` or
+ `effective`. The default value is `effective`.
diff --git a/...ation/ProjectFile/prj/processes/process/TH2M/initial_stress/i_initial_stress.md b/...ation/ProjectFile/prj/processes/process/TH2M/initial_stress/i_initial_stress.md
@@ -0,0 +1,3 @@
+An optional tag to input the parameter of initial stress in porous media. It is
+ used in processes `LIE`, `HydroMechanics`, `ThermoHydroMechanics`,
+ `RichardsMechanics`, `ThermoRichardsMechanics`, and `TH2M`.
diff --git a/Documentation/ProjectFile/prj/processes/process/TH2M/t_initial_stress.md b/Documentation/ProjectFile/prj/processes/process/TH2M/t_initial_stress.md
diff --git a/ProcessLib/TH2M/CreateTH2MProcess.cpp b/ProcessLib/TH2M/CreateTH2MProcess.cpp
@@ -17,6 +17,7 @@
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
 #include "ParameterLib/Utils.h"
+#include "ProcessLib/Common/HydroMechanics/CreateInitialStress.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/TH2M/PhaseTransitionModels/NoPhaseTransition.h"
 #include "ProcessLib/TH2M/PhaseTransitionModels/PhaseTransition.h"
@@ -204,12 +205,8 @@ std::unique_ptr<Process> createTH2MProcess(
     }
 
     // Initial stress conditions
-    auto const initial_stress = ParameterLib::findOptionalTagParameter<double>(
-        //! \ogs_file_param_special{prj__processes__process__TH2M__initial_stress}
-        config, "initial_stress", parameters,
-        // Symmetric tensor size, 4 or 6, not a Kelvin vector.
-        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim),
-        &mesh);
+    auto initial_stress = ProcessLib::createInitialStress<DisplacementDim>(
+        config, parameters, mesh);
 
     auto const mass_lumping =
         //! \ogs_file_param{prj__processes__process__TH2M__mass_lumping}
@@ -232,7 +229,7 @@ std::unique_ptr<Process> createTH2MProcess(
         std::move(solid_constitutive_relations),
         std::move(phase_transition_model),
         reference_temperature,
-        initial_stress,
+        std::move(initial_stress),
         specific_body_force,
         mass_lumping,
         use_TaylorHood_elements};

diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -12,6 +12,7 @@
 
 #include <Eigen/LU>
 
+#include "ConstitutiveRelations/ConstitutiveData.h"
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Property.h"
 #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
@@ -21,6 +22,7 @@
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/Reflection/ReflectionSetIPData.h"
 #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
+#include "TH2MProcessData.h"
 
 namespace ProcessLib
 {
@@ -835,13 +837,13 @@ std::size_t TH2MLocalAssembler<
             this->element_.getID());
     }
 
-    if (name == "sigma" && this->process_data_.initial_stress != nullptr)
+    if (name == "sigma" && this->process_data_.initial_stress.value)
     {
         OGS_FATAL(
             "Setting initial conditions for stress from integration "
             "point data and from a parameter '{:s}' is not possible "
             "simultaneously.",
-            this->process_data_.initial_stress->name);
+            this->process_data_.initial_stress.value->name);
     }
 
     if (name.starts_with("material_state_variable_"))
@@ -899,6 +901,9 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         local_x.template segment<capillary_pressure_size>(
             capillary_pressure_index);
 
+    auto const p_GR =
+        local_x.template segment<gas_pressure_size>(gas_pressure_index);
+
     auto const temperature =
         local_x.template segment<temperature_size>(temperature_index);
 
@@ -956,10 +961,11 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         // Set volumetric strain rate for the general case without swelling.
         vars.volumetric_strain = Invariants::trace(eps);
 
-        this->prev_states_[ip].S_L_data->S_L =
+        double const S_L =
             medium.property(MPL::PropertyType::saturation)
                 .template value<double>(
                     vars, pos, t, std::numeric_limits<double>::quiet_NaN());
+        this->prev_states_[ip].S_L_data->S_L = S_L;
 
         // TODO (naumov) Double computation of C_el might be avoided if
         // updateConstitutiveVariables is called before. But it might interfere
@@ -977,6 +983,24 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate)
                 ? eps + C_el.inverse() * sigma_sw
                 : eps;
+
+        if (this->process_data_.initial_stress.isTotalStress())
+        {
+            auto const alpha_b =
+                medium.property(MPL::PropertyType::biot_coefficient)
+                    .template value<double>(vars, pos, t, 0.0 /*dt*/);
+
+            vars.liquid_saturation = S_L;
+            double const bishop =
+                medium.property(MPL::PropertyType::bishops_effective_stress)
+                    .template value<double>(vars, pos, t, 0.0 /*dt*/);
+
+            this->current_states_[ip].eff_stress_data.sigma.noalias() +=
+                alpha_b * Np.dot(p_GR - bishop * capillary_pressure) *
+                Invariants::identity2;
+            this->prev_states_[ip].eff_stress_data =
+                this->current_states_[ip].eff_stress_data;
+        }
     }
 
     // local_x_prev equal to local_x s.t. the local_x_dot is zero.

diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
@@ -132,14 +132,15 @@ class TH2MLocalAssembler : public LocalAssemblerInterface<DisplacementDim>
                                                                 ip_data.N_u))};
 
             /// Set initial stress from parameter.
-            if (this->process_data_.initial_stress != nullptr)
+            if (this->process_data_.initial_stress.value)
             {
                 this->current_states_[ip].eff_stress_data.sigma.noalias() =
                     MathLib::KelvinVector::symmetricTensorToKelvinVector<
-                        DisplacementDim>((*this->process_data_.initial_stress)(
-                        std::numeric_limits<
-                            double>::quiet_NaN() /* time independent */,
-                        x_position));
+                        DisplacementDim>(
+                        (*this->process_data_.initial_stress.value)(
+                            std::numeric_limits<
+                                double>::quiet_NaN() /* time independent */,
+                            x_position));
             }
 
             double const t = 0;  // TODO (naumov) pass t from top

diff --git a/ProcessLib/TH2M/TH2MProcessData.h b/ProcessLib/TH2M/TH2MProcessData.h
@@ -18,6 +18,7 @@
 #include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "ParameterLib/Parameter.h"
 #include "PhaseTransitionModels/PhaseTransitionModel.h"
+#include "ProcessLib/Common/HydroMechanics/InitialStress.h"
 
 namespace ProcessLib
 {
@@ -39,9 +40,7 @@ struct TH2MProcessData
 
     ParameterLib::Parameter<double> const& reference_temperature;
 
-    /// Optional, initial stress field. A symmetric tensor, short vector
-    /// representation of length 4 or 6, ParameterLib::Parameter<double>.
-    ParameterLib::Parameter<double> const* const initial_stress;
+    InitialStress const initial_stress;
 
     /// Specific body forces applied to solid and both fluid phases.
     /// It is usually used to apply gravitational forces.

diff --git a/ProcessLib/TH2M/Tests.cmake b/ProcessLib/TH2M/Tests.cmake
@@ -435,6 +435,24 @@ AddTest(
     results_heatpipe_strict_ts_23_t_40000_000000_2.vtu results_heatpipe_strict_ts_23_t_40000_000000_2.vtu saturation saturation 1e-9 1e-8
 )
 
+AddTest(
+    NAME TH2MTotalInitialStress
+    PATH TH2M/TotalInitialStress
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS total_initial_stress_H2M.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    RUNTIME 1
+    DIFF_DATA
+    total_initial_stress_H2M_ts_0_t_0.000000.vtu total_initial_stress_H2M_ts_0_t_0.000000.vtu capillary_pressure capillary_pressure 1.0e-10 1.e-10
+    total_initial_stress_H2M_ts_0_t_0.000000.vtu total_initial_stress_H2M_ts_0_t_0.000000.vtu p_GR p_GR 1.0e-10 1.e-10
+    total_initial_stress_H2M_ts_0_t_0.000000.vtu total_initial_stress_H2M_ts_0_t_0.000000.vtu sigma sigma 1.0e-8 1.e-10
+    total_initial_stress_H2M_ts_1_t_1000000.000000.vtu total_initial_stress_H2M_ts_1_t_1000000.000000.vtu capillary_pressure capillary_pressure 1.0e-10 1.e-10
+    total_initial_stress_H2M_ts_1_t_1000000.000000.vtu total_initial_stress_H2M_ts_1_t_1000000.000000.vtu p_GR p_GR 1.e-10 1.e-10
+    total_initial_stress_H2M_ts_1_t_1000000.000000.vtu total_initial_stress_H2M_ts_1_t_1000000.000000.vtu sigma sigma 1.0e-8 1.e-10
+)
+
 if(NOT OGS_USE_PETSC)
     NotebookTest(NOTEBOOKFILE TH2M/TH2/heatpipe/heatpipe.ipynb RUNTIME 140)
 endif()
diff --git a/Tests/Data/TH2M/TotalInitialStress/cube_1x1x1.gml b/Tests/Data/TH2M/TotalInitialStress/cube_1x1x1.gml
@@ -0,0 +1,44 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<?xml-stylesheet type="text/xsl" href="OpenGeoSysGLI.xsl"?>
+
+<OpenGeoSysGLI xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ogs="http://www.opengeosys.org">
+    <name>cube_1x1x1_geometry</name>
+    <points>
+        <point id="0" x="0" y="0" z="0"/>
+        <point id="1" x="0" y="0" z="1"/>
+        <point id="2" x="0" y="1" z="1"/>
+        <point id="3" x="0" y="1" z="0"/>
+
+        <point id="4" x="1" y="0" z="0"/>
+        <point id="5" x="1" y="0" z="1"/>
+        <point id="6" x="1" y="1" z="1"/>
+        <point id="7" x="1" y="1" z="0"/>
+    </points>
+
+    <surfaces>
+        <surface id="0" name="left"><!-- x=0 -->
+            <element p1="0" p2="1" p3="2"/>
+            <element p1="0" p2="3" p3="2"/>
+        </surface>
+        <surface id="1" name="right"><!-- x=1 -->
+            <element p1="4" p2="6" p3="5"/>
+            <element p1="4" p2="6" p3="7"/>
+        </surface>
+        <surface id="2" name="top"><!-- z=1 -->
+            <element p1="1" p2="2" p3="5"/>
+            <element p1="5" p2="2" p3="6"/>
+        </surface>
+        <surface id="3" name="bottom"><!-- z=0 -->
+            <element p1="0" p2="3" p3="4"/>
+            <element p1="4" p2="3" p3="7"/>
+        </surface>
+        <surface id="4" name="front"><!-- y=0 -->
+            <element p1="0" p2="1" p3="4"/>
+            <element p1="4" p2="1" p3="5"/>
+        </surface>
+        <surface id="5" name="back"><!-- y=1 -->
+            <element p1="2" p2="3" p3="6"/>
+            <element p1="6" p2="3" p3="7"/>
+        </surface>
+    </surfaces>
+</OpenGeoSysGLI>
diff --git a/Tests/Data/TH2M/TotalInitialStress/cube_1x1x1_hex_1e0.vtu b/Tests/Data/TH2M/TotalInitialStress/cube_1x1x1_hex_1e0.vtu
@@ -0,0 +1,24 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
+  <UnstructuredGrid>
+    <Piece NumberOfPoints="8"                    NumberOfCells="1"                   >
+      <PointData>
+        <DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="0"                    RangeMax="7"                    offset="0"                   />
+      </PointData>
+      <CellData>
+        <DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="0"                    RangeMax="0"                    offset="96"                  />
+      </CellData>
+      <Points>
+        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.7320508076"         offset="120"                 />
+      </Points>
+      <Cells>
+        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="388"                 />
+        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="484"                 />
+        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="508"                 />
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+  <AppendedData encoding="base64">
+   _QAAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAGAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAAAAAAAAAAAA==wAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADwPwAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAPA/AAAAAAAAAAAAAAAAAADwPwAAAAAAAPA/AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA8D8AAAAAAADwPwAAAAAAAAAAAAAAAAAA8D8AAAAAAAAAAAAAAAAAAPA/AAAAAAAA8D8AAAAAAADwPwAAAAAAAPA/AAAAAAAA8D8=QAAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAwAAAAAAAAACAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAHAAAAAAAAAAYAAAAAAAAACAAAAAAAAAAIAAAAAAAAAA==AQAAAAAAAAAM
+  </AppendedData>
+</VTKFile>