OGS
ProcessLib::ComponentTransport::ComponentTransportProcess Class Referencefinal

# ComponentTransport process

## Governing equations

The flow process is described by

$\phi \frac{\partial \rho}{\partial p} \frac{\partial p}{\partial t} + \phi \frac{\partial \rho}{\partial C} \frac{\partial C}{\partial t} - \nabla \cdot \left[\frac{\kappa}{\mu(C)} \rho \nabla \left( p + \rho g z \right)\right] + Q_p = 0,$

where the storage $$S$$ has been substituted by $$\phi \frac{\partial \rho}{\partial p}$$, $$\phi$$ is the porosity, $$C$$ is the concentration, $$p$$ is the pressure, $$\kappa$$ is permeability, $$\mu$$ is viscosity of the fluid, $$\rho$$ is the density of the fluid, and $$g$$ is the gravitational acceleration.

The mass transport process is described by

$\phi R C \frac{\partial \rho}{\partial p} \frac{\partial p}{\partial t} + \phi R \left(\rho + C \frac{\partial \rho}{\partial C}\right) \frac{\partial C}{\partial t} - \nabla \cdot \left[\frac{\kappa}{\mu(C)} \rho C \nabla \left( p + \rho g z \right) + \rho D \nabla C\right] + Q_C + R \vartheta \phi \rho C = 0,$

where $$R$$ is the retardation factor, $$\vec{q} = -\frac{\kappa}{\mu(C)} \nabla \left( p + \rho g z \right)$$ is the Darcy velocity, $$D$$ is the hydrodynamic dispersion tensor, $$\vartheta$$ is the decay rate.

For the hydrodynamic dispersion tensor the relation

$D = (\phi D_d + \beta_T \|\vec{q}\|) I + (\beta_L - \beta_T) \frac{\vec{q} \vec{q}^T}{\|\vec{q}\|}$

is implemented, where $$D_d$$ is the molecular diffusion coefficient, $$\beta_L$$ the longitudinal dispersivity of chemical species, and $$\beta_T$$ the transverse dispersivity of chemical species.

The implementation uses a monolithic approach, i.e., both processes are assembled within one global system of equations.

## Process Coupling

The advective term of the concentration equation is given by the confined groundwater flow process, i.e., the concentration distribution depends on Darcy velocity of the groundwater flow process. On the other hand the concentration dependencies of the viscosity and density in the groundwater flow couples the H process to the C process.

Note
At the moment there is not any coupling by source or sink terms, i.e., the coupling is implemented only by density and viscosity changes due to concentration changes as well as by the temporal derivatives of each variable.

Definition at line 94 of file ComponentTransportProcess.h.

#include <ComponentTransportProcess.h>

Inheritance diagram for ProcessLib::ComponentTransport::ComponentTransportProcess:
[legend]
Collaboration diagram for ProcessLib::ComponentTransport::ComponentTransportProcess:
[legend]

## Public Member Functions

ComponentTransportProcess (std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > &&process_variables, ComponentTransportProcessData &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme, std::unique_ptr< ProcessLib::SurfaceFluxData > &&surfaceflux, std::unique_ptr< ChemistryLib::ChemicalSolverInterface > &&chemical_solver_interface)

Eigen::Vector3d getFlux (std::size_t const element_id, MathLib::Point3d const &p, double const t, std::vector< GlobalVector * > const &x) const override

void solveReactionEquation (std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, double const t, double const dt, NumLib::EquationSystem &ode_sys, int const process_id) override

void computeSecondaryVariableConcrete (double const, double const, std::vector< GlobalVector * > const &x, GlobalVector const &, int const) override

void postTimestepConcreteProcess (std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) override

ODESystem interface
bool isLinear () const override

Public Member Functions inherited from ProcessLib::Process
Process (std::string name_, MeshLib::Mesh &mesh, std::unique_ptr< AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > &&process_variables, SecondaryVariableCollection &&secondary_variables, const bool use_monolithic_scheme=true)

