OGS
ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim > Class Template Reference

Detailed Description

template<typename ShapeFunction, int GlobalDim>
class ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >

Definition at line 205 of file ComponentTransportFEM.h.

#include <ComponentTransportFEM.h>

Inheritance diagram for ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >:
[legend]
Collaboration diagram for ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >:
[legend]

Public Member Functions

 LocalAssemblerData (MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, ComponentTransportProcessData const &process_data, std::vector< std::reference_wrapper< ProcessVariable > > const &transport_process_variables)
void setChemicalSystemID (std::size_t const) override
void initializeChemicalSystemConcrete (Eigen::VectorXd const &local_x, double const t) override
void setChemicalSystemConcrete (Eigen::VectorXd const &local_x, double const t, double dt) override
void postSpeciationCalculation (std::size_t const ele_id, double const t, double const dt) override
void assemble (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
void assembleBlockMatrices (GlobalDimVectorType const &b, int const component_id, double const t, double const dt, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< LocalBlockMatrixType > KCC, Eigen::Ref< LocalBlockMatrixType > MCC, Eigen::Ref< LocalBlockMatrixType > MCp, Eigen::Ref< LocalBlockMatrixType > MpC, Eigen::Ref< LocalBlockMatrixType > Kpp, Eigen::Ref< LocalBlockMatrixType > Mpp, Eigen::Ref< LocalSegmentVectorType > Bp)
void assembleKCmCn (int const component_id, double const t, double const dt, Eigen::Ref< LocalBlockMatrixType > KCmCn, double const stoichiometric_coefficient, double const kinetic_prefactor)
void assembleForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
void assembleHydraulicEquation (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
void assembleHeatTransportEquation (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &)
void assembleComponentTransportEquation (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &, int const transport_process_id)
void assembleWithJacobianForStaggeredScheme (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
void assembleWithJacobianHydraulicEquation (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void assembleWithJacobianComponentTransportEquation (double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, int const component_id)
void assembleReactionEquationConcrete (double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id) override
std::vector< double > const & getIntPtLiquidDensity (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 & calculateIntPtLiquidDensity (const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &T_nodal_values, std::vector< double > &cache) const
std::vector< double > const & getIntPtDarcyVelocity (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 & calculateIntPtDarcyVelocity (const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &T_nodal_values, std::vector< double > &cache) const
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix (const unsigned integration_point) const override
 Provides the shape matrix at the given integration point.
Eigen::Vector3d getFlux (MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
void computeSecondaryVariableConcrete (double const t, double const, Eigen::VectorXd const &local_x, Eigen::VectorXd const &) override
void computeReactionRelatedSecondaryVariable (std::size_t const ele_id) override
std::vector< double > const & getIntPtMolarFlux (const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< double > &cache, int const component_id) const override
void postTimestepConcrete (Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
Public Member Functions inherited from ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface
 ComponentTransportLocalAssemblerInterface ()=default
void initializeChemicalSystem (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t)
void setChemicalSystem (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
void assembleReactionEquation (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, int const process_id)
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
virtual ~LocalAssemblerInterface ()=default
virtual void setInitialConditions (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, int const process_id)
virtual void initialize (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
virtual void preAssemble (double const, double const, std::vector< double > const &)
virtual void assembleWithJacobian (double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
virtual void computeSecondaryVariable (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
virtual void preTimestep (std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table, GlobalVector const &x, double const t, double const delta_t)
virtual void postTimestep (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
void postNonLinearSolver (std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, int const process_id)
virtual Eigen::Vector3d getFlux (MathLib::Point3d const &, double const, std::vector< std::vector< double > > const &) const
 Fits to staggered scheme.
virtual int getNumberOfVectorElementsForDeformation () const
Public Member Functions inherited from NumLib::ExtrapolatableElement
virtual ~ExtrapolatableElement ()=default

Private Types

using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices
using LocalBlockMatrixType
using LocalSegmentVectorType
using LocalMatrixType
using LocalVectorType = Eigen::Matrix<double, Eigen::Dynamic, 1>
using NodalVectorType = typename ShapeMatricesType::NodalVectorType
using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType
using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
using GlobalDimNodalMatrixType
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType

Private Member Functions

void updateAveragePorosity (std::size_t const ele_id)
double getHeatEnergyCoefficient (MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
GlobalDimMatrixType getThermalConductivityDispersivity (MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
NodalVectorType getLocalTemperature (double const t, Eigen::VectorXd const &local_x) const

Private Attributes

const int temperature_index = -1
const int first_concentration_index = -1
MeshLib::Element const & _element
ComponentTransportProcessData const & _process_data
NumLib::GenericIntegrationMethod const & _integration_method
std::vector< std::reference_wrapper< ProcessVariable > > const _transport_process_variables
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data

Static Private Attributes

static const int pressure_index = 0
static const int pressure_size = ShapeFunction::NPOINTS
static const int temperature_size = ShapeFunction::NPOINTS
static const int concentration_size

Member Typedef Documentation

◆ GlobalDimMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType
private

Definition at line 237 of file ComponentTransportFEM.h.

◆ GlobalDimNodalMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::GlobalDimNodalMatrixType
private
Initial value:
typename ShapeMatricesType::GlobalDimNodalMatrixType
EigenFixedShapeMatrixPolicy::GlobalDimNodalMatrixType
MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType
Definition ShapeMatrixPolicy.h:75

Definition at line 235 of file ComponentTransportFEM.h.

◆ GlobalDimVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType
private

Definition at line 234 of file ComponentTransportFEM.h.

◆ LocalBlockMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::LocalBlockMatrixType
private
Initial value:
typename ShapeMatricesType::template MatrixType<pressure_size,
pressure_size>

Definition at line 221 of file ComponentTransportFEM.h.

◆ LocalMatrixType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::LocalMatrixType
private
Initial value:
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>

Definition at line 227 of file ComponentTransportFEM.h.

◆ LocalSegmentVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::LocalSegmentVectorType
private
Initial value:
typename ShapeMatricesType::template VectorType<pressure_size>

Definition at line 224 of file ComponentTransportFEM.h.

◆ LocalVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::LocalVectorType = Eigen::Matrix<double, Eigen::Dynamic, 1>
private

Definition at line 229 of file ComponentTransportFEM.h.

◆ NodalRowVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType
private

Definition at line 232 of file ComponentTransportFEM.h.

◆ NodalVectorType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::NodalVectorType = typename ShapeMatricesType::NodalVectorType
private

Definition at line 231 of file ComponentTransportFEM.h.

◆ ShapeMatrices

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::ShapeMatrices = typename ShapeMatricesType::ShapeMatrices
private

Definition at line 219 of file ComponentTransportFEM.h.

◆ ShapeMatricesType

template<typename ShapeFunction, int GlobalDim>
using ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>
private

Definition at line 218 of file ComponentTransportFEM.h.

Constructor & Destructor Documentation

◆ LocalAssemblerData()

template<typename ShapeFunction, int GlobalDim>
ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::LocalAssemblerData ( MeshLib::Element const & element,
std::size_t const local_matrix_size,
NumLib::GenericIntegrationMethod const & integration_method,
bool is_axially_symmetric,
ComponentTransportProcessData const & process_data,
std::vector< std::reference_wrapper< ProcessVariable > > const & transport_process_variables )
inline

Definition at line 240 of file ComponentTransportFEM.h.

248 : temperature_index(process_data.isothermal ? -1
249 : ShapeFunction::NPOINTS),
250 first_concentration_index(process_data.isothermal
251 ? ShapeFunction::NPOINTS
252 : 2 * ShapeFunction::NPOINTS),
253 _element(element),
254 _process_data(process_data),
255 _integration_method(integration_method),
256 _transport_process_variables(transport_process_variables)
257 {
258 (void)local_matrix_size;
259
260 unsigned const n_integration_points =
261 _integration_method.getNumberOfPoints();
262 _ip_data.reserve(n_integration_points);
263
264 ParameterLib::SpatialPosition const element_pos({}, _element.getID(),
265 {});
266
267 double const aperture_size =
268 _process_data.aperture_size(0.0, element_pos)[0];
269
270 auto const shape_matrices =
271 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
272 GlobalDim>(element, is_axially_symmetric,
273 _integration_method);
274 auto const& Ns =
275 _process_data.shape_matrix_cache
276 .NsHigherOrder<typename ShapeFunction::MeshElement>();
277 auto const& medium =
278 *_process_data.media_map.getMedium(_element.getID());
279 for (unsigned ip = 0; ip < n_integration_points; ip++)
280 {
281 ParameterLib::SpatialPosition const pos(
282 {}, _element.getID(),
283 MathLib::Point3d(
284 NumLib::interpolateCoordinates<ShapeFunction,
285 ShapeMatricesType>(_element,
286 Ns[ip])));
287
288 _ip_data.emplace_back(
289 shape_matrices[ip].dNdx,
290 _integration_method.getWeightedPoint(ip).getWeight() *
291 shape_matrices[ip].integralMeasure *
292 shape_matrices[ip].detJ * aperture_size);
293
294 _ip_data[ip].porosity =
295 medium[MaterialPropertyLib::PropertyType::porosity]
296 .template initialValue<double>(
297 pos, std::numeric_limits<double>::quiet_NaN() /*t*/);
298
299 _ip_data[ip].pushBackState();
300 }
301 }
ProcessLib::ComponentTransport::LocalAssemblerData::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition ComponentTransportFEM.h:2313
ProcessLib::ComponentTransport::LocalAssemblerData::_ip_data
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data
Definition ComponentTransportFEM.h:2317
ProcessLib::ComponentTransport::LocalAssemblerData::_process_data
ComponentTransportProcessData const & _process_data
Definition ComponentTransportFEM.h:2311
ProcessLib::ComponentTransport::LocalAssemblerData::ShapeMatricesType
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
Definition ComponentTransportFEM.h:218
ProcessLib::ComponentTransport::LocalAssemblerData::_transport_process_variables
std::vector< std::reference_wrapper< ProcessVariable > > const _transport_process_variables
Definition ComponentTransportFEM.h:2315

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, first_concentration_index, NumLib::initShapeMatrices(), NumLib::interpolateCoordinates(), MaterialPropertyLib::porosity, and temperature_index.

Member Function Documentation

◆ assemble()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assemble ( double const t,
double const dt,
std::vector< double > const & local_x,
std::vector< double > const & ,
std::vector< double > & local_M_data,
std::vector< double > & local_K_data,
std::vector< double > & local_b_data )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 470 of file ComponentTransportFEM.h.

476 {
477 auto const local_matrix_size = local_x.size();
478 // Nodal DOFs include pressure
479 int const num_nodal_dof = 1 + _transport_process_variables.size();
480 // This assertion is valid only if all nodal d.o.f. use the same shape
481 // matrices.
482 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
483
484 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
485 local_M_data, local_matrix_size, local_matrix_size);
486 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
487 local_K_data, local_matrix_size, local_matrix_size);
488 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
489 local_b_data, local_matrix_size);
490
491 // Get block matrices
492 auto Kpp = local_K.template block<pressure_size, pressure_size>(
493 pressure_index, pressure_index);
494 auto Mpp = local_M.template block<pressure_size, pressure_size>(
495 pressure_index, pressure_index);
496 auto Bp = local_b.template segment<pressure_size>(pressure_index);
497
498 auto local_p = Eigen::Map<const NodalVectorType>(
499 &local_x[pressure_index], pressure_size);
500
501 auto const& b =
502 _process_data
503 .projected_specific_body_force_vectors[_element.getID()];
504
505 auto const number_of_components = num_nodal_dof - 1;
506 for (int component_id = 0; component_id < number_of_components;
507 ++component_id)
508 {
509 /* Partitioned assembler matrix
510 * | pp | pc1 | pc2 | pc3 |
511 * |-----|-----|-----|-----|
512 * | c1p | c1c1| 0 | 0 |
513 * |-----|-----|-----|-----|
514 * | c2p | 0 | c2c2| 0 |
515 * |-----|-----|-----|-----|
516 * | c3p | 0 | 0 | c3c3|
517 */
518 auto concentration_index =
519 pressure_size + component_id * concentration_size;
520
521 auto KCC =
522 local_K.template block<concentration_size, concentration_size>(
523 concentration_index, concentration_index);
524 auto MCC =
525 local_M.template block<concentration_size, concentration_size>(
526 concentration_index, concentration_index);
527 auto MCp =
528 local_M.template block<concentration_size, pressure_size>(
529 concentration_index, pressure_index);
530 auto MpC =
531 local_M.template block<pressure_size, concentration_size>(
532 pressure_index, concentration_index);
533
534 auto local_C = Eigen::Map<const NodalVectorType>(
535 &local_x[concentration_index], concentration_size);
536
537 assembleBlockMatrices(b, component_id, t, dt, local_C, local_p, KCC,
538 MCC, MCp, MpC, Kpp, Mpp, Bp);
539
540 if (_process_data.chemical_solver_interface)
541 {
542 auto const stoichiometric_matrix =
543 _process_data.chemical_solver_interface
544 ->getStoichiometricMatrix();
545
546 assert(stoichiometric_matrix);
547
548 for (Eigen::SparseMatrix<double>::InnerIterator it(
549 *stoichiometric_matrix, component_id);
550 it;
551 ++it)
552 {
553 auto const stoichiometric_coefficient = it.value();
554 auto const coupled_component_id = it.row();
555 auto const kinetic_prefactor =
556 _process_data.chemical_solver_interface
557 ->getKineticPrefactor(coupled_component_id);
558
559 auto const concentration_index =
560 pressure_size + component_id * concentration_size;
561 auto const coupled_concentration_index =
562 pressure_size +
563 coupled_component_id * concentration_size;
564 auto KCmCn = local_K.template block<concentration_size,
565 concentration_size>(
566 concentration_index, coupled_concentration_index);
567
568 // account for the coupling between components
569 assembleKCmCn(component_id, t, dt, KCmCn,
570 stoichiometric_coefficient,
571 kinetic_prefactor);
572 }
573 }
574 }
575 }
ProcessLib::ComponentTransport::LocalAssemblerData::assembleBlockMatrices
void assembleBlockMatrices(GlobalDimVectorType const &b, int const component_id, double const t, double const dt, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< LocalBlockMatrixType > KCC, Eigen::Ref< LocalBlockMatrixType > MCC, Eigen::Ref< LocalBlockMatrixType > MCp, Eigen::Ref< LocalBlockMatrixType > MpC, Eigen::Ref< LocalBlockMatrixType > Kpp, Eigen::Ref< LocalBlockMatrixType > Mpp, Eigen::Ref< LocalSegmentVectorType > Bp)
Definition ComponentTransportFEM.h:577
ProcessLib::ComponentTransport::LocalAssemblerData::assembleKCmCn
void assembleKCmCn(int const component_id, double const t, double const dt, Eigen::Ref< LocalBlockMatrixType > KCmCn, double const stoichiometric_coefficient, double const kinetic_prefactor)
Definition ComponentTransportFEM.h:776
MathLib::createZeroedVector
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Definition EigenMapTools.h:146
MathLib::createZeroedMatrix
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:25

References _element, _process_data, _transport_process_variables, assembleBlockMatrices(), assembleKCmCn(), concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), pressure_index, and pressure_size.

◆ assembleBlockMatrices()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleBlockMatrices ( GlobalDimVectorType const & b,
int const component_id,
double const t,
double const dt,
Eigen::Ref< const NodalVectorType > const & C_nodal_values,
Eigen::Ref< const NodalVectorType > const & p_nodal_values,
Eigen::Ref< LocalBlockMatrixType > KCC,
Eigen::Ref< LocalBlockMatrixType > MCC,
Eigen::Ref< LocalBlockMatrixType > MCp,
Eigen::Ref< LocalBlockMatrixType > MpC,
Eigen::Ref< LocalBlockMatrixType > Kpp,
Eigen::Ref< LocalBlockMatrixType > Mpp,
Eigen::Ref< LocalSegmentVectorType > Bp )
inline

Definition at line 577 of file ComponentTransportFEM.h.

589 {
590 unsigned const n_integration_points =
591 _integration_method.getNumberOfPoints();
592
593 MaterialPropertyLib::VariableArray vars;
594
595 // Get material properties
596 auto const& medium =
597 *_process_data.media_map.getMedium(_element.getID());
598 // Select the only valid for component transport liquid phase.
599 auto const& phase =
600 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
601
602 // Assume that the component name is the same as the process variable
603 // name. Components are shifted by one because the first one is always
604 // pressure.
605 auto const& component = phase.component(
606 _transport_process_variables[component_id].get().getName());
607
608 LocalBlockMatrixType KCC_Laplacian =
609 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
610
611 std::vector<GlobalDimVectorType> ip_flux_vector;
612 double average_velocity_norm = 0.0;
613 if (!_process_data.non_advective_form)
614 {
615 ip_flux_vector.reserve(n_integration_points);
616 }
617
618 auto const& Ns =
619 _process_data.shape_matrix_cache
620 .NsHigherOrder<typename ShapeFunction::MeshElement>();
621
622 for (unsigned ip(0); ip < n_integration_points; ++ip)
623 {
624 auto& ip_data = _ip_data[ip];
625 auto const& dNdx = ip_data.dNdx;
626 auto const& N = Ns[ip];
627 auto const& w = ip_data.integration_weight;
628 auto& porosity = ip_data.porosity;
629
630 ParameterLib::SpatialPosition const pos(
631 {}, _element.getID(),
632 MathLib::Point3d(
633 NumLib::interpolateCoordinates<ShapeFunction,
634 ShapeMatricesType>(_element,
635 N)));
636
637 double C_int_pt = 0.0;
638 double p_int_pt = 0.0;
639
640 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
641 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
642
643 vars.concentration = C_int_pt;
644 vars.liquid_phase_pressure = p_int_pt;
645
646 // update according to a particular porosity model
647 porosity = medium[MaterialPropertyLib::PropertyType::porosity]
648 .template value<double>(vars, pos, t, dt);
649 vars.porosity = porosity;
650
651 auto const& retardation_factor =
652 component[MaterialPropertyLib::PropertyType::retardation_factor]
653 .template value<double>(vars, pos, t, dt);
654
655 auto const& solute_dispersivity_transverse = medium.template value<
656 double>(
657 MaterialPropertyLib::PropertyType::transversal_dispersivity);
658
659 auto const& solute_dispersivity_longitudinal =
660 medium.template value<double>(
661 MaterialPropertyLib::PropertyType::
662 longitudinal_dispersivity);
663
664 // Use the fluid density model to compute the density
665 // TODO (renchao): concentration of which component as the argument
666 // for calculation of fluid density
667 auto const density =
668 phase[MaterialPropertyLib::PropertyType::density]
669 .template value<double>(vars, pos, t, dt);
670
671 auto const decay_rate =
672 component[MaterialPropertyLib::PropertyType::decay_rate]
673 .template value<double>(vars, pos, t, dt);
674
675 auto const& pore_diffusion_coefficient =
676 MaterialPropertyLib::formEigenTensor<GlobalDim>(
677 component[MaterialPropertyLib::PropertyType::pore_diffusion]
678 .value(vars, pos, t, dt));
679
680 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
681 medium[MaterialPropertyLib::PropertyType::permeability].value(
682 vars, pos, t, dt));
683
684 // Use the viscosity model to compute the viscosity
685 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
686 .template value<double>(vars, pos, t, dt);
687
688 double storage = 0;
689 if (medium.hasProperty(MaterialPropertyLib::PropertyType::storage))
690 {
691 storage = medium[MaterialPropertyLib::PropertyType::storage]
692 .template value<double>(vars, pos, t, dt);
693 }
694
695 GlobalDimMatrixType const K_over_mu = K / mu;
696 GlobalDimVectorType const velocity =
697 _process_data.has_gravity
698 ? GlobalDimVectorType(-K_over_mu *
699 (dNdx * p_nodal_values - density * b))
700 : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
701
702 const double drho_dp =
703 phase[MaterialPropertyLib::PropertyType::density]
704 .template dValue<double>(
705 vars,
706 MaterialPropertyLib::Variable::liquid_phase_pressure,
707 pos, t, dt);
708
709 const double drho_dC =
710 phase[MaterialPropertyLib::PropertyType::density]
711 .template dValue<double>(
712 vars, MaterialPropertyLib::Variable::concentration, pos,
713 t, dt);
714
715 GlobalDimMatrixType const hydrodynamic_dispersion =
716 NumLib::computeHydrodynamicDispersion(
717 _process_data.stabilizer, _element.getID(),
718 pore_diffusion_coefficient, velocity, porosity,
719 solute_dispersivity_transverse,
720 solute_dispersivity_longitudinal);
721
722 const double R_times_phi(retardation_factor * porosity);
723 GlobalDimVectorType const mass_density_flow = velocity * density;
724 auto const N_t_N = (N.transpose() * N).eval();
725
726 if (_process_data.non_advective_form)
727 {
728 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
729 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
730 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
731 }
732 else
733 {
734 ip_flux_vector.emplace_back(mass_density_flow);
735 average_velocity_norm += velocity.norm();
736 }
737 MCC.noalias() += N_t_N * (R_times_phi * density * w);
738 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
739 KCC_Laplacian.noalias() +=
740 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
741
742 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
743
744 // Calculate Mpp, Kpp, and bp in the first loop over components
745 if (component_id == 0)
746 {
747 Mpp.noalias() +=
748 N_t_N * (porosity * drho_dp * w + density * storage * w);
749 Kpp.noalias() +=
750 dNdx.transpose() * K_over_mu * dNdx * (density * w);
751
752 if (_process_data.has_gravity)
753 {
754 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
755 (density * density * w);
756 }
757 }
758 }
759
760 if (!_process_data.non_advective_form)
761 {
762 NumLib::assembleAdvectionMatrix<
763 typename ShapeFunction::MeshElement>(
764 _process_data.stabilizer,
765 _ip_data,
766 _process_data.shape_matrix_cache,
767 ip_flux_vector,
768 average_velocity_norm /
769 static_cast<double>(n_integration_points),
770 KCC_Laplacian);
771 }
772
773 KCC.noalias() += KCC_Laplacian;
774 }
getName
std::string getName(std::string const &line)
Returns the name/title from the "Zone"-description.
Definition TecPlotTools.cpp:56
ProcessLib::ComponentTransport::LocalAssemblerData::GlobalDimVectorType
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
Definition ComponentTransportFEM.h:234
ProcessLib::ComponentTransport::LocalAssemblerData::LocalBlockMatrixType
typename ShapeMatricesType::template MatrixType< pressure_size, pressure_size > LocalBlockMatrixType
Definition ComponentTransportFEM.h:221
ProcessLib::ComponentTransport::LocalAssemblerData::GlobalDimMatrixType
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
Definition ComponentTransportFEM.h:237
MaterialPropertyLib::formEigenTensor
constexpr Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Definition FormEigenTensor.h:163
NumLib::detail::shapeFunctionInterpolate
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition Interpolation.h:20
NumLib::computeHydrodynamicDispersion
Eigen::MatrixXd computeHydrodynamicDispersion(NumericalStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
Definition HydrodynamicDispersion.h:65

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::detail::assembleAdvectionMatrix(), NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, MaterialPropertyLib::formEigenTensor(), getName(), NumLib::interpolateCoordinates(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, MaterialPropertyLib::retardation_factor, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::storage, MaterialPropertyLib::transversal_dispersivity, and MaterialPropertyLib::viscosity.

Referenced by assemble().

◆ assembleComponentTransportEquation()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleComponentTransportEquation ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
std::vector< double > & local_M_data,
std::vector< double > & local_K_data,
std::vector< double > & ,
int const transport_process_id )
inline

Definition at line 1134 of file ComponentTransportFEM.h.

1139 {
1140 assert(static_cast<int>(local_x.size()) ==
1141 pressure_size +
1142 concentration_size *
1143 static_cast<int>(_transport_process_variables.size()) +
1144 (_process_data.isothermal ? 0 : temperature_size));
1145
1146 auto const local_p =
1147 local_x.template segment<pressure_size>(pressure_index);
1148
1149 NodalVectorType local_T = getLocalTemperature(t, local_x);
1150
1151 auto const local_C = local_x.template segment<concentration_size>(
1152 first_concentration_index +
1153 (transport_process_id - (_process_data.isothermal ? 1 : 2)) *
1154 concentration_size);
1155 auto const local_p_prev =
1156 local_x_prev.segment<pressure_size>(pressure_index);
1157
1158 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1159 local_M_data, concentration_size, concentration_size);
1160 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1161 local_K_data, concentration_size, concentration_size);
1162
1163 LocalBlockMatrixType KCC_Laplacian =
1164 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1165
1166 unsigned const n_integration_points =
1167 _integration_method.getNumberOfPoints();
1168
1169 std::vector<GlobalDimVectorType> ip_flux_vector;
1170 double average_velocity_norm = 0.0;
1171 if (!_process_data.non_advective_form)
1172 {
1173 ip_flux_vector.reserve(n_integration_points);
1174 }
1175
1176 auto const& b =
1177 _process_data
1178 .projected_specific_body_force_vectors[_element.getID()];
1179
1180 MaterialPropertyLib::VariableArray vars;
1181 MaterialPropertyLib::VariableArray vars_prev;
1182
1183 auto const& medium =
1184 *_process_data.media_map.getMedium(_element.getID());
1185 auto const& phase =
1186 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
1187 auto const component_id =
1188 transport_process_id - (_process_data.isothermal ? 1 : 2);
1189 auto const& component = phase.component(
1190 _transport_process_variables[component_id].get().getName());
1191
1192 auto const& Ns =
1193 _process_data.shape_matrix_cache
1194 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1195
1196 for (unsigned ip(0); ip < n_integration_points; ++ip)
1197 {
1198 auto& ip_data = _ip_data[ip];
1199 auto const& dNdx = ip_data.dNdx;
1200 auto const& w = ip_data.integration_weight;
1201 auto const& N = Ns[ip];
1202 auto& porosity = ip_data.porosity;
1203 auto const& porosity_prev = ip_data.porosity_prev;
1204
1205 ParameterLib::SpatialPosition const pos(
1206 {}, _element.getID(),
1207 MathLib::Point3d(
1208 NumLib::interpolateCoordinates<ShapeFunction,
1209 ShapeMatricesType>(_element,
1210 N)));
1211
1212 double const C_int_pt = N.dot(local_C);
1213 double const p_int_pt = N.dot(local_p);
1214 double const T_int_pt = N.dot(local_T);
1215
1216 vars.concentration = C_int_pt;
1217 vars.liquid_phase_pressure = p_int_pt;
1218 vars.temperature = T_int_pt;
1219
1220 // porosity
1221 {
1222 vars_prev.porosity = porosity_prev;
1223
1224 porosity =
1225 _process_data.chemically_induced_porosity_change
1226 ? porosity_prev
1227 : medium[MaterialPropertyLib::PropertyType::porosity]
1228 .template value<double>(vars, vars_prev, pos, t,
1229 dt);
1230
1231 vars.porosity = porosity;
1232 }
1233
1234 auto const& retardation_factor =
1235 component[MaterialPropertyLib::PropertyType::retardation_factor]
1236 .template value<double>(vars, pos, t, dt);
1237
1238 auto const& solute_dispersivity_transverse = medium.template value<
1239 double>(
1240 MaterialPropertyLib::PropertyType::transversal_dispersivity);
1241 auto const& solute_dispersivity_longitudinal =
1242 medium.template value<double>(
1243 MaterialPropertyLib::PropertyType::
1244 longitudinal_dispersivity);
1245
1246 // Use the fluid density model to compute the density
1247 auto const density =
1248 phase[MaterialPropertyLib::PropertyType::density]
1249 .template value<double>(vars, pos, t, dt);
1250 auto const decay_rate =
1251 component[MaterialPropertyLib::PropertyType::decay_rate]
1252 .template value<double>(vars, pos, t, dt);
1253
1254 auto const& pore_diffusion_coefficient =
1255 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1256 component[MaterialPropertyLib::PropertyType::pore_diffusion]
1257 .value(vars, pos, t, dt));
1258
1259 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1260 medium[MaterialPropertyLib::PropertyType::permeability].value(
1261 vars, pos, t, dt));
1262 // Use the viscosity model to compute the viscosity
1263 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1264 .template value<double>(vars, pos, t, dt);
1265
1266 GlobalDimMatrixType const K_over_mu = K / mu;
1267 GlobalDimVectorType const velocity =
1268 _process_data.has_gravity
1269 ? GlobalDimVectorType(-K_over_mu *
1270 (dNdx * local_p - density * b))
1271 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1272
1273 GlobalDimMatrixType const hydrodynamic_dispersion =
1274 NumLib::computeHydrodynamicDispersion(
1275 _process_data.stabilizer, _element.getID(),
1276 pore_diffusion_coefficient, velocity, porosity,
1277 solute_dispersivity_transverse,
1278 solute_dispersivity_longitudinal);
1279
1280 double const R_times_phi = retardation_factor * porosity;
1281 auto const N_t_N = (N.transpose() * N).eval();
1282
1283 if (_process_data.non_advective_form)
1284 {
1285 const double drho_dC =
1286 phase[MaterialPropertyLib::PropertyType::density]
1287 .template dValue<double>(
1288 vars, MaterialPropertyLib::Variable::concentration,
1289 pos, t, dt);
1290 local_M.noalias() +=
1291 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1292 }
1293
1294 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1295
1296 // coupling term
1297 if (_process_data.non_advective_form)
1298 {
1299 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1300
1301 const double drho_dp =
1302 phase[MaterialPropertyLib::PropertyType::density]
1303 .template dValue<double>(vars,
1304 MaterialPropertyLib::Variable::
1305 liquid_phase_pressure,
1306 pos, t, dt);
1307
1308 local_K.noalias() +=
1309 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1310 dNdx.transpose() * velocity * N * (density * w);
1311 }
1312 else
1313 {
1314 ip_flux_vector.emplace_back(velocity * density);
1315 average_velocity_norm += velocity.norm();
1316 }
1317 local_K.noalias() +=
1318 N_t_N * (decay_rate * R_times_phi * density * w);
1319
1320 KCC_Laplacian.noalias() += dNdx.transpose() *
1321 hydrodynamic_dispersion * dNdx *
1322 (density * w);
1323 }
1324
1325 if (!_process_data.non_advective_form)
1326 {
1327 NumLib::assembleAdvectionMatrix<
1328 typename ShapeFunction::MeshElement>(
1329 _process_data.stabilizer, _ip_data,
1330 _process_data.shape_matrix_cache, ip_flux_vector,
1331 average_velocity_norm /
1332 static_cast<double>(n_integration_points),
1333 KCC_Laplacian);
1334 }
1335 local_K.noalias() += KCC_Laplacian;
1336 }
ProcessLib::ComponentTransport::LocalAssemblerData::getLocalTemperature
NodalVectorType getLocalTemperature(double const t, Eigen::VectorXd const &local_x) const
Definition ComponentTransportFEM.h:2385
ProcessLib::ComponentTransport::LocalAssemblerData::NodalVectorType
typename ShapeMatricesType::NodalVectorType NodalVectorType
Definition ComponentTransportFEM.h:231

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::detail::assembleAdvectionMatrix(), NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getLocalTemperature(), getName(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, MaterialPropertyLib::retardation_factor, MaterialPropertyLib::VariableArray::temperature, temperature_size, MaterialPropertyLib::transversal_dispersivity, and MaterialPropertyLib::viscosity.

Referenced by assembleForStaggeredScheme().

◆ assembleForStaggeredScheme()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleForStaggeredScheme ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
int const process_id,
std::vector< double > & local_M_data,
std::vector< double > & local_K_data,
std::vector< double > & local_b_data )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 828 of file ComponentTransportFEM.h.

835 {
836 if (process_id == _process_data.hydraulic_process_id)
837 {
838 assembleHydraulicEquation(t, dt, local_x, local_x_prev,
839 local_M_data, local_K_data, local_b_data);
840 }
841 else if (process_id == _process_data.thermal_process_id)
842 {
843 assembleHeatTransportEquation(t, dt, local_x, local_x_prev,
844 local_M_data, local_K_data,
845 local_b_data);
846 }
847 else
848 {
849 // Go for assembling in an order of transport process id.
850 assembleComponentTransportEquation(t, dt, local_x, local_x_prev,
851 local_M_data, local_K_data,
852 local_b_data, process_id);
853 }
854 }
ProcessLib::ComponentTransport::LocalAssemblerData::assembleHeatTransportEquation
void assembleHeatTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &)
Definition ComponentTransportFEM.h:996
ProcessLib::ComponentTransport::LocalAssemblerData::assembleHydraulicEquation
void assembleHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
Definition ComponentTransportFEM.h:856
ProcessLib::ComponentTransport::LocalAssemblerData::assembleComponentTransportEquation
void assembleComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &, int const transport_process_id)
Definition ComponentTransportFEM.h:1134

References _process_data, assembleComponentTransportEquation(), assembleHeatTransportEquation(), and assembleHydraulicEquation().

◆ assembleHeatTransportEquation()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleHeatTransportEquation ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & ,
std::vector< double > & local_M_data,
std::vector< double > & local_K_data,
std::vector< double > &  )
inline

Definition at line 996 of file ComponentTransportFEM.h.

1002 {
1003 // In the staggered HTC process, number of components might be non-zero.
1004 assert(local_x.size() ==
1005 pressure_size + temperature_size +
1006 concentration_size *
1007 static_cast<int>(_transport_process_variables.size()));
1008
1009 auto const local_p =
1010 local_x.template segment<pressure_size>(pressure_index);
1011 auto const local_T = getLocalTemperature(t, local_x);
1012 auto const local_C = local_x.template segment<concentration_size>(
1013 first_concentration_index);
1014
1015 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1016 local_M_data, temperature_size, temperature_size);
1017 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1018 local_K_data, temperature_size, temperature_size);
1019
1020 auto const& process_data = this->_process_data;
1021 auto const& medium =
1022 *process_data.media_map.getMedium(this->_element.getID());
1023 auto const& liquid_phase =
1024 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
1025
1026 auto const& b =
1027 _process_data
1028 .projected_specific_body_force_vectors[_element.getID()];
1029
1030 MaterialPropertyLib::VariableArray vars;
1031
1032 unsigned const n_integration_points =
1033 this->_integration_method.getNumberOfPoints();
1034
1035 std::vector<GlobalDimVectorType> ip_flux_vector;
1036 double average_velocity_norm = 0.0;
1037 ip_flux_vector.reserve(n_integration_points);
1038
1039 auto const& Ns =
1040 _process_data.shape_matrix_cache
1041 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1042
1043 for (unsigned ip(0); ip < n_integration_points; ip++)
1044 {
1045 auto const& ip_data = this->_ip_data[ip];
1046 auto const& dNdx = ip_data.dNdx;
1047 auto const& w = ip_data.integration_weight;
1048 auto const& N = Ns[ip];
1049
1050 ParameterLib::SpatialPosition const pos(
1051 {}, this->_element.getID(),
1052 MathLib::Point3d(
1053 NumLib::interpolateCoordinates<ShapeFunction,
1054 ShapeMatricesType>(_element,
1055 N)));
1056
1057 double p_at_xi = 0.;
1058 NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi);
1059 double T_at_xi = 0.;
1060 NumLib::shapeFunctionInterpolate(local_T, N, T_at_xi);
1061 double const C_int_pt = N.dot(local_C);
1062
1063 vars.temperature = T_at_xi;
1064 vars.liquid_phase_pressure = p_at_xi;
1065
1066 vars.liquid_saturation = 1.0;
1067
1068 auto const porosity =
1069 medium.property(MaterialPropertyLib::PropertyType::porosity)
1070 .template value<double>(vars, pos, t, dt);
1071 vars.porosity = porosity;
1072 vars.concentration = C_int_pt;
1073
1074 // Use the fluid density model to compute the density
1075 auto const fluid_density =
1076 liquid_phase
1077 .property(MaterialPropertyLib::PropertyType::density)
1078 .template value<double>(vars, pos, t, dt);
1079 vars.density = fluid_density;
1080 auto const specific_heat_capacity_fluid =
1081 liquid_phase
1082 .property(MaterialPropertyLib::specific_heat_capacity)
1083 .template value<double>(vars, pos, t, dt);
1084
1085 // Assemble mass matrix
1086 local_M.noalias() +=
1087 N.transpose() * N *
1088 (this->getHeatEnergyCoefficient(vars, porosity, fluid_density,
1089 specific_heat_capacity_fluid,
1090 pos, t, dt) *
1091 w);
1092
1093 // Assemble Laplace matrix
1094 auto const viscosity =
1095 liquid_phase
1096 .property(MaterialPropertyLib::PropertyType::viscosity)
1097 .template value<double>(vars, pos, t, dt);
1098
1099 auto const intrinsic_permeability =
1100 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1101 medium
1102 .property(
1103 MaterialPropertyLib::PropertyType::permeability)
1104 .value(vars, pos, t, dt));
1105
1106 GlobalDimMatrixType const K_over_mu =
1107 intrinsic_permeability / viscosity;
1108 GlobalDimVectorType const velocity =
1109 process_data.has_gravity
1110 ? GlobalDimVectorType(-K_over_mu *
1111 (dNdx * local_p - fluid_density * b))
1112 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1113
1114 GlobalDimMatrixType const thermal_conductivity_dispersivity =
1115 this->getThermalConductivityDispersivity(
1116 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1117 pos, t, dt);
1118
1119 local_K.noalias() +=
1120 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1121
1122 ip_flux_vector.emplace_back(velocity * fluid_density *
1123 specific_heat_capacity_fluid);
1124 average_velocity_norm += velocity.norm();
1125 }
1126
1127 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1128 process_data.stabilizer, this->_ip_data,
1129 _process_data.shape_matrix_cache, ip_flux_vector,
1130 average_velocity_norm / static_cast<double>(n_integration_points),
1131 local_K);
1132 }
ProcessLib::ComponentTransport::LocalAssemblerData::getThermalConductivityDispersivity
GlobalDimMatrixType getThermalConductivityDispersivity(MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
Definition ComponentTransportFEM.h:2344
ProcessLib::ComponentTransport::LocalAssemblerData::getHeatEnergyCoefficient
double getHeatEnergyCoefficient(MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
Definition ComponentTransportFEM.h:2319
NumLib::detail::assembleAdvectionMatrix
void assembleAdvectionMatrix(IPData const &ip_data_vector, NumLib::ShapeMatrixCache const &shape_matrix_cache, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
Definition AdvectionMatrixAssembler.h:22

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::detail::assembleAdvectionMatrix(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getHeatEnergyCoefficient(), getLocalTemperature(), getThermalConductivityDispersivity(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::specific_heat_capacity, MaterialPropertyLib::VariableArray::temperature, temperature_size, and MaterialPropertyLib::viscosity.

Referenced by assembleForStaggeredScheme().

◆ assembleHydraulicEquation()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleHydraulicEquation ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
std::vector< double > & local_M_data,
std::vector< double > & local_K_data,
std::vector< double > & local_b_data )
inline

Definition at line 856 of file ComponentTransportFEM.h.

863 {
864 auto const local_p =
865 local_x.template segment<pressure_size>(pressure_index);
866 auto const local_C = local_x.template segment<concentration_size>(
867 first_concentration_index);
868 auto const local_C_prev =
869 local_x_prev.segment<concentration_size>(first_concentration_index);
870
871 NodalVectorType local_T = getLocalTemperature(t, local_x);
872
873 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
874 local_M_data, pressure_size, pressure_size);
875 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
876 local_K_data, pressure_size, pressure_size);
877 auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
878 local_b_data, pressure_size);
879
880 unsigned const n_integration_points =
881 _integration_method.getNumberOfPoints();
882
883 auto const& b =
884 _process_data
885 .projected_specific_body_force_vectors[_element.getID()];
886
887 auto const& medium =
888 *_process_data.media_map.getMedium(_element.getID());
889 auto const& phase =
890 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
891
892 MaterialPropertyLib::VariableArray vars;
893 MaterialPropertyLib::VariableArray vars_prev;
894
895 auto const& Ns =
896 _process_data.shape_matrix_cache
897 .NsHigherOrder<typename ShapeFunction::MeshElement>();
898
899 for (unsigned ip(0); ip < n_integration_points; ++ip)
900 {
901 auto& ip_data = _ip_data[ip];
902 auto const& dNdx = ip_data.dNdx;
903 auto const& w = ip_data.integration_weight;
904 auto const& N = Ns[ip];
905 auto& porosity = ip_data.porosity;
906 auto const& porosity_prev = ip_data.porosity_prev;
907
908 ParameterLib::SpatialPosition const pos(
909 {}, _element.getID(),
910 MathLib::Point3d(
911 NumLib::interpolateCoordinates<ShapeFunction,
912 ShapeMatricesType>(_element,
913 N)));
914
915 double const C_int_pt = N.dot(local_C);
916 double const p_int_pt = N.dot(local_p);
917 double const T_int_pt = N.dot(local_T);
918
919 vars.concentration = C_int_pt;
920 vars.liquid_phase_pressure = p_int_pt;
921 vars.temperature = T_int_pt;
922
923 // porosity
924 {
925 vars_prev.porosity = porosity_prev;
926
927 porosity =
928 _process_data.chemically_induced_porosity_change
929 ? porosity_prev
930 : medium[MaterialPropertyLib::PropertyType::porosity]
931 .template value<double>(vars, vars_prev, pos, t,
932 dt);
933
934 vars.porosity = porosity;
935 }
936
937 // Use the fluid density model to compute the density
938 // TODO: Concentration of which component as one of arguments for
939 // calculation of fluid density
940 auto const density =
941 phase[MaterialPropertyLib::PropertyType::density]
942 .template value<double>(vars, pos, t, dt);
943
944 double storage = 0;
945 if (medium.hasProperty(MaterialPropertyLib::PropertyType::storage))
946 {
947 storage = medium[MaterialPropertyLib::PropertyType::storage]
948 .template value<double>(vars, pos, t, dt);
949 }
950
951 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
952 medium[MaterialPropertyLib::PropertyType::permeability].value(
953 vars, pos, t, dt));
954
955 // Use the viscosity model to compute the viscosity
956 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
957 .template value<double>(vars, pos, t, dt);
958
959 GlobalDimMatrixType const K_over_mu = K / mu;
960
961 const double drho_dp =
962 phase[MaterialPropertyLib::PropertyType::density]
963 .template dValue<double>(
964 vars,
965 MaterialPropertyLib::Variable::liquid_phase_pressure,
966 pos, t, dt);
967 const double drho_dC =
968 phase[MaterialPropertyLib::PropertyType::density]
969 .template dValue<double>(
970 vars, MaterialPropertyLib::Variable::concentration, pos,
971 t, dt);
972
973 // matrix assembly
974 local_M.noalias() +=
975 N.transpose() * N *
976 (porosity * drho_dp * w + density * storage * w);
977 local_K.noalias() +=
978 w * dNdx.transpose() * density * K_over_mu * dNdx;
979
980 if (_process_data.has_gravity)
981 {
982 local_b.noalias() +=
983 w * density * density * dNdx.transpose() * K_over_mu * b;
984 }
985
986 // coupling term
987 {
988 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
989
990 local_b.noalias() -=
991 N.transpose() * (porosity * drho_dC * C_dot * w);
992 }
993 }
994 }

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getLocalTemperature(), NumLib::interpolateCoordinates(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, MaterialPropertyLib::storage, MaterialPropertyLib::VariableArray::temperature, and MaterialPropertyLib::viscosity.

Referenced by assembleForStaggeredScheme().

◆ assembleKCmCn()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleKCmCn ( int const component_id,
double const t,
double const dt,
Eigen::Ref< LocalBlockMatrixType > KCmCn,
double const stoichiometric_coefficient,
double const kinetic_prefactor )
inline

Definition at line 776 of file ComponentTransportFEM.h.

780 {
781 unsigned const n_integration_points =
782 _integration_method.getNumberOfPoints();
783
784 MaterialPropertyLib::VariableArray vars;
785
786 auto const& medium =
787 *_process_data.media_map.getMedium(_element.getID());
788 auto const& phase =
789 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
790 auto const& component = phase.component(
791 _transport_process_variables[component_id].get().getName());
792
793 auto const& Ns =
794 _process_data.shape_matrix_cache
795 .NsHigherOrder<typename ShapeFunction::MeshElement>();
796
797 for (unsigned ip(0); ip < n_integration_points; ++ip)
798 {
799 auto& ip_data = _ip_data[ip];
800 auto const& w = ip_data.integration_weight;
801 auto const& N = Ns[ip];
802 auto& porosity = ip_data.porosity;
803
804 ParameterLib::SpatialPosition const pos(
805 {}, _element.getID(),
806 MathLib::Point3d(
807 NumLib::interpolateCoordinates<ShapeFunction,
808 ShapeMatricesType>(_element,
809 N)));
810
811 auto const retardation_factor =
812 component[MaterialPropertyLib::PropertyType::retardation_factor]
813 .template value<double>(vars, pos, t, dt);
814
815 porosity = medium[MaterialPropertyLib::PropertyType::porosity]
816 .template value<double>(vars, pos, t, dt);
817
818 auto const density =
819 phase[MaterialPropertyLib::PropertyType::density]
820 .template value<double>(vars, pos, t, dt);
821
822 KCmCn.noalias() -= N.transpose() * N *
823 (stoichiometric_coefficient * kinetic_prefactor *
824 retardation_factor * porosity * density * w);
825 }
826 }

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::density, getName(), NumLib::interpolateCoordinates(), MaterialPropertyLib::porosity, and MaterialPropertyLib::retardation_factor.

Referenced by assemble().

◆ assembleReactionEquationConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleReactionEquationConcrete ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
std::vector< double > & local_M_data,
std::vector< double > & local_K_data,
std::vector< double > & local_b_data,
int const transport_process_id )
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 1617 of file ComponentTransportFEM.h.

