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

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 pos;
265 pos.setElementID(_element.getID());
266
267 double const aperture_size = _process_data.aperture_size(0.0, pos)[0];
268
269 auto const shape_matrices =
270 NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
271 GlobalDim>(element, is_axially_symmetric,
272 _integration_method);
273 auto const& medium =
274 *_process_data.media_map.getMedium(_element.getID());
275 for (unsigned ip = 0; ip < n_integration_points; ip++)
276 {
277 _ip_data.emplace_back(
278 shape_matrices[ip].dNdx,
279 _integration_method.getWeightedPoint(ip).getWeight() *
280 shape_matrices[ip].integralMeasure *
281 shape_matrices[ip].detJ * aperture_size);
282
283 _ip_data[ip].porosity =
284 medium[MaterialPropertyLib::PropertyType::porosity]
285 .template initialValue<double>(
286 pos, std::numeric_limits<double>::quiet_NaN() /*t*/);
287
288 _ip_data[ip].pushBackState();
289 }
290 }
ProcessLib::ComponentTransport::LocalAssemblerData::_integration_method
NumLib::GenericIntegrationMethod const & _integration_method
Definition ComponentTransportFEM.h:2214
ProcessLib::ComponentTransport::LocalAssemblerData::_ip_data
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data
Definition ComponentTransportFEM.h:2218
ProcessLib::ComponentTransport::LocalAssemblerData::_process_data
ComponentTransportProcessData const & _process_data
Definition ComponentTransportFEM.h:2212
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:2216

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, first_concentration_index, NumLib::initShapeMatrices(), MaterialPropertyLib::porosity, ParameterLib::SpatialPosition::setElementID(), 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 457 of file ComponentTransportFEM.h.

