OGS
ProcessLib::ComponentTransport::ComponentTransportProcess Class Referencefinal

Detailed Description

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 95 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, bool const is_linear, bool const ls_compute_only_upon_timestep_change)
 
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, std::vector< GlobalVector * > const &x_prev, const double t, const double dt, int const process_id) override
 
bool shouldLinearSolverComputeOnlyUponTimestepChange () const 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.
 
void postTimestep (std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, const double delta_t, int const process_id)
 Postprocessing after a complete timestep.
 
void postNonLinearSolver (std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, 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_prev, int const process_id)
 compute secondary variables for the coupled equations or for output.
 
NumLib::IterationResult postIteration (GlobalVector const &x) final
 
void initialize (std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
 
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 updateDeactivatedSubdomains (double const time, const int process_id)
 
virtual bool isMonolithicSchemeUsed () 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 &x_prev, 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 &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) final
 
void preOutput (const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id)
 
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
 
std::vector< std::size_t > const & getActiveElementIDs () const
 
SecondaryVariableCollection const & getSecondaryVariables () const
 
std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > const & getIntegrationPointWriters () const
 
- Public Member Functions inherited from ProcessLib::SubmeshAssemblySupport
virtual std::vector< std::string > initializeAssemblyOnSubmeshes (std::vector< std::reference_wrapper< MeshLib::Mesh > > const &meshes)
 
virtual ~SubmeshAssemblySupport ()=default
 

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().
 
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 &x_prev, 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 &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
 
void preOutputConcreteProcess (const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id) 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
 
std::vector< MeshLib::PropertyVector< double > * > _residua
 
AssembledMatrixCache _asm_mat_cache
 
bool const _ls_compute_only_upon_timestep_change
 

Additional Inherited Members

- Public Attributes inherited from ProcessLib::Process
std::string const name
 
- Static Public Attributes inherited from ProcessLib::Process
static PROCESSLIB_EXPORT const std::string constant_one_parameter_name = "constant_one"
 
- 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, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
 
virtual void constructDofTable ()
 
void constructMonolithicProcessDofTable ()
 
void constructDofTableOfSpecifiedProcessStaggeredScheme (const int specified_process_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
 
std::unique_ptr< ProcessLib::AbstractJacobianAssembler_jacobian_assembler
 
VectorMatrixAssembler _global_assembler
 
const bool _use_monolithic_scheme
 
unsigned const _integration_order
 
std::vector< std::unique_ptr< MeshLib::IntegrationPointWriter > > _integration_point_writer
 
GlobalSparsityPattern _sparsity_pattern
 
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > _process_variables
 
std::vector< BoundaryConditionCollection_boundary_conditions
 

Constructor & Destructor Documentation

◆ 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,
bool const is_linear,
bool const ls_compute_only_upon_timestep_change )

Definition at line 39 of file ComponentTransportProcess.cpp.

