From 04dcb4523c507cac782c51fcbb7c5bca45addac8 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon, 20 Feb 2023 17:28:55 +0100
Subject: [PATCH] [HC] Refactored the implementation of the isotropic diffusion
stability
---
.../ComponentTransportFEM.h | 113 +++++++++++-------
1 file changed, 71 insertions(+), 42 deletions(-)
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 7eae80a42c6..1e7bd7ad69d 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -12,6 +12,7 @@
#include <Eigen/Core>
#include <Eigen/Sparse>
+#include <functional>
#include <numeric>
#include <vector>
@@ -237,6 +238,10 @@ class LocalAssemblerData : public ComponentTransportLocalAssemblerInterface
using GlobalDimNodalMatrixType =
typename ShapeMatricesType::GlobalDimNodalMatrixType;
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
+ using LaplaceCoefficientFunction = std::function<GlobalDimMatrixType(
+ const LocalAssemblerData<ShapeFunction, GlobalDim>&,
+ GlobalDimMatrixType const&, GlobalDimMatrixType const&,
+ GlobalDimVectorType const&, double const, double const, double const)>;
public:
LocalAssemblerData(
@@ -286,6 +291,12 @@ public:
_ip_data[ip].pushBackState();
}
+ _laplace_coefficient_function =
+ _process_data.stabilizer
+ ? &LocalAssemblerData<ShapeFunction, GlobalDim>::
+ getHydrodynamicDispersionWithArtificialDiffusion
+ : &LocalAssemblerData<ShapeFunction,
+ GlobalDim>::getHydrodynamicDispersion;
}
void setChemicalSystemID(std::size_t const /*mesh_item_id*/) override
@@ -543,38 +554,6 @@ public:
}
}
- GlobalDimMatrixType getHydrodynamicDispersion(
- GlobalDimMatrixType const& I,
- GlobalDimMatrixType const& pore_diffusion_coefficient,
- GlobalDimVectorType const& velocity,
- double const porosity,
- double const solute_dispersivity_transverse,
- double const solute_dispersivity_longitudinal)
- {
- double const velocity_magnitude = velocity.norm();
-
- if (velocity_magnitude == 0.0)
- {
- return porosity * pore_diffusion_coefficient;
- }
-
- GlobalDimMatrixType const D_c = GlobalDimMatrixType(
- porosity * pore_diffusion_coefficient +
- solute_dispersivity_transverse * velocity_magnitude * I +
- (solute_dispersivity_longitudinal -
- solute_dispersivity_transverse) /
- velocity_magnitude * velocity * velocity.transpose());
-
- if (!_process_data.stabilizer)
- {
- return D_c;
- }
-
- return D_c + _process_data.stabilizer->computeArtificialDiffusion(
- _element.getID(), velocity_magnitude) *
- I;
- }
-
void assembleBlockMatrices(
GlobalDimVectorType const& b, int const component_id, double const t,
double const dt,
@@ -692,10 +671,10 @@ public:
t, dt);
GlobalDimMatrixType const hydrodynamic_dispersion =
- getHydrodynamicDispersion(I, pore_diffusion_coefficient,
- velocity, porosity,
- solute_dispersivity_transverse,
- solute_dispersivity_longitudinal);
+ _laplace_coefficient_function(
+ *this, I, pore_diffusion_coefficient, velocity, porosity,
+ solute_dispersivity_transverse,
+ solute_dispersivity_longitudinal);
const double R_times_phi(retardation_factor * porosity);
GlobalDimVectorType const mass_density_flow = velocity * density;
@@ -1053,10 +1032,10 @@ public:
: GlobalDimVectorType(-K_over_mu * dNdx * local_p);
GlobalDimMatrixType const hydrodynamic_dispersion =
- getHydrodynamicDispersion(I, pore_diffusion_coefficient,
- velocity, porosity,
- solute_dispersivity_transverse,
- solute_dispersivity_longitudinal);
+ _laplace_coefficient_function(
+ *this, I, pore_diffusion_coefficient, velocity, porosity,
+ solute_dispersivity_transverse,
+ solute_dispersivity_longitudinal);
double const R_times_phi = retardation_factor * porosity;
auto const N_t_N = (N.transpose() * N).eval();
@@ -1341,8 +1320,8 @@ public:
? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
: GlobalDimVectorType(-k / mu * dNdx * p);
- GlobalDimMatrixType const D =
- getHydrodynamicDispersion(I, Dp, q, phi, alpha_T, alpha_L);
+ GlobalDimMatrixType const D = _laplace_coefficient_function(
+ *this, I, Dp, q, phi, alpha_T, alpha_L);
// matrix assembly
local_Jac.noalias() +=
@@ -1820,6 +1799,56 @@ private:
Eigen::aligned_allocator<
IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType>>>
_ip_data;
+
+ LaplaceCoefficientFunction _laplace_coefficient_function;
+
+ GlobalDimMatrixType getHydrodynamicDispersion(
+ GlobalDimMatrixType const& I,
+ GlobalDimMatrixType const& pore_diffusion_coefficient,
+ GlobalDimVectorType const& velocity,
+ double const porosity,
+ double const solute_dispersivity_transverse,
+ double const solute_dispersivity_longitudinal) const
+ {
+ double const velocity_magnitude = velocity.norm();
+ if (velocity_magnitude == 0.0)
+ {
+ return porosity * pore_diffusion_coefficient;
+ }
+
+ return GlobalDimMatrixType(
+ porosity * pore_diffusion_coefficient +
+ solute_dispersivity_transverse * velocity_magnitude * I +
+ (solute_dispersivity_longitudinal -
+ solute_dispersivity_transverse) /
+ velocity_magnitude * velocity * velocity.transpose());
+ }
+
+ GlobalDimMatrixType getHydrodynamicDispersionWithArtificialDiffusion(
+ GlobalDimMatrixType const& I,
+ GlobalDimMatrixType const& pore_diffusion_coefficient,
+ GlobalDimVectorType const& velocity,
+ double const porosity,
+ double const solute_dispersivity_transverse,
+ double const solute_dispersivity_longitudinal) const
+ {
+ double const velocity_magnitude = velocity.norm();
+ if (velocity_magnitude == 0.0)
+ {
+ return porosity * pore_diffusion_coefficient;
+ }
+
+ GlobalDimMatrixType const D_c = GlobalDimMatrixType(
+ porosity * pore_diffusion_coefficient +
+ solute_dispersivity_transverse * velocity_magnitude * I +
+ (solute_dispersivity_longitudinal -
+ solute_dispersivity_transverse) /
+ velocity_magnitude * velocity * velocity.transpose());
+
+ return D_c + _process_data.stabilizer->computeArtificialDiffusion(
+ _element.getID(), velocity_magnitude) *
+ I;
+ }
};
} // namespace ComponentTransport
--
GitLab