463 {
464 auto const local_matrix_size = local_x.size();
465 // Nodal DOFs include pressure
466 int const num_nodal_dof = 1 + _transport_process_variables.size();
467 // This assertion is valid only if all nodal d.o.f. use the same shape
468 // matrices.
469 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
470
471 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
472 local_M_data, local_matrix_size, local_matrix_size);
473 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
474 local_K_data, local_matrix_size, local_matrix_size);
475 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
476 local_b_data, local_matrix_size);
477
478 // Get block matrices
479 auto Kpp = local_K.template block<pressure_size, pressure_size>(
480 pressure_index, pressure_index);
481 auto Mpp = local_M.template block<pressure_size, pressure_size>(
482 pressure_index, pressure_index);
483 auto Bp = local_b.template segment<pressure_size>(pressure_index);
484
485 auto local_p = Eigen::Map<const NodalVectorType>(
486 &local_x[pressure_index], pressure_size);
487
488 auto const& b =
489 _process_data
490 .projected_specific_body_force_vectors[_element.getID()];
491
492 auto const number_of_components = num_nodal_dof - 1;
493 for (int component_id = 0; component_id < number_of_components;
494 ++component_id)
495 {
496 /* Partitioned assembler matrix
497 * | pp | pc1 | pc2 | pc3 |
498 * |-----|-----|-----|-----|
499 * | c1p | c1c1| 0 | 0 |
500 * |-----|-----|-----|-----|
501 * | c2p | 0 | c2c2| 0 |
502 * |-----|-----|-----|-----|
503 * | c3p | 0 | 0 | c3c3|
504 */
505 auto concentration_index =
506 pressure_size + component_id * concentration_size;
507
508 auto KCC =
509 local_K.template block<concentration_size, concentration_size>(
510 concentration_index, concentration_index);
511 auto MCC =
512 local_M.template block<concentration_size, concentration_size>(
513 concentration_index, concentration_index);
514 auto MCp =
515 local_M.template block<concentration_size, pressure_size>(
516 concentration_index, pressure_index);
517 auto MpC =
518 local_M.template block<pressure_size, concentration_size>(
519 pressure_index, concentration_index);
520
521 auto local_C = Eigen::Map<const NodalVectorType>(
522 &local_x[concentration_index], concentration_size);
523
524 assembleBlockMatrices(b, component_id, t, dt, local_C, local_p, KCC,
525 MCC, MCp, MpC, Kpp, Mpp, Bp);
526
527 if (_process_data.chemical_solver_interface)
528 {
529 auto const stoichiometric_matrix =
530 _process_data.chemical_solver_interface
531 ->getStoichiometricMatrix();
532
533 assert(stoichiometric_matrix);
534
535 for (Eigen::SparseMatrix<double>::InnerIterator it(
536 *stoichiometric_matrix, component_id);
537 it;
538 ++it)
539 {
540 auto const stoichiometric_coefficient = it.value();
541 auto const coupled_component_id = it.row();
542 auto const kinetic_prefactor =
543 _process_data.chemical_solver_interface
544 ->getKineticPrefactor(coupled_component_id);
545
546 auto const concentration_index =
547 pressure_size + component_id * concentration_size;
548 auto const coupled_concentration_index =
549 pressure_size +
550 coupled_component_id * concentration_size;
551 auto KCmCn = local_K.template block<concentration_size,
552 concentration_size>(
553 concentration_index, coupled_concentration_index);
554
555 // account for the coupling between components
556 assembleKCmCn(component_id, t, dt, KCmCn,
557 stoichiometric_coefficient,
558 kinetic_prefactor);
559 }
560 }
561 }
562 }
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:564
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:765
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 564 of file ComponentTransportFEM.h.

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

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, 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(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, MaterialPropertyLib::retardation_factor, ParameterLib::SpatialPosition::setElementID(), 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 819 of file ComponentTransportFEM.h.

826 {
827 if (process_id == _process_data.hydraulic_process_id)
828 {
829 assembleHydraulicEquation(t, dt, local_x, local_x_prev,
830 local_M_data, local_K_data, local_b_data);
831 }
832 else if (process_id == _process_data.thermal_process_id)
833 {
834 assembleHeatTransportEquation(t, dt, local_x, local_x_prev,
835 local_M_data, local_K_data,
836 local_b_data);
837 }
838 else
839 {
840 // Go for assembling in an order of transport process id.
841 assembleComponentTransportEquation(t, dt, local_x, local_x_prev,
842 local_M_data, local_K_data,
843 local_b_data, process_id);
844 }
845 }
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:982
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:847
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:1118

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 982 of file ComponentTransportFEM.h.

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

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

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getLocalTemperature(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, ParameterLib::SpatialPosition::setElementID(), 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 765 of file ComponentTransportFEM.h.

769 {
770 unsigned const n_integration_points =
771 _integration_method.getNumberOfPoints();
772
773 ParameterLib::SpatialPosition pos;
774 pos.setElementID(_element.getID());
775
776 MaterialPropertyLib::VariableArray vars;
777
778 auto const& medium =
779 *_process_data.media_map.getMedium(_element.getID());
780 auto const& phase = medium.phase("AqueousLiquid");
781 auto const& component = phase.component(
782 _transport_process_variables[component_id].get().getName());
783
784 auto const& Ns =
785 _process_data.shape_matrix_cache
786 .NsHigherOrder<typename ShapeFunction::MeshElement>();
787
788 for (unsigned ip(0); ip < n_integration_points; ++ip)
789 {
790 auto& ip_data = _ip_data[ip];
791 auto const& w = ip_data.integration_weight;
792 auto const& N = Ns[ip];
793 auto& porosity = ip_data.porosity;
794
795 // set position with N as the shape matrix at the current
796 // integration point
797 pos.setCoordinates(MathLib::Point3d(
798 NumLib::interpolateCoordinates<ShapeFunction,
799 ShapeMatricesType>(_element,
800 N)));
801
802 auto const retardation_factor =
803 component[MaterialPropertyLib::PropertyType::retardation_factor]
804 .template value<double>(vars, pos, t, dt);
805
806 porosity = medium[MaterialPropertyLib::PropertyType::porosity]
807 .template value<double>(vars, pos, t, dt);
808
809 auto const density =
810 phase[MaterialPropertyLib::PropertyType::density]
811 .template value<double>(vars, pos, t, dt);
812
813 KCmCn.noalias() -= N.transpose() * N *
814 (stoichiometric_coefficient * kinetic_prefactor *
815 retardation_factor * porosity * density * w);
816 }
817 }

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::density, getName(), NumLib::interpolateCoordinates(), MaterialPropertyLib::porosity, MaterialPropertyLib::retardation_factor, ParameterLib::SpatialPosition::setCoordinates(), and ParameterLib::SpatialPosition::setElementID().

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 1598 of file ComponentTransportFEM.h.

1603 {
1604 auto const local_C = local_x.template segment<concentration_size>(
1605 first_concentration_index +
1606 (transport_process_id - 1) * concentration_size);
1607
1608 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1609 local_M_data, concentration_size, concentration_size);
1610 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1611 local_K_data, concentration_size, concentration_size);
1612 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1613 local_b_data, concentration_size);
1614
1615 unsigned const n_integration_points =
1616 _integration_method.getNumberOfPoints();
1617
1618 ParameterLib::SpatialPosition pos;
1619 pos.setElementID(_element.getID());
1620
1621 MaterialPropertyLib::VariableArray vars;
1622 MaterialPropertyLib::VariableArray vars_prev;
1623
1624 auto const& medium =
1625 *_process_data.media_map.getMedium(_element.getID());
1626 auto const component_id = transport_process_id - 1;
1627
1628 auto const& Ns =
1629 _process_data.shape_matrix_cache
1630 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1631
1632 for (unsigned ip(0); ip < n_integration_points; ++ip)
1633 {
1634 auto& ip_data = _ip_data[ip];
1635 auto const w = ip_data.integration_weight;
1636 auto const& N = Ns[ip];
1637 auto& porosity = ip_data.porosity;
1638 auto const& porosity_prev = ip_data.porosity_prev;
1639 auto const chemical_system_id = ip_data.chemical_system_id;
1640
1641 double C_int_pt = 0.0;
1642 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
1643
1644 vars.concentration = C_int_pt;
1645
1646 auto const porosity_dot = (porosity - porosity_prev) / dt;
1647
1648 // porosity
1649 {
1650 vars_prev.porosity = porosity_prev;
1651
1652 porosity =
1653 _process_data.chemically_induced_porosity_change
1654 ? porosity_prev
1655 : medium[MaterialPropertyLib::PropertyType::porosity]
1656 .template value<double>(vars, vars_prev, pos, t,
1657 dt);
1658 }
1659
1660 local_M.noalias() += w * N.transpose() * porosity * N;
1661
1662 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1663
1664 if (chemical_system_id == -1)
1665 {
1666 continue;
1667 }
1668
1669 auto const C_post_int_pt =
1670 _process_data.chemical_solver_interface->getConcentration(
1671 component_id, chemical_system_id);
1672
1673 local_b.noalias() += N.transpose() * ((C_post_int_pt - C_int_pt) /
1674 dt * porosity * w);
1675 }
1676 }

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), first_concentration_index, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), 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 1450 of file ComponentTransportFEM.h.

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, NumLib::detail::assembleAdvectionMatrix(), NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getName(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::longitudinal_dispersivity, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, MaterialPropertyLib::retardation_factor, ParameterLib::SpatialPosition::setElementID(), 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 1322 of file ComponentTransportFEM.h.

1327 {
1328 if (process_id == _process_data.hydraulic_process_id)
1329 {
1330 assembleWithJacobianHydraulicEquation(t, dt, local_x, local_x_prev,
1331 local_b_data, local_Jac_data);
1332 }
1333 else
1334 {
1335 int const component_id = process_id - 1;
1336 assembleWithJacobianComponentTransportEquation(
1337 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1338 component_id);
1339 }
1340 }
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:1342
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:1450

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 1342 of file ComponentTransportFEM.h.

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

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, ParameterLib::SpatialPosition::setElementID(), 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 1883 of file ComponentTransportFEM.h.

1889 {
1890 auto const n_integration_points =
1891 _integration_method.getNumberOfPoints();
1892
1893 cache.clear();
1894 auto cache_mat = MathLib::createZeroedMatrix<
1895 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1896 cache, GlobalDim, n_integration_points);
1897
1898 ParameterLib::SpatialPosition pos;
1899 pos.setElementID(_element.getID());
1900
1901 auto const& b =
1902 _process_data
1903 .projected_specific_body_force_vectors[_element.getID()];
1904
1905 MaterialPropertyLib::VariableArray vars;
1906
1907 auto const& medium =
1908 *_process_data.media_map.getMedium(_element.getID());
1909 auto const& phase = medium.phase("AqueousLiquid");
1910
1911 auto const& Ns =
1912 _process_data.shape_matrix_cache
1913 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1914
1915 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1916 {
1917 auto const& ip_data = _ip_data[ip];
1918 auto const& dNdx = ip_data.dNdx;
1919 auto const& N = Ns[ip];
1920 auto const& porosity = ip_data.porosity;
1921
1922 double C_int_pt = 0.0;
1923 double p_int_pt = 0.0;
1924 double T_int_pt = 0.0;
1925
1926 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1927 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1928 NumLib::shapeFunctionInterpolate(T_nodal_values, N, T_int_pt);
1929
1930 vars.concentration = C_int_pt;
1931 vars.liquid_phase_pressure = p_int_pt;
1932 vars.porosity = porosity;
1933 vars.temperature = T_int_pt;
1934
1935 // TODO (naumov) Temporary value not used by current material
1936 // models. Need extension of secondary variables interface.
1937 double const dt = std::numeric_limits<double>::quiet_NaN();
1938 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1939 medium[MaterialPropertyLib::PropertyType::permeability].value(
1940 vars, pos, t, dt));
1941 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1942 .template value<double>(vars, pos, t, dt);
1943 GlobalDimMatrixType const K_over_mu = K / mu;
1944
1945 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1946 if (_process_data.has_gravity)
1947 {
1948 auto const rho_w =
1949 phase[MaterialPropertyLib::PropertyType::density]
1950 .template value<double>(vars, pos, t, dt);
1951 // here it is assumed that the vector b is directed 'downwards'
1952 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1953 }
1954 }
1955
1956 return cache;
1957 }

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::VariableArray::concentration, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), 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 1753 of file ComponentTransportFEM.h.