1622 {
1623 auto const local_C = local_x.template segment<concentration_size>(
1624 first_concentration_index +
1625 (transport_process_id - 1) * concentration_size);
1626
1627 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1628 local_M_data, concentration_size, concentration_size);
1629 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1630 local_K_data, concentration_size, concentration_size);
1631 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1632 local_b_data, concentration_size);
1633
1634 unsigned const n_integration_points =
1635 _integration_method.getNumberOfPoints();
1636
1637 MaterialPropertyLib::VariableArray vars;
1638 MaterialPropertyLib::VariableArray vars_prev;
1639
1640 auto const& medium =
1641 *_process_data.media_map.getMedium(_element.getID());
1642 auto const component_id = transport_process_id - 1;
1643
1644 auto const& Ns =
1645 _process_data.shape_matrix_cache
1646 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1647
1648 for (unsigned ip(0); ip < n_integration_points; ++ip)
1649 {
1650 auto& ip_data = _ip_data[ip];
1651 auto const w = ip_data.integration_weight;
1652 auto const& N = Ns[ip];
1653 auto& porosity = ip_data.porosity;
1654 auto const& porosity_prev = ip_data.porosity_prev;
1655 auto const chemical_system_id = ip_data.chemical_system_id;
1656
1657 ParameterLib::SpatialPosition const pos(
1658 {}, _element.getID(),
1659 MathLib::Point3d(
1660 NumLib::interpolateCoordinates<ShapeFunction,
1661 ShapeMatricesType>(_element,
1662 N)));
1663
1664 double C_int_pt = 0.0;
1665 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
1666
1667 vars.concentration = C_int_pt;
1668
1669 auto const porosity_dot = (porosity - porosity_prev) / dt;
1670
1671 // porosity
1672 {
1673 vars_prev.porosity = porosity_prev;
1674
1675 porosity =
1676 _process_data.chemically_induced_porosity_change
1677 ? porosity_prev
1678 : medium[MaterialPropertyLib::PropertyType::porosity]
1679 .template value<double>(vars, vars_prev, pos, t,
1680 dt);
1681 }
1682
1683 local_M.noalias() += w * N.transpose() * porosity * N;
1684
1685 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1686
1687 if (chemical_system_id == -1)
1688 {
1689 continue;
1690 }
1691
1692 auto const C_post_int_pt =
1693 _process_data.chemical_solver_interface->getConcentration(
1694 component_id, chemical_system_id);
1695
1696 local_b.noalias() += N.transpose() * ((C_post_int_pt - C_int_pt) /
1697 dt * porosity * w);
1698 }
1699 }

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), first_concentration_index, NumLib::interpolateCoordinates(), MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, and NumLib::detail::shapeFunctionInterpolate().

