Merge branch 'THM_FixPrevIceVolumeFraction' into 'master'

[PL/THM] Explicitly set phi_fr_prev See merge request ogs/ogs!4744
ufz · Oct 6, 2023 · b8e243a · b8e243a
2 parents d846de4 + 69ca768
commit b8e243a
Show file tree

Hide file tree

Showing 9 changed files with 673 additions and 4 deletions.
diff --git a/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h b/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
@@ -71,15 +71,14 @@ struct IntegrationPointData final
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
 
-    double phi_fr = 0;
-    double phi_fr_prev;
+    double phi_fr = std::numeric_limits<double>::quiet_NaN();
+    double phi_fr_prev = std::numeric_limits<double>::quiet_NaN();
     double integration_weight;
 
     double porosity;
 
     void pushBackState()
     {
-        phi_fr_prev = phi_fr;
         eps0_prev = eps0;
         eps_prev = eps;
         eps_m_prev = eps_m;

diff --git a/ProcessLib/ThermoHydroMechanics/Tests.cmake b/ProcessLib/ThermoHydroMechanics/Tests.cmake
@@ -7,6 +7,7 @@ if (NOT OGS_USE_MPI)
     OgsTest(PROJECTFILE ThermoHydroMechanics/ColumnDeformationFreezing/TM.prj RUNTIME 13)
     OgsTest(PROJECTFILE ThermoHydroMechanics/HeatingHomogeneousDomain/hex_THM.prj RUNTIME 30)
     OgsTest(PROJECTFILE ThermoHydroMechanics/9percentWaterFreezingExpansion/UnitSquare.prj RUNTIME 1)
+    OgsTest(PROJECTFILE ThermoHydroMechanics/Linear/HeatTransportInStationaryFlow/WithFreezingPhase.prj RUNTIME 1)
 endif()
 
 AddTest(

diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -389,7 +389,15 @@ ConstitutiveRelationsValues<DisplacementDim> ThermoHydroMechanicsLocalAssembler<
                     vars, MaterialPropertyLib::Variable::temperature,
                     x_position, t, dt);
 
-        double const& phi_fr_prev = ip_data.phi_fr_prev;
+        double const phi_fr_prev = [&]()
+        {
+            MaterialPropertyLib::VariableArray vars_prev;
+            vars_prev.temperature = T_prev_int_pt;
+            return (*medium)[MaterialPropertyLib::PropertyType::volume_fraction]
+                .template value<double>(vars_prev, x_position, t, dt);
+        }();
+        ip_data.phi_fr_prev = phi_fr_prev;
+
         // alpha_T^I
         MathLib::KelvinVector::KelvinVectorType<
             DisplacementDim> const ice_linear_thermal_expansion_coefficient =