From 9cbdb4f504adf96f2a6a8800d64d52f507c67b2c Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Fri, 5 Apr 2024 13:26:42 +0200
Subject: [PATCH] [PL/TH2M] Move internal energies to states
This also eliminates the IntegrationPointData::pushBackState().
---
.../ConstitutiveRelations/ConstitutiveData.h | 4 ++++
.../ConstitutiveRelations/InternalEnergy.h | 22 +++++++++++++++++++
ProcessLib/TH2M/IntegrationPointData.h | 6 -----
ProcessLib/TH2M/TH2MFEM-impl.h | 18 ++++++++-------
ProcessLib/TH2M/TH2MFEM.h | 2 --
5 files changed, 36 insertions(+), 16 deletions(-)
create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 00e55b8cbab..f56b9186b2a 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -16,6 +16,7 @@
#include "Enthalpy.h"
#include "EquivalentPlasticStrainData.h"
#include "FluidDensity.h"
+#include "InternalEnergy.h"
#include "MassMoleFractions.h"
#include "MechanicalStrain.h"
#include "PermeabilityData.h"
@@ -52,6 +53,7 @@ struct StatefulData
MechanicalStrainData<DisplacementDim> mechanical_strain_data;
PureLiquidDensityData rho_W_LR;
ConstituentDensityData constituent_density_data;
+ InternalEnergyData internal_energy_data;
static auto reflect()
{
@@ -72,6 +74,7 @@ struct StatefulDataPrev
PrevState<MechanicalStrainData<DisplacementDim>> mechanical_strain_data;
PrevState<PureLiquidDensityData> rho_W_LR;
PrevState<ConstituentDensityData> constituent_density_data;
+ PrevState<InternalEnergyData> internal_energy_data;
StatefulDataPrev<DisplacementDim>& operator=(
StatefulData<DisplacementDim> const& state)
@@ -82,6 +85,7 @@ struct StatefulDataPrev
mechanical_strain_data = state.mechanical_strain_data;
rho_W_LR = state.rho_W_LR;
constituent_density_data = state.constituent_density_data;
+ internal_energy_data = state.internal_energy_data;
return *this;
}
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h
new file mode 100644
index 00000000000..f5f393704bc
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h
@@ -0,0 +1,22 @@
+/**
+ * \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 "BaseLib/StrongType.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+using InternalEnergyData =
+ BaseLib::StrongType<double, struct InternalEnergyTag>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index 8e1e7da6631..5a5c1b0d9e3 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -41,10 +41,6 @@ struct IntegrationPointData final
double rho_L_h_L = std::numeric_limits<double>::quiet_NaN();
double rho_S_h_S = std::numeric_limits<double>::quiet_NaN();
- // internal energies
- double rho_u_eff = std::numeric_limits<double>::quiet_NaN();
- double rho_u_eff_prev = std::numeric_limits<double>::quiet_NaN();
-
GlobalDimVectorType d_CG;
GlobalDimVectorType d_WG;
GlobalDimVectorType d_CL;
@@ -55,8 +51,6 @@ struct IntegrationPointData final
double integration_weight = std::numeric_limits<double>::quiet_NaN();
- void pushBackState() { rho_u_eff_prev = rho_u_eff; }
-
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index d00cee689e1..41abf2db569 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -265,9 +265,10 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_out.enthalpy_data.h_S = ip_cv.solid_heat_capacity_data() * T;
auto const u_S = ip_out.enthalpy_data.h_S;
- 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 * ip_out.solid_density_data.rho_SR * u_S;
+ current_state.internal_energy_data() =
+ phi_G * ip_out.fluid_density_data.rho_GR * c.uG +
+ phi_L * ip_out.fluid_density_data.rho_LR * c.uL +
+ 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;
@@ -898,9 +899,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- auto& ip_data = _ip_data[ip];
- ip_data.pushBackState();
-
this->material_states_[ip].pushBackState();
this->prev_states_[ip] = this->current_states_[ip];
}
@@ -1131,7 +1129,9 @@ void TH2MLocalAssembler<
auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
- auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
+ auto const rho_u_eff_dot = (current_state.internal_energy_data() -
+ **prev_state.internal_energy_data) /
+ dt;
GlobalDimMatrixType const k_over_mu_G =
ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
@@ -1623,7 +1623,9 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
auto const rho_h_eff = ip.rho_G_h_G + ip.rho_L_h_L + ip.rho_S_h_S;
- auto const rho_u_eff_dot = (ip.rho_u_eff - ip.rho_u_eff_prev) / dt;
+ auto const rho_u_eff_dot = (current_state.internal_energy_data() -
+ **prev_state.internal_energy_data) /
+ dt;
GlobalDimMatrixType const k_over_mu_G =
ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index b2936f9dc63..88713a38d24 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -147,7 +147,6 @@ private:
this->solid_material_.initializeInternalStateVariables(
t, x_position, *material_state.material_state_variables);
- ip_data.pushBackState();
material_state.pushBackState();
}
@@ -167,7 +166,6 @@ private:
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
- _ip_data[ip].pushBackState();
this->material_states_[ip].pushBackState();
}
--
GitLab