◆ assembleWithJacobianComponentTransportEquation()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleWithJacobianComponentTransportEquation ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
std::vector< double > & local_b_data,
std::vector< double > & local_Jac_data,
int const component_id )
inline

Definition at line 1472 of file ComponentTransportFEM.h.

1476 {
1477 auto const concentration_index =
1478 first_concentration_index + component_id * concentration_size;
1479
1480 auto const p = local_x.template segment<pressure_size>(pressure_index);
1481 auto const c =
1482 local_x.template segment<concentration_size>(concentration_index);
1483 auto const c_prev =
1484 local_x_prev.segment<concentration_size>(concentration_index);
1485
1486 auto const local_T = getLocalTemperature(t, local_x);
1487
1488 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1489 local_Jac_data, concentration_size, concentration_size);
1490 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1491 local_b_data, concentration_size);
1492
1493 LocalBlockMatrixType KCC_Laplacian =
1494 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1495
1496 unsigned const n_integration_points =
1497 _integration_method.getNumberOfPoints();
1498
1499 std::vector<GlobalDimVectorType> ip_flux_vector;
1500 double average_velocity_norm = 0.0;
1501 ip_flux_vector.reserve(n_integration_points);
1502
1503 auto const& b =
1504 _process_data
1505 .projected_specific_body_force_vectors[_element.getID()];
1506
1507 MaterialPropertyLib::VariableArray vars;
1508 MaterialPropertyLib::VariableArray vars_prev;
1509
1510 auto const& medium =
1511 *_process_data.media_map.getMedium(_element.getID());
1512 auto const& phase =
1513 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
1514 auto const& component = phase.component(
1515 _transport_process_variables[component_id].get().getName());
1516
1517 auto const& Ns =
1518 _process_data.shape_matrix_cache
1519 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1520
1521 for (unsigned ip(0); ip < n_integration_points; ++ip)
1522 {
1523 auto& ip_data = _ip_data[ip];
1524 auto const& dNdx = ip_data.dNdx;
1525 auto const& w = ip_data.integration_weight;
1526 auto const& N = Ns[ip];
1527 auto& phi = ip_data.porosity;
1528 auto const& phi_prev = ip_data.porosity_prev;
1529
1530 ParameterLib::SpatialPosition const pos(
1531 {}, _element.getID(),
1532 MathLib::Point3d(
1533 NumLib::interpolateCoordinates<ShapeFunction,
1534 ShapeMatricesType>(_element,
1535 N)));
1536
1537 double const p_ip = N.dot(p);
1538 double const c_ip = N.dot(c);
1539
1540 vars.liquid_phase_pressure = p_ip;
1541 vars.concentration = c_ip;
1542 vars.temperature = N.dot(local_T);
1543
1544 // porosity
1545 {
1546 vars_prev.porosity = phi_prev;
1547
1548 phi = _process_data.chemically_induced_porosity_change
1549 ? phi_prev
1550 : medium[MaterialPropertyLib::PropertyType::porosity]
1551 .template value<double>(vars, vars_prev, pos, t,
1552 dt);
1553
1554 vars.porosity = phi;
1555 }
1556
1557 auto const R =
1558 component[MaterialPropertyLib::PropertyType::retardation_factor]
1559 .template value<double>(vars, pos, t, dt);
1560
1561 auto const alpha_T = medium.template value<double>(
1562 MaterialPropertyLib::PropertyType::transversal_dispersivity);
1563 auto const alpha_L = medium.template value<double>(
1564 MaterialPropertyLib::PropertyType::longitudinal_dispersivity);
1565
1566 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1567 .template value<double>(vars, pos, t, dt);
1568 // first-order decay constant
1569 auto const alpha =
1570 component[MaterialPropertyLib::PropertyType::decay_rate]
1571 .template value<double>(vars, pos, t, dt);
1572
1573 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1574 component[MaterialPropertyLib::PropertyType::pore_diffusion]
1575 .value(vars, pos, t, dt));
1576
1577 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1578 medium[MaterialPropertyLib::PropertyType::permeability].value(
1579 vars, pos, t, dt));
1580 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1581 .template value<double>(vars, pos, t, dt);
1582 // Darcy flux
1583 GlobalDimVectorType const q =
1584 _process_data.has_gravity
1585 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1586 : GlobalDimVectorType(-k / mu * dNdx * p);
1587
1588 GlobalDimMatrixType const D = NumLib::computeHydrodynamicDispersion(
1589 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1590 alpha_L);
1591
1592 // matrix assembly
1593 local_Jac.noalias() +=
1594 N.transpose() * N * (rho * phi * R * (alpha + 1 / dt) * w);
1595
1596 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1597
1598 auto const cdot = (c - c_prev) / dt;
1599 local_rhs.noalias() -=
1600 N.transpose() * N * (cdot + alpha * c) * (rho * phi * R * w);
1601
1602 ip_flux_vector.emplace_back(q * rho);
1603 average_velocity_norm += q.norm();
1604 }
1605
1606 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1607 _process_data.stabilizer, _ip_data,
1608 _process_data.shape_matrix_cache, ip_flux_vector,
1609 average_velocity_norm / static_cast<double>(n_integration_points),
1610 KCC_Laplacian);
1611
1612 local_rhs.noalias() -= KCC_Laplacian * c;
1613
1614 local_Jac.noalias() += KCC_Laplacian;
1615 }

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::detail::assembleAdvectionMatrix(), NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getLocalTemperature(), getName(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::longitudinal_dispersivity, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, MaterialPropertyLib::retardation_factor, MaterialPropertyLib::VariableArray::temperature, MaterialPropertyLib::transversal_dispersivity, and MaterialPropertyLib::viscosity.

Referenced by assembleWithJacobianForStaggeredScheme().

◆ assembleWithJacobianForStaggeredScheme()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleWithJacobianForStaggeredScheme ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
int const process_id,
std::vector< double > & local_b_data,
std::vector< double > & local_Jac_data )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1338 of file ComponentTransportFEM.h.

