Merge branch 'total_sigma_0_RM' into 'master'

[RM] Initial total stress input

See merge request ogs/ogs!5120

Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/initial_stress

+../initial_stress
\ No newline at end of file

Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_initial_stress.md

-\copydoc ProcessLib::RichardsMechanics::RichardsMechanicsProcessData::initial_stress

ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp

@@ -19,6 +19,7 @@
 #include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "NumLib/CreateNewtonRaphsonSolverParameters.h"
 #include "ParameterLib/Utils.h"
+#include "ProcessLib/Common/HydroMechanics/CreateInitialStress.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
 #include "RichardsMechanicsProcess.h"
@@ -166,12 +167,9 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
     DBUG("Media properties verified.");
 
     // Initial stress conditions
-    auto const initial_stress = ParameterLib::findOptionalTagParameter<double>(
-        //! \ogs_file_param_special{prj__processes__process__RICHARDS_MECHANICS__initial_stress}
-        config, "initial_stress", parameters,
-        // Symmetric tensor size, 4 or 6, not a Kelvin vector.
-        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim),
-        &mesh);
+    auto const initial_stress =
+        ProcessLib::createInitialStress<DisplacementDim>(config, parameters,
+                                                         mesh);
 
     std::optional<MicroPorosityParameters> micro_porosity_parameters;
     if (auto const micro_porosity_config =

ProcessLib/RichardsMechanics/LocalAssemblerInterface.h

@@ -76,13 +76,13 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
                 element_.getID());
         }
 
-        if (name == "sigma" && process_data_.initial_stress != nullptr)
+        if (name == "sigma" && process_data_.initial_stress.value != nullptr)
         {
             OGS_FATAL(
                 "Setting initial conditions for stress from integration "
                 "point data and from a parameter '{:s}' is not possible "
                 "simultaneously.",
-                process_data_.initial_stress->name);
+                process_data_.initial_stress.value->name);
         }
 
         if (name.starts_with("material_state_variable_"))

ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h

@@ -219,6 +219,10 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
     auto const& solid_phase = medium->phase("Solid");
 
+    auto const& identity2 = MathLib::KelvinVector::Invariants<
+        MathLib::KelvinVector::kelvin_vector_dimensions(
+            DisplacementDim)>::identity2;
+
     unsigned const n_integration_points =
         this->integration_method_.getNumberOfPoints();
     for (unsigned ip = 0; ip < n_integration_points; ip++)
@@ -257,6 +261,38 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         S_L_prev = medium->property(MPL::PropertyType::saturation)
                        .template value<double>(variables, x_position, t, dt);
 
+        if (this->process_data_.initial_stress.isTotalStress())
+        {
+            auto const alpha_b =
+                medium->property(MPL::PropertyType::biot_coefficient)
+                    .template value<double>(variables, x_position, t, dt);
+
+            variables.liquid_saturation = S_L_prev;
+            double const chi_S_L =
+                medium->property(MPL::PropertyType::bishops_effective_stress)
+                    .template value<double>(variables, x_position, t, dt);
+
+            // Initial stresses are total stress, which were assigned to
+            // sigma_eff in
+            // RichardsMechanicsLocalAssembler::initializeConcrete().
+            auto& sigma_eff =
+                std::get<ProcessLib::ThermoRichardsMechanics::
+                             ConstitutiveStress_StrainTemperature::
+                                 EffectiveStressData<DisplacementDim>>(
+                    this->current_states_[ip]);
+
+            auto& sigma_eff_prev = std::get<
+                PrevState<ProcessLib::ThermoRichardsMechanics::
+                              ConstitutiveStress_StrainTemperature::
+                                  EffectiveStressData<DisplacementDim>>>(
+                this->prev_states_[ip]);
+
+            // Reset sigma_eff to effective stress
+            sigma_eff.sigma_eff.noalias() +=
+                chi_S_L * alpha_b * (-p_cap_ip) * identity2;
+            sigma_eff_prev->sigma_eff = sigma_eff.sigma_eff;
+        }
+
         if (medium->hasProperty(MPL::PropertyType::saturation_micro))
         {
             MPL::VariableArray vars;

ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h

@@ -131,17 +131,20 @@ public:
                                                                 ip_data.N_u))};
 
             /// Set initial stress from parameter.