1759 {
1760 auto const n_integration_points =
1761 _integration_method.getNumberOfPoints();
1762
1763 cache.clear();
1764 auto cache_vec = MathLib::createZeroedVector<Eigen::VectorXd>(
1765 cache, n_integration_points);
1766
1767 ParameterLib::SpatialPosition pos;
1768 pos.setElementID(_element.getID());
1769
1770 MaterialPropertyLib::VariableArray vars;
1771
1772 auto const& medium =
1773 *_process_data.media_map.getMedium(_element.getID());
1774 auto const& phase = medium.phase("AqueousLiquid");
1775
1776 auto const& Ns =
1777 _process_data.shape_matrix_cache
1778 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1779
1780 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1781 {
1782 auto const& N = Ns[ip];
1783
1784 double C_int_pt = 0.0;
1785 double p_int_pt = 0.0;
1786 double T_int_pt = 0.0;
1787
1788 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1789 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1790 NumLib::shapeFunctionInterpolate(T_nodal_values, N, T_int_pt);
1791
1792 vars.concentration = C_int_pt;
1793 vars.liquid_phase_pressure = p_int_pt;
1794 vars.temperature = T_int_pt;
1795
1796 // TODO (naumov) Temporary value not used by current material
1797 // models. Need extension of secondary variables interface.
1798 double const dt = std::numeric_limits<double>::quiet_NaN();
1799
1800 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
1801 .template value<double>(vars, pos, t, dt);
1802 cache_vec[ip] = rho_w;
1803 }
1804
1805 return cache;
1806 }