1343 {
1344 if (process_id == _process_data.hydraulic_process_id)
1345 {
1346 assembleWithJacobianHydraulicEquation(t, dt, local_x, local_x_prev,
1347 local_b_data, local_Jac_data);
1348 }
1349 else
1350 {
1351 int const component_id = process_id - 1;
1352 assembleWithJacobianComponentTransportEquation(
1353 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1354 component_id);
1355 }
1356 }
ProcessLib::ComponentTransport::LocalAssemblerData::assembleWithJacobianHydraulicEquation
void assembleWithJacobianHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
Definition ComponentTransportFEM.h:1358
ProcessLib::ComponentTransport::LocalAssemblerData::assembleWithJacobianComponentTransportEquation
void assembleWithJacobianComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, int const component_id)
Definition ComponentTransportFEM.h:1472

References _process_data, assembleWithJacobianComponentTransportEquation(), and assembleWithJacobianHydraulicEquation().

◆ assembleWithJacobianHydraulicEquation()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::assembleWithJacobianHydraulicEquation ( double const t,
double const dt,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const & local_x_prev,
std::vector< double > & local_b_data,
std::vector< double > & local_Jac_data )
inline

Definition at line 1358 of file ComponentTransportFEM.h.

1362 {
1363 auto const p = local_x.template segment<pressure_size>(pressure_index);
1364 auto const c = local_x.template segment<concentration_size>(
1365 first_concentration_index);
1366
1367 auto const p_prev = local_x_prev.segment<pressure_size>(pressure_index);
1368 auto const c_prev =
1369 local_x_prev.segment<concentration_size>(first_concentration_index);
1370
1371 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1372 local_Jac_data, pressure_size, pressure_size);
1373 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1374 local_b_data, pressure_size);
1375
1376 unsigned const n_integration_points =
1377 _integration_method.getNumberOfPoints();
1378
1379 auto const& b =
1380 _process_data
1381 .projected_specific_body_force_vectors[_element.getID()];
1382
1383 auto const& medium =
1384 *_process_data.media_map.getMedium(_element.getID());
1385 auto const& phase =
1386 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
1387
1388 MaterialPropertyLib::VariableArray vars;
1389 MaterialPropertyLib::VariableArray vars_prev;
1390
1391 auto const& Ns =
1392 _process_data.shape_matrix_cache
1393 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1394
1395 for (unsigned ip(0); ip < n_integration_points; ++ip)
1396 {
1397 auto& ip_data = _ip_data[ip];
1398 auto const& dNdx = ip_data.dNdx;
1399 auto const& w = ip_data.integration_weight;
1400 auto const& N = Ns[ip];
1401 auto& phi = ip_data.porosity;
1402 auto const& phi_prev = ip_data.porosity_prev;
1403
1404 ParameterLib::SpatialPosition const pos(
1405 {}, _element.getID(),
1406 MathLib::Point3d(
1407 NumLib::interpolateCoordinates<ShapeFunction,
1408 ShapeMatricesType>(_element,
1409 N)));
1410
1411 double const p_ip = N.dot(p);
1412 double const c_ip = N.dot(c);
1413
1414 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1415
1416 vars.liquid_phase_pressure = p_ip;
1417 vars.concentration = c_ip;
1418
1419 // porosity
1420 {
1421 vars_prev.porosity = phi_prev;
1422
1423 phi = _process_data.chemically_induced_porosity_change
1424 ? phi_prev
1425 : medium[MaterialPropertyLib::PropertyType::porosity]
1426 .template value<double>(vars, vars_prev, pos, t,
1427 dt);
1428
1429 vars.porosity = phi;
1430 }
1431
1432 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1433 .template value<double>(vars, pos, t, dt);
1434
1435 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1436 medium[MaterialPropertyLib::PropertyType::permeability].value(
1437 vars, pos, t, dt));
1438
1439 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1440 .template value<double>(vars, pos, t, dt);
1441
1442 auto const drho_dp =
1443 phase[MaterialPropertyLib::PropertyType::density]
1444 .template dValue<double>(
1445 vars,
1446 MaterialPropertyLib::Variable::liquid_phase_pressure,
1447 pos, t, dt);
1448 auto const drho_dc =
1449 phase[MaterialPropertyLib::PropertyType::density]
1450 .template dValue<double>(
1451 vars, MaterialPropertyLib::Variable::concentration, pos,
1452 t, dt);
1453
1454 // matrix assembly
1455 local_Jac.noalias() +=
1456 N.transpose() * N * (phi * drho_dp / dt * w) +
1457 w * dNdx.transpose() * rho * k / mu * dNdx;
1458
1459 local_rhs.noalias() -=
1460 N.transpose() * (drho_dp * N * p_prev + drho_dc * cdot_ip) *
1461 (phi * w) +
1462 dNdx.transpose() * k / mu * dNdx * p * (rho * w);
1463
1464 if (_process_data.has_gravity)
1465 {
1466 local_rhs.noalias() +=
1467 w * rho * dNdx.transpose() * k / mu * rho * b;
1468 }
1469 }
1470 }

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), NumLib::interpolateCoordinates(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, and MaterialPropertyLib::viscosity.

Referenced by assembleWithJacobianForStaggeredScheme().

◆ calculateIntPtDarcyVelocity()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::calculateIntPtDarcyVelocity ( const double t,
Eigen::Ref< const NodalVectorType > const & p_nodal_values,
Eigen::Ref< const NodalVectorType > const & C_nodal_values,
Eigen::Ref< const NodalVectorType > const & T_nodal_values,
std::vector< double > & cache ) const
inline

Definition at line 1907 of file ComponentTransportFEM.h.

1913 {
1914 auto const n_integration_points =
1915 _integration_method.getNumberOfPoints();
1916
1917 cache.clear();
1918 auto cache_mat = MathLib::createZeroedMatrix<
1919 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1920 cache, GlobalDim, n_integration_points);
1921
1922 auto const& b =
1923 _process_data
1924 .projected_specific_body_force_vectors[_element.getID()];
1925
1926 MaterialPropertyLib::VariableArray vars;
1927
1928 auto const& medium =
1929 *_process_data.media_map.getMedium(_element.getID());
1930 auto const& phase =
1931 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
1932
1933 auto const& Ns =
1934 _process_data.shape_matrix_cache
1935 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1936
1937 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1938 {
1939 auto const& ip_data = _ip_data[ip];
1940 auto const& dNdx = ip_data.dNdx;
1941 auto const& N = Ns[ip];
1942 auto const& porosity = ip_data.porosity;
1943
1944 ParameterLib::SpatialPosition const pos(
1945 {}, _element.getID(),
1946 MathLib::Point3d(
1947 NumLib::interpolateCoordinates<ShapeFunction,
1948 ShapeMatricesType>(_element,
1949 N)));
1950
1951 double C_int_pt = 0.0;
1952 double p_int_pt = 0.0;
1953 double T_int_pt = 0.0;
1954
1955 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1956 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1957 NumLib::shapeFunctionInterpolate(T_nodal_values, N, T_int_pt);
1958
1959 vars.concentration = C_int_pt;
1960 vars.liquid_phase_pressure = p_int_pt;
1961 vars.porosity = porosity;
1962 vars.temperature = T_int_pt;
1963
1964 // TODO (naumov) Temporary value not used by current material
1965 // models. Need extension of secondary variables interface.
1966 double const dt = std::numeric_limits<double>::quiet_NaN();
1967 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1968 medium[MaterialPropertyLib::PropertyType::permeability].value(
1969 vars, pos, t, dt));
1970 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1971 .template value<double>(vars, pos, t, dt);
1972 GlobalDimMatrixType const K_over_mu = K / mu;
1973
1974 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1975 if (_process_data.has_gravity)
1976 {
1977 auto const rho_w =
1978 phase[MaterialPropertyLib::PropertyType::density]
1979 .template value<double>(vars, pos, t, dt);
1980 // here it is assumed that the vector b is directed 'downwards'
1981 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1982 }
1983 }
1984
1985 return cache;
1986 }

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::VariableArray::concentration, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::formEigenTensor(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::VariableArray::porosity, NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::temperature, and MaterialPropertyLib::viscosity.

Referenced by computeSecondaryVariableConcrete(), and getIntPtDarcyVelocity().

◆ calculateIntPtLiquidDensity()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::calculateIntPtLiquidDensity ( const double t,
Eigen::Ref< const NodalVectorType > const & p_nodal_values,
Eigen::Ref< const NodalVectorType > const & C_nodal_values,
Eigen::Ref< const NodalVectorType > const & T_nodal_values,
std::vector< double > & cache ) const
inline

Definition at line 1776 of file ComponentTransportFEM.h.

1782 {
1783 auto const n_integration_points =
1784 _integration_method.getNumberOfPoints();
1785
1786 cache.clear();
1787 auto cache_vec = MathLib::createZeroedVector<Eigen::VectorXd>(
1788 cache, n_integration_points);
1789
1790 ParameterLib::SpatialPosition pos;
1791 pos.setElementID(_element.getID());
1792
1793 MaterialPropertyLib::VariableArray vars;
1794
1795 auto const& medium =
1796 *_process_data.media_map.getMedium(_element.getID());
1797 auto const& phase =
1798 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
1799
1800 auto const& Ns =
1801 _process_data.shape_matrix_cache
1802 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1803
1804 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1805 {
1806 auto const& N = Ns[ip];
1807
1808 double C_int_pt = 0.0;
1809 double p_int_pt = 0.0;
1810 double T_int_pt = 0.0;
1811
1812 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1813 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1814 NumLib::shapeFunctionInterpolate(T_nodal_values, N, T_int_pt);
1815
1816 vars.concentration = C_int_pt;
1817 vars.liquid_phase_pressure = p_int_pt;
1818 vars.temperature = T_int_pt;
1819
1820 // TODO (naumov) Temporary value not used by current material
1821 // models. Need extension of secondary variables interface.
1822 double const dt = std::numeric_limits<double>::quiet_NaN();
1823
1824 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
1825 .template value<double>(vars, pos, t, dt);
1826 cache_vec[ip] = rho_w;
1827 }
1828
1829 return cache;
1830 }

References _element, _integration_method, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::VariableArray::concentration, MathLib::createZeroedVector(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::liquid_phase_pressure, ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), and MaterialPropertyLib::VariableArray::temperature.

Referenced by getIntPtLiquidDensity().

◆ computeReactionRelatedSecondaryVariable()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::computeReactionRelatedSecondaryVariable ( std::size_t const ele_id)
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 2143 of file ComponentTransportFEM.h.

2145 {
2146 auto const n_integration_points =
2147 _integration_method.getNumberOfPoints();
2148
2149 if (_process_data.chemically_induced_porosity_change)
2150 {
2151 auto const& medium = *_process_data.media_map.getMedium(ele_id);
2152
2153 for (auto& ip_data : _ip_data)
2154 {
2155 ip_data.porosity = ip_data.porosity_prev;
2156
2157 _process_data.chemical_solver_interface
2158 ->updatePorosityPostReaction(ip_data.chemical_system_id,
2159 medium, ip_data.porosity);
2160 }
2161
2162 updateAveragePorosity(ele_id);
2163 }
2164
2165 std::vector<GlobalIndexType> chemical_system_indices;
2166 chemical_system_indices.reserve(n_integration_points);
2167 std::transform(_ip_data.begin(), _ip_data.end(),
2168 std::back_inserter(chemical_system_indices),
2169 [](auto const& ip_data)
2170 { return ip_data.chemical_system_id; });
2171
2172 _process_data.chemical_solver_interface->computeSecondaryVariable(
2173 ele_id, chemical_system_indices);
2174 }

References _integration_method, _ip_data, _process_data, and updateAveragePorosity().

◆ computeSecondaryVariableConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::computeSecondaryVariableConcrete ( double const t,
double const ,
Eigen::VectorXd const & local_x,
Eigen::VectorXd const &  )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 2065 of file ComponentTransportFEM.h.

2070 {
2071 auto const local_p =
2072 local_x.template segment<pressure_size>(pressure_index);
2073 auto const local_C = local_x.template segment<concentration_size>(
2074 first_concentration_index);
2075 auto const local_T = getLocalTemperature(t, local_x);
2076
2077 std::vector<double> ele_velocity;
2078 calculateIntPtDarcyVelocity(t, local_p, local_C, local_T, ele_velocity);
2079
2080 auto const n_integration_points =
2081 _integration_method.getNumberOfPoints();
2082 auto const ele_velocity_mat =
2083 MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
2084
2085 auto const ele_id = _element.getID();
2086 Eigen::Map<LocalVectorType>(
2087 &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
2088 GlobalDim) =
2089 ele_velocity_mat.rowwise().sum() / n_integration_points;
2090
2091 ParameterLib::SpatialPosition pos;
2092 pos.setElementID(ele_id);
2093
2094 MaterialPropertyLib::VariableArray vars;
2095 auto const& medium = *_process_data.media_map.getMedium(ele_id);
2096 auto const& Ns =
2097 _process_data.shape_matrix_cache
2098 .NsHigherOrder<typename ShapeFunction::MeshElement>();
2099
2100 double permeability_avg = 0.0;
2101 for (unsigned ip = 0; ip < n_integration_points; ++ip)
2102 {
2103 auto const& ip_data = _ip_data[ip];
2104 auto const& N = Ns[ip];
2105
2106 double C_int_pt = 0.0;
2107 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
2108
2109 double p_int_pt = 0.0;
2110 NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
2111
2112 double T_int_pt = 0.0;
2113 NumLib::shapeFunctionInterpolate(local_T, N, T_int_pt);
2114
2115 vars.concentration = C_int_pt;
2116 vars.liquid_phase_pressure = p_int_pt;
2117 vars.porosity = ip_data.porosity;
2118 vars.temperature = T_int_pt;
2119
2120 pos.setCoordinates(MathLib::Point3d(
2121 NumLib::interpolateCoordinates<ShapeFunction,
2122 ShapeMatricesType>(_element,
2123 N)));
2124
2125 // TODO (naumov) Temporary value not used by current material
2126 // models. Need extension of secondary variables interface.
2127 double const dt = std::numeric_limits<double>::quiet_NaN();
2128 auto const permeability_tensor =
2129 MaterialPropertyLib::formEigenTensor<GlobalDim>(
2130 medium[MaterialPropertyLib::PropertyType::permeability]
2131 .value(vars, pos, t, dt));
2132 permeability_avg += permeability_tensor.trace() / GlobalDim;
2133 }
2134 (*_process_data.mesh_prop_permeability)[ele_id] =
2135 permeability_avg / n_integration_points;
2136
2137 if (!_process_data.chemically_induced_porosity_change)
2138 {
2139 updateAveragePorosity(ele_id);
2140 }
2141 }
ProcessLib::ComponentTransport::LocalAssemblerData::calculateIntPtDarcyVelocity
std::vector< double > const & calculateIntPtDarcyVelocity(const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &T_nodal_values, std::vector< double > &cache) const
Definition ComponentTransportFEM.h:1907
MathLib::toMatrix
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition EigenMapTools.h:65

References _element, _integration_method, _ip_data, _process_data, calculateIntPtDarcyVelocity(), MaterialPropertyLib::VariableArray::concentration, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getLocalTemperature(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::VariableArray::porosity, pressure_index, ParameterLib::SpatialPosition::setCoordinates(), ParameterLib::SpatialPosition::setElementID(), NumLib::detail::shapeFunctionInterpolate(), MaterialPropertyLib::VariableArray::temperature, MathLib::toMatrix(), and updateAveragePorosity().

◆ getFlux()

template<typename ShapeFunction, int GlobalDim>
Eigen::Vector3d ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getFlux ( MathLib::Point3d const & ,
double const ,
std::vector< double > const &  ) const
inlineoverridevirtual

Computes the flux in the point p_local_coords that is given in local coordinates using the values from local_x. Fits to monolithic scheme.

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 1998 of file ComponentTransportFEM.h.

2001 {
2002 auto const local_p = Eigen::Map<const NodalVectorType>(
2003 &local_x[pressure_index], pressure_size);
2004 auto const local_C = Eigen::Map<const NodalVectorType>(
2005 &local_x[first_concentration_index], concentration_size);
2006
2007 // Eval shape matrices at given point
2008 // Note: Axial symmetry is set to false here, because we only need dNdx
2009 // here, which is not affected by axial symmetry.
2010 auto const shape_matrices =
2011 NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
2012 GlobalDim>(
2013 _element, false /*is_axially_symmetric*/,
2014 std::array{pnt_local_coords})[0];
2015
2016 ParameterLib::SpatialPosition const pos(
2017 {}, _element.getID(),
2018 MathLib::Point3d(NumLib::interpolateCoordinates<ShapeFunction,
2019 ShapeMatricesType>(
2020 _element, shape_matrices.N)));
2021
2022 auto const& b =
2023 _process_data
2024 .projected_specific_body_force_vectors[_element.getID()];
2025
2026 MaterialPropertyLib::VariableArray vars;
2027
2028 auto const& medium =
2029 *_process_data.media_map.getMedium(_element.getID());
2030 auto const& phase =
2031 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
2032
2033 // local_x contains the local concentration and pressure values
2034 double c_int_pt;
2035 NumLib::shapeFunctionInterpolate(local_C, shape_matrices.N, c_int_pt);
2036 vars.concentration = c_int_pt;
2037
2038 double p_int_pt;
2039 NumLib::shapeFunctionInterpolate(local_p, shape_matrices.N, p_int_pt);
2040 vars.liquid_phase_pressure = p_int_pt;
2041
2042 // TODO (naumov) Temporary value not used by current material models.
2043 // Need extension of secondary variables interface.
2044 double const dt = std::numeric_limits<double>::quiet_NaN();
2045 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2046 medium[MaterialPropertyLib::PropertyType::permeability].value(
2047 vars, pos, t, dt));
2048
2049 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
2050 .template value<double>(vars, pos, t, dt);
2051 GlobalDimMatrixType const K_over_mu = K / mu;
2052
2053 GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
2054 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
2055 .template value<double>(vars, pos, t, dt);
2056 if (_process_data.has_gravity)
2057 {
2058 q += K_over_mu * rho_w * b;
2059 }
2060 Eigen::Vector3d flux(0.0, 0.0, 0.0);
2061 flux.head<GlobalDim>() = rho_w * q;
2062 return flux;
2063 }

References _element, _process_data, MaterialPropertyLib::AqueousLiquid, NumLib::computeShapeMatrices(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, pressure_index, pressure_size, NumLib::detail::shapeFunctionInterpolate(), and MaterialPropertyLib::viscosity.

◆ getHeatEnergyCoefficient()

template<typename ShapeFunction, int GlobalDim>
double ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getHeatEnergyCoefficient ( MaterialPropertyLib::VariableArray const & vars,
const double porosity,
const double fluid_density,
const double specific_heat_capacity_fluid,
ParameterLib::SpatialPosition const & pos,
double const t,
double const dt )
inlineprivate

Definition at line 2319 of file ComponentTransportFEM.h.

2324 {
2325 auto const& medium =
2326 *_process_data.media_map.getMedium(this->_element.getID());
2327 auto const& solid_phase =
2328 medium.phase(MaterialPropertyLib::PhaseName::Solid);
2329
2330 auto const specific_heat_capacity_solid =
2331 solid_phase
2332 .property(
2333 MaterialPropertyLib::PropertyType::specific_heat_capacity)
2334 .template value<double>(vars, pos, t, dt);
2335
2336 auto const solid_density =
2337 solid_phase.property(MaterialPropertyLib::PropertyType::density)
2338 .template value<double>(vars, pos, t, dt);
2339
2340 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2341 fluid_density * specific_heat_capacity_fluid * porosity;
2342 }

References _process_data, MaterialPropertyLib::density, MeshLib::Element::getID(), MaterialPropertyLib::Solid, and MaterialPropertyLib::specific_heat_capacity.

Referenced by assembleHeatTransportEquation().

◆ getIntPtDarcyVelocity()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getIntPtDarcyVelocity ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 1832 of file ComponentTransportFEM.h.

1837 {
1838 assert(x.size() == dof_table.size());
1839
1840 auto const n_processes = x.size();
1841 std::vector<std::vector<double>> local_x;
1842 local_x.reserve(n_processes);
1843
1844 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1845 {
1846 auto const indices =
1847 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1848 assert(!indices.empty());
1849 local_x.push_back(x[process_id]->get(indices));
1850 }
1851
1852 // only one process, must be monolithic.
1853 if (n_processes == 1)
1854 {
1855 auto const local_p = Eigen::Map<const NodalVectorType>(
1856 &local_x[0][pressure_index], pressure_size);
1857 auto const local_C = Eigen::Map<const NodalVectorType>(
1858 &local_x[0][first_concentration_index], concentration_size);
1859 NodalVectorType local_T;
1860 local_T.setConstant(ShapeFunction::NPOINTS,
1861 std::numeric_limits<double>::quiet_NaN());
1862 int const temperature_index =
1863 _process_data.isothermal ? -1 : ShapeFunction::NPOINTS;
1864 if (temperature_index != -1)
1865 {
1866 local_T = Eigen::Map<const NodalVectorType>(
1867 &local_x[0][temperature_index], temperature_size);
1868 }
1869 return calculateIntPtDarcyVelocity(t, local_p, local_C, local_T,
1870 cache);
1871 }
1872
1873 // multiple processes, must be staggered.
1874 {
1875 constexpr int pressure_process_id = 0;
1876 int concentration_process_id = 1;
1877 // Normally temperature process is not there,
1878 // hence set the default temperature index to -1
1879 int temperature_process_id = -1;
1880
1881 // check whether temperature process exists
1882 if (!_process_data.isothermal)
1883 {
1884 // if temperature process exists, its id is 1
1885 temperature_process_id = 1;
1886 // then the concentration index shifts to 2
1887 concentration_process_id = 2;
1888 }
1889
1890 auto const local_p = Eigen::Map<const NodalVectorType>(
1891 &local_x[pressure_process_id][0], pressure_size);
1892 auto const local_C = Eigen::Map<const NodalVectorType>(
1893 &local_x[concentration_process_id][0], concentration_size);
1894 NodalVectorType local_T;
1895 local_T.setConstant(ShapeFunction::NPOINTS,
1896 std::numeric_limits<double>::quiet_NaN());
1897 if (temperature_process_id != -1)
1898 {
1899 local_T = Eigen::Map<const NodalVectorType>(
1900 &local_x[temperature_process_id][0], temperature_size);
1901 }
1902 return calculateIntPtDarcyVelocity(t, local_p, local_C, local_T,
1903 cache);
1904 }
1905 }
NumLib::getIndices
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
Definition DOFTableUtil.cpp:105

References _element, _process_data, calculateIntPtDarcyVelocity(), concentration_size, first_concentration_index, NumLib::getIndices(), pressure_index, pressure_size, temperature_index, and temperature_size.

◆ getIntPtLiquidDensity()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getIntPtLiquidDensity ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_table,
std::vector< double > & cache ) const
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 1701 of file ComponentTransportFEM.h.

1706 {
1707 assert(x.size() == dof_table.size());
1708
1709 auto const n_processes = x.size();
1710 std::vector<std::vector<double>> local_x;
1711 local_x.reserve(n_processes);
1712
1713 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1714 {
1715 auto const indices =
1716 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1717 assert(!indices.empty());
1718 local_x.push_back(x[process_id]->get(indices));
1719 }
1720
1721 // only one process, must be monolithic.
1722 if (n_processes == 1)
1723 {
1724 auto const local_p = Eigen::Map<const NodalVectorType>(
1725 &local_x[0][pressure_index], pressure_size);
1726 auto const local_C = Eigen::Map<const NodalVectorType>(
1727 &local_x[0][first_concentration_index], concentration_size);
1728 NodalVectorType local_T;
1729 local_T.setConstant(ShapeFunction::NPOINTS,
1730 std::numeric_limits<double>::quiet_NaN());
1731 int const temperature_index =
1732 _process_data.isothermal ? -1 : ShapeFunction::NPOINTS;
1733 if (temperature_index != -1)
1734 {
1735 local_T = Eigen::Map<const NodalVectorType>(
1736 &local_x[0][temperature_index], temperature_size);
1737 }
1738 return calculateIntPtLiquidDensity(t, local_p, local_C, local_T,
1739 cache);
1740 }
1741
1742 // multiple processes, must be staggered.
1743 {
1744 constexpr int pressure_process_id = 0;
1745 int concentration_process_id = 1;
1746 // Normally temperature process is not there,
1747 // hence set the default temperature index to -1
1748 int temperature_process_id = -1;
1749
1750 // check whether temperature process exists
1751 if (!_process_data.isothermal)
1752 {
1753 // if temperature process exists, its id is 1
1754 temperature_process_id = 1;
1755 // then the concentration index shifts to 2
1756 concentration_process_id = 2;
1757 }
1758
1759 auto const local_p = Eigen::Map<const NodalVectorType>(
1760 &local_x[pressure_process_id][0], pressure_size);
1761 auto const local_C = Eigen::Map<const NodalVectorType>(
1762 &local_x[concentration_process_id][0], concentration_size);
1763 NodalVectorType local_T;
1764 local_T.setConstant(ShapeFunction::NPOINTS,
1765 std::numeric_limits<double>::quiet_NaN());
1766 if (temperature_process_id != -1)
1767 {
1768 local_T = Eigen::Map<const NodalVectorType>(
1769 &local_x[temperature_process_id][0], temperature_size);
1770 }
1771 return calculateIntPtLiquidDensity(t, local_p, local_C, local_T,
1772 cache);
1773 }
1774 }
ProcessLib::ComponentTransport::LocalAssemblerData::calculateIntPtLiquidDensity
std::vector< double > const & calculateIntPtLiquidDensity(const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &T_nodal_values, std::vector< double > &cache) const
Definition ComponentTransportFEM.h:1776

References _element, _process_data, calculateIntPtLiquidDensity(), concentration_size, first_concentration_index, NumLib::getIndices(), pressure_index, pressure_size, temperature_index, and temperature_size.

◆ getIntPtMolarFlux()

template<typename ShapeFunction, int GlobalDim>
std::vector< double > const & ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getIntPtMolarFlux ( const double t,
std::vector< GlobalVector * > const & x,
std::vector< NumLib::LocalToGlobalIndexMap const * > const & dof_tables,
std::vector< double > & cache,
int const component_id ) const
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 2176 of file ComponentTransportFEM.h.

2180 {
2181 std::vector<double> local_x_vec;
2182
2183 auto const n_processes = x.size();
2184 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
2185 {
2186 auto const indices =
2187 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
2188 assert(!indices.empty());
2189 auto const local_solution = x[process_id]->get(indices);
2190 local_x_vec.insert(std::end(local_x_vec),
2191 std::begin(local_solution),
2192 std::end(local_solution));
2193 }
2194 auto const local_x = MathLib::toVector(local_x_vec);
2195
2196 auto const p = local_x.template segment<pressure_size>(pressure_index);
2197 auto const c = local_x.template segment<concentration_size>(
2198 first_concentration_index + component_id * concentration_size);
2199
2200 auto const n_integration_points =
2201 _integration_method.getNumberOfPoints();
2202
2203 cache.clear();
2204 auto cache_mat = MathLib::createZeroedMatrix<
2205 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
2206 cache, GlobalDim, n_integration_points);
2207
2208 auto const& b =
2209 _process_data
2210 .projected_specific_body_force_vectors[_element.getID()];
2211
2212 MaterialPropertyLib::VariableArray vars;
2213
2214 auto const& medium =
2215 *_process_data.media_map.getMedium(_element.getID());
2216 auto const& phase =
2217 medium.phase(MaterialPropertyLib::PhaseName::AqueousLiquid);
2218
2219 auto const& component = phase.component(
2220 _transport_process_variables[component_id].get().getName());
2221
2222 auto const& Ns =
2223 _process_data.shape_matrix_cache
2224 .NsHigherOrder<typename ShapeFunction::MeshElement>();
2225
2226 for (unsigned ip = 0; ip < n_integration_points; ++ip)
2227 {
2228 auto const& ip_data = _ip_data[ip];
2229 auto const& dNdx = ip_data.dNdx;
2230 auto const& N = Ns[ip];
2231 auto const& phi = ip_data.porosity;
2232
2233 ParameterLib::SpatialPosition const pos(
2234 {}, _element.getID(),
2235 MathLib::Point3d(
2236 NumLib::interpolateCoordinates<ShapeFunction,
2237 ShapeMatricesType>(_element,
2238 N)));
2239
2240 double const p_ip = N.dot(p);
2241 double const c_ip = N.dot(c);
2242
2243 vars.concentration = c_ip;
2244 vars.liquid_phase_pressure = p_ip;
2245 vars.porosity = phi;
2246
2247 double const dt = std::numeric_limits<double>::quiet_NaN();
2248
2249 auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2250 medium[MaterialPropertyLib::PropertyType::permeability].value(
2251 vars, pos, t, dt));
2252 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
2253 .template value<double>(vars, pos, t, dt);
2254 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
2255 .template value<double>(vars, pos, t, dt);
2256
2257 // Darcy flux
2258 GlobalDimVectorType const q =
2259 _process_data.has_gravity
2260 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
2261 : GlobalDimVectorType(-k / mu * dNdx * p);
2262
2263 auto const alpha_T = medium.template value<double>(
2264 MaterialPropertyLib::PropertyType::transversal_dispersivity);
2265 auto const alpha_L = medium.template value<double>(
2266 MaterialPropertyLib::PropertyType::longitudinal_dispersivity);
2267 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2268 component[MaterialPropertyLib::PropertyType::pore_diffusion]
2269 .value(vars, pos, t, dt));
2270
2271 // Hydrodynamic dispersion
2272 GlobalDimMatrixType const D = NumLib::computeHydrodynamicDispersion(
2273 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
2274 alpha_L);
2275
2276 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
2277 }
2278
2279 return cache;
2280 }
MathLib::toVector
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Definition EigenMapTools.h:160

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), NumLib::getIndices(), getName(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::longitudinal_dispersivity, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::VariableArray::porosity, pressure_index, MathLib::toVector(), MaterialPropertyLib::transversal_dispersivity, and MaterialPropertyLib::viscosity.

◆ getLocalTemperature()

template<typename ShapeFunction, int GlobalDim>
NodalVectorType ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getLocalTemperature ( double const t,
Eigen::VectorXd const & local_x ) const
inlineprivate

Definition at line 2385 of file ComponentTransportFEM.h.

2387 {
2388 NodalVectorType local_T;
2389 if (_process_data.isothermal)
2390 {
2391 if (_process_data.temperature)
2392 {
2393 local_T = _process_data.temperature->getNodalValuesOnElement(
2394 _element, t);
2395 }
2396 else
2397 {
2398 local_T = NodalVectorType::Zero(temperature_size);
2399 }
2400 }
2401 else
2402 {
2403 local_T =
2404 local_x.template segment<temperature_size>(temperature_index);
2405 }
2406 return local_T;
2407 }

References _element, _process_data, temperature_index, and temperature_size.

Referenced by assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleWithJacobianComponentTransportEquation(), and computeSecondaryVariableConcrete().

◆ getShapeMatrix()

template<typename ShapeFunction, int GlobalDim>
Eigen::Map< const Eigen::RowVectorXd > ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getShapeMatrix ( const unsigned integration_point) const
inlineoverridevirtual

Provides the shape matrix at the given integration point.

Implements NumLib::ExtrapolatableElement.

Definition at line 1988 of file ComponentTransportFEM.h.

1990 {
1991 auto const& N = _process_data.shape_matrix_cache.NsHigherOrder<
1992 typename ShapeFunction::MeshElement>()[integration_point];
1993
1994 // assumes N is stored contiguously in memory
1995 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1996 }

References _process_data.

◆ getThermalConductivityDispersivity()

template<typename ShapeFunction, int GlobalDim>
GlobalDimMatrixType ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::getThermalConductivityDispersivity ( MaterialPropertyLib::VariableArray const & vars,
const double fluid_density,
const double specific_heat_capacity_fluid,
const GlobalDimVectorType & velocity,
ParameterLib::SpatialPosition const & pos,
double const t,
double const dt )
inlineprivate

Definition at line 2344 of file ComponentTransportFEM.h.

2350 {
2351 auto const& medium =
2352 *_process_data.media_map.getMedium(_element.getID());
2353
2354 auto thermal_conductivity =
2355 MaterialPropertyLib::formEigenTensor<GlobalDim>(
2356 medium
2357 .property(
2358 MaterialPropertyLib::PropertyType::thermal_conductivity)
2359 .value(vars, pos, t, dt));
2360
2361 auto const thermal_dispersivity_transversal =
2362 medium
2363 .property(MaterialPropertyLib::PropertyType::
2364 thermal_transversal_dispersivity)
2365 .template value<double>();
2366
2367 auto const thermal_dispersivity_longitudinal =
2368 medium
2369 .property(MaterialPropertyLib::PropertyType::
2370 thermal_longitudinal_dispersivity)
2371 .template value<double>();
2372
2373 // Thermal conductivity is moved outside and zero matrix is passed
2374 // instead due to multiplication with fluid's density times specific
2375 // heat capacity.
2376 return thermal_conductivity +
2377 fluid_density * specific_heat_capacity_fluid *
2378 NumLib::computeHydrodynamicDispersion(
2379 _process_data.stabilizer, _element.getID(),
2380 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2381 velocity, 0 /* phi */, thermal_dispersivity_transversal,
2382 thermal_dispersivity_longitudinal);
2383 }

References _element, _process_data, NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::formEigenTensor(), and MaterialPropertyLib::thermal_conductivity.

Referenced by assembleHeatTransportEquation().

◆ initializeChemicalSystemConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::initializeChemicalSystemConcrete ( Eigen::VectorXd const & local_x,
double const t )
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 323 of file ComponentTransportFEM.h.

325 {
326 assert(_process_data.chemical_solver_interface);
327
328 auto const& medium =
329 *_process_data.media_map.getMedium(_element.getID());
330
331 auto const& Ns =
332 _process_data.shape_matrix_cache
333 .NsHigherOrder<typename ShapeFunction::MeshElement>();
334
335 unsigned const n_integration_points =
336 _integration_method.getNumberOfPoints();
337
338 for (unsigned ip = 0; ip < n_integration_points; ip++)
339 {
340 auto& ip_data = _ip_data[ip];
341 auto const& N = Ns[ip];
342 auto const& chemical_system_id = ip_data.chemical_system_id;
343
344 ParameterLib::SpatialPosition const pos(
345 {}, _element.getID(),
346 MathLib::Point3d(
347 NumLib::interpolateCoordinates<ShapeFunction,
348 ShapeMatricesType>(_element,
349 N)));
350
351 auto const n_component = _transport_process_variables.size();
352 std::vector<double> C_int_pt(n_component);
353 for (unsigned component_id = 0; component_id < n_component;
354 ++component_id)
355 {
356 auto const concentration_index =
357 first_concentration_index +
358 component_id * concentration_size;
359 auto const local_C =
360 local_x.template segment<concentration_size>(
361 concentration_index);
362
363 NumLib::shapeFunctionInterpolate(local_C, N,
364 C_int_pt[component_id]);
365 }
366
367 _process_data.chemical_solver_interface
368 ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id,
369 medium, pos, t);
370 }
371 }

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, concentration_size, first_concentration_index, NumLib::interpolateCoordinates(), and NumLib::detail::shapeFunctionInterpolate().

