Merge branch 'TH2M_removeStrainState' into 'master'

TH2M; remove strain state See merge request ogs/ogs!4835
ufz · Jan 16, 2024 · 8bc2598 · 8bc2598
2 parents e86f565 + 4d55379
commit 8bc2598
Show file tree

Hide file tree

Showing 15 changed files with 165 additions and 162 deletions.
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
@@ -41,18 +41,16 @@ struct IntegrationPointData final
             MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
         sigma_eff.setZero(kelvin_vector_size);
         sigma_sw.setZero(kelvin_vector_size);
-        eps.setZero(kelvin_vector_size);
+        eps.resize(
+            kelvin_vector_size);  // Later initialization from displacement
         eps_m.setZero(kelvin_vector_size);
         eps_m_prev.resize(kelvin_vector_size);
-
-        // Previous time step values are not initialized and are set later.
-        eps_prev.resize(kelvin_vector_size);
         sigma_eff_prev.resize(kelvin_vector_size);
     }
 
     typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
     typename BMatricesType::KelvinVectorType sigma_sw, sigma_sw_prev;
-    typename BMatricesType::KelvinVectorType eps, eps_prev;
+    typename BMatricesType::KelvinVectorType eps;
     typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
 
     typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
@@ -183,7 +181,6 @@ struct IntegrationPointData final
 
     void pushBackState()
     {
-        eps_prev = eps;
         eps_m_prev = eps_m;
         sigma_eff_prev = sigma_eff;
         sigma_sw_prev = sigma_sw;

diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -111,6 +111,8 @@ std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
 
     auto const displacement =
         local_x.template segment<displacement_size>(displacement_index);
+    auto const displacement_prev =
+        local_x_prev.template segment<displacement_size>(displacement_index);
 
     auto const& medium = *_process_data.media_map.getMedium(_element.getID());
     auto const& gas_phase = medium.phase("Gas");
@@ -281,11 +283,11 @@ std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
         MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
             dthermal_strain = ip_data.alpha_T_SR * (T - T_prev);
 
-        auto& eps_prev = ip_data.eps_prev;
         auto& eps_m = ip_data.eps_m;
         auto& eps_m_prev = ip_data.eps_m_prev;
 
-        eps_m.noalias() = eps_m_prev + eps - eps_prev - dthermal_strain;
+        eps_m.noalias() =
+            eps_m_prev + eps - Bu * displacement_prev - dthermal_strain;
 
         if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
         {
@@ -896,7 +898,6 @@ std::size_t TH2MLocalAssembler<
     if (name.starts_with("material_state_variable_") && name.ends_with("_ip"))
     {
         std::string const variable_name = name.substr(24, name.size() - 24 - 3);
-        DBUG("Setting material state variable '{:s}'", variable_name);
 
         // Using first ip data for solid material. TODO (naumov) move solid
         // material into element, store only material state in IPs.
@@ -908,10 +909,16 @@ std::size_t TH2MLocalAssembler<
                              { return iv.name == variable_name; });
             iv != end(internal_variables))
         {
+            DBUG("Setting material state variable '{:s}'", variable_name);
             return ProcessLib::setIntegrationPointDataMaterialStateVariables(
                 values, _ip_data, &IpData::material_state_variables,
                 iv->reference);
         }
+
+        WARN(
+            "Could not find variable {:s} in solid material model's "
+            "internal variables.",
+            variable_name);
     }
     return 0;
 }
@@ -940,6 +947,9 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     auto const temperature =
         local_x.template segment<temperature_size>(temperature_index);
 
+    auto const displacement =
+        local_x.template segment<displacement_size>(displacement_index);
+
     constexpr double dt = std::numeric_limits<double>::quiet_NaN();
     auto const& medium = *_process_data.media_map.getMedium(_element.getID());
     auto const& solid_phase = medium.phase("Solid");
@@ -954,6 +964,12 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         auto& ip_data = _ip_data[ip];
         auto const& Np = ip_data.N_p;
         auto const& NT = Np;
