diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index f0463af81abd76246d937d546c258f1302e8a723..d55b03082197bd04cb94d4f28660c82114a8923b 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -27,6 +27,7 @@
#include "PureLiquidDensity.h"
#include "Saturation.h"
#include "SolidCompressibility.h"
+#include "SolidDensity.h"
#include "SolidMechanics.h"
#include "SolidThermalExpansion.h"
#include "Swelling.h"
@@ -95,6 +96,7 @@ struct OutputData
FluidDensityData fluid_density_data;
VapourPartialPressureData vapour_pressure_data;
PorosityData porosity_data;
+ SolidDensityData solid_density_data;
static auto reflect()
{
@@ -106,7 +108,8 @@ struct OutputData
&Self::mass_mole_fractions_data,
&Self::fluid_density_data,
&Self::vapour_pressure_data,
- &Self::porosity_data);
+ &Self::porosity_data,
+ &Self::solid_density_data);
}
};
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
new file mode 100644
index 0000000000000000000000000000000000000000..399f73851f0158ac2155ede0bd0ef5fc03ddc5a8
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h
@@ -0,0 +1,33 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include "Base.h"
+#include "ProcessLib/Reflection/ReflectionData.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+struct SolidDensityData
+{
+ double rho_SR = nan;
+
+ static auto reflect()
+ {
+ using Self = SolidDensityData;
+ namespace R = ProcessLib::Reflection;
+
+ return std::tuple{
+ R::makeReflectionData("solid_density", &Self::rho_SR)};
+ }
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index a23f6c4512eccf7ff77d751c2570b3802a83a58f..ec9de5fdb91b1fc2a3efbb11e14b5b8828cfd52c 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -36,9 +36,6 @@ struct IntegrationPointData final
typename ShapeMatricesTypePressure::NodalRowVectorType N_p;
typename ShapeMatricesTypePressure::GlobalDimNodalMatrixType dNdx_p;
- // phase intrinsic densities
- double rhoSR = std::numeric_limits<double>::quiet_NaN();
-
// phase enthalpies
double rho_G_h_G = std::numeric_limits<double>::quiet_NaN();
double rho_L_h_L = std::numeric_limits<double>::quiet_NaN();
diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h
index 39c28bffc2f03bf3c0aef41c5b403dff9444b8ec..9e8d86a0863b988bedb2a1a481c85bd38837a952 100644
--- a/ProcessLib/TH2M/LocalAssemblerInterface.h
+++ b/ProcessLib/TH2M/LocalAssemblerInterface.h
@@ -102,12 +102,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
- virtual std::vector<double> const& getIntPtSolidDensity(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
virtual std::vector<double> const& getIntPtEnthalpySolid(
const double t,
std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 334c939363dfd5b166aa3788d2fb5ac40a42ae99..b4a80dc60bd3bd8985c07e63f5d8489f18230553 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -261,11 +261,11 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
vars, pos, t, std::numeric_limits<double>::quiet_NaN());
#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- auto const rhoSR =
+ ip_out.solid_density_data.rho_SR =
rho_ref_SR * (1. - ip_cv.s_therm_exp_data.thermal_volume_strain +
(ip_cv.biot_data() - 1.) * div_u);
#else // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
- auto const rhoSR = rho_ref_SR;
+ ip_out.solid_density_data.rho_SR = rho_ref_SR;
#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION
auto const& c = ip_cv.phase_transition_data;
@@ -291,15 +291,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_data.rho_u_eff = phi_G * ip_out.fluid_density_data.rho_GR * c.uG +
phi_L * ip_out.fluid_density_data.rho_LR * c.uL +
- phi_S * rhoSR * u_S;
+ phi_S * ip_out.solid_density_data.rho_SR * u_S;
ip_data.rho_G_h_G =
phi_G * ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G;
ip_data.rho_L_h_L =
phi_L * ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L;
- ip_data.rho_S_h_S = phi_S * rhoSR * ip_data.h_S;
-
- ip_data.rhoSR = rhoSR;
+ ip_data.rho_S_h_S =
+ phi_S * ip_out.solid_density_data.rho_SR * ip_data.h_S;
// for variable output
auto const xmCL = 1. - ip_out.mass_mole_fractions_data.xmWL;
@@ -365,8 +364,10 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
phi_G * ip_out.fluid_density_data.rho_GR * c.du_G_dT +
phi_L * drho_LR_dT * c.uL +
phi_L * ip_out.fluid_density_data.rho_LR * c.du_L_dT +
- phi_S * drho_SR_dT * u_S + phi_S * rhoSR * cpS +
- ip_cv.porosity_d_data.dphi_S_dT * rhoSR * u_S;
+ phi_S * drho_SR_dT * u_S +
+ phi_S * ip_out.solid_density_data.rho_SR * cpS +
+ ip_cv.porosity_d_data.dphi_S_dT * ip_out.solid_density_data.rho_SR *
+ u_S;
// ds_L_dp_GR = 0;
// ds_G_dp_GR = -ds_L_dp_GR;
@@ -466,8 +467,10 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_out.enthalpy_data.h_L +
phi_L * drho_LR_dT * ip_out.enthalpy_data.h_L +
phi_L * ip_out.fluid_density_data.rho_LR * c.dh_L_dT +
- ip_cv.porosity_d_data.dphi_S_dT * rhoSR * ip_data.h_S +
- phi_S * drho_SR_dT * ip_data.h_S + phi_S * rhoSR * cpS;
+ ip_cv.porosity_d_data.dphi_S_dT * ip_out.solid_density_data.rho_SR *
+ ip_data.h_S +
+ phi_S * drho_SR_dT * ip_data.h_S +
+ phi_S * ip_out.solid_density_data.rho_SR * cpS;
ip_cv.drho_u_eff_dp_GR =
/*(dphi_G_dp_GR = 0) * ip_out.fluid_density_data.rho_GR * c.uG +*/
@@ -1184,14 +1187,12 @@ void TH2MLocalAssembler<
auto const phi_L = s_L * ip_out.porosity_data.phi;
auto const phi_S = 1. - ip_out.porosity_data.phi;
- // solid phase density
- auto& rho_SR = ip.rhoSR;
-
auto const rhoGR = ip_out.fluid_density_data.rho_GR;
auto const rhoLR = ip_out.fluid_density_data.rho_LR;
// effective density
- auto const rho = phi_G * rhoGR + phi_L * rhoLR + phi_S * rho_SR;
+ auto const rho = phi_G * rhoGR + phi_L * rhoLR +
+ phi_S * ip_out.solid_density_data.rho_SR;
// abbreviations
const double rho_C_FR =
@@ -1677,14 +1678,12 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
auto const phi_L = s_L * ip_out.porosity_data.phi;
auto const phi_S = 1. - ip_out.porosity_data.phi;
- // solid phase density
- auto& rho_SR = ip.rhoSR;
-
auto const rhoGR = ip_out.fluid_density_data.rho_GR;
auto const rhoLR = ip_out.fluid_density_data.rho_LR;
// effective density
- auto const rho = phi_G * rhoGR + phi_L * rhoLR + phi_S * rho_SR;
+ auto const rho = phi_G * rhoGR + phi_L * rhoLR +
+ phi_S * ip_out.solid_density_data.rho_SR;
// abbreviations
const double rho_C_FR =
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 701fd58cb14381050e2f36a7e80a9094b350435d..da9643ed1706c890f62afe42056776b839226845 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -235,16 +235,6 @@ private:
// be stored in some container to avoid defining one method for each
// variable.
- virtual std::vector<double> const& getIntPtSolidDensity(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const override
- {
- return ProcessLib::getIntegrationPointScalarData(_ip_data,
- &IpData::rhoSR, cache);
- }
-
virtual std::vector<double> const& getIntPtEnthalpySolid(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index 13744aa42e24915fc184a36c52e8f04bfaaddef9..34d2e0d1db4de1e5c312f543bbb073b56b3420b2 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -186,10 +186,6 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
&LocalAssemblerInterface<
DisplacementDim>::getIntPtDiffusionVelocityLiquidLiquid);
- add_secondary_variable(
- "solid_density", 1,
- &LocalAssemblerInterface<DisplacementDim>::getIntPtSolidDensity);
-
add_secondary_variable(
"enthalpy_solid", 1,
&LocalAssemblerInterface<DisplacementDim>::getIntPtEnthalpySolid);