14#include <range/v3/algorithm/copy.hpp>
15#include <range/v3/view/drop.hpp>
38namespace ComponentTransport
43 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
44 std::vector<std::unique_ptr<ParameterLib::ParameterBase>>
const& parameters,
45 unsigned const integration_order,
46 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
50 bool const use_monolithic_scheme,
51 std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux,
52 std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
53 chemical_solver_interface,
55 bool const ls_compute_only_upon_timestep_change)
56 :
Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
57 integration_order, std::move(process_variables),
58 std::move(secondary_variables), use_monolithic_scheme),
59 _process_data(std::move(process_data)),
60 _surfaceflux(std::move(surfaceflux)),
61 _chemical_solver_interface(std::move(chemical_solver_interface)),
62 _asm_mat_cache{is_linear, use_monolithic_scheme},
63 _ls_compute_only_upon_timestep_change{
64 ls_compute_only_upon_timestep_change}
66 if (ls_compute_only_upon_timestep_change)
72 "Using the linear solver compute() method only upon timestep "
73 "change only makes sense for linear model equations.");
77 "You specified that the ComponentTransport linear solver will do "
78 "the compute() step only upon timestep change. This is an expert "
79 "option. It is your responsibility to ensure that "
80 "the conditions for the correct use of this feature are met! "
81 "Otherwise OGS might compute garbage without being recognized. "
82 "There is no safety net!");
85 "You specified that the ComponentTransport linear solver will do "
86 "the compute() step only upon timestep change. This option will "
87 "only be used by the Picard non-linear solver. The Newton-Raphson "
88 "non-linear solver will silently ignore this setting, i.e., it "
89 "won't do any harm, there, but you won't observe the speedup you "
93 auto residuum_name = [&](
auto const& pv) -> std::string
95 std::string
const& pv_name = pv.getName();
96 if (pv_name ==
"pressure")
98 return "LiquidMassFlowRate";
100 if (pv_name ==
"temperature")
102 return "HeatFlowRate";
104 return pv_name +
"FlowRate";
109 int const process_id = 0;
112 _residua.push_back(MeshLib::getOrCreateMeshProperty<double>(
120 _residua.push_back(MeshLib::getOrCreateMeshProperty<double>(
130 unsigned const integration_order)
142 std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
143 transport_process_variables;
146 const int process_id = 0;
151 transport_process_variables.push_back(*pv_iter);
158 for (
auto const& pv :
162 transport_process_variables.push_back(pv[0]);
166 ProcessLib::createLocalAssemblers<LocalAssemblerData>(
194 std::vector<GlobalVector*>& x,
double const t,
int const process_id)
201 if (process_id !=
static_cast<int>(x.size() - 1))
208 [](
auto const process_solution)
209 { MathLib::LinAlg::setLocalAccessibleVector(*process_solution); });
211 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
212 dof_tables.reserve(x.size());
213 std::generate_n(std::back_inserter(dof_tables), x.size(),
214 [&]() { return _local_to_global_index_map.get(); });
223 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
224 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
227 DBUG(
"Assemble ComponentTransportProcess.");
231 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
238 std::generate_n(std::back_inserter(dof_tables),
240 [&]() { return _local_to_global_index_map.get(); });
257 matrix_specification, M, K, b);
261 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
262 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
265 DBUG(
"AssembleWithJacobian ComponentTransportProcess.");
270 "The AssembleWithJacobian() for ComponentTransportProcess is "
271 "implemented only for staggered scheme.");
276 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
279 [&]() { return _local_to_global_index_map.get(); });
285 x_prev, process_id, M, K, b, Jac);
291 std::negate<double>());
295 std::size_t
const element_id,
298 std::vector<GlobalVector*>
const& x)
const
300 std::vector<GlobalIndexType> indices_cache;
304 std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
305 x.size(), r_c_indices.rows};
306 auto const local_xs =
313 std::vector<GlobalVector*>& x, std::vector<GlobalVector*>
const& x_prev,
315 int const process_id)
320 INFO(
"Update process solutions via the look-up table approach");
334 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
335 dof_tables.reserve(x.size());
336 std::generate_n(std::back_inserter(dof_tables), x.size(),
337 [&]() { return _local_to_global_index_map.get(); });
344 dof_tables, x, t, dt);
347 time_phreeqc.
start();
353 INFO(
"[time] Phreeqc took {:g} s.", time_phreeqc.
elapsed());
357 postSpeciationCalculation,
366 std::size_t matrix_id = 0u;
368 matrix_specification, matrix_id);
370 matrix_specification, matrix_id);
390 matrix_specification, matrix_id);
410 auto& equation_system =
static_cast<EqSys&
>(ode_sys);
411 equation_system.applyKnownSolutionsPicard(A, rhs, *x[process_id]);
413 auto& linear_solver =
416 linear_solver.solve(A, rhs, *x[process_id]);
428 std::vector<GlobalVector*>
const& x,
430 int const process_id)
437 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
438 dof_tables.reserve(x.size());
439 std::generate_n(std::back_inserter(dof_tables), x.size(),
440 [&]() { return _local_to_global_index_map.get(); });
455 computeReactionRelatedSecondaryVariable,
460 std::vector<GlobalVector*>
const& x,
461 std::vector<GlobalVector*>
const& x_prev,
464 int const process_id)
471 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
472 dof_tables.reserve(x.size());
473 std::generate_n(std::back_inserter(dof_tables), x.size(),
474 [&]() { return _local_to_global_index_map.get(); });
493 std::vector<GlobalVector*>
const& x,
494 std::vector<GlobalVector*>
const& x_prev,
495 int const process_id)
500 matrix_specification);
502 matrix_specification);
504 matrix_specification);
512 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
519 std::generate_n(std::back_inserter(dof_tables),
521 [&]() { return _local_to_global_index_map.get(); });
525 time_residuum.
start();
529 auto const residuum =
531 for (std::size_t variable_id = 0; variable_id <
_residua.size();
534 transformVariableFromGlobalVector(
535 residuum, variable_id, *dof_tables[0], *
_residua[variable_id],
536 std::negate<double>());
542 *x_prev[process_id], *M, *K, *b);
543 transformVariableFromGlobalVector(residuum, 0, *dof_tables[process_id],
545 std::negate<double>());
548 INFO(
"[time] Computing residuum flow rates took {:g} s",
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the RunTime class.
double elapsed() const
Get the elapsed time in seconds.
void start()
Start the timer.
GlobalLinearSolver & linear_solver
Global vector based on Eigen vector.
bool isAxiallySymmetric() const
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
std::size_t getNumberOfNodes() const
Get the number of nodes.
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
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
virtual void setChemicalSystemID(std::size_t const)=0
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)
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) 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
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 initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
std::unique_ptr< ChemistryLib::ChemicalSolverInterface > _chemical_solver_interface
std::vector< MeshLib::PropertyVector< double > * > _residua
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
ComponentTransportProcessData _process_data
void computeSecondaryVariableConcrete(double const, double const, std::vector< GlobalVector * > const &x, GlobalVector const &, int 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 postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, const double dt, int const process_id) override
std::vector< std::unique_ptr< ComponentTransportLocalAssemblerInterface > > _local_assemblers
std::unique_ptr< ProcessLib::SurfaceFluxData > _surfaceflux
ComponentTransportProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const ¶meters, 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
AssembledMatrixCache _asm_mat_cache
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)
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)
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > _process_variables
SecondaryVariableCollection _secondary_variables
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
VectorMatrixAssembler _global_assembler
unsigned const _integration_order
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
NumLib::Extrapolator & getExtrapolator() const
const bool _use_monolithic_scheme
Handles configuration of several secondary variables from the project file.
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > 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)
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
void copy(PETScVector const &x, PETScVector &y)
void setLocalAccessibleVector(PETScVector const &x)
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
void aypx(PETScVector &y, PetscScalar const a, PETScVector const &x)
void finalizeVectorAssembly(VEC_T &)
General function to finalize the vector assembly.
bool finalizeMatrixAssembly(MAT_T &)
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)
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalIndexType > &indices)
GlobalVector computeResiduum(double const dt, GlobalVector const &x, GlobalVector const &x_prev, GlobalMatrix const &M, GlobalMatrix const &K, GlobalVector const &b)
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType > > const &indices)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
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< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, VectorMatrixAssembler &global_assembler, VectorOfLocalAssemblers const &local_assemblers, std::vector< std::size_t > const &active_element_ids)
NumLib::NumericalStabilization stabilizer
std::unique_ptr< LookupTable > lookup_table
ChemistryLib::ChemicalSolverInterface *const chemical_solver_interface
static const int hydraulic_process_id
MeshLib::PropertyVector< double > * mesh_prop_velocity
int const mesh_space_dimension
MeshLib::PropertyVector< double > * mesh_prop_porosity