HydroMechanics non-equilibrium initial states.

Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/nonequilibrium_state_variables/i_nonequilibrium_state_variables.md

@@ -0,0 +1,3 @@
+Non-equilibrium initial state variables which are not explicitly in equilibrium
+with any boundary conditions. They are not considered in the weak form but
+affect the constitutive behaviour of the materials.

Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/nonequilibrium_state_variables/t_pressure.md

@@ -0,0 +1 @@
+Non-equilibrium pressure parameter.

Documentation/ProjectFile/prj/processes/process/HYDRO_MECHANICS/nonequilibrium_state_variables/t_stress.md

@@ -0,0 +1 @@
+Non-equilibrium stress parameter.

ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp

@@ -188,13 +188,45 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
         config.getConfigParameter<double>(
             "reference_temperature", std::numeric_limits<double>::quiet_NaN());
 
+    // Non-equilibrium variables
+    ParameterLib::Parameter<double> const* nonequilibrium_stress = nullptr;
+    ParameterLib::Parameter<double> const* nonequilibrium_pressure = nullptr;
+    const auto& nonequilibrium_state_variables_config =
+        //! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__nonequilibrium_state_variables}
+        config.getConfigSubtreeOptional("nonequilibrium_state_variables");
+    if (nonequilibrium_state_variables_config)
+    {
+        auto const nonequilibrium_stress_parameter_name =
+            nonequilibrium_state_variables_config
+                //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__nonequilibrium_state_variables__stress}
+                ->getConfigParameterOptional<std::string>("stress");
+        if (nonequilibrium_stress_parameter_name)
+        {
+            nonequilibrium_stress = &ParameterLib::findParameter<double>(
+                *nonequilibrium_stress_parameter_name,
+                parameters,
+                MathLib::KelvinVector::KelvinVectorDimensions<
+                    DisplacementDim>::value);
+        }
+        auto const nonequilibrium_pressure_parameter_name =
+            nonequilibrium_state_variables_config
+                //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__nonequilibrium_state_variables__pressure}
+                ->getConfigParameterOptional<std::string>("pressure");
+        if (nonequilibrium_pressure_parameter_name)
+        {
+            nonequilibrium_pressure = &ParameterLib::findParameter<double>(
+                *nonequilibrium_pressure_parameter_name, parameters, 1);
+        }
+    }
+
     HydroMechanicsProcessData<DisplacementDim> process_data{
         materialIDs(mesh),      std::move(solid_constitutive_relations),
         intrinsic_permeability, specific_storage,
         fluid_viscosity,        fluid_density,
         biot_coefficient,       porosity,
         solid_density,          specific_body_force,
-        reference_temperature};
+        reference_temperature,  nonequilibrium_stress,
+        nonequilibrium_pressure};
 
     SecondaryVariableCollection secondary_variables;
 

ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h

@@ -59,26 +59,28 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
             _process_data.solid_materials, _process_data.material_ids,
             e.getID());
 
+    if (_process_data.nonequilibrium_pressure)
+    {
+        // TODO (naumov): Use proper time value to get the parameter value.
+        double const t = std::numeric_limits<double>::quiet_NaN();
+        p_neq = _process_data.nonequilibrium_pressure
+                    ->getNodalValuesOnElement(_element, t)
+                    .template topRows<
+                        ShapeFunctionPressure::MeshElement::n_all_nodes>();
+    }
+
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         _ip_data.emplace_back(solid_material);
         auto& ip_data = _ip_data[ip];
+
+        //
+        // The shape matrices
+        //
         auto const& sm_u = shape_matrices_u[ip];
         _ip_data[ip].integration_weight =
             _integration_method.getWeightedPoint(ip).getWeight() *
             sm_u.integralMeasure * sm_u.detJ;
-
-        // Initialize current time step values
-        static const int kelvin_vector_size =
-            MathLib::KelvinVector::KelvinVectorDimensions<
-                DisplacementDim>::value;
-        ip_data.sigma_eff.setZero(kelvin_vector_size);
-        ip_data.eps.setZero(kelvin_vector_size);
-
-        // Previous time step values are not initialized and are set later.
-        ip_data.eps_prev.resize(kelvin_vector_size);
-        ip_data.sigma_eff_prev.resize(kelvin_vector_size);
-
         ip_data.N_u_op = ShapeMatricesTypeDisplacement::template MatrixType<
             DisplacementDim, displacement_size>::Zero(DisplacementDim,
                                                       displacement_size);
@@ -97,6 +99,35 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         ip_data.dNdx_p = shape_matrices_p[ip].dNdx;
 
         _secondary_data.N_u[ip] = shape_matrices_u[ip].N;