◆ postSpeciationCalculation()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::postSpeciationCalculation ( std::size_t const ele_id,
double const t,
double const dt )
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 443 of file ComponentTransportFEM.h.

445 {
446 if (!_process_data.chemically_induced_porosity_change)
447 {
448 return;
449 }
450
451 auto const& medium = *_process_data.media_map.getMedium(ele_id);
452
453 ParameterLib::SpatialPosition pos;
454 pos.setElementID(ele_id);
455
456 for (auto& ip_data : _ip_data)
457 {
458 ip_data.porosity = ip_data.porosity_prev;
459
460 _process_data.chemical_solver_interface
461 ->updateVolumeFractionPostReaction(ip_data.chemical_system_id,
462 medium, pos,
463 ip_data.porosity, t, dt);
464
465 _process_data.chemical_solver_interface->updatePorosityPostReaction(
466 ip_data.chemical_system_id, medium, ip_data.porosity);
467 }
468 }

References _ip_data, _process_data, and ParameterLib::SpatialPosition::setElementID().

◆ postTimestepConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::postTimestepConcrete ( Eigen::VectorXd const & ,
Eigen::VectorXd const & ,
double const ,
double const ,
int const  )
inlineoverridevirtual

Reimplemented from ProcessLib::LocalAssemblerInterface.

