From d7dd37a6dab427a8e43d5ee36df7c3b95e83e02b Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Fri, 5 Apr 2024 16:31:14 +0200
Subject: [PATCH] [PL/TH2M] Move darcy velocity to output data
---
.../ConstitutiveRelations/ConstitutiveData.h | 5 +-
.../ConstitutiveRelations/DarcyVelocity.h | 36 ++++++
ProcessLib/TH2M/IntegrationPointData.h | 3 -
ProcessLib/TH2M/LocalAssemblerInterface.h | 12 --
ProcessLib/TH2M/TH2MFEM-impl.h | 108 ++++++------------
ProcessLib/TH2M/TH2MFEM.h | 15 ---
ProcessLib/TH2M/TH2MProcess.cpp | 7 --
7 files changed, 75 insertions(+), 111 deletions(-)
create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
index 847a62eb2d1..7c4611cea70 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h
@@ -12,6 +12,7 @@
#include "Biot.h"
#include "Bishops.h"
#include "ConstitutiveDensity.h"
+#include "DarcyVelocity.h"
#include "DiffusionVelocity.h"
#include "ElasticTangentStiffnessData.h"
#include "Enthalpy.h"
@@ -105,6 +106,7 @@ struct OutputData
PorosityData porosity_data;
SolidDensityData solid_density_data;
DiffusionVelocityData<DisplacementDim> diffusion_velocity_data;
+ DarcyVelocityData<DisplacementDim> darcy_velocity_data;
static auto reflect()
{
@@ -118,7 +120,8 @@ struct OutputData
&Self::vapour_pressure_data,
&Self::porosity_data,
&Self::solid_density_data,
- &Self::diffusion_velocity_data);
+ &Self::diffusion_velocity_data,
+ &Self::darcy_velocity_data);
}
};
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
new file mode 100644
index 00000000000..84e66d31d02
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
@@ -0,0 +1,36 @@
+/**
+ * \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
+{
+template <int DisplacementDim>
+struct DarcyVelocityData
+{
+ GlobalDimVector<DisplacementDim> w_GS;
+ GlobalDimVector<DisplacementDim> w_LS;
+
+ static auto reflect()
+ {
+ using Self = DarcyVelocityData<DisplacementDim>;
+ namespace R = ProcessLib::Reflection;
+
+ return std::tuple{
+ R::makeReflectionData("velocity_gas", &Self::w_GS),
+ R::makeReflectionData("velocity_liquid", &Self::w_LS)};
+ }
+};
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index 76ff44a96f0..dffc86438d0 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;
- GlobalDimVectorType w_GS;
- GlobalDimVectorType w_LS;
-
double integration_weight = std::numeric_limits<double>::quiet_NaN();
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h
index aafc2bddafa..30e624d0b54 100644
--- a/ProcessLib/TH2M/LocalAssemblerInterface.h
+++ b/ProcessLib/TH2M/LocalAssemblerInterface.h
@@ -66,18 +66,6 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
std::string_view name, double const* values,
int const integration_order) = 0;
- virtual std::vector<double> const& getIntPtDarcyVelocityGas(
- 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& getIntPtDarcyVelocityLiquid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const = 0;
-
// TODO move to NumLib::ExtrapolatableElement
unsigned getNumberOfIntegrationPoints() const
{
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index ff308665c59..2f14daaaeeb 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -428,14 +428,17 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
ip_cv.dS_L_dp_cap() /
ip_cv.viscosity_data.mu_LR;
- ip_data.w_GS = k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
- k_over_mu_G * gradpGR;
- ip_data.w_LS = k_over_mu_L * gradpCap +
- k_over_mu_L * ip_out.fluid_density_data.rho_LR * b -
- k_over_mu_L * gradpGR;
+ ip_out.darcy_velocity_data.w_GS =
+ k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
+ k_over_mu_G * gradpGR;
+ ip_out.darcy_velocity_data.w_LS =
+ k_over_mu_L * gradpCap +
+ k_over_mu_L * ip_out.fluid_density_data.rho_LR * b -
+ k_over_mu_L * gradpGR;
ip_cv.drho_GR_h_w_eff_dp_GR_Npart =
- c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G * ip_data.w_GS +
+ c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G *
+ ip_out.darcy_velocity_data.w_GS +
ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G *
k_over_mu_G * c.drho_GR_dp_GR * b;
ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart =
@@ -445,7 +448,8 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
k_over_mu_L;
ip_cv.drho_LR_h_w_eff_dp_cap_Npart =
- -drho_LR_dp_cap * ip_out.enthalpy_data.h_L * ip_data.w_LS -
+ -drho_LR_dp_cap * ip_out.enthalpy_data.h_L *
+ ip_out.darcy_velocity_data.w_LS -
ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L *
k_over_mu_L * drho_LR_dp_cap * b;
ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart =
@@ -454,10 +458,14 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
k_over_mu_L;
ip_cv.drho_GR_h_w_eff_dT =
- c.drho_GR_dT * ip_out.enthalpy_data.h_G * ip_data.w_GS +
- ip_out.fluid_density_data.rho_GR * c.dh_G_dT * ip_data.w_GS +
- c.drho_LR_dT * ip_out.enthalpy_data.h_L * ip_data.w_LS +
- ip_out.fluid_density_data.rho_LR * c.dh_L_dT * ip_data.w_LS;
+ c.drho_GR_dT * ip_out.enthalpy_data.h_G *
+ ip_out.darcy_velocity_data.w_GS +
+ ip_out.fluid_density_data.rho_GR * c.dh_G_dT *
+ ip_out.darcy_velocity_data.w_GS +
+ c.drho_LR_dT * ip_out.enthalpy_data.h_L *
+ ip_out.darcy_velocity_data.w_LS +
+ ip_out.fluid_density_data.rho_LR * c.dh_L_dT *
+ ip_out.darcy_velocity_data.w_LS;
// TODO (naumov) + k_over_mu_G * drho_GR_dT * b + k_over_mu_L *
// drho_LR_dT * b
@@ -1327,8 +1335,10 @@ void TH2MLocalAssembler<
fT.noalias() -= NTT * rho_u_eff_dot * w;
fT.noalias() += gradNTT *
- (rhoGR * ip_out.enthalpy_data.h_G * ip.w_GS +
- rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) *
+ (rhoGR * ip_out.enthalpy_data.h_G *
+ ip_out.darcy_velocity_data.w_GS +
+ rhoLR * ip_out.enthalpy_data.h_L *
+ ip_out.darcy_velocity_data.w_LS) *
w;
fT.noalias() += gradNTT *
@@ -1340,9 +1350,10 @@ void TH2MLocalAssembler<
ip_out.diffusion_velocity_data.d_WG) *
w;
- fT.noalias() +=
- NTT * (rhoGR * ip.w_GS.transpose() + rhoLR * ip.w_LS.transpose()) *
- b * w;
+ fT.noalias() += NTT *
+ (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() +
+ rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) *
+ b * w;
// ---------------------------------------------------------------------
// - displacement equation
@@ -2057,8 +2068,10 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// fT_2
fT.noalias() += gradNTT *
- (rhoGR * ip_out.enthalpy_data.h_G * ip.w_GS +
- rhoLR * ip_out.enthalpy_data.h_L * ip.w_LS) *
+ (rhoGR * ip_out.enthalpy_data.h_G *
+ ip_out.darcy_velocity_data.w_GS +
+ rhoLR * ip_out.enthalpy_data.h_L *
+ ip_out.darcy_velocity_data.w_LS) *
w;
// dfT_2/dp_GR
@@ -2088,9 +2101,10 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w;
// fT_3
- fT.noalias() +=
- NTT * (rhoGR * ip.w_GS.transpose() + rhoLR * ip.w_LS.transpose()) *
- b * w;
+ fT.noalias() += NTT *
+ (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() +
+ rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) *
+ b * w;
fT.noalias() += gradNTT *
(current_state.constituent_density_data.rho_C_GR *
@@ -2223,58 +2237,6 @@ void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.noalias() += KUpC;
}
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim>
-std::vector<double> const& TH2MLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
- getIntPtDarcyVelocityGas(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- this->integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = _ip_data[ip].w_GS;
- }
-
- return cache;
-}
-
-template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
- int DisplacementDim>
-std::vector<double> const& TH2MLocalAssembler<
- ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
- getIntPtDarcyVelocityLiquid(
- const double /*t*/,
- std::vector<GlobalVector*> const& /*x*/,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
- std::vector<double>& cache) const
-{
- unsigned const n_integration_points =
- this->integration_method_.getNumberOfPoints();
-
- cache.clear();
- auto cache_matrix = MathLib::createZeroedMatrix<Eigen::Matrix<
- double, DisplacementDim, Eigen::Dynamic, Eigen::RowMajor>>(
- cache, DisplacementDim, n_integration_points);
-
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- cache_matrix.col(ip).noalias() = _ip_data[ip].w_LS;
- }
-
- return cache;
-}
-
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
int DisplacementDim>
void TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h
index 288d78f32a6..ca7e25ccbba 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -188,21 +188,6 @@ private:
return Eigen::Map<const Eigen::RowVectorXd>(N_u.data(), N_u.size());
}
- // TODO: Here is some refactoring potential. All secondary variables could
- // be stored in some container to avoid defining one method for each
- // variable.
-
- std::vector<double> const& getIntPtDarcyVelocityGas(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
- std::vector<double> const& getIntPtDarcyVelocityLiquid(
- const double t,
- std::vector<GlobalVector*> const& x,
- std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
- std::vector<double>& cache) const override;
-
std::tuple<
std::vector<ConstitutiveRelations::ConstitutiveData<DisplacementDim>>,
std::vector<
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index 51cd4a76d80..b3289ca9d14 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -163,13 +163,6 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
LocalAssemblerInterface<DisplacementDim>::getReflectionDataForOutput(),
_secondary_variables, getExtrapolator(), local_assemblers_);
- add_secondary_variable(
- "velocity_gas", mesh.getDimension(),
- &LocalAssemblerInterface<DisplacementDim>::getIntPtDarcyVelocityGas);
- add_secondary_variable(
- "velocity_liquid", mesh.getDimension(),
- &LocalAssemblerInterface<DisplacementDim>::getIntPtDarcyVelocityLiquid);
-
//
// enable output of internal variables defined by material models
//
--
GitLab