+
+        //
+        // Initialize current time step values
+        //
+        if (_process_data.nonequilibrium_stress)
+        {
+            // Computation of non-equilibrium stress.
+            ParameterLib::SpatialPosition x_position;
+            x_position.setCoordinates(MathLib::Point3d(
+                interpolateCoordinates<ShapeFunctionDisplacement,
+                                       ShapeMatricesTypeDisplacement>(
+                    e, ip_data.N_u)));
+            std::vector<double> sigma_neq_data = (*_process_data
+                                                       .nonequilibrium_stress)(
+                std::numeric_limits<double>::quiet_NaN() /* time independent */,
+                x_position);
+            ip_data.sigma_neq =
+                Eigen::Map<typename BMatricesType::KelvinVectorType>(
+                    sigma_neq_data.data(), KelvinVectorSize, 1);
+        }
+        // Initialization from non-equilibrium sigma, which is zero by
+        // default, or is set to some value.
+        ip_data.sigma_eff = ip_data.sigma_neq;
+
+        ip_data.eps.setZero(KelvinVectorSize);
+
+        // Previous time step values are not initialized and are set later.
+        ip_data.eps_prev.resize(KelvinVectorSize);
+        ip_data.sigma_eff_prev.resize(KelvinVectorSize);
     }
 }
 
@@ -196,6 +227,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
         auto& eps = _ip_data[ip].eps;
         auto const& sigma_eff = _ip_data[ip].sigma_eff;
+        auto const& sigma_neq = _ip_data[ip].sigma_neq;
 
         double const S = _process_data.specific_storage(t, x_position)[0];
         double const K_over_mu =
@@ -225,8 +257,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
         double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
         local_rhs.template segment<displacement_size>(displacement_index)
-            .noalias() -=
-            (B.transpose() * sigma_eff - N_u_op.transpose() * rho * b) * w;
+            .noalias() -= (B.transpose() * (sigma_eff - sigma_neq) -
+                           N_u_op.transpose() * rho * b) *
+                          w;
 
         //
         // displacement equation, pressure part
@@ -269,11 +302,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
     // pressure equation
     local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
-        laplace_p * p + storage_p * p_dot + Kpu * u_dot;
+        laplace_p * (p - p_neq) + storage_p * p_dot + Kpu * u_dot;
 
     // displacement equation
     local_rhs.template segment<displacement_size>(displacement_index)
-        .noalias() += Kup * p;
+        .noalias() += Kup * (p - p_neq);
 }
 
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,

ProcessLib/HydroMechanics/HydroMechanicsFEM.h

@@ -50,6 +50,12 @@ struct IntegrationPointData final
     typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
     typename BMatricesType::KelvinVectorType eps, eps_prev;
 
+    /// Non-equilibrium initial stress.
+    typename BMatricesType::KelvinVectorType sigma_neq =
+        BMatricesType::KelvinVectorType::Zero(
+            MathLib::KelvinVector::KelvinVectorDimensions<
+                DisplacementDim>::value);
+
     typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
     typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
 
@@ -415,6 +421,11 @@ class HydroMechanicsLocalAssembler : public LocalAssemblerInterface
                              ShapeFunctionDisplacement::NPOINTS>;
     std::vector<IpData, Eigen::aligned_allocator<IpData>> _ip_data;
 
+    // Non-equilibrium pressure values of the element, initialized by default
+    // to zero unless non-equilibrium pressure parameter is present.
+    typename ShapeMatricesTypePressure::NodalVectorType p_neq =
+        ShapeMatricesTypePressure::NodalVectorType::Zero(pressure_size);
+
     IntegrationMethod _integration_method;
     MeshLib::Element const& _element;
     bool const _is_axially_symmetric;

ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp

@@ -270,6 +270,61 @@ void HydroMechanicsProcess<DisplacementDim>::initializeBoundaryConditions()
         *_local_to_global_index_map, mechanical_process_id);
 }
 
