From 99e739256ce20c155eaf88ea64b792b97c24f6cd Mon Sep 17 00:00:00 2001
From: Norbert Grunwald <norbert.grunwald@ufz.de>
Date: Tue, 21 Jun 2022 08:21:42 +0200
Subject: [PATCH] [PL/TH2M] Output of diffusion velocities.
---
ProcessLib/TH2M/IntegrationPointData.h | 2 +
ProcessLib/TH2M/LocalAssemblerInterface.h | 24 +++++
ProcessLib/TH2M/TH2MFEM-impl.h | 123 ++++++++++++++++++++++
ProcessLib/TH2M/TH2MFEM.h | 20 ++++
ProcessLib/TH2M/TH2MProcess.cpp | 13 +++
5 files changed, 182 insertions(+)
diff --git a/ProcessLib/TH2M/IntegrationPointData.h b/ProcessLib/TH2M/IntegrationPointData.h
index 3c11736a84e..131c8386000 100644
--- a/ProcessLib/TH2M/IntegrationPointData.h
+++ b/ProcessLib/TH2M/IntegrationPointData.h
@@ -158,6 +158,8 @@ struct IntegrationPointData final
GlobalDimMatrixType lambda;
GlobalDimVectorType d_CG;
GlobalDimVectorType d_WG;
+ GlobalDimVectorType d_CL;
+ GlobalDimVectorType d_WL;
GlobalDimVectorType w_GS;
GlobalDimVectorType w_LS;
diff --git a/ProcessLib/TH2M/LocalAssemblerInterface.h b/ProcessLib/TH2M/LocalAssemblerInterface.h
index 53a5b52df1d..d584d9c734d 100644
--- a/ProcessLib/TH2M/LocalAssemblerInterface.h
+++ b/ProcessLib/TH2M/LocalAssemblerInterface.h
@@ -60,6 +60,30 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
+ virtual std::vector<double> const& getIntPtDiffusionVelocityVapourGas(
+ 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& getIntPtDiffusionVelocityGasGas(
+ 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& getIntPtDiffusionVelocitySoluteLiquid(
+ 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& getIntPtDiffusionVelocityLiquidLiquid(
+ 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& getIntPtLiquidDensity(
const double t,
std::vector<GlobalVector*> const& x,
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index 1cf040c06e4..c2f2ce76d64 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -396,6 +396,7 @@ std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
ip_data.xmCG = 1. - c.xmWG;
ip_data.xmWG = c.xmWG;
ip_data.xmWL = c.xmWL;
+ auto const xmCL = 1. - c.xmWL;
ip_data.diffusion_coefficient_vapour = c.diffusion_coefficient_vapour;
ip_data.diffusion_coefficient_solute = c.diffusion_coefficient_solute;
@@ -411,6 +412,11 @@ std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
ip_data.dxmWG_dT * gradT;
const GlobalDimVectorType gradxmCG = -gradxmWG;
+ const GlobalDimVectorType gradxmWL = ip_data.dxmWL_dpGR * gradpGR +
+ ip_data.dxmWL_dpCap * gradpCap +
+ ip_data.dxmWL_dT * gradT;
+ const GlobalDimVectorType gradxmCL = -gradxmWL;
+
// Todo: factor -phiAlpha / xmZetaAlpha * DZetaAlpha can be evaluated in
// the respective phase transition model, here only the multiplication
// with the gradient of the mass fractions should take place.
@@ -429,6 +435,19 @@ std::vector<ConstitutiveVariables<DisplacementDim>> TH2MLocalAssembler<
ip_data.diffusion_coefficient_vapour *
gradxmWG;
+ ip_data.d_CL = xmCL == 0. ? 0. * gradxmCL // Keep d_CL's dimension and
+ // prevent division by zero
+ : -phi_L / xmCL *
+ ip_data.diffusion_coefficient_solute *
+ gradxmCL;
+
+ ip_data.d_WL = ip_data.xmWL == 0.
+ ? 0. * gradxmWL // Keep d_WG's dimension and prevent
+ // division by zero
+ : -phi_L / ip_data.xmWL *
+ ip_data.diffusion_coefficient_solute *
+ gradxmWL;
+
// ---------------------------------------------------------------------
// Derivatives for Jacobian
// ---------------------------------------------------------------------
@@ -2154,6 +2173,110 @@ std::vector<double> const& TH2MLocalAssembler<
return cache;
}
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ int DisplacementDim>
+std::vector<double> const& TH2MLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
+ getIntPtDiffusionVelocityVapourGas(
+ 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 =
+ _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].d_WG;
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ int DisplacementDim>
+std::vector<double> const& TH2MLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
+ getIntPtDiffusionVelocityGasGas(
+ 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 =
+ _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].d_CG;
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ int DisplacementDim>
+std::vector<double> const& TH2MLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
+ getIntPtDiffusionVelocitySoluteLiquid(
+ 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 =
+ _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].d_CL;
+ }
+
+ return cache;
+}
+
+template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
+ int DisplacementDim>
+std::vector<double> const& TH2MLocalAssembler<
+ ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
+ getIntPtDiffusionVelocityLiquidLiquid(
+ 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 =
+ _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].d_WL;
+ }
+
+ 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 8691f2e48fc..4535cd2f03d 100644
--- a/ProcessLib/TH2M/TH2MFEM.h
+++ b/ProcessLib/TH2M/TH2MFEM.h
@@ -179,6 +179,26 @@ private:
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const override;
+ std::vector<double> const& getIntPtDiffusionVelocityVapourGas(
+ 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& getIntPtDiffusionVelocityGasGas(
+ 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& getIntPtDiffusionVelocitySoluteLiquid(
+ 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& getIntPtDiffusionVelocityLiquidLiquid(
+ const double t,
+ std::vector<GlobalVector*> const& x,
+ std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+ std::vector<double>& cache) const override;
std::vector<ConstitutiveVariables<DisplacementDim>>
updateConstitutiveVariables(Eigen::VectorXd const& local_x,
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index bf2a5004746..2ba3ddf98d6 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -203,6 +203,19 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
add_secondary_variable(
"velocity_liquid", mesh.getDimension(),
&LocalAssemblerInterface::getIntPtDarcyVelocityLiquid);
+ add_secondary_variable(
+ "diffusion_velocity_vapour_gas", mesh.getDimension(),
+ &LocalAssemblerInterface::getIntPtDiffusionVelocityVapourGas);
+ add_secondary_variable(
+ "diffusion_velocity_gas_gas", mesh.getDimension(),
+ &LocalAssemblerInterface::getIntPtDiffusionVelocityGasGas);
+ add_secondary_variable(
+ "diffusion_velocity_solute_liquid", mesh.getDimension(),
+ &LocalAssemblerInterface::getIntPtDiffusionVelocitySoluteLiquid);
+ add_secondary_variable(
+ "diffusion_velocity_liquid_liquid", mesh.getDimension(),
+ &LocalAssemblerInterface::getIntPtDiffusionVelocityLiquidLiquid);
+
add_secondary_variable("saturation", 1,
&LocalAssemblerInterface::getIntPtSaturation);
add_secondary_variable("vapour_pressure", 1,
--
GitLab