55 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
56 integration_order, std::move(process_variables),
57 std::move(secondary_variables), use_monolithic_scheme),
58 _process_data(std::move(process_data)),
59 _surfaceflux(std::move(surfaceflux)),
60 _chemical_solver_interface(std::move(chemical_solver_interface)),
61 _asm_mat_cache{is_linear, use_monolithic_scheme},
63 ls_compute_only_upon_timestep_change}
64{
65 if (ls_compute_only_upon_timestep_change)
66 {
67 // TODO move this feature to some common location for all processes.
68 if (!is_linear)
69 {
71 "Using the linear solver compute() method only upon timestep "
72 "change only makes sense for linear model equations.");
73 }
74
75 WARN(
76 "You specified that the ComponentTransport linear solver will do "
77 "the compute() step only upon timestep change. This is an expert "
78 "option. It is your responsibility to ensure that "
79 "the conditions for the correct use of this feature are met! "
80 "Otherwise OGS might compute garbage without being recognized. "
81 "There is no safety net!");
82
83 WARN(
84 "You specified that the ComponentTransport linear solver will do "
85 "the compute() step only upon timestep change. This option will "
86 "only be used by the Picard non-linear solver. The Newton-Raphson "
87 "non-linear solver will silently ignore this setting, i.e., it "
88 "won't do any harm, there, but you won't observe the speedup you "
89 "probably expect.");
90 }
91
92 auto residuum_name = [&](auto const& pv) -> std::string
93 {
94 std::string const& pv_name = pv.getName();
95 if (pv_name == "pressure")
96 {
97 return "LiquidMassFlowRate";
98 }
99 if (pv_name == "temperature")
100 {
101 return "HeatFlowRate";
102 }
103 return pv_name + "FlowRate";
104 };
105
107 {
108 int const process_id = 0;
109 for (auto const& pv : _process_variables[process_id])
110 {
111 _residua.push_back(MeshLib::getOrCreateMeshProperty<double>(
112 mesh, residuum_name(pv.get()), MeshLib::MeshItemType::Node, 1));
113 }
114 }
115 else
116 {
117 for (auto const& pv : _process_variables)
118 {
119 _residua.push_back(MeshLib::getOrCreateMeshProperty<double>(
120 mesh, residuum_name(pv[0].get()), MeshLib::MeshItemType::Node,
121 1));
122 }
123 }
124}
#define OGS_FATAL(...)
Definition Error.h:26
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
std::unique_ptr< ChemistryLib::ChemicalSolverInterface > _chemical_solver_interface
std::vector< MeshLib::PropertyVector< double > * > _residua
std::unique_ptr< ProcessLib::SurfaceFluxData > _surfaceflux
std::string const name
Definition Process.h:351
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:28
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > _process_variables
Definition Process.h:387
const bool _use_monolithic_scheme
Definition Process.h:366
auto & get(Tuples &... ts)
Definition Get.h:67

References ProcessLib::Process::_process_variables, _residua, ProcessLib::Process::_use_monolithic_scheme, MeshLib::Node, OGS_FATAL, and WARN().

Member Function Documentation

◆ assembleConcreteProcess()

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

Implements ProcessLib::Process.

Definition at line 221 of file ComponentTransportProcess.cpp.

225{
226 DBUG("Assemble ComponentTransportProcess.");
227
228 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
229
230 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
231 dof_tables;
233 {
234 dof_tables.push_back(std::ref(*_local_to_global_index_map));
235 }
236 else
237 {
238 std::generate_n(
239 std::back_inserter(dof_tables), _process_variables.size(),
240 [&]() { return std::ref(*_local_to_global_index_map); });
241 }
242
243 _asm_mat_cache.assemble(t, dt, x, x_prev, process_id, M, K, b, dof_tables,
246
247 // TODO (naumov) What about temperature? A test with FCT would reveal any
248 // problems.
249 if (process_id == _process_data.hydraulic_process_id)
250 {
251 return;
252 }
253 auto const matrix_specification = getMatrixSpecifications(process_id);
254
256 x_prev, process_id,
257 matrix_specification, M, K, b);
258}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
std::vector< std::unique_ptr< ComponentTransportLocalAssemblerInterface > > _local_assemblers
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition Process.h:155
VectorMatrixAssembler _global_assembler
Definition Process.h:364
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition Process.h:357
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
Definition Process.cpp:194
void computeFluxCorrectedTransport(NumericalStabilization const &stabilizer, const double t, const double dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, const MathLib::MatrixSpecifications &matrix_specification, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
void assemble(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > const &dof_tables, VectorMatrixAssembler &global_assembler, VectorOfLocalAssemblers const &local_assemblers, std::vector< std::size_t > const &active_element_ids)

References _asm_mat_cache, ProcessLib::Process::_global_assembler, _local_assemblers, ProcessLib::Process::_local_to_global_index_map, _process_data, ProcessLib::Process::_process_variables, ProcessLib::Process::_use_monolithic_scheme, ProcessLib::AssembledMatrixCache::assemble(), NumLib::computeFluxCorrectedTransport(), DBUG(), ProcessLib::ProcessVariable::getActiveElementIDs(), ProcessLib::Process::getMatrixSpecifications(), ProcessLib::Process::getProcessVariables(), ProcessLib::ComponentTransport::ComponentTransportProcessData::hydraulic_process_id, and ProcessLib::ComponentTransport::ComponentTransportProcessData::stabilizer.

Referenced by preOutputConcreteProcess().

◆ assembleWithJacobianConcreteProcess()

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

Implements ProcessLib::Process.

Definition at line 260 of file ComponentTransportProcess.cpp.

264{
265 DBUG("AssembleWithJacobian ComponentTransportProcess.");
266
268 {
269 OGS_FATAL(
270 "The AssembleWithJacobian() for ComponentTransportProcess is "
271 "implemented only for staggered scheme.");
272 }
273
274 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
275
276 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
277 dof_tables;
278
279 std::generate_n(std::back_inserter(dof_tables), _process_variables.size(),
280 [&]() { return std::ref(*_local_to_global_index_map); });
281
282 // Call global assembler for each local assembly item.
285 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x,
286 x_prev, process_id, M, K, b, Jac);
287
288 // b is the negated residumm used in the Newton's method.
289 // Here negating b is to recover the primitive residuum.
291 *_residua[process_id],
292 std::negate<double>());
293}
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 &x_prev, 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 map_function)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)

