Merge branch 'initialStress-PhaseField' into 'master'

Restart simulation capability for phase field process

See merge request ogs/ogs!5071

ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp

@@ -17,6 +17,7 @@
 #include "ParameterLib/Utils.h"
 #include "PhaseFieldProcess.h"
 #include "PhaseFieldProcessData.h"
+#include "ProcessLib/Common/HydroMechanics/CreateInitialStress.h"
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
 
@@ -210,6 +211,10 @@ std::unique_ptr<Process> createPhaseFieldProcess(
             phasefield_model_string.c_str());
     }();
 
+    // Initial stress conditions
+    auto initial_stress = ProcessLib::createInitialStress<DisplacementDim>(
+        config, parameters, mesh);
+
     auto const softening_curve = [&]
     {
         auto const softening_curve_string =
@@ -284,6 +289,7 @@ std::unique_ptr<Process> createPhaseFieldProcess(
         crack_resistance,
         crack_length_scale,
         solid_density,
+        initial_stress,
         specific_body_force,
         pressurized_crack,
         propagating_pressurized_crack,

ProcessLib/PhaseField/LocalAssemblerInterface.h

@@ -23,6 +23,14 @@ struct PhaseFieldLocalAssemblerInterface
     : public ProcessLib::LocalAssemblerInterface,
       public NumLib::ExtrapolatableElement
 {
+    virtual std::size_t setIPDataInitialConditions(
+        std::string_view const name, double const* values,
+        int const integration_order) = 0;
+
+    virtual std::vector<double> getSigma() const = 0;
+
+    virtual std::vector<double> getEpsilon() const = 0;
+
     virtual std::vector<double> const& getIntPtSigma(
         const double t,
         std::vector<GlobalVector*> const& x,

ProcessLib/PhaseField/PhaseFieldFEM-impl.h

@@ -415,5 +415,73 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::computeEnergy(
     surface_energy += element_surface_energy;
     pressure_work += element_pressure_work;
 }
+
+template <typename ShapeFunctionDisplacement, int DisplacementDim>
+std::size_t PhaseFieldLocalAssembler<
+    ShapeFunctionDisplacement,
+    DisplacementDim>::setIPDataInitialConditions(std::string_view const name,
+                                                 double const* values,
+                                                 int const integration_order)
+{
+    if (integration_order !=
+        static_cast<int>(_integration_method.getIntegrationOrder()))
+    {
+        OGS_FATAL(
+            "Setting integration point initial conditions; The integration "
+            "order of the local assembler for element {:d} is different from "
+            "the integration order in the initial condition.",
+            _element.getID());
+    }
+
+    if (name == "sigma")
+    {
+        if (_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.value->name);
+        }
+
+        return ProcessLib::setIntegrationPointKelvinVectorData<DisplacementDim>(
+            values, _ip_data, &IpData::sigma);
+    }
+
+    return 0;
+}
+
+template <typename ShapeFunctionDisplacement, int DisplacementDim>
+std::vector<double> PhaseFieldLocalAssembler<ShapeFunctionDisplacement,
+                                             DisplacementDim>::getSigma() const
+{
+    return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
+        _ip_data, &IpData::sigma);
+}
+
+template <typename ShapeFunctionDisplacement, int DisplacementDim>
+std::vector<double> PhaseFieldLocalAssembler<
+    ShapeFunctionDisplacement, DisplacementDim>::getEpsilon() const
+{
+    auto const kelvin_vector_size =
+        MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+
+    std::vector<double> ip_epsilon_values;
+    auto cache_mat = MathLib::createZeroedMatrix<Eigen::Matrix<
+        double, Eigen::Dynamic, kelvin_vector_size, Eigen::RowMajor>>(
+        ip_epsilon_values, n_integration_points, kelvin_vector_size);
+
+    for (unsigned ip = 0; ip < n_integration_points; ++ip)
+    {
+        auto const& eps = _ip_data[ip].eps;
+        cache_mat.row(ip) =
+            MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps);
+    }
+
+    return ip_epsilon_values;
+}
+
 }  // namespace PhaseField
 }  // namespace ProcessLib

ProcessLib/PhaseField/PhaseFieldFEM.h

@@ -262,6 +262,12 @@ class PhaseFieldLocalAssembler : public PhaseFieldLocalAssemblerInterface
         }
     }
 
+    /// Returns number of read integration points.
+    std::size_t setIPDataInitialConditions(
+        std::string_view const name,
+        double const* values,
+        int const integration_order) override;
+
     void assemble(double const /*t*/, double const /*dt*/,
                   std::vector<double> const& /*local_x*/,
                   std::vector<double> const& /*local_x_prev*/,
@@ -287,7 +293,27 @@ class PhaseFieldLocalAssembler : public PhaseFieldLocalAssemblerInterface
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
-            _ip_data[ip].pushBackState();
+            auto& ip_data = _ip_data[ip];
+
+            ParameterLib::SpatialPosition const x_position{
+                std::nullopt, _element.getID(), ip,
+                MathLib::Point3d(
+                    NumLib::interpolateCoordinates<ShapeFunction,
+                                                   ShapeMatricesType>(
+                        _element, ip_data.N))};
+
+            /// Set initial stress from parameter.
+            if (_process_data.initial_stress.value)
+            {
+                ip_data.sigma =
+                    MathLib::KelvinVector::symmetricTensorToKelvinVector<
+                        DisplacementDim>((*_process_data.initial_stress.value)(
+                        std::numeric_limits<
+                            double>::quiet_NaN() /* time independent */,
+                        x_position));
+            }
+
+            ip_data.pushBackState();
         }
     }
 