References _element, _integration_method, _process_data, 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 2058 of file ComponentTransportFEM.h.

2060 {
2061 auto const n_integration_points =
2062 _integration_method.getNumberOfPoints();
2063
2064 if (_process_data.chemically_induced_porosity_change)
2065 {
2066 auto const& medium = *_process_data.media_map.getMedium(ele_id);
2067
2068 for (auto& ip_data : _ip_data)
2069 {
2070 ip_data.porosity = ip_data.porosity_prev;
2071
2072 _process_data.chemical_solver_interface
2073 ->updatePorosityPostReaction(ip_data.chemical_system_id,
2074 medium, ip_data.porosity);
2075 }
2076
2077 (*_process_data.mesh_prop_porosity)[ele_id] =
2078 std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
2079 [](double const s, auto const& ip)
2080 { return s + ip.porosity; }) /
2081 n_integration_points;
2082 }
2083
2084 std::vector<GlobalIndexType> chemical_system_indices;
2085 chemical_system_indices.reserve(n_integration_points);
2086 std::transform(_ip_data.begin(), _ip_data.end(),
2087 std::back_inserter(chemical_system_indices),
2088 [](auto const& ip_data)
2089 { return ip_data.chemical_system_id; });
2090
2091 _process_data.chemical_solver_interface->computeSecondaryVariable(
2092 ele_id, chemical_system_indices);
2093 }

References _integration_method, _ip_data, and _process_data.

◆ 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 2031 of file ComponentTransportFEM.h.

2036 {
2037 auto const local_p =
2038 local_x.template segment<pressure_size>(pressure_index);
2039 auto const local_C = local_x.template segment<concentration_size>(
2040 first_concentration_index);
2041 auto const local_T = getLocalTemperature(t, local_x);
2042
2043 std::vector<double> ele_velocity;
2044 calculateIntPtDarcyVelocity(t, local_p, local_C, local_T, ele_velocity);
2045
2046 auto const n_integration_points =
2047 _integration_method.getNumberOfPoints();
2048 auto const ele_velocity_mat =
2049 MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
2050
2051 auto const ele_id = _element.getID();
2052 Eigen::Map<LocalVectorType>(
2053 &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
2054 GlobalDim) =
2055 ele_velocity_mat.rowwise().sum() / n_integration_points;
2056 }
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:1883
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, _process_data, calculateIntPtDarcyVelocity(), first_concentration_index, getLocalTemperature(), pressure_index, and MathLib::toMatrix().

◆ 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 1969 of file ComponentTransportFEM.h.