Definition at line 2282 of file ComponentTransportFEM.h.

2286 {
2287 unsigned const n_integration_points =
2288 _integration_method.getNumberOfPoints();
2289
2290 for (unsigned ip = 0; ip < n_integration_points; ip++)
2291 {
2292 _ip_data[ip].pushBackState();
2293 }
2294 }

References _integration_method, and _ip_data.

◆ setChemicalSystemConcrete()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::setChemicalSystemConcrete ( Eigen::VectorXd const & local_x,
double const t,
double dt )
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 373 of file ComponentTransportFEM.h.

375 {
376 assert(_process_data.chemical_solver_interface);
377
378 auto const& medium =
379 _process_data.media_map.getMedium(_element.getID());
380
381 MaterialPropertyLib::VariableArray vars;
382 MaterialPropertyLib::VariableArray vars_prev;
383
384 auto const& Ns =
385 _process_data.shape_matrix_cache
386 .NsHigherOrder<typename ShapeFunction::MeshElement>();
387
388 unsigned const n_integration_points =
389 _integration_method.getNumberOfPoints();
390
391 for (unsigned ip = 0; ip < n_integration_points; ip++)
392 {
393 auto& ip_data = _ip_data[ip];
394 auto const& N = Ns[ip];
395 auto& porosity = ip_data.porosity;
396 auto const& porosity_prev = ip_data.porosity_prev;
397 auto const& chemical_system_id = ip_data.chemical_system_id;
398
399 ParameterLib::SpatialPosition const pos(
400 {}, _element.getID(),
401 MathLib::Point3d(
402 NumLib::interpolateCoordinates<ShapeFunction,
403 ShapeMatricesType>(_element,
404 N)));
405
406 auto const n_component = _transport_process_variables.size();
407 std::vector<double> C_int_pt(n_component);
408
409 for (unsigned component_id = 0; component_id < n_component;
410 ++component_id)
411 {
412 auto const concentration_index =
413 first_concentration_index +
414 component_id * concentration_size;
415 auto const local_C =
416 local_x.template segment<concentration_size>(
417 concentration_index);
418
419 NumLib::shapeFunctionInterpolate(local_C, N,
420 C_int_pt[component_id]);
421 }
422
423 {
424 vars_prev.porosity = porosity_prev;
425
426 porosity =
427 _process_data.chemically_induced_porosity_change
428 ? porosity_prev
429 : medium
430 ->property(
431 MaterialPropertyLib::PropertyType::porosity)
432 .template value<double>(vars, vars_prev, pos, t,
433 dt);
434
435 vars.porosity = porosity;
436 }
437
438 _process_data.chemical_solver_interface->setChemicalSystemConcrete(
439 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
440 }
441 }

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, concentration_size, first_concentration_index, NumLib::interpolateCoordinates(), MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, and NumLib::detail::shapeFunctionInterpolate().

◆ setChemicalSystemID()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::setChemicalSystemID ( std::size_t const )
inlineoverridevirtual

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 303 of file ComponentTransportFEM.h.

304 {
305 assert(_process_data.chemical_solver_interface);
306 // chemical system index map
307 auto& chemical_system_index_map =
308 _process_data.chemical_solver_interface->chemical_system_index_map;
309
310 unsigned const n_integration_points =
311 _integration_method.getNumberOfPoints();
312 for (unsigned ip = 0; ip < n_integration_points; ip++)
313 {
314 _ip_data[ip].chemical_system_id =
315 chemical_system_index_map.empty()
316 ? 0
317 : chemical_system_index_map.back() + 1;
318 chemical_system_index_map.push_back(
319 _ip_data[ip].chemical_system_id);
320 }
321 }

References _integration_method, _ip_data, and _process_data.

◆ updateAveragePorosity()

template<typename ShapeFunction, int GlobalDim>
void ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::updateAveragePorosity ( std::size_t const ele_id)
inlineprivate

Averages ip_data.porosity over the element's integration points and writes it into the porosity_avg mesh output for ele_id.

Definition at line 2299 of file ComponentTransportFEM.h.

2300 {
2301 auto const n_integration_points =
2302 _integration_method.getNumberOfPoints();
2303 (*_process_data.mesh_prop_porosity)[ele_id] =
2304 std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
2305 [](double const s, auto const& ip)
2306 { return s + ip.porosity; }) /
2307 n_integration_points;
2308 }