@@ -327,6 +353,13 @@ class PhaseFieldLocalAssembler : public PhaseFieldLocalAssemblerInterface
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
     }
 
+    // TODO (naumov) This method is same as getIntPtSigma but for arguments and
+    // the ordering of the cache_mat.
+    // There should be only one.
+    std::vector<double> getSigma() const override;
+
+    std::vector<double> getEpsilon() const override;
+
 private:
     std::vector<double> const& getIntPtSigma(
         const double /*t*/,

ProcessLib/PhaseField/PhaseFieldProcess.cpp

@@ -17,6 +17,7 @@
 #include "PhaseFieldFEM.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/SmallDeformation/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/SetIPDataInitialConditions.h"
 
 namespace ProcessLib
 {
@@ -47,6 +48,13 @@ PhaseFieldProcess<DisplacementDim>::PhaseFieldProcess(
 
     _nodal_forces = MeshLib::getOrCreateMeshProperty<double>(
         mesh, "NodalForces", MeshLib::MeshItemType::Node, DisplacementDim);
+
+    _integration_point_writer.emplace_back(
+        std::make_unique<MeshLib::IntegrationPointWriter>(
+            "sigma_ip",
+            static_cast<int>(mesh.getDimension() == 2 ? 4 : 6) /*n components*/,
+            integration_order, _local_assemblers,
+            &LocalAssemblerInterface::getSigma));
 }
 
 template <int DisplacementDim>
@@ -159,6 +167,9 @@ void PhaseFieldProcess<DisplacementDim>::initializeConcreteProcess(
                          getExtrapolator(), _local_assemblers,
                          &LocalAssemblerInterface::getIntPtEpsilonTensile));
 
+    setIPDataInitialConditions(_integration_point_writer, mesh.getProperties(),
+                               _local_assemblers);
+
     // Initialize local assemblers after all variables have been set.
     GlobalExecutor::executeMemberOnDereferenced(
         &LocalAssemblerInterface::initialize, _local_assemblers,

ProcessLib/PhaseField/PhaseFieldProcessData.h

@@ -16,6 +16,7 @@
 
 #include "MeshLib/PropertyVector.h"
 #include "ParameterLib/Parameter.h"
+#include "ProcessLib/Common/HydroMechanics/InitialStress.h"
 
 namespace MaterialLib
 {
@@ -209,6 +210,10 @@ struct PhaseFieldProcessData
     ParameterLib::Parameter<double> const& crack_resistance;
     ParameterLib::Parameter<double> const& crack_length_scale;
     ParameterLib::Parameter<double> const& solid_density;
+    /// Optional, initial stress field. A symmetric tensor, short vector
+    /// representation of length 4 or 6, ParameterLib::Parameter<double>.
+    InitialStress const initial_stress;
+
     Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
     bool pressurized_crack = false;
     bool propagating_pressurized_crack = false;

ProcessLib/PhaseField/Tests.cmake

@@ -97,6 +97,21 @@ AddTest(
         expected_surfing_ts_20_t_1_000000_0.vtu surfing_ts_20_t_1.000000.vtu phasefield phasefield 1e-6 0
 )
 
+AddTest(
+    NAME PhaseField_2D_surfing_AT1_vd_restart
+    PATH PhaseField/surfing
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS surfing_restart.xml
+    WRAPPER mpirun
+    WRAPPER_ARGS -np 1
+    TESTER vtkdiff
+    REQUIREMENTS OGS_USE_MPI
+    RUNTIME 18
+    DIFF_DATA
+        expected_surfing_ts_20_t_1_000000_0.vtu surfing_restart_ts_10_t_1.000000.vtu displacement displacement 1e-10 0
+        expected_surfing_ts_20_t_1_000000_0.vtu surfing_restart_ts_10_t_1.000000.vtu phasefield phasefield 1e-10 0
+)
+
 AddTest(
     NAME PhaseField_2D_K_regime_HF_2cores
     PATH PhaseField/k_regime_HF

Tests/Data/PhaseField/surfing/surfing.prj

@@ -108,10 +108,6 @@
         </processes>
         <output>
             <variables>
-                <variable>displacement</variable>
-                <variable>phasefield</variable>
-                <variable>sigma</variable>
-                <variable>epsilon</variable>
             </variables>
             <type>VTK</type>
             <prefix>surfing</prefix>

Tests/Data/PhaseField/surfing/surfing_restart.xml

@@ -0,0 +1,18 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProjectDiff base_file="surfing.prj">
+    <!-- Replace the mesh file -->
+    <replace sel="/*/meshes/mesh[1]/text()">surfing_ts_10_t_0.500000.vtu</replace>
+    <!-- Modify the displacement0 parameter -->
+    <replace sel="/*/parameters/parameter[name='displacement0']/type/text()">MeshNode</replace>
+    <remove sel="/*/parameters/parameter[name='displacement0']/values"/>
+    <add sel="/*/parameters/parameter[name='displacement0']">
+        <field_name>displacement</field_name>
+    </add>
+    <!-- Modify the phasefield_ic parameter -->
+    <replace sel="/*/parameters/parameter[name='phasefield_ic']/field_name/text()">phasefield</replace>
+    <!-- Modify the t_initial times:0 to 0.5 -->
+    <replace sel="/*/time_loop/processes/process[1]/time_stepping/t_initial/text()">0.5</replace>
+    <replace sel="/*/time_loop/processes/process[2]/time_stepping/t_initial/text()">0.5</replace>
+    <!-- Modify the output prefix -->
+    <replace sel="/*/time_loop/output/prefix/text()">surfing_restart</replace>
+</OpenGeoSysProjectDiff>