1972 {
1973 auto const local_p = Eigen::Map<const NodalVectorType>(
1974 &local_x[pressure_index], pressure_size);
1975 auto const local_C = Eigen::Map<const NodalVectorType>(
1976 &local_x[first_concentration_index], concentration_size);
1977
1978 // Eval shape matrices at given point
1979 // Note: Axial symmetry is set to false here, because we only need dNdx
1980 // here, which is not affected by axial symmetry.
1981 auto const shape_matrices =
1982 NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
1983 GlobalDim>(
1984 _element, false /*is_axially_symmetric*/,
1985 std::array{pnt_local_coords})[0];
1986
1987 ParameterLib::SpatialPosition pos;
1988 pos.setElementID(_element.getID());
1989 auto const& b =
1990 _process_data
1991 .projected_specific_body_force_vectors[_element.getID()];
1992
1993 MaterialPropertyLib::VariableArray vars;
1994
1995 auto const& medium =
1996 *_process_data.media_map.getMedium(_element.getID());
1997 auto const& phase = medium.phase("AqueousLiquid");
1998
1999 // local_x contains the local concentration and pressure values
2000 double c_int_pt;
2001 NumLib::shapeFunctionInterpolate(local_C, shape_matrices.N, c_int_pt);
2002 vars.concentration = c_int_pt;
2003
2004 double p_int_pt;
2005 NumLib::shapeFunctionInterpolate(local_p, shape_matrices.N, p_int_pt);
2006 vars.liquid_phase_pressure = p_int_pt;
2007
2008 // TODO (naumov) Temporary value not used by current material models.
2009 // Need extension of secondary variables interface.
2010 double const dt = std::numeric_limits<double>::quiet_NaN();
2011 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2012 medium[MaterialPropertyLib::PropertyType::permeability].value(
2013 vars, pos, t, dt));
2014
2015 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
2016 .template value<double>(vars, pos, t, dt);
2017 GlobalDimMatrixType const K_over_mu = K / mu;
2018
2019 GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
2020 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
2021 .template value<double>(vars, pos, t, dt);
2022 if (_process_data.has_gravity)
2023 {
2024 q += K_over_mu * rho_w * b;
2025 }
2026 Eigen::Vector3d flux(0.0, 0.0, 0.0);
2027 flux.head<GlobalDim>() = rho_w * q;
2028 return flux;
2029 }