void preTimestep (std::vector< GlobalVector * > const &x, const double t, const double delta_t, const int process_id)
Preprocessing before starting assembly for new timestep. More...

void postTimestep (std::vector< GlobalVector * > const &x, const double t, const double delta_t, int const process_id)
Postprocessing after a complete timestep. More...

void postNonLinearSolver (GlobalVector const &x, GlobalVector const &xdot, const double t, double const dt, int const process_id)

void preIteration (const unsigned iter, GlobalVector const &x) final

void computeSecondaryVariable (double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id)
compute secondary variables for the coupled equations or for output. More...

NumLib::IterationResult postIteration (GlobalVector const &x) final

void initialize ()

void setInitialConditions (std::vector< GlobalVector * > &process_solutions, std::vector< GlobalVector * > const &process_solutions_prev, double const t, int const process_id)

MathLib::MatrixSpecifications getMatrixSpecifications (const int process_id) const override

void setCoupledSolutionsForStaggeredScheme (CoupledSolutionsForStaggeredScheme *const coupled_solutions)

void updateDeactivatedSubdomains (double const time, const int process_id)

bool isMonolithicSchemeUsed () const

virtual void setCoupledTermForTheStaggeredSchemeToLocalAssemblers (int const)

virtual void extrapolateIntegrationPointValuesToNodes (const double, std::vector< GlobalVector * > const &, std::vector< GlobalVector * > &)

void preAssemble (const double t, double const dt, GlobalVector const &x) final

void assemble (const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) final

void assembleWithJacobian (const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) final

std::vector< NumLib::IndexValueVector< GlobalIndexType > > const * getKnownSolutions (double const t, GlobalVector const &x, int const process_id) const final

virtual NumLib::LocalToGlobalIndexMap const & getDOFTable (const int) const

MeshLib::MeshgetMesh () const

std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables (const int process_id) const

SecondaryVariableCollection const & getSecondaryVariables () const

std::vector< std::unique_ptr< IntegrationPointWriter > > const & getIntegrationPointWriters () const

virtual Eigen::Vector3d getFlux (std::size_t, MathLib::Point3d const &, double const, std::vector< GlobalVector * > const &) const

virtual void solveReactionEquation (std::vector< GlobalVector * > &, std::vector< GlobalVector * > const &, double const, double const, NumLib::EquationSystem &, int const)

## Private Member Functions

