From c85c38becacccd12f14e55643c1ce63b71e09c16 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Tue, 16 Apr 2024 18:57:38 +0200
Subject: [PATCH] [PL/TH2M] Extract gravity model
computes volumetric body force.
---
.../ConstitutiveRelations/ConstitutiveData.h | 2 +
.../ConstitutiveModels.h | 2 +
.../TH2M/ConstitutiveRelations/Gravity.cpp | 39 ++++++++++++++
.../TH2M/ConstitutiveRelations/Gravity.h | 42 +++++++++++++++
ProcessLib/TH2M/TH2MFEM-impl.h | 51 +++++++------------
5 files changed, 103 insertions(+), 33 deletions(-)
create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp
create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/Gravity.h
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 3171263e249..68ebf146b5a 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -20,6 +20,7 @@
#include "Enthalpy.h"
#include "EquivalentPlasticStrainData.h"
#include "FluidDensity.h"
+#include "Gravity.h"
#include "InternalEnergy.h"
#include "MassMoleFractions.h"
#include "MechanicalStrain.h"
@@ -157,6 +158,7 @@ struct ConstitutiveTempData
ThermalConductivityData<DisplacementDim> thermal_conductivity_data;
EffectiveVolumetricEnthalpy effective_volumetric_enthalpy_data;
AdvectionData<DisplacementDim> advection_data;
+ VolumetricBodyForce<DisplacementDim> volumetric_body_force;
FC1Data<DisplacementDim> fC_1;
FC2aData fC_2a;
FC3aData fC_3a;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index 1fb329654da..df7bee1cd8a 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -15,6 +15,7 @@
#include "CEquation.h"
#include "ElasticTangentStiffnessModel.h"
#include "Enthalpy.h"
+#include "Gravity.h"
#include "InternalEnergy.h"
#include "MechanicalStrain.h"
#include "PermeabilityModel.h"
@@ -83,6 +84,7 @@ struct ConstitutiveModels
AdvectionModel<DisplacementDim> advection_model;
InternalEnergyModel internal_energy_model;
EffectiveVolumetricEnthalpyModel effective_volumetric_enthalpy_model;
+ GravityModel<DisplacementDim> gravity_model;
FC1Model<DisplacementDim> fC_1_model;
FC2aModel fC_2a_model;
FC3aModel fC_3a_model;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp
new file mode 100644
index 00000000000..8a4e994a876
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp
@@ -0,0 +1,39 @@
+/**
+ * \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
+ */
+
+#include "Gravity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void GravityModel<DisplacementDim>::eval(
+ FluidDensityData const& fluid_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ VolumetricBodyForce<DisplacementDim>& volumetric_body_force) const
+{
+ auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi;
+ auto const phi_L = S_L_data.S_L * porosity_data.phi;
+ auto const phi_S = 1. - porosity_data.phi;
+
+ auto const rhoGR = fluid_density_data.rho_GR;
+ auto const rhoLR = fluid_density_data.rho_LR;
+ auto const rho =
+ phi_G * rhoGR + phi_L * rhoLR + phi_S * solid_density_data.rho_SR;
+ *volumetric_body_force = rho * specific_body_force();
+}
+
+template struct GravityModel<2>;
+template struct GravityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h
new file mode 100644
index 00000000000..dc1aac245e6
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h
@@ -0,0 +1,42 @@
+/**
+ * \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 "FluidDensity.h"
+#include "Porosity.h"
+#include "Saturation.h"
+#include "SolidDensity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+using VolumetricBodyForce =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct VolumetricBodyForceTag>;
+
+template <int DisplacementDim>
+struct GravityModel
+{
+ void eval(
+ FluidDensityData const& fluid_density_data,
+ PorosityData const& porosity_data,
+ SaturationData const& S_L_data,
+ SolidDensityData const& solid_density_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ VolumetricBodyForce<DisplacementDim>& volumetric_body_force) const;
+};
+
+extern template struct GravityModel<2>;
+extern template struct GravityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 42610ef3881..cc118746435 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -259,6 +259,15 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_cv.viscosity_data,
ip_cv.advection_data);
+ models.gravity_model.eval(
+ ip_out.fluid_density_data,
+ ip_out.porosity_data,
+ current_state.S_L_data,
+ ip_out.solid_density_data,
+ ConstitutiveRelations::SpecificBodyForceData<DisplacementDim>{
+ this->process_data_.specific_body_force},
+ ip_cv.volumetric_body_force);
+
auto const& c = ip_cv.phase_transition_data;
auto const phi_L =
@@ -1099,25 +1108,12 @@ void TH2MLocalAssembler<
double const pCap_prev = Np.dot(capillary_pressure_prev);
auto const s_L = current_state.S_L_data.S_L;
- auto const s_G = 1. - s_L;
auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
auto& alpha_B = ip_cv.biot_data();
auto const& b = this->process_data_.specific_body_force;
- // volume fraction
- auto const phi_G = s_G * ip_out.porosity_data.phi;
- auto const phi_L = s_L * ip_out.porosity_data.phi;
- auto const phi_S = 1. - ip_out.porosity_data.phi;
-
- 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 * ip_out.solid_density_data.rho_SR;
-
// ---------------------------------------------------------------------
// C-component equation
// ---------------------------------------------------------------------
@@ -1214,9 +1210,10 @@ void TH2MLocalAssembler<
KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w;
- fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
- N_u_op(Nu).transpose() * rho * b) *
- w;
+ fU.noalias() -=
+ (BuT * current_state.eff_stress_data.sigma -
+ N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
+ w;
if (this->process_data_.apply_mass_lumping)
{
@@ -1435,26 +1432,13 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
double const pCap_prev = Np.dot(capillary_pressure_prev);
double const T_prev = NT.dot(temperature_prev);
- auto& s_L = current_state.S_L_data.S_L;
- auto const s_G = 1. - s_L;
+ auto const& s_L = current_state.S_L_data.S_L;
auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt;
auto const alpha_B = ip_cv.biot_data();
auto const& b = this->process_data_.specific_body_force;
- // volume fraction
- auto const phi_G = s_G * ip_out.porosity_data.phi;
- auto const phi_L = s_L * ip_out.porosity_data.phi;
- auto const phi_S = 1. - ip_out.porosity_data.phi;
-
- 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 * ip_out.solid_density_data.rho_SR;
-
// ---------------------------------------------------------------------
// C-component equation
// ---------------------------------------------------------------------
@@ -1817,9 +1801,10 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.noalias() += BuT * ip_cd.s_mech_data.stiffness_tensor * Bu * w;
// fU_1
- fU.noalias() -= (BuT * current_state.eff_stress_data.sigma -
- N_u_op(Nu).transpose() * rho * b) *
- w;
+ fU.noalias() -=
+ (BuT * current_state.eff_stress_data.sigma -
+ N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) *
+ w;
// KuT
local_Jac
--
GitLab