References _element, _process_data, NumLib::computeShapeMatrices(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, pressure_index, pressure_size, ParameterLib::SpatialPosition::setElementID(), 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 2220 of file ComponentTransportFEM.h.

2225 {
2226 auto const& medium =
2227 *_process_data.media_map.getMedium(this->_element.getID());
2228 auto const& solid_phase = medium.phase("Solid");
2229
2230 auto const specific_heat_capacity_solid =
2231 solid_phase
2232 .property(
2233 MaterialPropertyLib::PropertyType::specific_heat_capacity)
2234 .template value<double>(vars, pos, t, dt);
2235
2236 auto const solid_density =
2237 solid_phase.property(MaterialPropertyLib::PropertyType::density)
2238 .template value<double>(vars, pos, t, dt);
2239
2240 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2241 fluid_density * specific_heat_capacity_fluid * porosity;
2242 }

References _process_data, MaterialPropertyLib::density, MeshLib::Element::getID(), 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 1808 of file ComponentTransportFEM.h.

1813 {
1814 assert(x.size() == dof_table.size());
1815
1816 auto const n_processes = x.size();
1817 std::vector<std::vector<double>> local_x;
1818 local_x.reserve(n_processes);
1819
1820 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1821 {
1822 auto const indices =
1823 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1824 assert(!indices.empty());
1825 local_x.push_back(x[process_id]->get(indices));
1826 }
1827
1828 // only one process, must be monolithic.
1829 if (n_processes == 1)
1830 {
1831 auto const local_p = Eigen::Map<const NodalVectorType>(
1832 &local_x[0][pressure_index], pressure_size);
1833 auto const local_C = Eigen::Map<const NodalVectorType>(
1834 &local_x[0][first_concentration_index], concentration_size);
1835 NodalVectorType local_T;
1836 local_T.setConstant(ShapeFunction::NPOINTS,
1837 std::numeric_limits<double>::quiet_NaN());
1838 int const temperature_index =
1839 _process_data.isothermal ? -1 : ShapeFunction::NPOINTS;
1840 if (temperature_index != -1)
1841 {
1842 local_T = Eigen::Map<const NodalVectorType>(
1843 &local_x[0][temperature_index], temperature_size);
1844 }
1845 return calculateIntPtDarcyVelocity(t, local_p, local_C, local_T,
1846 cache);
1847 }
1848
1849 // multiple processes, must be staggered.
1850 {
1851 constexpr int pressure_process_id = 0;
1852 int concentration_process_id = 1;
1853 // Normally temperature process is not there,
1854 // hence set the default temperature index to -1
1855 int temperature_process_id = -1;
1856
1857 // check whether temperature process exists
1858 if (!_process_data.isothermal)
1859 {
1860 // if temperature process exists, its id is 1
1861 temperature_process_id = 1;
1862 // then the concentration index shifts to 2
1863 concentration_process_id = 2;
1864 }
1865
1866 auto const local_p = Eigen::Map<const NodalVectorType>(
1867 &local_x[pressure_process_id][0], pressure_size);
1868 auto const local_C = Eigen::Map<const NodalVectorType>(
1869 &local_x[concentration_process_id][0], concentration_size);
1870 NodalVectorType local_T;
1871 local_T.setConstant(ShapeFunction::NPOINTS,
1872 std::numeric_limits<double>::quiet_NaN());
1873 if (temperature_process_id != -1)
1874 {
1875 local_T = Eigen::Map<const NodalVectorType>(
1876 &local_x[temperature_process_id][0], temperature_size);
1877 }
1878 return calculateIntPtDarcyVelocity(t, local_p, local_C, local_T,
1879 cache);
1880 }
1881 }
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 1678 of file ComponentTransportFEM.h.

1683 {
1684 assert(x.size() == dof_table.size());
1685
1686 auto const n_processes = x.size();
1687 std::vector<std::vector<double>> local_x;
1688 local_x.reserve(n_processes);
1689
1690 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1691 {
1692 auto const indices =
1693 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1694 assert(!indices.empty());
1695 local_x.push_back(x[process_id]->get(indices));
1696 }
1697
1698 // only one process, must be monolithic.
1699 if (n_processes == 1)
1700 {
1701 auto const local_p = Eigen::Map<const NodalVectorType>(
1702 &local_x[0][pressure_index], pressure_size);
1703 auto const local_C = Eigen::Map<const NodalVectorType>(
1704 &local_x[0][first_concentration_index], concentration_size);
1705 NodalVectorType local_T;
1706 local_T.setConstant(ShapeFunction::NPOINTS,
1707 std::numeric_limits<double>::quiet_NaN());
1708 int const temperature_index =
1709 _process_data.isothermal ? -1 : ShapeFunction::NPOINTS;
1710 if (temperature_index != -1)
1711 {
1712 local_T = Eigen::Map<const NodalVectorType>(
1713 &local_x[0][temperature_index], temperature_size);
1714 }
1715 return calculateIntPtLiquidDensity(t, local_p, local_C, local_T,
1716 cache);
1717 }
1718
1719 // multiple processes, must be staggered.
1720 {
1721 constexpr int pressure_process_id = 0;
1722 int concentration_process_id = 1;
1723 // Normally temperature process is not there,
1724 // hence set the default temperature index to -1
1725 int temperature_process_id = -1;
1726
1727 // check whether temperature process exists
1728 if (!_process_data.isothermal)
1729 {
1730 // if temperature process exists, its id is 1
1731 temperature_process_id = 1;
1732 // then the concentration index shifts to 2
1733 concentration_process_id = 2;
1734 }
1735
1736 auto const local_p = Eigen::Map<const NodalVectorType>(
1737 &local_x[pressure_process_id][0], pressure_size);
1738 auto const local_C = Eigen::Map<const NodalVectorType>(
1739 &local_x[concentration_process_id][0], concentration_size);
1740 NodalVectorType local_T;
1741 local_T.setConstant(ShapeFunction::NPOINTS,
1742 std::numeric_limits<double>::quiet_NaN());
1743 if (temperature_process_id != -1)
1744 {
1745 local_T = Eigen::Map<const NodalVectorType>(
1746 &local_x[temperature_process_id][0], temperature_size);
1747 }
1748 return calculateIntPtLiquidDensity(t, local_p, local_C, local_T,
1749 cache);
1750 }
1751 }
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:1753

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 2095 of file ComponentTransportFEM.h.

2099 {
2100 std::vector<double> local_x_vec;
2101
2102 auto const n_processes = x.size();
2103 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
2104 {
2105 auto const indices =
2106 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
2107 assert(!indices.empty());
2108 auto const local_solution = x[process_id]->get(indices);
2109 local_x_vec.insert(std::end(local_x_vec),
2110 std::begin(local_solution),
2111 std::end(local_solution));
2112 }
2113 auto const local_x = MathLib::toVector(local_x_vec);
2114
2115 auto const p = local_x.template segment<pressure_size>(pressure_index);
2116 auto const c = local_x.template segment<concentration_size>(
2117 first_concentration_index + component_id * concentration_size);
2118
2119 auto const n_integration_points =
2120 _integration_method.getNumberOfPoints();
2121
2122 cache.clear();
2123 auto cache_mat = MathLib::createZeroedMatrix<
2124 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
2125 cache, GlobalDim, n_integration_points);
2126
2127 ParameterLib::SpatialPosition pos;
2128 pos.setElementID(_element.getID());
2129
2130 auto const& b =
2131 _process_data
2132 .projected_specific_body_force_vectors[_element.getID()];
2133
2134 MaterialPropertyLib::VariableArray vars;
2135
2136 auto const& medium =
2137 *_process_data.media_map.getMedium(_element.getID());
2138 auto const& phase = medium.phase("AqueousLiquid");
2139
2140 auto const& component = phase.component(
2141 _transport_process_variables[component_id].get().getName());
2142
2143 auto const& Ns =
2144 _process_data.shape_matrix_cache
2145 .NsHigherOrder<typename ShapeFunction::MeshElement>();
2146
2147 for (unsigned ip = 0; ip < n_integration_points; ++ip)
2148 {
2149 auto const& ip_data = _ip_data[ip];
2150 auto const& dNdx = ip_data.dNdx;
2151 auto const& N = Ns[ip];
2152 auto const& phi = ip_data.porosity;
2153
2154 double const p_ip = N.dot(p);
2155 double const c_ip = N.dot(c);
2156
2157 vars.concentration = c_ip;
2158 vars.liquid_phase_pressure = p_ip;
2159 vars.porosity = phi;
2160
2161 double const dt = std::numeric_limits<double>::quiet_NaN();
2162
2163 auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2164 medium[MaterialPropertyLib::PropertyType::permeability].value(
2165 vars, pos, t, dt));
2166 auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
2167 .template value<double>(vars, pos, t, dt);
2168 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
2169 .template value<double>(vars, pos, t, dt);
2170
2171 // Darcy flux
2172 GlobalDimVectorType const q =
2173 _process_data.has_gravity
2174 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
2175 : GlobalDimVectorType(-k / mu * dNdx * p);
2176
2177 auto const alpha_T = medium.template value<double>(
2178 MaterialPropertyLib::PropertyType::transversal_dispersivity);
2179 auto const alpha_L = medium.template value<double>(
2180 MaterialPropertyLib::PropertyType::longitudinal_dispersivity);
2181 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
2182 component[MaterialPropertyLib::PropertyType::pore_diffusion]
2183 .value(vars, pos, t, dt));
2184
2185 // Hydrodynamic dispersion
2186 GlobalDimMatrixType const D = NumLib::computeHydrodynamicDispersion(
2187 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
2188 alpha_L);
2189
2190 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
2191 }
2192
2193 return cache;
2194 }
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, NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), NumLib::getIndices(), getName(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::longitudinal_dispersivity, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::VariableArray::porosity, pressure_index, ParameterLib::SpatialPosition::setElementID(), 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 2285 of file ComponentTransportFEM.h.

2287 {
2288 NodalVectorType local_T;
2289 if (_process_data.isothermal)
2290 {
2291 if (_process_data.temperature)
2292 {
2293 local_T = _process_data.temperature->getNodalValuesOnElement(
2294 _element, t);
2295 }
2296 else
2297 {
2298 local_T = NodalVectorType::Zero(temperature_size);
2299 }
2300 }
2301 else
2302 {
2303 local_T =
2304 local_x.template segment<temperature_size>(temperature_index);
2305 }
2306 return local_T;
2307 }

References _element, _process_data, temperature_index, and temperature_size.

Referenced by assembleComponentTransportEquation(), assembleHeatTransportEquation(), assembleHydraulicEquation(), 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 1959 of file ComponentTransportFEM.h.

1961 {
1962 auto const& N = _process_data.shape_matrix_cache.NsHigherOrder<
1963 typename ShapeFunction::MeshElement>()[integration_point];
1964
1965 // assumes N is stored contiguously in memory
1966 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1967 }

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 2244 of file ComponentTransportFEM.h.

2250 {
2251 auto const& medium =
2252 *_process_data.media_map.getMedium(_element.getID());
2253
2254 auto thermal_conductivity =
2255 MaterialPropertyLib::formEigenTensor<GlobalDim>(
2256 medium
2257 .property(
2258 MaterialPropertyLib::PropertyType::thermal_conductivity)
2259 .value(vars, pos, t, dt));
2260
2261 auto const thermal_dispersivity_transversal =
2262 medium
2263 .property(MaterialPropertyLib::PropertyType::
2264 thermal_transversal_dispersivity)
2265 .template value<double>();
2266
2267 auto const thermal_dispersivity_longitudinal =
2268 medium
2269 .property(MaterialPropertyLib::PropertyType::
2270 thermal_longitudinal_dispersivity)
2271 .template value<double>();
2272
2273 // Thermal conductivity is moved outside and zero matrix is passed
2274 // instead due to multiplication with fluid's density times specific
2275 // heat capacity.
2276 return thermal_conductivity +
2277 fluid_density * specific_heat_capacity_fluid *
2278 NumLib::computeHydrodynamicDispersion(
2279 _process_data.stabilizer, _element.getID(),
2280 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2281 velocity, 0 /* phi */, thermal_dispersivity_transversal,
2282 thermal_dispersivity_longitudinal);
2283 }

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 312 of file ComponentTransportFEM.h.

314 {
315 assert(_process_data.chemical_solver_interface);
316
317 auto const& medium =
318 *_process_data.media_map.getMedium(_element.getID());
319
320 ParameterLib::SpatialPosition pos;
321 pos.setElementID(_element.getID());
322
323 auto const& Ns =
324 _process_data.shape_matrix_cache
325 .NsHigherOrder<typename ShapeFunction::MeshElement>();
326
327 unsigned const n_integration_points =
328 _integration_method.getNumberOfPoints();
329
330 for (unsigned ip = 0; ip < n_integration_points; ip++)
331 {
332 auto& ip_data = _ip_data[ip];
333 auto const& N = Ns[ip];
334 auto const& chemical_system_id = ip_data.chemical_system_id;
335
336 // set position with N as the shape matrix at the current
337 // integration point
338 pos.setCoordinates(MathLib::Point3d(
339 NumLib::interpolateCoordinates<ShapeFunction,
340 ShapeMatricesType>(_element,
341 N)));
342
343 auto const n_component = _transport_process_variables.size();
344 std::vector<double> C_int_pt(n_component);
345 for (unsigned component_id = 0; component_id < n_component;
346 ++component_id)
347 {
348 auto const concentration_index =
349 first_concentration_index +
350 component_id * concentration_size;
351 auto const local_C =
352 local_x.template segment<concentration_size>(
353 concentration_index);
354
355 NumLib::shapeFunctionInterpolate(local_C, N,
356 C_int_pt[component_id]);
357 }
358
359 _process_data.chemical_solver_interface
360 ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id,
361 medium, pos, t);
362 }
363 }

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, concentration_size, first_concentration_index, NumLib::interpolateCoordinates(), ParameterLib::SpatialPosition::setCoordinates(), ParameterLib::SpatialPosition::setElementID(), 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 430 of file ComponentTransportFEM.h.