void initializeConcreteProcess (NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize(). More...

void setInitialConditionsConcreteProcess (std::vector< GlobalVector * > &x, double const t, int const process_id) override

void assembleConcreteProcess (const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override

void assembleWithJacobianConcreteProcess (const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override

## Private Attributes

ComponentTransportProcessData _process_data

std::vector< std::unique_ptr< ComponentTransportLocalAssemblerInterface > > _local_assemblers

std::unique_ptr< ProcessLib::SurfaceFluxData_surfaceflux

std::unique_ptr< ChemistryLib::ChemicalSolverInterface_chemical_solver_interface

MeshLib::PropertyVector< double > * _molar_flow_rate = nullptr

Public Types inherited from ProcessLib::Process
using NonlinearSolver = NumLib::NonlinearSolverBase

using TimeDiscretization = NumLib::TimeDiscretization

Public Attributes inherited from ProcessLib::Process
std::string const name

Protected Member Functions inherited from ProcessLib::Process
NumLib::ExtrapolatorgetExtrapolator () const

NumLib::LocalToGlobalIndexMap const & getSingleComponentDOFTable () const

void initializeProcessBoundaryConditionsAndSourceTerms (const NumLib::LocalToGlobalIndexMap &dof_table, const int process_id)

virtual void constructDofTable ()

void constructMonolithicProcessDofTable ()

void constructDofTableOfSpecifiedProcessStaggeredScheme (const int specified_prosess_id)

virtual std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData () const

std::vector< GlobalIndexTypegetIndicesOfResiduumWithoutInitialCompensation () const override

Protected Attributes inherited from ProcessLib::Process
MeshLib::Mesh_mesh

std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes

std::unique_ptr< NumLib::LocalToGlobalIndexMap_local_to_global_index_map

SecondaryVariableCollection _secondary_variables

VectorMatrixAssembler _global_assembler

const bool _use_monolithic_scheme

CoupledSolutionsForStaggeredScheme_coupled_solutions

unsigned const _integration_order

std::vector< std::unique_ptr< IntegrationPointWriter > > _integration_point_writer

GlobalSparsityPattern _sparsity_pattern

std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > _process_variables

std::vector< BoundaryConditionCollection_boundary_conditions

## ◆ ComponentTransportProcess()

 ProcessLib::ComponentTransport::ComponentTransportProcess::ComponentTransportProcess ( std::string name, MeshLib::Mesh & mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > && jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const & parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > && process_variables, ComponentTransportProcessData && process_data, SecondaryVariableCollection && secondary_variables, bool const use_monolithic_scheme, std::unique_ptr< ProcessLib::SurfaceFluxData > && surfaceflux, std::unique_ptr< ChemistryLib::ChemicalSolverInterface > && chemical_solver_interface )

Definition at line 30 of file ComponentTransportProcess.cpp.

44 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
45 integration_order, std::move(process_variables),
46 std::move(secondary_variables), use_monolithic_scheme),
47 _process_data(std::move(process_data)),
48 _surfaceflux(std::move(surfaceflux)),
49 _chemical_solver_interface(std::move(chemical_solver_interface))
50{
51 // Todo(renchao): Need to adapt for the multi-component case.
52 _molar_flow_rate = MeshLib::getOrCreateMeshProperty<double>(
53 mesh, "MolarFlowRate", MeshLib::MeshItemType::Node, 1);
54}
std::unique_ptr< ChemistryLib::ChemicalSolverInterface > _chemical_solver_interface
std::unique_ptr< ProcessLib::SurfaceFluxData > _surfaceflux
std::string const name
Definition: Process.h:325
Process(std::string name_, MeshLib::Mesh &mesh, std::unique_ptr< AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > &&process_variables, SecondaryVariableCollection &&secondary_variables, const bool use_monolithic_scheme=true)
Definition: Process.cpp:22

References _molar_flow_rate, and MeshLib::Node.

## ◆ assembleConcreteProcess()

 void ProcessLib::ComponentTransport::ComponentTransportProcess::assembleConcreteProcess ( const double t, double const dt, std::vector< GlobalVector * > const & x, std::vector< GlobalVector * > const & xdot, int const process_id, GlobalMatrix & M, GlobalMatrix & K, GlobalVector & b )
overrideprivatevirtual

Implements ProcessLib::Process.

Definition at line 147 of file ComponentTransportProcess.cpp.

151{
152 DBUG("Assemble ComponentTransportProcess.");
153
154 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
155
156 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
157 dof_tables;
159 {
160 dof_tables.push_back(std::ref(*_local_to_global_index_map));
161 }
162 else
163 {
164 std::generate_n(
165 std::back_inserter(dof_tables), _process_variables.size(),
166 [&]() { return std::ref(*_local_to_global_index_map); });
167 }
168 // Call global assembler for each local assembly item.
171 pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
172 b);
173}
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:29
std::vector< std::unique_ptr< ComponentTransportLocalAssemblerInterface > > _local_assemblers
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > _process_variables
Definition: Process.h:362
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:146
VectorMatrixAssembler _global_assembler
Definition: Process.h:335
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:331
const bool _use_monolithic_scheme
Definition: Process.h:337
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
static const double t
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)

## ◆ assembleWithJacobianConcreteProcess()

 void ProcessLib::ComponentTransport::ComponentTransportProcess::assembleWithJacobianConcreteProcess ( const double t, double const dt, std::vector< GlobalVector * > const & x, std::vector< GlobalVector * > const & xdot, int const process_id, GlobalMatrix & M, GlobalMatrix & K, GlobalVector & b, GlobalMatrix & Jac )