+template <int DisplacementDim>
+void HydroMechanicsProcess<
+    DisplacementDim>::setInitialConditionsConcreteProcess(GlobalVector& x,
+                                                          double const t)
+{
+    //
+    // Add non-equilibrium pressure to the process variable.
+    //
+    if (_process_data.nonequilibrium_pressure == nullptr)
+    {
+        return;
+    }
+    auto const& p_neq = *_process_data.nonequilibrium_pressure;
+
+    // Get the pressure nodes.
+    if (_use_monolithic_scheme)
+    {
+        constexpr int pressure_variable_id = 0;
+        constexpr int pressure_component_id = 0;
+
+        auto const pressure_mesh_subset =
+            _local_to_global_index_map->getMeshSubset(pressure_variable_id,
+                                                      pressure_component_id);
+        auto const mesh_id = pressure_mesh_subset.getMeshID();
+
+        ParameterLib::SpatialPosition pos;
+        for (auto const* node : pressure_mesh_subset.getNodes())
+        {
+            MeshLib::Location const l(mesh_id, MeshLib::MeshItemType::Node,
+                                      node->getID());
+            pos.setNodeID(node->getID());
+            double const p_neq_value = p_neq(t, pos)[pressure_component_id];
+
+            auto global_index =
+                std::abs(_local_to_global_index_map->getGlobalIndex(
+                    l, pressure_variable_id, pressure_component_id));
+#ifdef USE_PETSC
+            // The global indices of the ghost entries of the global
+            // matrix or the global vectors need to be set as negative
+            // values for equation assembly, however the global indices
+            // start from zero. Therefore, any ghost entry with zero
+            // index is assigned an negative value of the vector size
+            // or the matrix dimension. To assign the initial value for
+            // the ghost entries, the negative indices of the ghost
+            // entries are restored to zero.
+            if (global_index == x.size())
+            {
+                global_index = 0;
+            }
+#endif
+            x.add(global_index, p_neq_value);
+        }
+    }
+}
+
 template <int DisplacementDim>
 void HydroMechanicsProcess<DisplacementDim>::assembleConcreteProcess(
     const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,

ProcessLib/HydroMechanics/HydroMechanicsProcess.h

@@ -71,6 +71,9 @@ class HydroMechanicsProcess final : public Process
 
     void initializeBoundaryConditions() override;
 
+    void setInitialConditionsConcreteProcess(GlobalVector& x,
+                                             double const t) override;
+
     void assembleConcreteProcess(const double t, GlobalVector const& x,
                                  GlobalMatrix& M, GlobalMatrix& K,
                                  GlobalVector& b) override;

ProcessLib/HydroMechanics/HydroMechanicsProcessData.h

@@ -46,7 +46,9 @@ struct HydroMechanicsProcessData
         ParameterLib::Parameter<double> const& solid_density_,
         Eigen::Matrix<double, DisplacementDim, 1>
             specific_body_force_,
-        double const reference_temperature_)
+        double const reference_temperature_,
+        ParameterLib::Parameter<double> const* const nonequilibrium_stress_,
+        ParameterLib::Parameter<double> const* const nonequilibrium_pressure_)
         : material_ids(material_ids_),
           solid_materials{std::move(solid_materials_)},
           intrinsic_permeability(intrinsic_permeability_),
@@ -56,6 +58,8 @@ struct HydroMechanicsProcessData
           biot_coefficient(biot_coefficient_),
           porosity(porosity_),
           solid_density(solid_density_),
+          nonequilibrium_stress(nonequilibrium_stress_),
+          nonequilibrium_pressure(nonequilibrium_pressure_),
           specific_body_force(std::move(specific_body_force_)),
           reference_temperature(reference_temperature_)
     {
@@ -97,6 +101,12 @@ struct HydroMechanicsProcessData
     ParameterLib::Parameter<double> const& porosity;
     /// Solid's density. A scalar quantity, ParameterLib::Parameter<double>.
     ParameterLib::Parameter<double> const& solid_density;
+
+    /// Optional, non-equilibrium stress parameter.
+    ParameterLib::Parameter<double> const* const nonequilibrium_stress;
+    /// Optional, non-equilibrium pressure parameter.
+    ParameterLib::Parameter<double> const* const nonequilibrium_pressure;
+
     /// Specific body forces applied to solid and fluid.
     /// It is usually used to apply gravitational forces.
     /// A vector of displacement dimension's length.

ProcessLib/HydroMechanics/Tests.cmake

@@ -1,3 +1,10 @@
+if (NOT OGS_USE_MPI)
+    OgsTest(PROJECTFILE HydroMechanics/InitialStates/displacement/into_initial_state.prj)
+    OgsTest(PROJECTFILE HydroMechanics/InitialStates/displacement/non_equilibrium_initial_state.prj)
+    OgsTest(PROJECTFILE HydroMechanics/InitialStates/pressure/into_initial_state.prj)
+    OgsTest(PROJECTFILE HydroMechanics/InitialStates/pressure/non_equilibrium_initial_state.prj)
+endif()
+
 # HydroMechanics; Small deformations, linear poroelastic (HML)
 
 ### With monolithic scheme

ProcessLib/Process.h

@@ -178,7 +178,7 @@ class Process
     /// processes. It is called by initialize().
     virtual void initializeBoundaryConditions();
 
-    virtual void setInitialConditionsConcreteProcess(GlobalVector const& /*x*/,
+    virtual void setInitialConditionsConcreteProcess(GlobalVector& /*x*/,
                                                      double const /*t*/)
     {
     }