-            if (this->process_data_.initial_stress != nullptr)
+            if (this->process_data_.initial_stress.value)
             {
                 std::get<ProcessLib::ThermoRichardsMechanics::
                              ConstitutiveStress_StrainTemperature::
                                  EffectiveStressData<DisplacementDim>>(SD)
                     .sigma_eff =
                     MathLib::KelvinVector::symmetricTensorToKelvinVector<
-                        DisplacementDim>((*this->process_data_.initial_stress)(
-                        std::numeric_limits<
-                            double>::quiet_NaN() /* time independent */,
-                        x_position));
+                        DisplacementDim>(
+                        // The data in process_data_.initial_stress.value can
+                        // be total stress or effective stress.
+                        (*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

ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h

@@ -17,6 +17,7 @@
 #include "ComputeMicroPorosity.h"
 #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 #include "ParameterLib/Parameter.h"
+#include "ProcessLib/Common/HydroMechanics/InitialStress.h"
 
 namespace MaterialLib
 {
@@ -42,9 +43,7 @@ struct RichardsMechanicsProcessData
                       MaterialLib::Solids::MechanicsBase<DisplacementDim>>>
         solid_materials;
 
-    /// 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 fluid.
     /// It is usually used to apply gravitational forces.

ProcessLib/RichardsMechanics/Tests.cmake

@@ -68,3 +68,38 @@ AddTest(
     confined_compression_fully_saturated_ts_120_t_1000.000000.vtu confined_compression_fully_saturated_restart_ts_100_t_1000.000000.vtu velocity velocity 1e-16 0
     confined_compression_fully_saturated_ts_420_t_4000.000000.vtu confined_compression_fully_saturated_restart_ts_400_t_4000.000000.vtu velocity velocity 1e-16 0
 )
+
+AddTest(
+    NAME RichardsMechanics_A2_total_initial_stress
+    PATH RichardsMechanics
+    EXECUTABLE ogs
+    RUNTIME 15
+    EXECUTABLE_ARGS A2_total_stress0.xml
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    A2_ts_3_t_4320.000000.vtu A2_total_stess0_test_ts_3_t_4320.000000.vtu displacement displacement 1e-16 0
+    A2_ts_42_t_20736.000000.vtu A2_total_stess0_test_ts_42_t_20736.000000.vtu displacement displacement 1e-16 0
+    A2_ts_76_t_2764800.000000.vtu A2_total_stess0_test_ts_76_t_2764800.000000.vtu displacement displacement 1e-16 0
+
+    A2_ts_3_t_4320.000000.vtu A2_total_stess0_test_ts_3_t_4320.000000.vtu pressure pressure 1e-16 1e-12
+    A2_ts_42_t_20736.000000.vtu A2_total_stess0_test_ts_42_t_20736.000000.vtu pressure pressure 1e-16 1e-12
+    A2_ts_76_t_2764800.000000.vtu A2_total_stess0_test_ts_76_t_2764800.000000.vtu pressure pressure 1e-16 1e-12
+
+    A2_ts_3_t_4320.000000.vtu A2_total_stess0_test_ts_3_t_4320.000000.vtu sigma sigma 5e-8 0
+    A2_ts_42_t_20736.000000.vtu A2_total_stess0_test_ts_42_t_20736.000000.vtu sigma sigma 5e-8 0
+    A2_ts_76_t_2764800.000000.vtu A2_total_stess0_test_ts_76_t_2764800.000000.vtu sigma sigma 5e-8 0
+
+    A2_ts_3_t_4320.000000.vtu A2_total_stess0_test_ts_3_t_4320.000000.vtu epsilon epsilon 5e-14 0
+    A2_ts_42_t_20736.000000.vtu A2_total_stess0_test_ts_42_t_20736.000000.vtu epsilon epsilon 5e-14 0
+    A2_ts_76_t_2764800.000000.vtu A2_total_stess0_test_ts_76_t_2764800.000000.vtu epsilon epsilon 5e-14 0
+
+    A2_ts_3_t_4320.000000.vtu A2_total_stess0_test_ts_3_t_4320.000000.vtu saturation saturation 4e-15 0
+    A2_ts_42_t_20736.000000.vtu A2_total_stess0_test_ts_42_t_20736.000000.vtu saturation saturation 4e-15 0
+    A2_ts_76_t_2764800.000000.vtu A2_total_stess0_test_ts_76_t_2764800.000000.vtu saturation saturation 4e-15 0
+
+    A2_ts_3_t_4320.000000.vtu A2_total_stess0_test_ts_3_t_4320.000000.vtu velocity velocity 1e-16 0
+    A2_ts_42_t_20736.000000.vtu A2_total_stess0_test_ts_42_t_20736.000000.vtu velocity velocity 1e-16 0
+    A2_ts_76_t_2764800.000000.vtu A2_total_stess0_test_ts_76_t_2764800.000000.vtu velocity velocity 1e-16 0
+)

Tests/Data/RichardsMechanics/A2_total_stress0.xml

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="A2.prj">
+    <remove sel="/*/processes/process[1]/initial_stress"/>
+    <add sel="/*/processes/process[1]">
+        <initial_stress type = "total">total_stress0</initial_stress>
+    </add>
+    <add sel="/*/parameters">
+        <parameter>
+            <name>total_stress0</name>
+            <type>Constant</type>
+            <values>-10e6 -10e6 -12e6 0.0 0.0 0.0</values>
+        </parameter>
+    </add>
+    <replace msel="/*/time_loop/output/prefix/text()">
+        A2_total_stess0_test
+    </replace>
+    <remove sel="/*/test_definition"/>
+
+</OpenGeoSysProjectDiff>