overrideprivatevirtual

Implements ProcessLib::Process.

Definition at line 175 of file ComponentTransportProcess.cpp.

179{
180 DBUG("AssembleWithJacobian ComponentTransportProcess.");
181
182 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
183
184 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
185 dof_tables;
186
187 assert(!_use_monolithic_scheme);
188 std::generate_n(std::back_inserter(dof_tables), _process_variables.size(),
189 [&]() { return std::ref(*_local_to_global_index_map); });
190
191 // Call global assembler for each local assembly item.
194 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
195 process_id, M, K, b, Jac);
196
197 if (process_id == _process_data.hydraulic_process_id)
198 {
199 return;
200 }
201
202 // b is the negated residumm used in the Newton's method.
203 // Here negating b is to recover the primitive residuum.
204 transformVariableFromGlobalVector(b, 0 /*variable id*/,
206 *_molar_flow_rate, std::negate<double>());
207}
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)
void transformVariableFromGlobalVector(GlobalVector const &input_vector, int const variable_id, NumLib::LocalToGlobalIndexMap const &local_to_global_index_map, MeshLib::PropertyVector< double > &output_vector, Functor mapFunction)
Definition: DOFTableUtil.h:59

## ◆ computeSecondaryVariableConcrete()

 void ProcessLib::ComponentTransport::ComponentTransportProcess::computeSecondaryVariableConcrete ( double const t, double const dt, std::vector< GlobalVector * > const & x, GlobalVector const & x_dot, int const process_id )
overridevirtual

Reimplemented from ProcessLib::Process.

Definition at line 340 of file ComponentTransportProcess.cpp.

346{
347 if (process_id != 0)
348 {
349 return;
350 }
351
352 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
353 dof_tables.reserve(x.size());
354 std::generate_n(std::back_inserter(dof_tables), x.size(),
355 [&]() { return _local_to_global_index_map.get(); });
356
357 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
360 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x,
361 x_dot, process_id);
362
364 {
365 return;
366 }
367
369 &ComponentTransportLocalAssemblerInterface::
370 computeReactionRelatedSecondaryVariable,
372}
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_dot, int const process_id)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)

## ◆ getFlux()

 Eigen::Vector3d ProcessLib::ComponentTransport::ComponentTransportProcess::getFlux ( std::size_t const element_id, MathLib::Point3d const & p, double const t, std::vector< GlobalVector * > const & x ) const
overridevirtual

Reimplemented from ProcessLib::Process.

Definition at line 209 of file ComponentTransportProcess.cpp.

214{
215 std::vector<GlobalIndexType> indices_cache;
216 auto const r_c_indices = NumLib::getRowColumnIndices(
217 element_id, *_local_to_global_index_map, indices_cache);
218
219 std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
220 x.size(), r_c_indices.rows};
221 auto const local_xs =
222 getCoupledLocalSolutions(x, indices_of_all_coupled_processes);
223
224 return _local_assemblers[element_id]->getFlux(p, t, local_xs);
225}
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalIndexType > &indices)
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType > > const &indices)

## ◆ initializeConcreteProcess()

 void ProcessLib::ComponentTransport::ComponentTransportProcess::initializeConcreteProcess ( NumLib::LocalToGlobalIndexMap const & dof_table, MeshLib::Mesh const & mesh, unsigned const integration_order )
overrideprivatevirtual

Process specific initialization called by initialize().

Implements ProcessLib::Process.

Definition at line 56 of file ComponentTransportProcess.cpp.