432 {
433 if (!_process_data.chemically_induced_porosity_change)
434 {
435 return;
436 }
437
438 auto const& medium = *_process_data.media_map.getMedium(ele_id);
439
440 ParameterLib::SpatialPosition pos;
441 pos.setElementID(ele_id);
442
443 for (auto& ip_data : _ip_data)
444 {
445 ip_data.porosity = ip_data.porosity_prev;
446
447 _process_data.chemical_solver_interface
448 ->updateVolumeFractionPostReaction(ip_data.chemical_system_id,
449 medium, pos,
450 ip_data.porosity, t, dt);
451
452 _process_data.chemical_solver_interface->updatePorosityPostReaction(
453 ip_data.chemical_system_id, medium, ip_data.porosity);
454 }
455 }

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 2196 of file ComponentTransportFEM.h.

2200 {
2201 unsigned const n_integration_points =
2202 _integration_method.getNumberOfPoints();
2203
2204 for (unsigned ip = 0; ip < n_integration_points; ip++)
2205 {
2206 _ip_data[ip].pushBackState();
2207 }
2208 }

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 365 of file ComponentTransportFEM.h.

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, concentration_size, first_concentration_index, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, ParameterLib::SpatialPosition::setElementID(), 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 292 of file ComponentTransportFEM.h.

293 {
294 assert(_process_data.chemical_solver_interface);
295 // chemical system index map
296 auto& chemical_system_index_map =
297 _process_data.chemical_solver_interface->chemical_system_index_map;
298
299 unsigned const n_integration_points =
300 _integration_method.getNumberOfPoints();
301 for (unsigned ip = 0; ip < n_integration_points; ip++)
302 {
303 _ip_data[ip].chemical_system_id =
304 chemical_system_index_map.empty()
305 ? 0
306 : chemical_system_index_map.back() + 1;
307 chemical_system_index_map.push_back(
308 _ip_data[ip].chemical_system_id);
309 }
310 }

References _integration_method, _ip_data, and _process_data.

Member Data Documentation

◆ _element

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

Definition at line 2211 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 2214 of file ComponentTransportFEM.h.

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

◆ _ip_data

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

Definition at line 2218 of file ComponentTransportFEM.h.

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

◆ _process_data

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

Definition at line 2212 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(), and setChemicalSystemID().

◆ _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 2216 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(), 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: