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:2229
ProcessLib::ComponentTransport::LocalAssemblerData::_ip_data
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data
Definition ComponentTransportFEM.h:2233
ProcessLib::ComponentTransport::LocalAssemblerData::_process_data
ComponentTransportProcessData const & _process_data
Definition ComponentTransportFEM.h:2227
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:2231

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::detail::assembleAdvectionMatrix(), NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, MaterialPropertyLib::formEigenTensor(), getName(), NumLib::interpolateCoordinates(), MaterialPropertyLib::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::permeability, MaterialPropertyLib::pore_diffusion, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, MaterialPropertyLib::retardation_factor, 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 1126 of file ComponentTransportFEM.h.

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::detail::assembleAdvectionMatrix(), NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getLocalTemperature(), getName(), 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 821 of file ComponentTransportFEM.h.

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

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

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::detail::assembleAdvectionMatrix(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, MaterialPropertyLib::VariableArray::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getHeatEnergyCoefficient(), getLocalTemperature(), getThermalConductivityDispersivity(), NumLib::interpolateCoordinates(), MaterialPropertyLib::VariableArray::liquid_phase_pressure, MaterialPropertyLib::VariableArray::liquid_saturation, MaterialPropertyLib::permeability, MaterialPropertyLib::porosity, MaterialPropertyLib::VariableArray::porosity, pressure_index, pressure_size, 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 849 of file ComponentTransportFEM.h.

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

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), getLocalTemperature(), 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 766 of file ComponentTransportFEM.h.

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, 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 1609 of file ComponentTransportFEM.h.

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

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

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::detail::assembleAdvectionMatrix(), NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::decay_rate, MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), 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 1331 of file ComponentTransportFEM.h.

1336 {
1337 if (process_id == _process_data.hydraulic_process_id)
1338 {
1339 assembleWithJacobianHydraulicEquation(t, dt, local_x, local_x_prev,
1340 local_b_data, local_Jac_data);
1341 }
1342 else
1343 {
1344 int const component_id = process_id - 1;
1345 assembleWithJacobianComponentTransportEquation(
1346 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1347 component_id);
1348 }
1349 }
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:1351
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:1460

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

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

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::concentration, MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MathLib::createZeroedVector(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), 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 1895 of file ComponentTransportFEM.h.

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

References _element, _integration_method, _ip_data, _process_data, MaterialPropertyLib::AqueousLiquid, 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 1764 of file ComponentTransportFEM.h.

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

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

Referenced by getIntPtLiquidDensity().

◆ computeReactionRelatedSecondaryVariable()

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

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 2072 of file ComponentTransportFEM.h.

2074 {
2075 auto const n_integration_points =
2076 _integration_method.getNumberOfPoints();
2077
2078 if (_process_data.chemically_induced_porosity_change)
2079 {
2080 auto const& medium = *_process_data.media_map.getMedium(ele_id);
2081
2082 for (auto& ip_data : _ip_data)
2083 {
2084 ip_data.porosity = ip_data.porosity_prev;
2085
2086 _process_data.chemical_solver_interface
2087 ->updatePorosityPostReaction(ip_data.chemical_system_id,
2088 medium, ip_data.porosity);
2089 }
2090
2091 (*_process_data.mesh_prop_porosity)[ele_id] =
2092 std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
2093 [](double const s, auto const& ip)
2094 { return s + ip.porosity; }) /
2095 n_integration_points;
2096 }
2097
2098 std::vector<GlobalIndexType> chemical_system_indices;
2099 chemical_system_indices.reserve(n_integration_points);
2100 std::transform(_ip_data.begin(), _ip_data.end(),
2101 std::back_inserter(chemical_system_indices),
2102 [](auto const& ip_data)
2103 { return ip_data.chemical_system_id; });
2104
2105 _process_data.chemical_solver_interface->computeSecondaryVariable(
2106 ele_id, chemical_system_indices);
2107 }

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

2050 {
2051 auto const local_p =
2052 local_x.template segment<pressure_size>(pressure_index);
2053 auto const local_C = local_x.template segment<concentration_size>(
2054 first_concentration_index);
2055 auto const local_T = getLocalTemperature(t, local_x);
2056
2057 std::vector<double> ele_velocity;
2058 calculateIntPtDarcyVelocity(t, local_p, local_C, local_T, ele_velocity);
2059
2060 auto const n_integration_points =
2061 _integration_method.getNumberOfPoints();
2062 auto const ele_velocity_mat =
2063 MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
2064
2065 auto const ele_id = _element.getID();
2066 Eigen::Map<LocalVectorType>(
2067 &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
2068 GlobalDim) =
2069 ele_velocity_mat.rowwise().sum() / n_integration_points;
2070 }
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:1895
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 1982 of file ComponentTransportFEM.h.

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

References _element, _process_data, MaterialPropertyLib::AqueousLiquid, 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 2235 of file ComponentTransportFEM.h.

2240 {
2241 auto const& medium =
2242 *_process_data.media_map.getMedium(this->_element.getID());
2243 auto const& solid_phase =
2244 medium.phase(MaterialPropertyLib::PhaseName::Solid);
2245
2246 auto const specific_heat_capacity_solid =
2247 solid_phase
2248 .property(
2249 MaterialPropertyLib::PropertyType::specific_heat_capacity)
2250 .template value<double>(vars, pos, t, dt);
2251
2252 auto const solid_density =
2253 solid_phase.property(MaterialPropertyLib::PropertyType::density)
2254 .template value<double>(vars, pos, t, dt);
2255
2256 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2257 fluid_density * specific_heat_capacity_fluid * porosity;
2258 }

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

Referenced by assembleHeatTransportEquation().

◆ getIntPtDarcyVelocity()

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

Implements ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface.

Definition at line 1820 of file ComponentTransportFEM.h.

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

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

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

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

References _element, _integration_method, _ip_data, _process_data, _transport_process_variables, MaterialPropertyLib::AqueousLiquid, NumLib::computeHydrodynamicDispersion(), MaterialPropertyLib::VariableArray::concentration, concentration_size, MathLib::createZeroedMatrix(), MaterialPropertyLib::density, first_concentration_index, MaterialPropertyLib::formEigenTensor(), NumLib::getIndices(), getName(), 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 2301 of file ComponentTransportFEM.h.

2303 {
2304 NodalVectorType local_T;
2305 if (_process_data.isothermal)
2306 {
2307 if (_process_data.temperature)
2308 {
2309 local_T = _process_data.temperature->getNodalValuesOnElement(
2310 _element, t);
2311 }
2312 else
2313 {
2314 local_T = NodalVectorType::Zero(temperature_size);
2315 }
2316 }
2317 else
2318 {
2319 local_T =
2320 local_x.template segment<temperature_size>(temperature_index);
2321 }
2322 return local_T;
2323 }

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

1974 {
1975 auto const& N = _process_data.shape_matrix_cache.NsHigherOrder<
1976 typename ShapeFunction::MeshElement>()[integration_point];
1977
1978 // assumes N is stored contiguously in memory
1979 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1980 }

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

2266 {
2267 auto const& medium =
2268 *_process_data.media_map.getMedium(_element.getID());
2269
2270 auto thermal_conductivity =
2271 MaterialPropertyLib::formEigenTensor<GlobalDim>(
2272 medium
2273 .property(
2274 MaterialPropertyLib::PropertyType::thermal_conductivity)
2275 .value(vars, pos, t, dt));
2276
2277 auto const thermal_dispersivity_transversal =
2278 medium
2279 .property(MaterialPropertyLib::PropertyType::
2280 thermal_transversal_dispersivity)
2281 .template value<double>();
2282
2283 auto const thermal_dispersivity_longitudinal =
2284 medium
2285 .property(MaterialPropertyLib::PropertyType::
2286 thermal_longitudinal_dispersivity)
2287 .template value<double>();
2288
2289 // Thermal conductivity is moved outside and zero matrix is passed
2290 // instead due to multiplication with fluid's density times specific
2291 // heat capacity.
2292 return thermal_conductivity +
2293 fluid_density * specific_heat_capacity_fluid *
2294 NumLib::computeHydrodynamicDispersion(
2295 _process_data.stabilizer, _element.getID(),
2296 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2297 velocity, 0 /* phi */, thermal_dispersivity_transversal,
2298 thermal_dispersivity_longitudinal);
2299 }

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

2215 {
2216 unsigned const n_integration_points =
2217 _integration_method.getNumberOfPoints();
2218
2219 for (unsigned ip = 0; ip < n_integration_points; ip++)
2220 {
2221 _ip_data[ip].pushBackState();
2222 }
2223 }

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 2226 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 2229 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 2233 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 2227 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 2231 of file ComponentTransportFEM.h.

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

◆ concentration_size

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

Definition at line 215 of file ComponentTransportFEM.h.

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

◆ first_concentration_index

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

Definition at line 211 of file ComponentTransportFEM.h.

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

◆ pressure_index

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

Definition at line 209 of file ComponentTransportFEM.h.

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

◆ pressure_size

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

Definition at line 213 of file ComponentTransportFEM.h.

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

◆ temperature_index

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

Definition at line 210 of file ComponentTransportFEM.h.

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

◆ temperature_size

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

Definition at line 214 of file ComponentTransportFEM.h.

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

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