References ProcessLib::Process::_global_assembler, _local_assemblers, ProcessLib::Process::_local_to_global_index_map, ProcessLib::Process::_process_variables, _residua, ProcessLib::Process::_use_monolithic_scheme, ProcessLib::VectorMatrixAssembler::assembleWithJacobian(), DBUG(), NumLib::SerialExecutor::executeSelectedMemberDereferenced(), ProcessLib::ProcessVariable::getActiveElementIDs(), ProcessLib::Process::getProcessVariables(), and OGS_FATAL.

◆ computeSecondaryVariableConcrete()

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

Reimplemented from ProcessLib::Process.

Definition at line 426 of file ComponentTransportProcess.cpp.

432{
433 if (process_id != 0)
434 {
435 return;
436 }
437
438 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
439 dof_tables.reserve(x.size());
440 std::generate_n(std::back_inserter(dof_tables), x.size(),
441 [&]() { return _local_to_global_index_map.get(); });
442
443 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
446 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x,
447 x_prev, process_id);
448
450 {
451 return;
452 }
453
455 &ComponentTransportLocalAssemblerInterface::
456 computeReactionRelatedSecondaryVariable,
458}
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)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)

References _chemical_solver_interface, _local_assemblers, ProcessLib::LocalAssemblerInterface::computeSecondaryVariable(), NumLib::SerialExecutor::executeSelectedMemberOnDereferenced(), ProcessLib::ProcessVariable::getActiveElementIDs(), and ProcessLib::Process::getProcessVariables().

◆ 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 295 of file ComponentTransportProcess.cpp.

300{
301 std::vector<GlobalIndexType> indices_cache;
302 auto const r_c_indices = NumLib::getRowColumnIndices(
303 element_id, *_local_to_global_index_map, indices_cache);
304
305 std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
306 x.size(), r_c_indices.rows};
307 auto const local_xs =
308 getCoupledLocalSolutions(x, indices_of_all_coupled_processes);
309
310 return _local_assemblers[element_id]->getFlux(p, t, local_xs);
311}
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)

References _local_assemblers, ProcessLib::Process::_local_to_global_index_map, ProcessLib::getCoupledLocalSolutions(), and NumLib::getRowColumnIndices().

◆ 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 126 of file ComponentTransportProcess.cpp.

130{
131 int const mesh_space_dimension = _process_data.mesh_space_dimension;
132
133 _process_data.mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
134 const_cast<MeshLib::Mesh&>(mesh), "velocity",
135 MeshLib::MeshItemType::Cell, mesh_space_dimension);
136
137 _process_data.mesh_prop_porosity = MeshLib::getOrCreateMeshProperty<double>(
138 const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
140
141 std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
142 transport_process_variables;
144 {
145 const int process_id = 0;
146 for (auto pv_iter = std::next(_process_variables[process_id].begin());
147 pv_iter != _process_variables[process_id].end();
148 ++pv_iter)
149 {
150 transport_process_variables.push_back(*pv_iter);
151 }
152 }
153 else
154 {
155 // All process variables but the pressure and optionally the temperature
156 // are transport variables.
157 for (auto const& pv :
159 ranges::views::drop(_process_data.isothermal ? 1 : 2))
160 {
161 transport_process_variables.push_back(pv[0]);
162 }
163 }
164
165 ProcessLib::createLocalAssemblers<LocalAssemblerData>(
166 mesh_space_dimension, mesh.getElements(), dof_table, _local_assemblers,
167 NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
168 _process_data, transport_process_variables);
169
171 {
175
176 _chemical_solver_interface->initialize();
177 }
178
180 "darcy_velocity",
182 mesh_space_dimension, getExtrapolator(), _local_assemblers,
184
186 "molar_flux",
188 mesh_space_dimension, getExtrapolator(), _local_assemblers,
190}
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:359
NumLib::Extrapolator & getExtrapolator() const
Definition Process.h:196
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)

References _chemical_solver_interface, _local_assemblers, _process_data, ProcessLib::Process::_process_variables, ProcessLib::Process::_secondary_variables, ProcessLib::Process::_use_monolithic_scheme, ProcessLib::SecondaryVariableCollection::addSecondaryVariable(), MeshLib::Cell, NumLib::SerialExecutor::executeSelectedMemberOnDereferenced(), MeshLib::Mesh::getElements(), ProcessLib::Process::getExtrapolator(), ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::getIntPtDarcyVelocity(), ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::getIntPtMolarFlux(), MeshLib::Mesh::isAxiallySymmetric(), ProcessLib::ComponentTransport::ComponentTransportProcessData::isothermal, ProcessLib::makeExtrapolator(), ProcessLib::ComponentTransport::ComponentTransportProcessData::mesh_prop_porosity, ProcessLib::ComponentTransport::ComponentTransportProcessData::mesh_prop_velocity, ProcessLib::ComponentTransport::ComponentTransportProcessData::mesh_space_dimension, and ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::setChemicalSystemID().

◆ isLinear()

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

◆ postTimestepConcreteProcess()

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

Reimplemented from ProcessLib::Process.

Definition at line 460 of file ComponentTransportProcess.cpp.

466{
467 if (process_id != 0)
468 {
469 return;
470 }
471
472 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
473 dof_tables.reserve(x.size());
474 std::generate_n(std::back_inserter(dof_tables), x.size(),
475 [&]() { return _local_to_global_index_map.get(); });
476
477 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
480 _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, x_prev, t,
481 dt, process_id);
482
483 if (!_surfaceflux) // computing the surfaceflux is optional
484 {
485 return;
486 }
487 _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
489}
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)
MeshLib::Mesh & _mesh
Definition Process.h:354
unsigned const _integration_order
Definition Process.h:371

References ProcessLib::Process::_integration_order, _local_assemblers, ProcessLib::Process::_mesh, _surfaceflux, NumLib::SerialExecutor::executeSelectedMemberOnDereferenced(), ProcessLib::ProcessVariable::getActiveElementIDs(), ProcessLib::Process::getProcessVariables(), and ProcessLib::LocalAssemblerInterface::postTimestep().

◆ preOutputConcreteProcess()

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

Reimplemented from ProcessLib::Process.

Definition at line 491 of file ComponentTransportProcess.cpp.

497{
498 auto const matrix_specification = getMatrixSpecifications(process_id);
499
501 matrix_specification);
503 matrix_specification);
505 matrix_specification);
506
507 M->setZero();
508 K->setZero();
509 b->setZero();
510
511 assembleConcreteProcess(t, dt, x, x_prev, process_id, *M, *K, *b);
512
513 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
514 dof_tables;
516 {
517 dof_tables.push_back(std::ref(*_local_to_global_index_map));
518 }
519 else
520 {
521 std::generate_n(
522 std::back_inserter(dof_tables), _process_variables.size(),
523 [&]() { return std::ref(*_local_to_global_index_map); });
524 }
525
526 BaseLib::RunTime time_residuum;
527 time_residuum.start();
528
530 {
531 auto const residuum =
532 computeResiduum(dt, *x[0], *x_prev[0], *M, *K, *b);
533 for (std::size_t variable_id = 0; variable_id < _residua.size();
534 ++variable_id)
535 {
537 residuum, variable_id, dof_tables[0], *_residua[variable_id],
538 std::negate<double>());
539 }
540 }
541 else
542 {
543 auto const residuum = computeResiduum(dt, *x[process_id],
544 *x_prev[process_id], *M, *K, *b);
545 transformVariableFromGlobalVector(residuum, 0, dof_tables[process_id],
546 *_residua[process_id],
547 std::negate<double>());
548 }
549
550 INFO("[time] Computing residuum flow rates took {:g} s",
551 time_residuum.elapsed());
552}
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
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
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
GlobalVector computeResiduum(double const dt, GlobalVector const &x, GlobalVector const &x_prev, GlobalMatrix const &M, GlobalMatrix const &K, GlobalVector const &b)

References ProcessLib::Process::_local_to_global_index_map, ProcessLib::Process::_process_variables, _residua, ProcessLib::Process::_use_monolithic_scheme, assembleConcreteProcess(), ProcessLib::computeResiduum(), BaseLib::RunTime::elapsed(), ProcessLib::Process::getMatrixSpecifications(), INFO(), and BaseLib::RunTime::start().

◆ 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 192 of file ComponentTransportProcess.cpp.

194{
196 {
197 return;
198 }
199
200 if (process_id != static_cast<int>(x.size() - 1))
201 {
202 return;
203 }
204
205 std::for_each(
206 x.begin(), x.end(),
207 [](auto const process_solution)
208 { MathLib::LinAlg::setLocalAccessibleVector(*process_solution); });
209
210 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
211 dof_tables.reserve(x.size());
212 std::generate_n(std::back_inserter(dof_tables), x.size(),
213 [&]() { return _local_to_global_index_map.get(); });
214
218 dof_tables, x, t);
219}
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)

References _chemical_solver_interface, _local_assemblers, NumLib::SerialExecutor::executeSelectedMemberOnDereferenced(), and ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::initializeChemicalSystem().

◆ shouldLinearSolverComputeOnlyUponTimestepChange()

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

◆ 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 313 of file ComponentTransportProcess.cpp.

317{
318 // todo (renchao): move chemical calculation to elsewhere.
319 if (_process_data.lookup_table && process_id == 0)
320 {
321 INFO("Update process solutions via the look-up table approach");
322 _process_data.lookup_table->lookup(x, x_prev, _mesh.getNumberOfNodes());
323
324 return;
325 }
326
328 {
329 return;
330 }
331
332 // Sequential non-iterative approach applied here to split the reactive
333 // transport process into the transport stage followed by the reaction
334 // stage.
335 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
336 dof_tables.reserve(x.size());
337 std::generate_n(std::back_inserter(dof_tables), x.size(),
338 [&]() { return _local_to_global_index_map.get(); });
339
340 if (process_id == 0)
341 {
345 dof_tables, x, t, dt);
346
347 BaseLib::RunTime time_phreeqc;
348 time_phreeqc.start();
349
350 _chemical_solver_interface->setAqueousSolutionsPrevFromDumpFile();
351
352 _chemical_solver_interface->executeSpeciationCalculation(dt);
353
354 INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
355
357 &ComponentTransportLocalAssemblerInterface::
358 postSpeciationCalculation,
360 dt);
361
362 return;
363 }
364
365 auto const matrix_specification = getMatrixSpecifications(process_id);
366
367 std::size_t matrix_id = 0u;
369 matrix_specification, matrix_id);
371 matrix_specification, matrix_id);
372 auto& b =
374
375 M.setZero();
376 K.setZero();
377 b.setZero();
378
379 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
382 _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt, M, K,
383 b, process_id);
384
385 // todo (renchao): incorporate Neumann boundary condition
389
391 matrix_specification, matrix_id);
392 auto& rhs =
394
395 A.setZero();
396 rhs.setZero();
397
400
401 // compute A
403 MathLib::LinAlg::aypx(A, 1.0 / dt, K);
404
405 // compute rhs
406 MathLib::LinAlg::matMult(M, *x[process_id], rhs);
407 MathLib::LinAlg::aypx(rhs, 1.0 / dt, b);
408
409 using Tag = NumLib::NonlinearSolverTag;
411 auto& equation_system = static_cast<EqSys&>(ode_sys);
412 equation_system.applyKnownSolutionsPicard(A, rhs, *x[process_id]);
413
414 auto& linear_solver =
417 linear_solver.solve(A, rhs, *x[process_id]);
418
424}
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition Mesh.h:100
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)
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:149
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.
bool finalizeMatrixAssembly(MAT_T &)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider

References _chemical_solver_interface, _local_assemblers, ProcessLib::Process::_mesh, _process_data, ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::assembleReactionEquation(), MathLib::LinAlg::aypx(), ProcessLib::ComponentTransport::ComponentTransportProcessData::chemical_solver_interface, MathLib::LinAlg::copy(), BaseLib::RunTime::elapsed(), NumLib::SerialExecutor::executeSelectedMemberOnDereferenced(), MathLib::finalizeMatrixAssembly(), MathLib::finalizeVectorAssembly(), ProcessLib::ProcessVariable::getActiveElementIDs(), NumLib::MatrixProvider::getMatrix(), ProcessLib::Process::getMatrixSpecifications(), MeshLib::Mesh::getNumberOfNodes(), ProcessLib::Process::getProcessVariables(), NumLib::VectorProvider::getVector(), INFO(), ChemistryLib::ChemicalSolverInterface::linear_solver, ProcessLib::ComponentTransport::ComponentTransportProcessData::lookup_table, MathLib::LinAlg::matMult(), NumLib::GlobalVectorProvider::provider, NumLib::GlobalMatrixProvider::provider, NumLib::MatrixProvider::releaseMatrix(), NumLib::VectorProvider::releaseVector(), ProcessLib::ComponentTransport::ComponentTransportLocalAssemblerInterface::setChemicalSystem(), MathLib::LinAlg::setLocalAccessibleVector(), MathLib::EigenVector::setZero(), and BaseLib::RunTime::start().

Member Data Documentation

◆ _asm_mat_cache

AssembledMatrixCache ProcessLib::ComponentTransport::ComponentTransportProcess::_asm_mat_cache
private

Definition at line 189 of file ComponentTransportProcess.h.

Referenced by assembleConcreteProcess(), and isLinear().

◆ _chemical_solver_interface

std::unique_ptr<ChemistryLib::ChemicalSolverInterface> ProcessLib::ComponentTransport::ComponentTransportProcess::_chemical_solver_interface
private

◆ _local_assemblers

◆ _ls_compute_only_upon_timestep_change

bool const ProcessLib::ComponentTransport::ComponentTransportProcess::_ls_compute_only_upon_timestep_change
private

◆ _process_data

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

◆ _residua

std::vector<MeshLib::PropertyVector<double>*> ProcessLib::ComponentTransport::ComponentTransportProcess::_residua
private

◆ _surfaceflux

std::unique_ptr<ProcessLib::SurfaceFluxData> ProcessLib::ComponentTransport::ComponentTransportProcess::_surfaceflux
private

Definition at line 182 of file ComponentTransportProcess.h.

Referenced by postTimestepConcreteProcess().


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