Skip to content
Snippets Groups Projects
Commit 19df81f1 authored by Tom Fischer's avatar Tom Fischer
Browse files

[PL] RichardsComponentTransport: Update secondary var. computation.

parent 0e3d135d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h
+ 90
− 83
View file @ 19df81f1
...@@ -16,6 +16,7 @@ ...@@ -16,6 +16,7 @@
#include "RichardsComponentTransportProcessData.h" #include "RichardsComponentTransportProcessData.h"
#include "MaterialLib/Fluid/FluidProperties/FluidProperties.h" #include "MaterialLib/Fluid/FluidProperties/FluidProperties.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h" #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "NumLib/Extrapolation/ExtrapolatableElement.h" #include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h" #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h" #include "NumLib/Fem/ShapeMatrixPolicy.h"
...@@ -64,23 +65,11 @@ public: ...@@ -64,23 +65,11 @@ public:
NumLib::LocalToGlobalIndexMap const& /*dof_table*/, NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0; std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityX( virtual std::vector<double> const& getIntPtDarcyVelocity(
const double /*t*/, const double t,
GlobalVector const& /*current_solution*/, GlobalVector const& current_solution,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/, NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<double>& /*cache*/) const = 0; std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityY(
const double /*t*/,
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
virtual std::vector<double> const& getIntPtDarcyVelocityZ(
const double /*t*/,
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
}; };
template <typename ShapeFunction, typename IntegrationMethod, template <typename ShapeFunction, typename IntegrationMethod,
...@@ -116,12 +105,7 @@ public: ...@@ -116,12 +105,7 @@ public:
RichardsComponentTransportProcessData const& process_data) RichardsComponentTransportProcessData const& process_data)
: _element(element), : _element(element),
_process_data(process_data), _process_data(process_data),
_integration_method(integration_order), _integration_method(integration_order)
_saturation(
std::vector<double>(_integration_method.getNumberOfPoints())),
_darcy_velocities(
GlobalDim,
std::vector<double>(_integration_method.getNumberOfPoints()))
{ {
// This assertion is valid only if all nodal d.o.f. use the same shape // This assertion is valid only if all nodal d.o.f. use the same shape
// matrices. // matrices.
...@@ -219,7 +203,6 @@ public: ...@@ -219,7 +203,6 @@ public:
.getCapillaryPressureSaturationModel(t, pos) .getCapillaryPressureSaturationModel(t, pos)
.getSaturation(pc_int_pt) .getSaturation(pc_int_pt)
: 1.0; : 1.0;
_saturation[ip] = Sw;
double const dSw_dpc = double const dSw_dpc =
(pc_int_pt > 0) (pc_int_pt > 0)
...@@ -228,6 +211,7 @@ public: ...@@ -228,6 +211,7 @@ public:
.getCapillaryPressureSaturationModel(t, pos) .getCapillaryPressureSaturationModel(t, pos)
.getdPcdS(Sw) .getdPcdS(Sw)
: 0.; : 0.;
// \todo the first argument has to be changed for non constant // \todo the first argument has to be changed for non constant
// porosity model // porosity model
auto const porosity = auto const porosity =
...@@ -262,7 +246,6 @@ public: ...@@ -262,7 +246,6 @@ public:
// Use the viscosity model to compute the viscosity // Use the viscosity model to compute the viscosity
auto const mu = _process_data.fluid_properties->getValue( auto const mu = _process_data.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars); MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
auto const K_times_k_rel_over_mu = K * (k_rel/mu); auto const K_times_k_rel_over_mu = K * (k_rel/mu);
GlobalDimVectorType const velocity = GlobalDimVectorType const velocity =
...@@ -315,24 +298,29 @@ public: ...@@ -315,24 +298,29 @@ public:
} }
} }
void computeSecondaryVariableConcrete( std::vector<double> const& getIntPtDarcyVelocity(
double const t, std::vector<double> const& local_x) override const double t,
GlobalVector const& current_solution,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<double>& cache) const override
{ {
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
auto const indices = NumLib::getIndices(_element.getID(), dof_table);
assert(!indices.empty());
auto const local_x = current_solution.get(indices);
cache.clear();
auto cache_mat = MathLib::createZeroedMatrix<
Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
cache, GlobalDim, n_integration_points);
SpatialPosition pos; SpatialPosition pos;
pos.setElementID(_element.getID()); pos.setElementID(_element.getID());
auto const& K =
_process_data.porous_media_properties.getIntrinsicPermeability(t,
pos);
MaterialLib::Fluid::FluidProperty::ArrayType vars; MaterialLib::Fluid::FluidProperty::ArrayType vars;
auto const mu = _process_data.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
GlobalDimMatrixType const K_over_mu = K / mu;
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
auto const p_nodal_values = Eigen::Map<const NodalVectorType>( auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
&local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS); &local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
...@@ -342,13 +330,34 @@ public: ...@@ -342,13 +330,34 @@ public:
auto const& N = ip_data.N; auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx; auto const& dNdx = ip_data.dNdx;
GlobalDimVectorType velocity = -K_over_mu * dNdx * p_nodal_values; pos.setIntegrationPoint(ip);
auto const& K =
_process_data.porous_media_properties.getIntrinsicPermeability(
t, pos);
auto const mu = _process_data.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
double C_int_pt = 0.0;
double p_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
// saturation
double const pc_int_pt = -p_int_pt;
double const Sw =
(pc_int_pt > 0)
? _process_data.porous_media_properties
.getCapillaryPressureSaturationModel(t, pos)
.getSaturation(pc_int_pt)
: 1.0;
auto const& k_rel = _process_data.porous_media_properties
.getRelativePermeability(t, pos)
.getValue(Sw);
cache_mat.col(ip).noalias() = -dNdx * p_nodal_values;
if (_process_data.has_gravity) if (_process_data.has_gravity)
{ {
double C_int_pt = 0.0;
double p_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt,
p_int_pt);
vars[static_cast<int>( vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt; MaterialLib::Fluid::PropertyVariableType::C)] = C_int_pt;
vars[static_cast<int>( vars[static_cast<int>(
...@@ -358,14 +367,11 @@ public: ...@@ -358,14 +367,11 @@ public:
MaterialLib::Fluid::FluidPropertyType::Density, vars); MaterialLib::Fluid::FluidPropertyType::Density, vars);
auto const b = _process_data.specific_body_force; auto const b = _process_data.specific_body_force;
// here it is assumed that the vector b is directed 'downwards' // here it is assumed that the vector b is directed 'downwards'
velocity += K_over_mu * rho_w * b; cache_mat.col(ip).noalias() += rho_w * b;
}
for (unsigned d = 0; d < GlobalDim; ++d)
{
_darcy_velocities[d][ip] = velocity[d];
} }
cache_mat.col(ip).noalias() = k_rel/mu * (K * cache_mat.col(ip));
} }
return cache;
} }
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix( Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
...@@ -378,43 +384,46 @@ public: ...@@ -378,43 +384,46 @@ public:
} }
std::vector<double> const& getIntPtSaturation( std::vector<double> const& getIntPtSaturation(
const double /*t*/, const double t,
GlobalVector const& /*current_solution*/, GlobalVector const& current_solution,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/, NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<double>& /*cache*/) const override std::vector<double>& cache) const override
{ {
assert(!_saturation.empty()); SpatialPosition pos;
return _saturation; pos.setElementID(_element.getID());
}
std::vector<double> const& getIntPtDarcyVelocityX( unsigned const n_integration_points =
const double /*t*/, _integration_method.getNumberOfPoints();
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 0);
return _darcy_velocities[0];
}
std::vector<double> const& getIntPtDarcyVelocityY( auto const indices = NumLib::getIndices(_element.getID(), dof_table);
const double /*t*/, assert(!indices.empty());
GlobalVector const& /*current_solution*/, auto const local_x = current_solution.get(indices);
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override
{
assert(_darcy_velocities.size() > 1);
return _darcy_velocities[1];
}
std::vector<double> const& getIntPtDarcyVelocityZ( cache.clear();
const double /*t*/, auto cache_vec = MathLib::createZeroedVector<LocalVectorType>(
GlobalVector const& /*current_solution*/, cache, n_integration_points);
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const override for (unsigned ip = 0; ip < n_integration_points; ++ip)
{ {
assert(_darcy_velocities.size() > 2); auto const& ip_data = _ip_data[ip];
return _darcy_velocities[2]; auto const& N = ip_data.N;
double C_int_pt = 0.0;
double p_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_x, N, C_int_pt, p_int_pt);
// saturation
double const pc_int_pt = -p_int_pt;
double const Sw =
(pc_int_pt > 0)
? _process_data.porous_media_properties
.getCapillaryPressureSaturationModel(t, pos)
.getSaturation(pc_int_pt)
: 1.0;
cache_vec[ip] = Sw;
}
return cache;
} }
private: private:
...@@ -428,8 +437,6 @@ private: ...@@ -428,8 +437,6 @@ private:
Eigen::aligned_allocator<IntegrationPointData< Eigen::aligned_allocator<IntegrationPointData<
NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>>> NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>>>
_ip_data; _ip_data;
std::vector<double> _saturation;
std::vector<std::vector<double>> _darcy_velocities;
}; };
} // namespace RichardsComponentTransport } // namespace RichardsComponentTransport
......
ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.cpp
+ 4
− 32
View file @ 19df81f1
...@@ -45,29 +45,12 @@ void RichardsComponentTransportProcess::initializeConcreteProcess( ...@@ -45,29 +45,12 @@ void RichardsComponentTransportProcess::initializeConcreteProcess(
mesh.isAxiallySymmetric(), integration_order, _process_data); mesh.isAxiallySymmetric(), integration_order, _process_data);
_secondary_variables.addSecondaryVariable( _secondary_variables.addSecondaryVariable(
"darcy_velocity_x", "darcy_velocity",
makeExtrapolator(1, getExtrapolator(), _local_assemblers, makeExtrapolator(mesh.getDimension(), getExtrapolator(),
_local_assemblers,
&RichardsComponentTransportLocalAssemblerInterface:: &RichardsComponentTransportLocalAssemblerInterface::
getIntPtDarcyVelocityX)); getIntPtDarcyVelocity));
if (mesh.getDimension() > 1)
{
_secondary_variables.addSecondaryVariable(
"darcy_velocity_y",
makeExtrapolator(
1, getExtrapolator(), _local_assemblers,
&RichardsComponentTransportLocalAssemblerInterface::
getIntPtDarcyVelocityY));
}
if (mesh.getDimension() > 2)
{
_secondary_variables.addSecondaryVariable(
"darcy_velocity_z",
makeExtrapolator(
1, getExtrapolator(), _local_assemblers,
&RichardsComponentTransportLocalAssemblerInterface::
getIntPtDarcyVelocityZ));
}
_secondary_variables.addSecondaryVariable( _secondary_variables.addSecondaryVariable(
"saturation", "saturation",
makeExtrapolator(1, getExtrapolator(), _local_assemblers, makeExtrapolator(1, getExtrapolator(), _local_assemblers,
...@@ -104,17 +87,6 @@ void RichardsComponentTransportProcess::assembleWithJacobianConcreteProcess( ...@@ -104,17 +87,6 @@ void RichardsComponentTransportProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac, _coupling_term); dx_dx, M, K, b, Jac, _coupling_term);
} }
void RichardsComponentTransportProcess::computeSecondaryVariableConcrete(
double const t, GlobalVector const& x,
StaggeredCouplingTerm const& coupling_term)
{
DBUG("Compute the Darcy velocity for RichardsComponentTransportProcess.");
GlobalExecutor::executeMemberOnDereferenced(
&RichardsComponentTransportLocalAssemblerInterface::
computeSecondaryVariable,
_local_assemblers, *_local_to_global_index_map, t, x, coupling_term);
}
} // namespace RichardsComponentTransport } // namespace RichardsComponentTransport
} // namespace ProcessLib } // namespace ProcessLib
ProcessLib/RichardsComponentTransport/RichardsComponentTransportProcess.h
+ 0
− 4
View file @ 19df81f1
...@@ -58,10 +58,6 @@ public: ...@@ -58,10 +58,6 @@ public:
bool isLinear() const override { return false; } bool isLinear() const override { return false; }
//! @} //! @}
void computeSecondaryVariableConcrete(
double const t, GlobalVector const& x,
StaggeredCouplingTerm const& coupling_term) override;
private: private:
void initializeConcreteProcess( void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table, NumLib::LocalToGlobalIndexMap const& dof_table,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment