OGS
ComponentTransportProcess.cpp
Go to the documentation of this file.
1 
12 
13 #include <cassert>
14 
15 #include "BaseLib/RunTime.h"
20 #include "MathLib/LinAlg/LinAlg.h"
25 
26 namespace ProcessLib
27 {
28 namespace ComponentTransport
29 {
31  std::string name,
32  MeshLib::Mesh& mesh,
33  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
34  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
35  unsigned const integration_order,
36  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
37  process_variables,
38  ComponentTransportProcessData&& process_data,
39  SecondaryVariableCollection&& secondary_variables,
40  bool const use_monolithic_scheme,
41  std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux,
42  std::unique_ptr<ChemistryLib::ChemicalSolverInterface>&&
43  chemical_solver_interface)
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 }
52 
54  NumLib::LocalToGlobalIndexMap const& dof_table,
55  MeshLib::Mesh const& mesh,
56  unsigned const integration_order)
57 {
58  _process_data.mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
59  const_cast<MeshLib::Mesh&>(mesh), "velocity",
61 
62  _process_data.mesh_prop_porosity = MeshLib::getOrCreateMeshProperty<double>(
63  const_cast<MeshLib::Mesh&>(mesh), "porosity_avg",
65 
66  const int process_id = 0;
67  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
68 
69  std::vector<std::reference_wrapper<ProcessLib::ProcessVariable>>
70  transport_process_variables;
72  {
73  for (auto pv_iter = std::next(_process_variables[process_id].begin());
74  pv_iter != _process_variables[process_id].end();
75  ++pv_iter)
76  {
77  transport_process_variables.push_back(*pv_iter);
78  }
79  }
80  else
81  {
82  for (auto pv_iter = std::next(_process_variables.begin());
83  pv_iter != _process_variables.end();
84  ++pv_iter)
85  {
86  transport_process_variables.push_back((*pv_iter)[0]);
87  }
88  }
89 
90  ProcessLib::createLocalAssemblers<LocalAssemblerData>(
91  mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
92  mesh.isAxiallySymmetric(), integration_order, _process_data,
93  transport_process_variables);
94 
96  {
100 
101  _chemical_solver_interface->initialize();
102  }
103 
105  "darcy_velocity",
109 }
110 
112  std::vector<GlobalVector*>& x, double const t, int const process_id)
113 {
115  {
116  return;
117  }
118 
119  if (process_id != static_cast<int>(x.size() - 1))
120  {
121  return;
122  }
123 
124  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
125 
126  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
127  dof_tables.reserve(x.size());
128  std::generate_n(std::back_inserter(dof_tables), x.size(),
129  [&]() { return _local_to_global_index_map.get(); });
130 
133  _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t);
134 
135  BaseLib::RunTime time_phreeqc;
136  time_phreeqc.start();
137 
138  _chemical_solver_interface->executeSpeciationCalculation(0. /*dt*/);
139 
141  t, _chemical_solver_interface->getIntPtProcessSolutions(), x);
142 
143  INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
144 }
145 
147  const double t, double const dt, std::vector<GlobalVector*> const& x,
148  std::vector<GlobalVector*> const& xdot, int const process_id,
150 {
151  DBUG("Assemble ComponentTransportProcess.");
152 
153  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
154 
155  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
156  dof_tables;
158  {
159  dof_tables.push_back(std::ref(*_local_to_global_index_map));
160  }
161  else
162  {
163  std::generate_n(
164  std::back_inserter(dof_tables), _process_variables.size(),
165  [&]() { return std::ref(*_local_to_global_index_map); });
166  }
167  // Call global assembler for each local assembly item.
170  pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
171  b);
172 }
173 
175  const double t, double const dt, std::vector<GlobalVector*> const& x,
176  std::vector<GlobalVector*> const& xdot, int const process_id,
178 {
179  DBUG("AssembleWithJacobian ComponentTransportProcess.");
180 
181  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
182  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
183  dof_table = {std::ref(*_local_to_global_index_map)};
184  // Call global assembler for each local assembly item.
187  _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
188  process_id, M, K, b, Jac);
189 }
190 
192  std::size_t const element_id,
193  MathLib::Point3d const& p,
194  double const t,
195  std::vector<GlobalVector*> const& x) const
196 {
197  std::vector<GlobalIndexType> indices_cache;
198  auto const r_c_indices = NumLib::getRowColumnIndices(
199  element_id, *_local_to_global_index_map, indices_cache);
200 
201  std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
202  x.size(), r_c_indices.rows};
203  auto const local_xs =
204  getCoupledLocalSolutions(x, indices_of_all_coupled_processes);
205 
206  return _local_assemblers[element_id]->getFlux(p, t, local_xs);
207 }
208 
211 {
212  DBUG("Set the coupled term for the staggered scheme to local assembers.");
213 
214  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
217  setStaggeredCoupledSolutions,
219 }
220 
222  std::vector<GlobalVector*>& x, std::vector<GlobalVector*> const& x_prev,
223  double const t, double const dt, NumLib::EquationSystem& ode_sys,
224  int const process_id)
225 {
226  // todo (renchao): move chemical calculation to elsewhere.
227  if (_process_data.lookup_table && process_id == 0)
228  {
229  INFO("Update process solutions via the look-up table approach");
230  _process_data.lookup_table->lookup(x, x_prev, _mesh.getNumberOfNodes());
231 
232  return;
233  }
234 
236  {
237  return;
238  }
239 
240  // Sequential non-iterative approach applied here to split the reactive
241  // transport process into the transport stage followed by the reaction
242  // stage.
243  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
244 
245  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
246  dof_tables.reserve(x.size());
247  std::generate_n(std::back_inserter(dof_tables), x.size(),
248  [&]() { return _local_to_global_index_map.get(); });
249 
250  if (process_id == 0)
251  {
254  _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt);
255 
256  BaseLib::RunTime time_phreeqc;
257  time_phreeqc.start();
258 
259  _chemical_solver_interface->setAqueousSolutionsPrevFromDumpFile();
260 
261  _chemical_solver_interface->executeSpeciationCalculation(dt);
262 
263  INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
264 
267  postSpeciationCalculation,
269 
270  return;
271  }
272 
273  auto const matrix_specification =
274  ode_sys.getMatrixSpecifications(process_id);
275 
276  std::size_t matrix_id = 0u;
278  matrix_specification, matrix_id);
280  matrix_specification, matrix_id);
281  auto& b =
283 
284  M.setZero();
285  K.setZero();
286  b.setZero();
287 
290  _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt, M, K,
291  b, process_id);
292 
293  // todo (renchao): incorporate Neumann boundary condition
297 
299  matrix_specification, matrix_id);
300  auto& rhs =
302 
303  A.setZero();
304  rhs.setZero();
305 
306  // compute A
307  MathLib::LinAlg::copy(M, A);
308  MathLib::LinAlg::aypx(A, 1.0 / dt, K);
309 
310  // compute rhs
311  MathLib::LinAlg::matMult(M, *x[process_id], rhs);
312  MathLib::LinAlg::aypx(rhs, 1.0 / dt, b);
313 
314  using Tag = NumLib::NonlinearSolverTag;
316  auto& equation_system = static_cast<EqSys&>(ode_sys);
317  equation_system.applyKnownSolutionsPicard(A, rhs, *x[process_id]);
318 
319  auto& linear_solver =
321  linear_solver.solve(A, rhs, *x[process_id]);
322 
327 }
328 
330  const double t,
331  std::vector<GlobalVector*> const& integration_point_values_vectors,
332  std::vector<GlobalVector*>& nodal_values_vectors)
333 {
334  auto& extrapolator = getExtrapolator();
335  auto const extrapolatables =
338  getInterpolatedLocalSolution);
339 
340  for (unsigned transport_process_id = 0;
341  transport_process_id < integration_point_values_vectors.size();
342  ++transport_process_id)
343  {
344  auto const& pv = _process_variables[transport_process_id + 1][0].get();
345  auto const& int_pt_C =
346  integration_point_values_vectors[transport_process_id];
347 
348  extrapolator.extrapolate(pv.getNumberOfGlobalComponents(),
349  extrapolatables, t, {int_pt_C},
350  {_local_to_global_index_map.get()});
351 
352  auto const& nodal_values = extrapolator.getNodalValues();
353  MathLib::LinAlg::copy(nodal_values,
354  *nodal_values_vectors[transport_process_id + 1]);
356  *nodal_values_vectors[transport_process_id + 1]);
357  }
358 }
359 
361  double const t,
362  double const dt,
363  std::vector<GlobalVector*> const& x,
364  GlobalVector const& x_dot,
365  int const process_id)
366 {
367  if (process_id != 0)
368  {
369  return;
370  }
371 
372  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
373  dof_tables.reserve(x.size());
374  std::generate_n(std::back_inserter(dof_tables), x.size(),
375  [&]() { return _local_to_global_index_map.get(); });
376 
377  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
380  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x,
381  x_dot, process_id);
382 }
383 
385  std::vector<GlobalVector*> const& x,
386  const double t,
387  const double dt,
388  int const process_id)
389 {
390  if (process_id != 0)
391  {
392  return;
393  }
394 
395  std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
396  dof_tables.reserve(x.size());
397  std::generate_n(std::back_inserter(dof_tables), x.size(),
398  [&]() { return _local_to_global_index_map.get(); });
399 
400  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
403  _local_assemblers, pv.getActiveElementIDs(), dof_tables, x, t, dt);
404 
405  if (!_surfaceflux) // computing the surfaceflux is optional
406  {
407  return;
408  }
409  _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
410  pv.getActiveElementIDs());
411 }
412 
413 } // namespace ComponentTransport
414 } // namespace ProcessLib
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
Definition of the RunTime class.
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
bool solve(EigenMatrix &A, EigenVector &b, EigenVector &x)
Global vector based on Eigen vector.
Definition: EigenVector.h:28
bool isAxiallySymmetric() const
Definition: Mesh.h:126
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:71
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:98
std::size_t getNumberOfNodes() const
Get the number of nodes.
Definition: Mesh.h:89
virtual MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const =0
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=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 & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
void setInitialConditionsConcreteProcess(std::vector< GlobalVector * > &x, double const t, int const process_id) 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
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
void extrapolateIntegrationPointValuesToNodes(const double t, std::vector< GlobalVector * > const &integration_point_values_vectors, std::vector< GlobalVector * > &nodal_values_vectors) 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 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, const double t, const double dt, int const process_id) override
void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(int const process_id) override
std::vector< std::unique_ptr< ComponentTransportLocalAssemblerInterface > > _local_assemblers
std::unique_ptr< ProcessLib::SurfaceFluxData > _surfaceflux
Eigen::Vector3d getFlux(std::size_t const element_id, MathLib::Point3d const &p, double const t, std::vector< GlobalVector * > const &x) const override
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)
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)
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)
std::vector< std::size_t > const & getActiveElementIDs() const
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:180
MeshLib::Mesh & _mesh
Definition: Process.h:321
std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > _process_variables
Definition: Process.h:355
SecondaryVariableCollection _secondary_variables
Definition: Process.h:326
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:334
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:144
VectorMatrixAssembler _global_assembler
Definition: Process.h:328
unsigned const _integration_order
Definition: Process.h:339
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:324
const bool _use_monolithic_scheme
Definition: Process.h:330
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< 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 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)
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
Definition: Types.h:20
void finalizeAssembly(PETScMatrix &A)
Definition: LinAlg.cpp:162
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:141
void aypx(PETScVector &y, double const a, PETScVector const &x)
Definition: LinAlg.cpp:50
void finalizeVectorAssembly(VEC_T &)
General function to finalize the vector assembly.
static const double u
static const double t
bool finalizeMatrixAssembly(MAT_T &)
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalIndexType > &indices)
ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection > makeExtrapolatable(LocalAssemblerCollection const &local_assemblers, IntegrationPointValuesMethod integration_point_values_method)
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)
ChemistryLib::ChemicalSolverInterface *const chemical_solver_interface