60{
61 _process_data.mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
62 const_cast<MeshLib::Mesh&>(mesh), "velocity",
63 MeshLib::MeshItemType::Cell, mesh.getDimension());
64
65 _process_data.mesh_prop_porosity = MeshLib::getOrCreateMeshProperty<double>(
66 const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
68
69 std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
70 transport_process_variables;
72 {
73 const int process_id = 0;
74 for (auto pv_iter = std::next(_process_variables[process_id].begin());
75 pv_iter != _process_variables[process_id].end();
76 ++pv_iter)
77 {
78 transport_process_variables.push_back(*pv_iter);
79 }
80 }
81 else
82 {
83 for (auto pv_iter = std::next(_process_variables.begin());
84 pv_iter != _process_variables.end();
85 ++pv_iter)
86 {
87 transport_process_variables.push_back((*pv_iter)[0]);
88 }
89 }
90
91 ProcessLib::createLocalAssemblers<LocalAssemblerData>(
92 mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
93 mesh.isAxiallySymmetric(), integration_order, _process_data,
94 transport_process_variables);
95
97 {
101
102 _chemical_solver_interface->initialize();
103 }
104
106 "darcy_velocity",
108 mesh.getDimension(), getExtrapolator(), _local_assemblers,
110
112 "molar_flux",
114 mesh.getDimension(), getExtrapolator(), _local_assemblers,
116}
virtual 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 =0
virtual std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
SecondaryVariableCollection _secondary_variables
Definition: Process.h:333
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:182
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)

## ◆ isLinear()

 bool ProcessLib::ComponentTransport::ComponentTransportProcess::isLinear ( ) const
inlineoverride

Definition at line 117 of file ComponentTransportProcess.h.

117{ return false; }

## ◆ postTimestepConcreteProcess()

 void ProcessLib::ComponentTransport::ComponentTransportProcess::postTimestepConcreteProcess ( std::vector< GlobalVector * > const & x, const double t, const double dt, int const process_id )
overridevirtual

Reimplemented from ProcessLib::Process.

Definition at line 374 of file ComponentTransportProcess.cpp.

379{
380 if (process_id != 0)
381 {
382 return;
383 }
384
385 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
386 dof_tables.reserve(x.size());
387 std::generate_n(std::back_inserter(dof_tables), x.size(),
388 [&]() { return _local_to_global_index_map.get(); });
389
390 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
393 _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt);
394
395 if (!_surfaceflux) // computing the surfaceflux is optional
396 {
397 return;
398 }
399 _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
401}
virtual void postTimestep(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)
MeshLib::Mesh & _mesh
Definition: Process.h:328
unsigned const _integration_order
Definition: Process.h:346

## ◆ setInitialConditionsConcreteProcess()

 void ProcessLib::ComponentTransport::ComponentTransportProcess::setInitialConditionsConcreteProcess ( std::vector< GlobalVector * > & x, double const t, int const process_id )
overrideprivatevirtual

Reimplemented from ProcessLib::Process.

Definition at line 118 of file ComponentTransportProcess.cpp.

120{
122 {
123 return;
124 }
125
126 if (process_id != static_cast<int>(x.size() - 1))
127 {
128 return;
129 }
130
131 std::for_each(
132 x.begin(), x.end(),
133 [](auto const process_solution)
134 { MathLib::LinAlg::setLocalAccessibleVector(*process_solution); });
135
136 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
137 dof_tables.reserve(x.size());
138 std::generate_n(std::back_inserter(dof_tables), x.size(),
139 [&]() { return _local_to_global_index_map.get(); });
140
144 dof_tables, x, t);
145}
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)

## ◆ solveReactionEquation()

 void ProcessLib::ComponentTransport::ComponentTransportProcess::solveReactionEquation ( std::vector< GlobalVector * > & x, std::vector< GlobalVector * > const & x_prev, double const t, double const dt, NumLib::EquationSystem & ode_sys, int const process_id )
overridevirtual

Reimplemented from ProcessLib::Process.

Definition at line 227 of file ComponentTransportProcess.cpp.

231{
232 // todo (renchao): move chemical calculation to elsewhere.
233 if (_process_data.lookup_table && process_id == 0)
234 {
235 INFO("Update process solutions via the look-up table approach");
236 _process_data.lookup_table->lookup(x, x_prev, _mesh.getNumberOfNodes());
237
238 return;
239 }
240
242 {
243 return;
244 }
245
246 // Sequential non-iterative approach applied here to split the reactive
247 // transport process into the transport stage followed by the reaction
248 // stage.
249 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
250 dof_tables.reserve(x.size());
251 std::generate_n(std::back_inserter(dof_tables), x.size(),
252 [&]() { return _local_to_global_index_map.get(); });
253
254 if (process_id == 0)
255 {
259 dof_tables, x, t, dt);
260
261 BaseLib::RunTime time_phreeqc;
262 time_phreeqc.start();
263
264 _chemical_solver_interface->setAqueousSolutionsPrevFromDumpFile();
265
266 _chemical_solver_interface->executeSpeciationCalculation(dt);
267
268 INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
269
271 &ComponentTransportLocalAssemblerInterface::
272 postSpeciationCalculation,
274 dt);
275
276 return;
277 }
278
279 auto const matrix_specification = getMatrixSpecifications(process_id);
280
281 std::size_t matrix_id = 0u;
283 matrix_specification, matrix_id);
285 matrix_specification, matrix_id);
286 auto& b =
288
289 M.setZero();
290 K.setZero();
291 b.setZero();
292
293 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
296 _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt, M, K,
297 b, process_id);
298
299 // todo (renchao): incorporate Neumann boundary condition
303
305 matrix_specification, matrix_id);
306 auto& rhs =
308
309 A.setZero();
310 rhs.setZero();
311
314
315 // compute A
317 MathLib::LinAlg::aypx(A, 1.0 / dt, K);
318
319 // compute rhs
320 MathLib::LinAlg::matMult(M, *x[process_id], rhs);
321 MathLib::LinAlg::aypx(rhs, 1.0 / dt, b);
322
323 using Tag = NumLib::NonlinearSolverTag;
325 auto& equation_system = static_cast<EqSys&>(ode_sys);
326 equation_system.applyKnownSolutionsPicard(A, rhs, *x[process_id]);
327
328 auto& linear_solver =
331 linear_solver.solve(A, rhs, *x[process_id]);
332
338}
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:34
Count the running time.
Definition: RunTime.h:29
double elapsed() const
Get the elapsed time in seconds.
Definition: RunTime.h:42
void start()
Start the timer.
Definition: RunTime.h:32
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:94
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
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)
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
Definition: Process.cpp:185
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
Definition: Types.h:20
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
void setLocalAccessibleVector(PETScVector const &x)
Definition: LinAlg.cpp:27
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:142
void aypx(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition: LinAlg.cpp:50
void finalizeVectorAssembly(VEC_T &)
General function to finalize the vector assembly.
static const double u
bool finalizeMatrixAssembly(MAT_T &)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
ChemistryLib::ChemicalSolverInterface *const chemical_solver_interface

## ◆ _chemical_solver_interface

 std::unique_ptr ProcessLib::ComponentTransport::ComponentTransportProcess::_chemical_solver_interface
private

Definition at line 171 of file ComponentTransportProcess.h.

## ◆ _local_assemblers

 std::vector > ProcessLib::ComponentTransport::ComponentTransportProcess::_local_assemblers
private

## ◆ _molar_flow_rate

 MeshLib::PropertyVector* ProcessLib::ComponentTransport::ComponentTransportProcess::_molar_flow_rate = nullptr
private

## ◆ _process_data

 ComponentTransportProcessData ProcessLib::ComponentTransport::ComponentTransportProcess::_process_data
private

Definition at line 163 of file ComponentTransportProcess.h.

## ◆ _surfaceflux

 std::unique_ptr ProcessLib::ComponentTransport::ComponentTransportProcess::_surfaceflux
private

Definition at line 168 of file ComponentTransportProcess.h.

Referenced by postTimestepConcreteProcess().

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