From 93edccf456470c6c959e79e31b5e304fc5cd578a Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 10 Jan 2023 16:47:06 +0100
Subject: [PATCH] [THM] Fixed a bug in using density dependent viscosity

---
 .../ThermoHydroMechanicsFEM-impl.h                  | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
index e1e695a9d18..910522822eb 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -715,9 +715,6 @@ std::vector<double> const& ThermoHydroMechanicsLocalAssembler<
         // TODO (naumov) Temporary value not used by current material models.
         // Need extension of secondary variables interface.
         double const dt = std::numeric_limits<double>::quiet_NaN();
-        auto const viscosity =
-            liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
-                .template value<double>(vars, x_position, t, dt);
 
         auto const alpha =
             medium
@@ -738,6 +735,13 @@ std::vector<double> const& ThermoHydroMechanicsLocalAssembler<
         vars.equivalent_plastic_strain =
             _ip_data[ip].material_state_variables->getEquivalentPlasticStrain();
 
+        auto const fluid_density =
+            liquid_phase.property(MaterialPropertyLib::PropertyType::density)
+                .template value<double>(vars, x_position, t, dt);
+        vars.density = fluid_density;
+        auto const viscosity =
+            liquid_phase.property(MaterialPropertyLib::PropertyType::viscosity)
+                .template value<double>(vars, x_position, t, dt);
         GlobalDimMatrixType K_over_mu =
             MaterialPropertyLib::formEigenTensor<DisplacementDim>(
                 medium
@@ -745,9 +749,6 @@ std::vector<double> const& ThermoHydroMechanicsLocalAssembler<
                     .value(vars, x_position, t, dt)) /
             viscosity;
 
-        auto const fluid_density =
-            liquid_phase.property(MaterialPropertyLib::PropertyType::density)
-                .template value<double>(vars, x_position, t, dt);
         auto const& b = _process_data.specific_body_force;
 
         auto const K_pT_thermal_osmosis =
-- 
GitLab