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>(
186 for (std::size_t component_id = 0;
187 component_id < transport_process_variables.size();
190 auto const& pv = transport_process_variables[component_id].get();
192 auto integration_point_values_accessor = [component_id](
195 std::vector<GlobalVector*>
const& x,
196 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
197 std::vector<double>& cache) ->
auto const&
204 pv.getName() +
"Flux",
207 integration_point_values_accessor));
212 std::vector<GlobalVector*>& x,
double const t,
int const process_id)
219 if (process_id !=
static_cast<int>(x.size() - 1))
226 [](
auto const process_solution)
227 { MathLib::LinAlg::setLocalAccessibleVector(*process_solution); });
229 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
230 dof_tables.reserve(x.size());
231 std::generate_n(std::back_inserter(dof_tables), x.size(),
232 [&]() { return _local_to_global_index_map.get(); });
241 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
242 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
245 DBUG(
"Assemble ComponentTransportProcess.");
247 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
254 std::generate_n(std::back_inserter(dof_tables),
256 [&]() { return _local_to_global_index_map.get(); });
273 matrix_specification, M, K, b);
277 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
278 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
281 DBUG(
"AssembleWithJacobian ComponentTransportProcess.");
286 "The AssembleWithJacobian() for ComponentTransportProcess is "
287 "implemented only for staggered scheme.");
290 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
293 [&]() { return _local_to_global_index_map.get(); });
299 process_id, &b, &Jac);
305 std::negate<double>());
309 std::size_t
const element_id,
312 std::vector<GlobalVector*>
const& x)
const
314 std::vector<GlobalIndexType> indices_cache;
318 std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
319 x.size(), r_c_indices.rows};
320 auto const local_xs =
327 std::vector<GlobalVector*>& x, std::vector<GlobalVector*>
const& x_prev,
329 int const process_id)
334 INFO(
"Update process solutions via the look-up table approach");
348 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
349 dof_tables.reserve(x.size());
350 std::generate_n(std::back_inserter(dof_tables), x.size(),
351 [&]() { return _local_to_global_index_map.get(); });
358 dof_tables, x, t, dt);
361 time_phreeqc.
start();
367 INFO(
"[time] Phreeqc took {:g} s.", time_phreeqc.
elapsed());
371 postSpeciationCalculation,
380 std::size_t matrix_id = 0u;
382 matrix_specification, matrix_id);
384 matrix_specification, matrix_id);
403 matrix_specification, matrix_id);
423 auto& equation_system =
static_cast<EqSys&
>(ode_sys);
424 equation_system.applyKnownSolutionsPicard(A, rhs, *x[process_id]);
426 auto& linear_solver =
429 linear_solver.solve(A, rhs, *x[process_id]);
441 std::vector<GlobalVector*>
const& x,
443 int const process_id)
450 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
451 dof_tables.reserve(x.size());
452 std::generate_n(std::back_inserter(dof_tables), x.size(),
453 [&]() { return _local_to_global_index_map.get(); });
467 computeReactionRelatedSecondaryVariable,
472 std::vector<GlobalVector*>
const& x,
473 std::vector<GlobalVector*>
const& x_prev,
476 int const process_id)
483 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
484 dof_tables.reserve(x.size());
485 std::generate_n(std::back_inserter(dof_tables), x.size(),
486 [&]() { return _local_to_global_index_map.get(); });
504 std::vector<GlobalVector*>
const& x,
505 std::vector<GlobalVector*>
const& x_prev,
506 int const process_id)
511 matrix_specification);
513 matrix_specification);
515 matrix_specification);
523 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
530 std::generate_n(std::back_inserter(dof_tables),
532 [&]() { return _local_to_global_index_map.get(); });
536 time_residuum.
start();
540 auto const residuum =
542 for (std::size_t variable_id = 0; variable_id <
_residua.size();
545 transformVariableFromGlobalVector(
546 residuum, variable_id, *dof_tables[0], *
_residua[variable_id],
547 std::negate<double>());
553 *x_prev[process_id], *M, *K, *b);
554 transformVariableFromGlobalVector(residuum, 0, *dof_tables[process_id],
556 std::negate<double>());
559 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 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)
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, int const component_id) const =0
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
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 assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector &b, GlobalMatrix &Jac) 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::vector< std::reference_wrapper< ProcessVariable > > > _process_variables
std::vector< std::size_t > const & getActiveElementIDs() const
SecondaryVariableCollection _secondary_variables
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, 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 *const M, GlobalMatrix *const K, GlobalVector *const 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