From 93b6074e973557a1ad6c41dc72fcaf1ff4c9ef0e Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 15 Feb 2017 17:51:14 +0100
Subject: [PATCH] [Coupling] Changed some variable or class member names and
made some other
minor changes according to Dima's comments.
---
.../HeatConduction/HeatConductionFEM-impl.h | 61 +++++----
ProcessLib/HeatConduction/HeatConductionFEM.h | 46 +++++--
.../HeatConduction/HeatConductionProcess.cpp | 4 +-
.../LiquidFlowLocalAssembler-impl.h | 129 ++++++++++--------
.../LiquidFlow/LiquidFlowLocalAssembler.h | 38 +++---
ProcessLib/LiquidFlow/LiquidFlowProcess.cpp | 5 +-
ProcessLib/LiquidFlow/LiquidFlowProcess.h | 2 +-
.../TestLiquidFlowMaterialProperties.cpp | 7 +-
8 files changed, 169 insertions(+), 123 deletions(-)
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM-impl.h b/ProcessLib/HeatConduction/HeatConductionFEM-impl.h
index a4593833d0f..598bc0d1d46 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM-impl.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM-impl.h
@@ -25,9 +25,10 @@ void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
assembleHeatTransportLiquidFlow(
double const t, int const material_id, SpatialPosition& pos,
int const gravitational_axis_id,
- double const gravitational_acceleration, Eigen::MatrixXd const& perm,
+ double const gravitational_acceleration,
+ Eigen::MatrixXd const& permeability,
ProcessLib::LiquidFlow::LiquidFlowMaterialProperties const&
- liquid_flow_prop,
+ liquid_flow_properties,
std::vector<double> const& local_x, std::vector<double> const& local_p,
std::vector<double>& local_M_data, std::vector<double>& local_K_data)
{
@@ -57,23 +58,28 @@ void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
double T = 0.;
NumLib::shapeFunctionInterpolate(local_x, sm.N, T);
- // Material parameters of solid phase
+ // Thermal conductivity of solid phase.
auto const k_s = _process_data.thermal_conductivity(t, pos)[0];
+ // Specific heat capacity of solid phase.
auto const cp_s = _process_data.heat_capacity(t, pos)[0];
+ // Density of solid phase.
auto const rho_s = _process_data.density(t, pos)[0];
- // Material parameters of fluid phase
- double const cp_f = liquid_flow_prop.getHeatCapacity(p, T);
- double const k_f = liquid_flow_prop.getThermalConductivity(p, T);
- double const rho_f = liquid_flow_prop.getLiquidDensity(p, T);
+ // Thermal conductivity of liquid.
+ double const k_f = liquid_flow_properties.getThermalConductivity(p, T);
+ // Specific heat capacity of liquid.
+ double const cp_f = liquid_flow_properties.getHeatCapacity(p, T);
+ // Density of liquid.
+ double const rho_f = liquid_flow_properties.getLiquidDensity(p, T);
- // Material parameter of porosity
- double const poro =
- liquid_flow_prop.getPorosity(material_id, porosity_variable, T);
+ // Porosity of porous media.
+ double const n = liquid_flow_properties.getPorosity(
+ material_id, porosity_variable, T);
- double const effective_cp =
- (1.0 - poro) * cp_s * rho_s + poro * cp_f * rho_f;
- double const effective_K = (1.0 - poro) * k_s + poro * k_f;
+ // Effective specific heat capacity.
+ double const effective_cp = (1.0 - n) * cp_s * rho_s + n * cp_f * rho_f;
+ // Effective thermal conductivity.
+ double const effective_K = (1.0 - n) * k_s + n * k_f;
double const integration_factor =
sm.detJ * wp.getWeight() * sm.integralMeasure;
@@ -83,12 +89,13 @@ void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
local_K.noalias() +=
sm.dNdx.transpose() * effective_K * sm.dNdx * integration_factor;
- // Compute fluid flow velocity
- double const mu = liquid_flow_prop.getViscosity(p, T); // Viscosity
+ // Compute fluid flow velocity:
+ double const mu =
+ liquid_flow_properties.getViscosity(p, T); // Viscosity
GlobalDimVectorType const& darcy_velocity =
LiquidFlowVelocityCalculator::calculateVelocity(
- local_p_vec, sm, perm, mu, rho_f * gravitational_acceleration,
- gravitational_axis_id);
+ local_p_vec, sm, permeability, mu,
+ rho_f * gravitational_acceleration, gravitational_axis_id);
// Add the advection term
local_K.noalias() += cp_f * rho_f * sm.N.transpose() *
(darcy_velocity.transpose() * sm.dNdx) *
@@ -105,32 +112,36 @@ void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
std::vector<double>& /*local_b_data*/,
LocalCouplingTerm const& coupled_term)
{
- for (auto const& coupled_process_map : coupled_term.coupled_processes)
+ for (auto const& coupled_process_pair : coupled_term.coupled_processes)
{
- if (coupled_process_map.first ==
+ if (coupled_process_pair.first ==
std::type_index(typeid(ProcessLib::LiquidFlow::LiquidFlowProcess)))
{
+ assert(
+ dynamic_cast<const ProcessLib::LiquidFlow::LiquidFlowProcess*>(
+ &(coupled_process_pair.second)) != nullptr);
+
ProcessLib::LiquidFlow::LiquidFlowProcess const& pcs =
static_cast<ProcessLib::LiquidFlow::LiquidFlowProcess const&>(
- coupled_process_map.second);
+ coupled_process_pair.second);
const auto liquid_flow_prop = pcs.getLiquidFlowMaterialProperties();
const auto local_p =
- coupled_term.local_coupled_xs.at(coupled_process_map.first);
+ coupled_term.local_coupled_xs.at(coupled_process_pair.first);
SpatialPosition pos;
pos.setElementID(_element.getID());
const int material_id = liquid_flow_prop->getMaterialID(pos);
- const Eigen::MatrixXd& perm = liquid_flow_prop->getPermeability(
+ const Eigen::MatrixXd& K = liquid_flow_prop->getPermeability(
material_id, t, pos, _element.getDimension());
- if (perm.size() == 1)
+ if (K.size() == 1)
{
assembleHeatTransportLiquidFlow<
IsotropicLiquidFlowVelocityCalculator>(
t, material_id, pos, pcs.getGravitationalAxisID(),
- pcs.getGravitationalacceleration(), perm, *liquid_flow_prop,
+ pcs.getGravitationalAcceleration(), K, *liquid_flow_prop,
local_x, local_p, local_M_data, local_K_data);
}
else
@@ -138,7 +149,7 @@ void LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
assembleHeatTransportLiquidFlow<
AnisotropicLiquidFlowVelocityCalculator>(
t, material_id, pos, pcs.getGravitationalAxisID(),
- pcs.getGravitationalacceleration(), perm, *liquid_flow_prop,
+ pcs.getGravitationalAcceleration(), K, *liquid_flow_prop,
local_x, local_p, local_M_data, local_K_data);
}
}
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h
index bf8a41ca964..3c65463b854 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM.h
@@ -199,27 +199,50 @@ private:
std::vector<std::vector<double>> _heat_fluxes;
+ /**
+ * @brief Assemble local matrices and vectors of the equations of
+ * heat transport process in porous media with full saturated liquid
+ * flow.
+ *
+ * @param t Time.
+ * @param material_id Material ID of the element.
+ * @param pos Position (t, x) of the element.
+ * @param gravitational_axis_id The ID of the axis, which indicates the
+ * direction of gravity portion of the
+ * Darcy's velocity.
+ * @param gravitational_acceleration Gravitational acceleration, 9.81 in the
+ * SI unit standard.
+ * @param permeability Intrinsic permeability of liquid
+ * in porous media.
+ * @param liquid_flow_properties Liquid flow properties.
+ * @param local_x Local vector of unknowns, e.g.
+ * nodal temperatures of the element.
+ * @param local_p Local vector of nodal pore pressures.
+ * @param local_M_data Local mass matrix.
+ * @param local_K_data Local Laplace matrix.
+ */
template <typename LiquidFlowVelocityCalculator>
- void assembleHeatTransportLiquidFlow( double const t,
- int const material_id, SpatialPosition& pos,
+ void assembleHeatTransportLiquidFlow(
+ double const t, int const material_id, SpatialPosition& pos,
int const gravitational_axis_id,
- double const gravitational_acceleration, Eigen::MatrixXd const& perm,
+ double const gravitational_acceleration,
+ Eigen::MatrixXd const& permeability,
ProcessLib::LiquidFlow::LiquidFlowMaterialProperties const&
- liquid_flow_prop,
+ liquid_flow_properties,
std::vector<double> const& local_x, std::vector<double> const& local_p,
std::vector<double>& local_M_data, std::vector<double>& local_K_data);
- /// Calculator of liquid fluid flow velocity for anisotropic permeability
+ /// Calculator of liquid flow velocity for isotropic permeability
/// tensor
struct IsotropicLiquidFlowVelocityCalculator
{
static GlobalDimVectorType calculateVelocity(
Eigen::Map<const NodalVectorType> const& local_p,
- ShapeMatrices const& sm, Eigen::MatrixXd const& perm,
+ ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
double const mu, double const rho_g,
int const gravitational_axis_id)
{
- const double K = perm(0, 0) / mu;
+ const double K = permeability(0, 0) / mu;
// Compute the velocity
GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p;
// gravity term
@@ -229,21 +252,22 @@ private:
}
};
- /// Calculator of liquid fluid flow velocity for isotropic permeability
+ /// Calculator of liquid flow velocity for anisotropic permeability
/// tensor
struct AnisotropicLiquidFlowVelocityCalculator
{
static GlobalDimVectorType calculateVelocity(
Eigen::Map<const NodalVectorType> const& local_p,
- ShapeMatrices const& sm, Eigen::MatrixXd const& perm,
+ ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
double const mu, double const rho_g,
int const gravitational_axis_id)
{
- GlobalDimVectorType darcy_velocity = -perm * sm.dNdx * local_p / mu;
+ GlobalDimVectorType darcy_velocity =
+ -permeability * sm.dNdx * local_p / mu;
if (gravitational_axis_id >= 0)
{
darcy_velocity.noalias() -=
- rho_g * perm.col(gravitational_axis_id) / mu;
+ rho_g * permeability.col(gravitational_axis_id) / mu;
}
return darcy_velocity;
}
diff --git a/ProcessLib/HeatConduction/HeatConductionProcess.cpp b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
index 012fb88bbc8..c32b9ec8c1c 100644
--- a/ProcessLib/HeatConduction/HeatConductionProcess.cpp
+++ b/ProcessLib/HeatConduction/HeatConductionProcess.cpp
@@ -44,8 +44,8 @@ void HeatConductionProcess::preTimestepConcreteProcess(GlobalVector const& x,
}
else
{
- auto x0 = _x_previous_timestep.get();
- MathLib::LinAlg::copy(x, *x0);
+ auto& x0 = *_x_previous_timestep;
+ MathLib::LinAlg::copy(x, x0);
}
}
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 4fe8431baee..e460bd7075c 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -35,20 +35,21 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
pos.setElementID(_element.getID());
const int material_id = _material_properties.getMaterialID(pos);
- const Eigen::MatrixXd& perm = _material_properties.getPermeability(
+ const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
material_id, t, pos, _element.getDimension());
// Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
- // the assert must be changed to perm.rows() == _element->getDimension()
- assert(perm.rows() == GlobalDim || perm.rows() == 1);
-
- if (perm.size() == 1) // isotropic or 1D problem.
- local_assemble<IsotropicCalculator>(material_id, t, local_x,
- local_M_data, local_K_data,
- local_b_data, pos, perm);
+ // the assert must be changed to permeability.rows() ==
+ // _element->getDimension()
+ assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
+
+ if (permeability.size() == 1) // isotropic or 1D problem.
+ assembleMatrixAndVector<IsotropicCalculator>(
+ material_id, t, local_x, local_M_data, local_K_data, local_b_data,
+ pos, permeability);
else
- local_assemble<AnisotropicCalculator>(material_id, t, local_x,
- local_M_data, local_K_data,
- local_b_data, pos, perm);
+ assembleMatrixAndVector<AnisotropicCalculator>(
+ material_id, t, local_x, local_M_data, local_K_data, local_b_data,
+ pos, permeability);
}
template <typename ShapeFunction, typename IntegrationMethod,
@@ -64,32 +65,35 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
pos.setElementID(_element.getID());
const int material_id = _material_properties.getMaterialID(pos);
- const Eigen::MatrixXd& perm = _material_properties.getPermeability(
+ const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
material_id, t, pos, _element.getDimension());
// Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
- // the assert must be changed to perm.rows() == _element->getDimension()
- assert(perm.rows() == GlobalDim || perm.rows() == 1);
+ // the assert must be changed to permeability.rows() ==
+ // _element->getDimension()
+ assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
const double dt = coupled_term.dt;
- for (auto const& coupled_process_map : coupled_term.coupled_processes)
+ for (auto const& coupled_process_pair : coupled_term.coupled_processes)
{
- if (coupled_process_map.first ==
+ if (coupled_process_pair.first ==
std::type_index(
typeid(ProcessLib::HeatConduction::HeatConductionProcess)))
{
const auto local_T0 =
- coupled_term.local_coupled_xs0.at(coupled_process_map.first);
+ coupled_term.local_coupled_xs0.at(coupled_process_pair.first);
const auto local_T1 =
- coupled_term.local_coupled_xs.at(coupled_process_map.first);
+ coupled_term.local_coupled_xs.at(coupled_process_pair.first);
- if (perm.size() == 1) // isotropic or 1D problem.
- local_assembleCoupledWithHeatTransport<IsotropicCalculator>(
+ if (permeability.size() == 1) // isotropic or 1D problem.
+ assembleWithCoupledWithHeatTransport<IsotropicCalculator>(
material_id, t, dt, local_x, local_T0, local_T1,
- local_M_data, local_K_data, local_b_data, pos, perm);
+ local_M_data, local_K_data, local_b_data, pos,
+ permeability);
else
- local_assembleCoupledWithHeatTransport<AnisotropicCalculator>(
+ assembleWithCoupledWithHeatTransport<AnisotropicCalculator>(
material_id, t, dt, local_x, local_T0, local_T1,
- local_M_data, local_K_data, local_b_data, pos, perm);
+ local_M_data, local_K_data, local_b_data, pos,
+ permeability);
}
else
{
@@ -104,12 +108,13 @@ template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
template <typename LaplacianGravityVelocityCalculator>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- local_assemble(const int material_id, double const t,
- std::vector<double> const& local_x,
- std::vector<double>& local_M_data,
- std::vector<double>& local_K_data,
- std::vector<double>& local_b_data,
- SpatialPosition const& pos, Eigen::MatrixXd const& perm)
+ assembleMatrixAndVector(const int material_id, double const t,
+ std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ SpatialPosition const& pos,
+ Eigen::MatrixXd const& permeability)
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -128,6 +133,8 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
// the integration loop for non-constant porosity and storage models.
double porosity_variable = 0.;
double storage_variable = 0.;
+ // Reference temperature for fluid property models because of no coupled
+ // heat transport process. TODO: Add an optional input for this.
const double temperature =
MaterialLib::PhysicalConstant::CelsiusZeroInKelvin + 18.0;
for (unsigned ip = 0; ip < n_integration_points; ip++)
@@ -156,7 +163,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
// Assemble Laplacian, K, and RHS by the gravitational term
LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
- local_K, local_b, sm, perm, integration_factor, mu, rho_g,
+ local_K, local_b, sm, permeability, integration_factor, mu, rho_g,
_gravitational_axis_id);
}
}
@@ -165,12 +172,12 @@ template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
template <typename LaplacianGravityVelocityCalculator>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- local_assembleCoupledWithHeatTransport(
+ assembleWithCoupledWithHeatTransport(
const int material_id, double const t, double const dt,
std::vector<double> const& local_x, std::vector<double> const& local_T0,
std::vector<double> const& local_T1, std::vector<double>& local_M_data,
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
- SpatialPosition const& pos, Eigen::MatrixXd const& perm)
+ SpatialPosition const& pos, Eigen::MatrixXd const& permeability)
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -219,7 +226,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
// Assemble Laplacian, K, and RHS by the gravitational term
LaplacianGravityVelocityCalculator::calculateLaplacianAndGravityTerm(
- local_K, local_b, sm, perm, integration_factor, mu, rho_g,
+ local_K, local_b, sm, permeability, integration_factor, mu, rho_g,
_gravitational_axis_id);
// Add the thermal expansion term
@@ -245,19 +252,20 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
SpatialPosition pos;
pos.setElementID(_element.getID());
const int material_id = _material_properties.getMaterialID(pos);
- const Eigen::MatrixXd& perm = _material_properties.getPermeability(
+ const Eigen::MatrixXd& permeability = _material_properties.getPermeability(
material_id, t, pos, _element.getDimension());
// Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
- // the assert must be changed to perm.rows() == _element->getDimension()
- assert(perm.rows() == GlobalDim || perm.rows() == 1);
+ // the assert must be changed to permeability.rows() ==
+ // _element->getDimension()
+ assert(permeability.rows() == GlobalDim || permeability.rows() == 1);
- if (perm.size() == 1) // isotropic or 1D problem.
+ if (permeability.size() == 1) // isotropic or 1D problem.
computeSecondaryVariableLocal<IsotropicCalculator>(t, local_x, pos,
- perm);
+ permeability);
else
computeSecondaryVariableLocal<AnisotropicCalculator>(t, local_x, pos,
- perm);
+ permeability);
}
template <typename ShapeFunction, typename IntegrationMethod,
@@ -267,7 +275,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
computeSecondaryVariableLocal(double const /*t*/,
std::vector<double> const& local_x,
SpatialPosition const& /*pos*/,
- Eigen::MatrixXd const& perm)
+ Eigen::MatrixXd const& permeability)
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -293,7 +301,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
const double mu = _material_properties.getViscosity(p, temperature);
LaplacianGravityVelocityCalculator::calculateVelocity(
- _darcy_velocities, local_p_vec, sm, perm, ip, mu, rho_g,
+ _darcy_velocities, local_p_vec, sm, permeability, ip, mu, rho_g,
_gravitational_axis_id);
}
}
@@ -307,7 +315,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
std::vector<double> const& local_x,
std::vector<double> const& local_T,
SpatialPosition const& /*pos*/,
- Eigen::MatrixXd const& perm)
+ Eigen::MatrixXd const& permeability)
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS);
@@ -332,7 +340,7 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
const double mu = _material_properties.getViscosity(p, T);
LaplacianGravityVelocityCalculator::calculateVelocity(
- _darcy_velocities, local_p_vec, sm, perm, ip, mu, rho_g,
+ _darcy_velocities, local_p_vec, sm, permeability, ip, mu, rho_g,
_gravitational_axis_id);
}
}
@@ -343,16 +351,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
IsotropicCalculator::calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
- Eigen::MatrixXd const& perm, double const integration_factor,
+ Eigen::MatrixXd const& permeability, double const integration_factor,
double const mu, double const rho_g, int const gravitational_axis_id)
{
- const double K = perm(0, 0) / mu;
+ const double K = permeability(0, 0) / mu;
const double fac = K * integration_factor;
local_K.noalias() += fac * sm.dNdx.transpose() * sm.dNdx;
if (gravitational_axis_id >= 0)
{
- // Note: Since perm, K, is a scalar number in this function,
+ // Note: Since permeability, K, is a scalar number in this function,
// (gradN)*K*V becomes K*(grad N)*V. With the gravity vector of V,
// the simplification of (grad N)*V is the gravitational_axis_id th
// column of the transpose of (grad N) multiplied with rho_g.
@@ -367,10 +375,11 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
IsotropicCalculator::calculateVelocity(
std::vector<std::vector<double>>& darcy_velocities,
Eigen::Map<const NodalVectorType> const& local_p,
- ShapeMatrices const& sm, Eigen::MatrixXd const& perm, unsigned const ip,
- double const mu, double const rho_g, int const gravitational_axis_id)
+ ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
+ unsigned const ip, double const mu, double const rho_g,
+ int const gravitational_axis_id)
{
- const double K = perm(0, 0) / mu;
+ const double K = permeability(0, 0) / mu;
// Compute the velocity
GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p;
// gravity term
@@ -388,18 +397,19 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
AnisotropicCalculator::calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
- Eigen::MatrixXd const& perm, double const integration_factor,
+ Eigen::MatrixXd const& permeability, double const integration_factor,
double const mu, double const rho_g, int const gravitational_axis_id)
{
const double fac = integration_factor / mu;
- local_K.noalias() += fac * sm.dNdx.transpose() * perm * sm.dNdx;
+ local_K.noalias() += fac * sm.dNdx.transpose() * permeability * sm.dNdx;
if (gravitational_axis_id >= 0)
{
- // Note: perm * gravity_vector = rho_g * perm.col(gravitational_axis_id)
+ // Note: permeability * gravity_vector = rho_g *
+ // permeability.col(gravitational_axis_id)
// i.e the result is the gravitational_axis_id th column of
- // perm multiplied with rho_g.
- local_b.noalias() -=
- fac * rho_g * sm.dNdx.transpose() * perm.col(gravitational_axis_id);
+ // permeability multiplied with rho_g.
+ local_b.noalias() -= fac * rho_g * sm.dNdx.transpose() *
+ permeability.col(gravitational_axis_id);
}
}
@@ -409,15 +419,16 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
AnisotropicCalculator::calculateVelocity(
std::vector<std::vector<double>>& darcy_velocities,
Eigen::Map<const NodalVectorType> const& local_p,
- ShapeMatrices const& sm, Eigen::MatrixXd const& perm, unsigned const ip,
- double const mu, double const rho_g, int const gravitational_axis_id)
+ ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
+ unsigned const ip, double const mu, double const rho_g,
+ int const gravitational_axis_id)
{
// Compute the velocity
- GlobalDimVectorType darcy_velocity = -perm * sm.dNdx * local_p / mu;
+ GlobalDimVectorType darcy_velocity = -permeability * sm.dNdx * local_p / mu;
if (gravitational_axis_id >= 0)
{
darcy_velocity.noalias() -=
- rho_g * perm.col(gravitational_axis_id) / mu;
+ rho_g * permeability.col(gravitational_axis_id) / mu;
}
for (unsigned d = 0; d < GlobalDim; ++d)
{
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 1e2a7a7ed64..6b2bf9e5ce3 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -146,14 +146,14 @@ private:
static void calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
- Eigen::MatrixXd const& perm, double const integration_factor,
- double const mu, double const rho_g,
- int const gravitational_axis_id);
+ Eigen::MatrixXd const& permeability,
+ double const integration_factor, double const mu,
+ double const rho_g, int const gravitational_axis_id);
static void calculateVelocity(
std::vector<std::vector<double>>& darcy_velocities,
Eigen::Map<const NodalVectorType> const& local_p,
- ShapeMatrices const& sm, Eigen::MatrixXd const& perm,
+ ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
unsigned const ip, double const mu, double const rho_g,
int const gravitational_axis_id);
};
@@ -167,40 +167,40 @@ private:
static void calculateLaplacianAndGravityTerm(
Eigen::Map<NodalMatrixType>& local_K,
Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
- Eigen::MatrixXd const& perm, double const integration_factor,
- double const mu, double const rho_g,
- int const gravitational_axis_id);
+ Eigen::MatrixXd const& permeability,
+ double const integration_factor, double const mu,
+ double const rho_g, int const gravitational_axis_id);
static void calculateVelocity(
std::vector<std::vector<double>>& darcy_velocities,
Eigen::Map<const NodalVectorType> const& local_p,
- ShapeMatrices const& sm, Eigen::MatrixXd const& perm,
+ ShapeMatrices const& sm, Eigen::MatrixXd const& permeability,
unsigned const ip, double const mu, double const rho_g,
int const gravitational_axis_id);
};
template <typename LaplacianGravityVelocityCalculator>
- void local_assemble(const int material_id, double const t,
- std::vector<double> const& local_x,
- std::vector<double>& local_M_data,
- std::vector<double>& local_K_data,
- std::vector<double>& local_b_data,
- SpatialPosition const& pos,
- Eigen::MatrixXd const& perm);
+ void assembleMatrixAndVector(const int material_id, double const t,
+ std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ SpatialPosition const& pos,
+ Eigen::MatrixXd const& permeability);
template <typename LaplacianGravityVelocityCalculator>
- void local_assembleCoupledWithHeatTransport(
+ void assembleWithCoupledWithHeatTransport(
const int material_id, double const t, double const dt,
std::vector<double> const& local_x, std::vector<double> const& local_T0,
std::vector<double> const& local_T1, std::vector<double>& local_M_data,
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
- SpatialPosition const& pos, Eigen::MatrixXd const& perm);
+ SpatialPosition const& pos, Eigen::MatrixXd const& permeability);
template <typename LaplacianGravityVelocityCalculator>
void computeSecondaryVariableLocal(double const /*t*/,
std::vector<double> const& local_x,
SpatialPosition const& pos,
- Eigen::MatrixXd const& perm);
+ Eigen::MatrixXd const& permeability);
template <typename LaplacianGravityVelocityCalculator>
void computeSecondaryVariableCoupledWithHeatTransportLocal(
@@ -208,7 +208,7 @@ private:
std::vector<double> const& local_x,
std::vector<double> const& local_T,
SpatialPosition const& pos,
- Eigen::MatrixXd const& perm);
+ Eigen::MatrixXd const& permeability);
const int _gravitational_axis_id;
const double _gravitational_acceleration;
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index ac347f517ff..4ab25b277b4 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -90,9 +90,8 @@ void LiquidFlowProcess::initializeConcreteProcess(
}
void LiquidFlowProcess::assembleConcreteProcess(
- const double t, GlobalVector const& x,
- GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
- StaggeredCouplingTerm const& coupling_term)
+ const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, StaggeredCouplingTerm const& coupling_term)
{
DBUG("Assemble LiquidFlowProcess.");
// Call global assembler for each local assembly item.
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index 5d110a8c162..fdb2092e790 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -77,7 +77,7 @@ public:
bool isLinear() const override { return true; }
int getGravitationalAxisID() const { return _gravitational_axis_id; }
- double getGravitationalacceleration() const
+ double getGravitationalAcceleration() const
{
return _gravitational_acceleration;
}
diff --git a/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp b/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp
index c6d445f99bd..de2cda71a2d 100644
--- a/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp
+++ b/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp
@@ -82,14 +82,15 @@ TEST(MaterialProcessLibLiquidFlow, checkLiquidFlowMaterialProperties)
const bool has_material_ids = false;
std::vector<std::unique_ptr<ProcessLib::ParameterBase>> parameters;
- const auto lprop = createLiquidFlowMaterialProperties(
+ const auto liquid_properties = createLiquidFlowMaterialProperties(
sub_config, parameters, has_material_ids, *dummy_property_vector);
ProcessLib::SpatialPosition pos;
pos.setElementID(0);
// Check permeability
- const Eigen::MatrixXd& perm = lprop->getPermeability(0, 0., pos, 1);
+ const Eigen::MatrixXd& perm =
+ liquid_properties->getPermeability(0, 0., pos, 1);
ASSERT_EQ(2.e-10, perm(0, 0));
ASSERT_EQ(0., perm(0, 1));
ASSERT_EQ(0., perm(0, 2));
@@ -103,6 +104,6 @@ TEST(MaterialProcessLibLiquidFlow, checkLiquidFlowMaterialProperties)
const double T = 273.15 + 60.0;
const double p = 1.e+6;
const double mass_coef =
- lprop->getMassCoefficient(0, 0., pos, p, T, 0., 0.);
+ liquid_properties->getMassCoefficient(0, 0., pos, p, T, 0., 0.);
ASSERT_NEAR(0.000100000093, mass_coef, 1.e-10);
}
--
GitLab