From eacfb8958d5db0f25850e4c0a98f84bba311d83d Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Mon, 8 Mar 2021 14:34:28 +0100
Subject: [PATCH] [PL/RM] Assemble. Compute p_SR before usage in rho

solid_grain_pressure used in some density models was
not set and caused an access error in the ExponentialProperty.
---
 .../RichardsMechanics/RichardsMechanicsFEM-impl.h     | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 04bd024fd98..df26fa1c646 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -372,10 +372,6 @@ void RichardsMechanicsLocalAssembler<
         auto const alpha =
             solid_phase.property(MPL::PropertyType::biot_coefficient)
                 .template value<double>(variables, x_position, t, dt);
-        auto const rho_SR =
-            solid_phase.property(MPL::PropertyType::density)
-                .template value<double>(variables, x_position, t, dt);
-
         auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
             t, x_position, dt, temperature);
 
@@ -546,6 +542,13 @@ void RichardsMechanicsLocalAssembler<
         _ip_data[ip].updateConstitutiveRelation(variables, t, x_position, dt,
                                                 temperature);
 
+        // p_SR
+        variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
+            p_FR - sigma_eff.dot(identity2) / (3 * (1 - phi));
+        auto const rho_SR =
+            solid_phase.property(MPL::PropertyType::density)
+                .template value<double>(variables, x_position, t, dt);
+
         //
         // displacement equation, displacement part
         //
-- 
GitLab