References _integration_method, _ip_data, and _process_data.

Referenced by computeReactionRelatedSecondaryVariable(), and computeSecondaryVariableConcrete().

Member Data Documentation

◆ _element

template<typename ShapeFunction, int GlobalDim>
MeshLib::Element const& ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::_element
private

Definition at line 2310 of file ComponentTransportFEM.h.

Referenced by LocalAssemblerData(), assemble(), assembleBlockMatrices(), assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleKCmCn(), assembleReactionEquationConcrete(), assembleWithJacobianComponentTransportEquation(), assembleWithJacobianHydraulicEquation(), calculateIntPtDarcyVelocity(), calculateIntPtLiquidDensity(), computeSecondaryVariableConcrete(), getFlux(), getIntPtDarcyVelocity(), getIntPtLiquidDensity(), getIntPtMolarFlux(), getLocalTemperature(), getThermalConductivityDispersivity(), initializeChemicalSystemConcrete(), and setChemicalSystemConcrete().

◆ _integration_method

template<typename ShapeFunction, int GlobalDim>
NumLib::GenericIntegrationMethod const& ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::_integration_method
private

Definition at line 2313 of file ComponentTransportFEM.h.

Referenced by LocalAssemblerData(), assembleBlockMatrices(), assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleKCmCn(), assembleReactionEquationConcrete(), assembleWithJacobianComponentTransportEquation(), assembleWithJacobianHydraulicEquation(), calculateIntPtDarcyVelocity(), calculateIntPtLiquidDensity(), computeReactionRelatedSecondaryVariable(), computeSecondaryVariableConcrete(), getIntPtMolarFlux(), initializeChemicalSystemConcrete(), postTimestepConcrete(), setChemicalSystemConcrete(), setChemicalSystemID(), and updateAveragePorosity().