+        auto const& Nu = ip_data.N_u;
+        auto const& gradNu = ip_data.dNdx_u;
+        auto const x_coord =
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, Nu);
         ParameterLib::SpatialPosition const pos{
             std::nullopt, _element.getID(), ip,
             MathLib::Point3d(
@@ -967,7 +983,14 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
         double const T = NT.dot(temperature);
         vars.temperature = T;
 
+        auto const Bu =
+            LinearBMatrix::computeBMatrix<DisplacementDim,
+                                          ShapeFunctionDisplacement::NPOINTS,
+                                          typename BMatricesType::BMatrixType>(
+                gradNu, Nu, x_coord, _is_axially_symmetric);
+
         auto& eps = ip_data.eps;
+        eps.noalias() = Bu * displacement;
 
         // Set volumetric strain rate for the general case without swelling.
         vars.volumetric_strain = Invariants::trace(eps);

diff --git a/ProcessLib/TH2M/Tests.cmake b/ProcessLib/TH2M/Tests.cmake
@@ -140,7 +140,7 @@ AddTest(
 
     result_TH2M_M_ts_2_t_2.000000.vtu result_TH2M_M_ts_2_t_2.000000.vtu velocity_liquid velocity_liquid 1e-8 1e-8
 
-    result_TH2M_M_ts_2_t_2.000000.vtu result_TH2M_M_ts_2_t_2.000000.vtu sigma sigma 6.8e-7 1e-8
+    result_TH2M_M_ts_2_t_2.000000.vtu result_TH2M_M_ts_2_t_2.000000.vtu sigma sigma 6.9e-7 1e-8
 
     result_TH2M_M_ts_2_t_2.000000.vtu result_TH2M_M_ts_2_t_2.000000.vtu epsilon epsilon 1e-8 1e-8
 
@@ -269,7 +269,7 @@ AddTest(
 
     result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu velocity_liquid velocity_liquid 1e-8 1e-8
 
-    result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu sigma sigma 3.1e-6 1e-8
+    result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu sigma sigma 3.3e-6 1e-8
 
     result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu result_1d_dirichlet_slab_ts_5_t_100000.000000.vtu epsilon epsilon 1e-8 1e-8
 

diff --git a/Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj b/Tests/Data/TH2M/H2M/OrthotropicSwelling/square.prj
@@ -488,7 +488,7 @@
         <vtkdiff>
             <regex>square_ts_.*.vtu</regex>
             <field>epsilon</field>
-            <absolute_tolerance>1e-15</absolute_tolerance>
+            <absolute_tolerance>2e-15</absolute_tolerance>
             <relative_tolerance>0</relative_tolerance>
         </vtkdiff>
         <vtkdiff>

diff --git a/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_gas.prj b/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_gas.prj
@@ -479,13 +479,13 @@
         <vtkdiff>
             <regex>HM_confined_compression_gas_ts_.*.vtu</regex>
             <field>sigma</field>
-            <absolute_tolerance>1e-15</absolute_tolerance>
+            <absolute_tolerance>7e-14</absolute_tolerance>
             <relative_tolerance>1e-15</relative_tolerance>
         </vtkdiff>
         <vtkdiff>
             <regex>HM_confined_compression_gas_ts_.*.vtu</regex>
             <field>epsilon</field>
-            <absolute_tolerance>1e-15</absolute_tolerance>
+            <absolute_tolerance>5e-15</absolute_tolerance>
             <relative_tolerance>1e-15</relative_tolerance>
         </vtkdiff>
         <vtkdiff>

diff --git a/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_liquid.prj b/Tests/Data/TH2M/HM/Confined_Compression/HM_confined_compression_liquid.prj
@@ -486,13 +486,13 @@
         <vtkdiff>
             <regex>HM_confined_compression_liquid_ts_.*.vtu</regex>
             <field>sigma</field>
-            <absolute_tolerance>1e-15</absolute_tolerance>
+            <absolute_tolerance>7e-14</absolute_tolerance>
             <relative_tolerance>1e-15</relative_tolerance>
         </vtkdiff>
         <vtkdiff>
             <regex>HM_confined_compression_liquid_ts_.*.vtu</regex>
             <field>epsilon</field>
-            <absolute_tolerance>1e-15</absolute_tolerance>
+            <absolute_tolerance>5e-15</absolute_tolerance>
             <relative_tolerance>1e-15</relative_tolerance>
         </vtkdiff>
         <vtkdiff>

diff --git a/Tests/Data/TH2M/M/MultiMaterialEhlers/square_1e1_2_matIDs_t_0.7000.vtu b/Tests/Data/TH2M/M/MultiMaterialEhlers/square_1e1_2_matIDs_t_0.7000.vtu
@@ -4,7 +4,6 @@
     <FieldData>
       <DataArray type="Int8" Name="IntegrationPointMetaData" NumberOfTuples="1002" format="appended" RangeMin="34"                   RangeMax="125"                  offset="0"                   />
       <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="25" format="appended" RangeMin="45"                   RangeMax="121"                  offset="312"                 />
-      <DataArray type="Float64" Name="epsilon_ip" NumberOfComponents="4" NumberOfTuples="48" format="appended" RangeMin="0.00049497474683"     RangeMax="0.00049497474683"     offset="400"                 />
       <DataArray type="Float64" Name="material_state_variable_ElasticStrain_ip" NumberOfComponents="4" NumberOfTuples="24" format="appended" RangeMin="0.00043074710448"     RangeMax="0.00043074710448"     offset="1000"                />
       <DataArray type="Float64" Name="material_state_variable_EquivalentPlasticStrain_ip" NumberOfTuples="24" format="appended" RangeMin="9.19908e-05"          RangeMax="9.19908e-05"          offset="1576"                />
       <DataArray type="Float64" Name="material_state_variable_damage.kappa_d_ip" NumberOfTuples="24" format="appended" RangeMin="0"                    RangeMax="0"                    offset="1728"                />

diff --git a/Tests/Data/TH2M/submesh_residuum_assembly/base.prj b/Tests/Data/TH2M/submesh_residuum_assembly/base.prj
@@ -221,7 +221,7 @@
                 <convergence_criterion>
                     <type>PerComponentDeltaX</type>
                     <norm_type>INFINITY_N</norm_type>
-                    <abstols>1e-3 1e-3 1e-5 1e-6 1.e-6</abstols>
+                    <abstols>1e-10 1e-7 1e-16 1e-16 1e-16</abstols>
                 </convergence_criterion>
                 <time_discretization>
                     <type>BackwardEuler</type>
@@ -342,7 +342,7 @@
         <process_variable>
             <name>displacement</name>
             <components>2</components>
-            <order>1</order>
+            <order>2</order>
             <boundary_conditions>
                 <boundary_condition>
                     <mesh>surface_y_min</mesh>

diff --git a/Tests/Data/TH2M/submesh_residuum_assembly/p_G.xml b/Tests/Data/TH2M/submesh_residuum_assembly/p_G.xml
@@ -105,7 +105,7 @@
             <vtkdiff>
                 <regex>p_G.*_ts_[0-9]*_.*\.vtu</regex>
                 <field>HeatFlowRate</field>
-                <absolute_tolerance>1e-15</absolute_tolerance>
+                <absolute_tolerance>2e-13</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
 

diff --git a/Tests/Data/TH2M/submesh_residuum_assembly/u.xml b/Tests/Data/TH2M/submesh_residuum_assembly/u.xml
@@ -7,7 +7,7 @@
     <add sel="/*/process_variables/process_variable[name=&quot;displacement&quot;]/boundary_conditions">
         <boundary_condition>
             <mesh>surface_x_max</mesh>
-            <type>Dirichlet</type>
+            <type>Neumann</type>
             <component>0</component>
             <parameter>u_right</parameter>
         </boundary_condition>
@@ -16,7 +16,7 @@
         <parameter>
             <name>u_right</name>
             <type>Constant</type>
-            <value>0.01</value>
+            <value>1e5</value>
         </parameter>
         <parameter>
             <name>u_initial</name>
@@ -32,14 +32,10 @@
         <parameter>
             <name>initial_stress</name>
             <type>Function</type>
-            <!-- sigma_ic = -lambda * eps_xx
-                 with lambda = E * nu / (1+nu) / (1-2*nu)
-                 for nu = 0.1, E = 1e9 Pa and eps_xx = 2.5e-3 we get
-                 sigma_ic = -2.84090909090909e5 Pa
-            -->
-            <expression>-2.84090909090909e5</expression>
-            <expression>-2.84090909090909e5</expression>
-            <expression>-2.84090909090909e5</expression>
+            <!-- sigma_ic corresponds to the force at the right boundary -->
+            <expression>1e5</expression>
+            <expression>0</expression>
+            <expression>0</expression>
             <expression>0</expression>
         </parameter>
     </add>
@@ -85,13 +81,13 @@
             <vtkdiff>
                 <regex>u.*_ts_[0-9]*_.*\.vtu</regex>
                 <field>temperature_interpolated</field>
-                <absolute_tolerance>1e-15</absolute_tolerance>
+                <absolute_tolerance>6e-7</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
             <vtkdiff>
                 <regex>u.*_ts_[0-9]*_.*\.vtu</regex>
                 <field>HeatFlowRate</field>
-                <absolute_tolerance>1e-15</absolute_tolerance>
+                <absolute_tolerance>5e-7</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
 
@@ -127,10 +123,16 @@
                 <absolute_tolerance>1e-15</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
+            <vtkdiff>
+                <regex>u.*_ts_[0-9]*_.*\.vtu</regex>
+                <field>sigma</field>
+                <absolute_tolerance>2e-7</absolute_tolerance>
+                <relative_tolerance>0</relative_tolerance>
+            </vtkdiff>
             <vtkdiff>
                 <regex>u.*_ts_[0-9]*_.*\.vtu</regex>
                 <field>NodalForces</field>
-                <absolute_tolerance>1.2e-8</absolute_tolerance>
+                <absolute_tolerance>1.4e-8</absolute_tolerance>
                 <relative_tolerance>0</relative_tolerance>
             </vtkdiff>
         </test_definition>