From 277410a1240a6520378ecd7a5a0962c9c0af27d2 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Tue, 26 Sep 2023 11:59:20 +0200
Subject: [PATCH] [PL/THM] Explicitly set phi_fr_prev, don't update

The phi_fr_prev variable initialization was not correct.
Setting it explicitly avoids any non-initialized states.
---
 ProcessLib/ThermoHydroMechanics/IntegrationPointData.h |  5 ++---
 .../ThermoHydroMechanicsFEM-impl.h                     | 10 +++++++++-
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h b/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
index 90d240ccd24..794beff556e 100644
--- 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/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
index 5993811df26..02274a65177 100644
--- 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 =
-- 
GitLab