◆ _ip_data

template<typename ShapeFunction, int GlobalDim>
std::vector<IntegrationPointData<GlobalDimNodalMatrixType> > ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::_ip_data
private

Definition at line 2317 of file ComponentTransportFEM.h.

Referenced by LocalAssemblerData(), assembleBlockMatrices(), assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleKCmCn(), assembleReactionEquationConcrete(), assembleWithJacobianComponentTransportEquation(), assembleWithJacobianHydraulicEquation(), calculateIntPtDarcyVelocity(), computeReactionRelatedSecondaryVariable(), computeSecondaryVariableConcrete(), getIntPtMolarFlux(), initializeChemicalSystemConcrete(), postSpeciationCalculation(), postTimestepConcrete(), setChemicalSystemConcrete(), setChemicalSystemID(), and updateAveragePorosity().

◆ _process_data

template<typename ShapeFunction, int GlobalDim>
ComponentTransportProcessData const& ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::_process_data
private

Definition at line 2311 of file ComponentTransportFEM.h.

Referenced by LocalAssemblerData(), assemble(), assembleBlockMatrices(), assembleComponentTransportEquation(), assembleForStaggeredScheme(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleKCmCn(), assembleReactionEquationConcrete(), assembleWithJacobianComponentTransportEquation(), assembleWithJacobianForStaggeredScheme(), assembleWithJacobianHydraulicEquation(), calculateIntPtDarcyVelocity(), calculateIntPtLiquidDensity(), computeReactionRelatedSecondaryVariable(), computeSecondaryVariableConcrete(), getFlux(), getHeatEnergyCoefficient(), getIntPtDarcyVelocity(), getIntPtLiquidDensity(), getIntPtMolarFlux(), getLocalTemperature(), getShapeMatrix(), getThermalConductivityDispersivity(), initializeChemicalSystemConcrete(), postSpeciationCalculation(), setChemicalSystemConcrete(), setChemicalSystemID(), and updateAveragePorosity().

◆ _transport_process_variables

template<typename ShapeFunction, int GlobalDim>
std::vector<std::reference_wrapper<ProcessVariable> > const ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::_transport_process_variables
private

Definition at line 2315 of file ComponentTransportFEM.h.

Referenced by LocalAssemblerData(), assemble(), assembleBlockMatrices(), assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleKCmCn(), assembleWithJacobianComponentTransportEquation(), getIntPtMolarFlux(), initializeChemicalSystemConcrete(), and setChemicalSystemConcrete().

◆ concentration_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::concentration_size
staticprivate
Initial value:
=
ShapeFunction::NPOINTS

Definition at line 215 of file ComponentTransportFEM.h.

Referenced by assemble(), assembleBlockMatrices(), assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleReactionEquationConcrete(), assembleWithJacobianComponentTransportEquation(), assembleWithJacobianHydraulicEquation(), getFlux(), getIntPtDarcyVelocity(), getIntPtLiquidDensity(), getIntPtMolarFlux(), initializeChemicalSystemConcrete(), and setChemicalSystemConcrete().

◆ first_concentration_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::first_concentration_index = -1
private

Definition at line 211 of file ComponentTransportFEM.h.

Referenced by LocalAssemblerData(), assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleReactionEquationConcrete(), assembleWithJacobianComponentTransportEquation(), assembleWithJacobianHydraulicEquation(), computeSecondaryVariableConcrete(), getFlux(), getIntPtDarcyVelocity(), getIntPtLiquidDensity(), getIntPtMolarFlux(), initializeChemicalSystemConcrete(), and setChemicalSystemConcrete().

◆ pressure_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::pressure_index = 0
staticprivate

Definition at line 209 of file ComponentTransportFEM.h.

Referenced by assemble(), assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleWithJacobianComponentTransportEquation(), assembleWithJacobianHydraulicEquation(), computeSecondaryVariableConcrete(), getFlux(), getIntPtDarcyVelocity(), getIntPtLiquidDensity(), and getIntPtMolarFlux().

◆ pressure_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::pressure_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 213 of file ComponentTransportFEM.h.

Referenced by assemble(), assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), assembleWithJacobianHydraulicEquation(), getFlux(), getIntPtDarcyVelocity(), and getIntPtLiquidDensity().

◆ temperature_index

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::temperature_index = -1
private

Definition at line 210 of file ComponentTransportFEM.h.

Referenced by LocalAssemblerData(), getIntPtDarcyVelocity(), getIntPtLiquidDensity(), and getLocalTemperature().

◆ temperature_size

template<typename ShapeFunction, int GlobalDim>
const int ProcessLib::ComponentTransport::LocalAssemblerData< ShapeFunction, GlobalDim >::temperature_size = ShapeFunction::NPOINTS
staticprivate

Definition at line 214 of file ComponentTransportFEM.h.

Referenced by assembleComponentTransportEquation(), assembleHeatTransportEquation(), getIntPtDarcyVelocity(), getIntPtLiquidDensity(), and getLocalTemperature().

The documentation for this class was generated from the following file: