OGS
ComponentTransportFEM.h
Go to the documentation of this file.
1 
11 #pragma once
12 
13 #include <Eigen/Core>
14 #include <Eigen/Sparse>
15 #include <numeric>
16 #include <vector>
17 
21 #include "MaterialLib/MPL/Medium.h"
32 #include "ParameterLib/Parameter.h"
36 
37 namespace ProcessLib
38 {
39 namespace ComponentTransport
40 {
41 template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
43 {
44  IntegrationPointData(NodalRowVectorType const& N_,
45  GlobalDimNodalMatrixType const& dNdx_,
46  double const& integration_weight_)
47  : N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
48  {
49  }
50 
52 
53  NodalRowVectorType const N;
54  GlobalDimNodalMatrixType const dNdx;
55  double const integration_weight;
56 
57  // -1 indicates that no chemical reaction takes place in the element to
58  // which the integration point belongs.
60 
61  double porosity = std::numeric_limits<double>::quiet_NaN();
62  double porosity_prev = std::numeric_limits<double>::quiet_NaN();
64 };
65 
69 {
70 public:
72 
73  virtual void setChemicalSystemID(std::size_t const /*mesh_item_id*/) = 0;
74 
76  std::size_t const mesh_item_id,
77  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
78  std::vector<GlobalVector*> const& x, double const t)
79  {
80  std::vector<double> local_x_vec;
81 
82  auto const n_processes = x.size();
83  for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
84  {
85  auto const indices =
86  NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
87  assert(!indices.empty());
88  auto const local_solution = x[process_id]->get(indices);
89  local_x_vec.insert(std::end(local_x_vec),
90  std::begin(local_solution),
91  std::end(local_solution));
92  }
93  auto const local_x = MathLib::toVector(local_x_vec);
94 
96  }
97 
99  std::size_t const mesh_item_id,
100  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
101  std::vector<GlobalVector*> const& x, double const t, double const dt)
102  {
103  std::vector<double> local_x_vec;
104 
105  auto const n_processes = x.size();
106  for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
107  {
108  auto const indices =
109  NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
110  assert(!indices.empty());
111  auto const local_solution = x[process_id]->get(indices);
112  local_x_vec.insert(std::end(local_x_vec),
113  std::begin(local_solution),
114  std::end(local_solution));
115  }
116  auto const local_x = MathLib::toVector(local_x_vec);
117 
118  setChemicalSystemConcrete(local_x, t, dt);
119  }
120 
122  std::size_t const mesh_item_id,
123  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
124  std::vector<GlobalVector*> const& x, double const t, double const dt,
125  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, int const process_id)
126  {
127  std::vector<double> local_x_vec;
128 
129  auto const n_processes = x.size();
130  for (std::size_t pcs_id = 0; pcs_id < n_processes; ++pcs_id)
131  {
132  auto const indices =
133  NumLib::getIndices(mesh_item_id, *dof_tables[pcs_id]);
134  assert(!indices.empty());
135  auto const local_solution = x[pcs_id]->get(indices);
136  local_x_vec.insert(std::end(local_x_vec),
137  std::begin(local_solution),
138  std::end(local_solution));
139  }
140  auto const local_x = MathLib::toVector(local_x_vec);
141 
142  auto const indices =
143  NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
144  auto const num_r_c = indices.size();
145 
146  std::vector<double> local_M_data;
147  local_M_data.reserve(num_r_c * num_r_c);
148  std::vector<double> local_K_data;
149  local_K_data.reserve(num_r_c * num_r_c);
150  std::vector<double> local_b_data;
151  local_b_data.reserve(num_r_c);
152 
153  assembleReactionEquationConcrete(t, dt, local_x, local_M_data,
154  local_K_data, local_b_data,
155  process_id);
156 
157  auto const r_c_indices =
159  if (!local_M_data.empty())
160  {
161  auto const local_M =
162  MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
163  M.add(r_c_indices, local_M);
164  }
165  if (!local_K_data.empty())
166  {
167  auto const local_K =
168  MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
169  K.add(r_c_indices, local_K);
170  }
171  if (!local_b_data.empty())
172  {
173  b.add(indices, local_b_data);
174  }
175  }
176 
177  virtual void postSpeciationCalculation(std::size_t const ele_id,
178  double const t, double const dt) = 0;
179 
181  std::size_t const ele_id) = 0;
182 
183  virtual std::vector<double> const& getIntPtDarcyVelocity(
184  const double t,
185  std::vector<GlobalVector*> const& x,
186  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
187  std::vector<double>& cache) const = 0;
188 
189  virtual std::vector<double> const& getIntPtMolarFlux(
190  const double t,
191  std::vector<GlobalVector*> const& x,
192  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
193  std::vector<double>& cache) const = 0;
194 
195 private:
197  Eigen::VectorXd const& /*local_x*/, double const /*t*/) = 0;
198 
199  virtual void setChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/,
200  double const /*t*/,
201  double const /*dt*/) = 0;
202 
204  double const t, double const dt, Eigen::VectorXd const& local_x,
205  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
206  std::vector<double>& local_b_data, int const transport_process_id) = 0;
207 };
208 
209 template <typename ShapeFunction, int GlobalDim>
211 {
212  // When monolithic scheme is adopted, nodal pressure and nodal concentration
213  // are accessed by vector index.
214  static const int pressure_index = 0;
215  static const int first_concentration_index = ShapeFunction::NPOINTS;
216 
217  static const int pressure_size = ShapeFunction::NPOINTS;
218  static const int concentration_size =
219  ShapeFunction::NPOINTS; // per component
220 
223 
225  typename ShapeMatricesType::template MatrixType<pressure_size,
226  pressure_size>;
228  typename ShapeMatricesType::template VectorType<pressure_size>;
229 
231  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
232  using LocalVectorType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
233 
236 
241 
242 public:
244  MeshLib::Element const& element,
245  std::size_t const local_matrix_size,
246  NumLib::GenericIntegrationMethod const& integration_method,
247  bool is_axially_symmetric,
248  ComponentTransportProcessData const& process_data,
249  std::vector<std::reference_wrapper<ProcessVariable>> const&
250  transport_process_variables)
251  : _element(element),
252  _process_data(process_data),
253  _integration_method(integration_method),
254  _transport_process_variables(transport_process_variables)
255  {
256  (void)local_matrix_size;
257 
258  unsigned const n_integration_points =
260  _ip_data.reserve(n_integration_points);
261 
263  pos.setElementID(_element.getID());
264 
265  auto const shape_matrices =
267  GlobalDim>(element, is_axially_symmetric,
269  auto const& medium =
270  *_process_data.media_map->getMedium(_element.getID());
271  for (unsigned ip = 0; ip < n_integration_points; ip++)
272  {
273  _ip_data.emplace_back(
274  shape_matrices[ip].N, shape_matrices[ip].dNdx,
276  shape_matrices[ip].integralMeasure *
277  shape_matrices[ip].detJ);
278 
279  pos.setIntegrationPoint(ip);
280 
281  _ip_data[ip].porosity =
283  .template initialValue<double>(
284  pos, std::numeric_limits<double>::quiet_NaN() /*t*/);
285 
286  _ip_data[ip].pushBackState();
287  }
288  }
289 
290  void setChemicalSystemID(std::size_t const /*mesh_item_id*/) override
291  {
293  // chemical system index map
294  auto& chemical_system_index_map =
296 
297  unsigned const n_integration_points =
299  for (unsigned ip = 0; ip < n_integration_points; ip++)
300  {
301  _ip_data[ip].chemical_system_id =
302  chemical_system_index_map.empty()
303  ? 0
304  : chemical_system_index_map.back() + 1;
305  chemical_system_index_map.push_back(
306  _ip_data[ip].chemical_system_id);
307  }
308  }
309 
310  void initializeChemicalSystemConcrete(Eigen::VectorXd const& local_x,
311  double const t) override
312  {
314 
315  auto const& medium =
316  *_process_data.media_map->getMedium(_element.getID());
317 
319  pos.setElementID(_element.getID());
320 
321  unsigned const n_integration_points =
323  for (unsigned ip = 0; ip < n_integration_points; ip++)
324  {
325  auto& ip_data = _ip_data[ip];
326  auto const& N = ip_data.N;
327  auto const& chemical_system_id = ip_data.chemical_system_id;
328 
329  auto const n_component = _transport_process_variables.size();
330  std::vector<double> C_int_pt(n_component);
331  for (unsigned component_id = 0; component_id < n_component;
332  ++component_id)
333  {
334  auto const concentration_index =
336  component_id * concentration_size;
337  auto const local_C =
338  local_x.template segment<concentration_size>(
339  concentration_index);
340 
342  C_int_pt[component_id]);
343  }
344 
346  ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id,
347  medium, pos, t);
348  }
349  }
350 
351  void setChemicalSystemConcrete(Eigen::VectorXd const& local_x,
352  double const t, double dt) override
353  {
355 
356  auto const& medium =
357  _process_data.media_map->getMedium(_element.getID());
358 
361 
363  pos.setElementID(_element.getID());
364 
365  unsigned const n_integration_points =
367  for (unsigned ip = 0; ip < n_integration_points; ip++)
368  {
369  auto& ip_data = _ip_data[ip];
370  auto const& N = ip_data.N;
371  auto& porosity = ip_data.porosity;
372  auto const& porosity_prev = ip_data.porosity_prev;
373  auto const& chemical_system_id = ip_data.chemical_system_id;
374 
375  auto const n_component = _transport_process_variables.size();
376  std::vector<double> C_int_pt(n_component);
377  for (unsigned component_id = 0; component_id < n_component;
378  ++component_id)
379  {
380  auto const concentration_index =
382  component_id * concentration_size;
383  auto const local_C =
384  local_x.template segment<concentration_size>(
385  concentration_index);
386 
388  C_int_pt[component_id]);
389  }
390 
391  {
392  vars_prev.porosity = porosity_prev;
393 
394  porosity =
396  ? porosity_prev
397  : medium
398  ->property(
400  .template value<double>(vars, vars_prev, pos, t,
401  dt);
402 
403  vars.porosity = porosity;
404  }
405 
407  C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
408  }
409  }
410 
411  void postSpeciationCalculation(std::size_t const ele_id, double const t,
412  double const dt) override
413  {
415  {
416  return;
417  }
418 
419  auto const& medium = *_process_data.media_map->getMedium(ele_id);
420 
422  pos.setElementID(ele_id);
423 
424  for (auto& ip_data : _ip_data)
425  {
426  ip_data.porosity = ip_data.porosity_prev;
427 
429  ->updateVolumeFractionPostReaction(ip_data.chemical_system_id,
430  medium, pos,
431  ip_data.porosity, t, dt);
432 
434  ip_data.chemical_system_id, medium, ip_data.porosity);
435  }
436  }
437 
438  void assemble(double const t, double const dt,
439  std::vector<double> const& local_x,
440  std::vector<double> const& /*local_xdot*/,
441  std::vector<double>& local_M_data,
442  std::vector<double>& local_K_data,
443  std::vector<double>& local_b_data) override
444  {
445  auto const local_matrix_size = local_x.size();
446  // Nodal DOFs include pressure
447  int const num_nodal_dof = 1 + _transport_process_variables.size();
448  // This assertion is valid only if all nodal d.o.f. use the same shape
449  // matrices.
450  assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
451 
452  auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
453  local_M_data, local_matrix_size, local_matrix_size);
454  auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
455  local_K_data, local_matrix_size, local_matrix_size);
456  auto local_b = MathLib::createZeroedVector<LocalVectorType>(
457  local_b_data, local_matrix_size);
458 
459  // Get block matrices
460  auto Kpp = local_K.template block<pressure_size, pressure_size>(
462  auto Mpp = local_M.template block<pressure_size, pressure_size>(
464  auto Bp = local_b.template segment<pressure_size>(pressure_index);
465 
466  auto local_p = Eigen::Map<const NodalVectorType>(
467  &local_x[pressure_index], pressure_size);
468 
469  auto const number_of_components = num_nodal_dof - 1;
470  for (int component_id = 0; component_id < number_of_components;
471  ++component_id)
472  {
473  /* Partitioned assembler matrix
474  * | pp | pc1 | pc2 | pc3 |
475  * |-----|-----|-----|-----|
476  * | c1p | c1c1| 0 | 0 |
477  * |-----|-----|-----|-----|
478  * | c2p | 0 | c2c2| 0 |
479  * |-----|-----|-----|-----|
480  * | c3p | 0 | 0 | c3c3|
481  */
482  auto concentration_index =
483  pressure_size + component_id * concentration_size;
484 
485  auto KCC =
486  local_K.template block<concentration_size, concentration_size>(
487  concentration_index, concentration_index);
488  auto MCC =
489  local_M.template block<concentration_size, concentration_size>(
490  concentration_index, concentration_index);
491  auto MCp =
492  local_M.template block<concentration_size, pressure_size>(
493  concentration_index, pressure_index);
494  auto MpC =
495  local_M.template block<pressure_size, concentration_size>(
496  pressure_index, concentration_index);
497 
498  auto local_C = Eigen::Map<const NodalVectorType>(
499  &local_x[concentration_index], concentration_size);
500 
501  assembleBlockMatrices(component_id, t, dt, local_C, local_p, KCC,
502  MCC, MCp, MpC, Kpp, Mpp, Bp);
503 
505  {
506  auto const stoichiometric_matrix =
509 
510  assert(stoichiometric_matrix);
511 
512  for (Eigen::SparseMatrix<double>::InnerIterator it(
513  *stoichiometric_matrix, component_id);
514  it;
515  ++it)
516  {
517  auto const stoichiometric_coefficient = it.value();
518  auto const coupled_component_id = it.row();
519  auto const kinetic_prefactor =
521  ->getKineticPrefactor(coupled_component_id);
522 
523  auto const concentration_index =
524  pressure_size + component_id * concentration_size;
525  auto const coupled_concentration_index =
526  pressure_size +
527  coupled_component_id * concentration_size;
528  auto KCmCn = local_K.template block<concentration_size,
530  concentration_index, coupled_concentration_index);
531 
532  // account for the coupling between components
533  assembleKCmCn(component_id, t, dt, KCmCn,
534  stoichiometric_coefficient,
535  kinetic_prefactor);
536  }
537  }
538  }
539  }
540 
542  GlobalDimMatrixType const& I,
543  GlobalDimMatrixType const& pore_diffusion_coefficient,
544  GlobalDimVectorType const& velocity,
545  double const porosity,
546  double const solute_dispersivity_transverse,
547  double const solute_dispersivity_longitudinal)
548  {
549  double const velocity_magnitude = velocity.norm();
550 
551  if (velocity_magnitude == 0.0)
552  {
553  return porosity * pore_diffusion_coefficient;
554  }
555 
557  porosity * pore_diffusion_coefficient +
558  solute_dispersivity_transverse * velocity_magnitude * I +
559  (solute_dispersivity_longitudinal -
560  solute_dispersivity_transverse) /
561  velocity_magnitude * velocity * velocity.transpose());
562 
564  {
565  return D_c;
566  }
567 
568  return D_c + _process_data.stabilizer->computeArtificialDiffusion(
569  _element.getID(), velocity_magnitude) *
570  I;
571  }
572 
574  int const component_id, double const t, double const dt,
575  Eigen::Ref<const NodalVectorType> const& C_nodal_values,
576  Eigen::Ref<const NodalVectorType> const& p_nodal_values,
577  Eigen::Ref<LocalBlockMatrixType> KCC,
578  Eigen::Ref<LocalBlockMatrixType> MCC,
579  Eigen::Ref<LocalBlockMatrixType> MCp,
580  Eigen::Ref<LocalBlockMatrixType> MpC,
581  Eigen::Ref<LocalBlockMatrixType> Kpp,
582  Eigen::Ref<LocalBlockMatrixType> Mpp,
583  Eigen::Ref<LocalSegmentVectorType> Bp)
584  {
585  unsigned const n_integration_points =
587 
589  pos.setElementID(_element.getID());
590 
591  auto const& b = _process_data.specific_body_force;
592 
594 
595  GlobalDimMatrixType const& I(
596  GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
597 
598  // Get material properties
599  auto const& medium =
600  *_process_data.media_map->getMedium(_element.getID());
601  // Select the only valid for component transport liquid phase.
602  auto const& phase = medium.phase("AqueousLiquid");
603 
604  // Assume that the component name is the same as the process variable
605  // name. Components are shifted by one because the first one is always
606  // pressure.
607  auto const& component = phase.component(
608  _transport_process_variables[component_id].get().getName());
609 
610  for (unsigned ip(0); ip < n_integration_points; ++ip)
611  {
612  pos.setIntegrationPoint(ip);
613 
614  auto& ip_data = _ip_data[ip];
615  auto const& N = ip_data.N;
616  auto const& dNdx = ip_data.dNdx;
617  auto const& w = ip_data.integration_weight;
618  auto& porosity = ip_data.porosity;
619 
620  double C_int_pt = 0.0;
621  double p_int_pt = 0.0;
622 
623  NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
624  NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
625 
626  vars.concentration = C_int_pt;
627  vars.phase_pressure = p_int_pt;
628 
629  // update according to a particular porosity model
631  .template value<double>(vars, pos, t, dt);
632  vars.porosity = porosity;
633 
634  auto const& retardation_factor =
636  .template value<double>(vars, pos, t, dt);
637 
638  auto const& solute_dispersivity_transverse = medium.template value<
639  double>(
641 
642  auto const& solute_dispersivity_longitudinal =
643  medium.template value<double>(
646 
647  // Use the fluid density model to compute the density
648  // TODO (renchao): concentration of which component as the argument
649  // for calculation of fluid density
650  auto const density =
652  .template value<double>(vars, pos, t, dt);
653 
654  auto const decay_rate =
656  .template value<double>(vars, pos, t, dt);
657 
658  auto const& pore_diffusion_coefficient =
659  MaterialPropertyLib::formEigenTensor<GlobalDim>(
661  .value(vars, pos, t, dt));
662 
663  auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
665  vars, pos, t, dt));
666 
667  // Use the viscosity model to compute the viscosity
668  auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
669  .template value<double>(vars, pos, t, dt);
670 
671  GlobalDimMatrixType const K_over_mu = K / mu;
672  GlobalDimVectorType const velocity =
674  ? GlobalDimVectorType(-K_over_mu *
675  (dNdx * p_nodal_values - density * b))
676  : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
677 
678  const double drho_dp =
680  .template dValue<double>(
682  pos, t, dt);
683 
684  const double drho_dC =
686  .template dValue<double>(
688  t, dt);
689 
690  GlobalDimMatrixType const hydrodynamic_dispersion =
691  getHydrodynamicDispersion(I, pore_diffusion_coefficient,
692  velocity, porosity,
693  solute_dispersivity_transverse,
694  solute_dispersivity_longitudinal);
695 
696  const double R_times_phi(retardation_factor * porosity);
697  GlobalDimVectorType const mass_density_flow = velocity * density;
698  auto const N_t_N = (N.transpose() * N).eval();
699 
701  {
702  MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
703  MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
704  KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
705  }
706  else
707  {
708  KCC.noalias() +=
709  N.transpose() * mass_density_flow.transpose() * dNdx * w;
710  }
711  MCC.noalias() += N_t_N * (R_times_phi * density * w);
712  KCC.noalias() += dNdx.transpose() * hydrodynamic_dispersion * dNdx *
713  (density * w) +
714  N_t_N * (decay_rate * R_times_phi * density * w);
715 
716  MpC.noalias() += N_t_N * (porosity * drho_dC * w);
717 
718  // Calculate Mpp, Kpp, and bp in the first loop over components
719  if (component_id == 0)
720  {
721  Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
722  Kpp.noalias() +=
723  dNdx.transpose() * K_over_mu * dNdx * (density * w);
724 
726  {
727  Bp.noalias() += dNdx.transpose() * K_over_mu * b *
728  (density * density * w);
729  }
730  }
731  }
732  }
733 
734  void assembleKCmCn(int const component_id, double const t, double const dt,
735  Eigen::Ref<LocalBlockMatrixType> KCmCn,
736  double const stoichiometric_coefficient,
737  double const kinetic_prefactor)
738  {
739  unsigned const n_integration_points =
741 
743  pos.setElementID(_element.getID());
744 
746 
747  auto const& medium =
748  *_process_data.media_map->getMedium(_element.getID());
749  auto const& phase = medium.phase("AqueousLiquid");
750  auto const& component = phase.component(
751  _transport_process_variables[component_id].get().getName());
752 
753  for (unsigned ip(0); ip < n_integration_points; ++ip)
754  {
755  auto& ip_data = _ip_data[ip];
756  auto const& N = ip_data.N;
757  auto const& w = ip_data.integration_weight;
758  auto& porosity = ip_data.porosity;
759 
760  auto const retardation_factor =
762  .template value<double>(vars, pos, t, dt);
763 
765  .template value<double>(vars, pos, t, dt);
766 
767  auto const density =
769  .template value<double>(vars, pos, t, dt);
770 
771  KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
772  kinetic_prefactor * retardation_factor *
773  porosity * density * N;
774  }
775  }
776 
777  void assembleForStaggeredScheme(double const t, double const dt,
778  Eigen::VectorXd const& local_x,
779  Eigen::VectorXd const& local_xdot,
780  int const process_id,
781  std::vector<double>& local_M_data,
782  std::vector<double>& local_K_data,
783  std::vector<double>& local_b_data) override
784  {
785  if (process_id == _process_data.hydraulic_process_id)
786  {
787  assembleHydraulicEquation(t, dt, local_x, local_xdot, local_M_data,
788  local_K_data, local_b_data);
789  }
790  else
791  {
792  // Go for assembling in an order of transport process id.
793  assembleComponentTransportEquation(t, dt, local_x, local_xdot,
794  local_M_data, local_K_data,
795  local_b_data, process_id);
796  }
797  }
798 
799  void assembleHydraulicEquation(double const t,
800  double const dt,
801  Eigen::VectorXd const& local_x,
802  Eigen::VectorXd const& local_xdot,
803  std::vector<double>& local_M_data,
804  std::vector<double>& local_K_data,
805  std::vector<double>& local_b_data)
806  {
807  auto const local_p =
808  local_x.template segment<pressure_size>(pressure_index);
809  auto const local_C = local_x.template segment<concentration_size>(
811  auto const local_Cdot =
812  local_xdot.segment<concentration_size>(first_concentration_index);
813 
814  auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
815  local_M_data, pressure_size, pressure_size);
816  auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
817  local_K_data, pressure_size, pressure_size);
818  auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
819  local_b_data, pressure_size);
820 
821  unsigned const n_integration_points =
823 
825  pos.setElementID(_element.getID());
826 
827  auto const& b = _process_data.specific_body_force;
828 
829  auto const& medium =
830  *_process_data.media_map->getMedium(_element.getID());
831  auto const& phase = medium.phase("AqueousLiquid");
832 
835 
836  for (unsigned ip(0); ip < n_integration_points; ++ip)
837  {
838  pos.setIntegrationPoint(ip);
839 
840  auto& ip_data = _ip_data[ip];
841  auto const& N = ip_data.N;
842  auto const& dNdx = ip_data.dNdx;
843  auto const& w = ip_data.integration_weight;
844  auto& porosity = ip_data.porosity;
845  auto const& porosity_prev = ip_data.porosity_prev;
846 
847  double C_int_pt = 0.0;
848  double p_int_pt = 0.0;
849 
850  NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
851  NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
852 
853  vars.concentration = C_int_pt;
854  vars.phase_pressure = p_int_pt;
855 
856  // porosity
857  {
858  vars_prev.porosity = porosity_prev;
859 
860  porosity =
862  ? porosity_prev
864  .template value<double>(vars, vars_prev, pos, t,
865  dt);
866 
867  vars.porosity = porosity;
868  }
869 
870  // Use the fluid density model to compute the density
871  // TODO: Concentration of which component as one of arguments for
872  // calculation of fluid density
873  auto const density =
875  .template value<double>(vars, pos, t, dt);
876 
877  auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
879  vars, pos, t, dt));
880 
881  // Use the viscosity model to compute the viscosity
882  auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
883  .template value<double>(vars, pos, t, dt);
884 
885  GlobalDimMatrixType const K_over_mu = K / mu;
886 
887  const double drho_dp =
889  .template dValue<double>(
891  pos, t, dt);
892  const double drho_dC =
894  .template dValue<double>(
896  t, dt);
897 
898  // matrix assembly
899  local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
900  local_K.noalias() +=
901  w * dNdx.transpose() * density * K_over_mu * dNdx;
902 
904  {
905  local_b.noalias() +=
906  w * density * density * dNdx.transpose() * K_over_mu * b;
907  }
908 
909  // coupling term
910  {
911  double dot_C_int_pt = 0.0;
912  NumLib::shapeFunctionInterpolate(local_Cdot, N, dot_C_int_pt);
913 
914  local_b.noalias() -=
915  w * N.transpose() * porosity * drho_dC * dot_C_int_pt;
916  }
917  }
918  }
919 
921  double const t, double const dt, Eigen::VectorXd const& local_x,
922  Eigen::VectorXd const& local_xdot, std::vector<double>& local_M_data,
923  std::vector<double>& local_K_data,
924  std::vector<double>& /*local_b_data*/, int const transport_process_id)
925  {
926  auto const local_p =
927  local_x.template segment<pressure_size>(pressure_index);
928  auto const local_C = local_x.template segment<concentration_size>(
930  (transport_process_id - 1) * concentration_size);
931  auto const local_pdot =
932  local_xdot.segment<pressure_size>(pressure_index);
933 
934  NodalVectorType local_T;
936  {
937  local_T =
939  }
940 
941  auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
942  local_M_data, concentration_size, concentration_size);
943  auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
944  local_K_data, concentration_size, concentration_size);
945 
946  unsigned const n_integration_points =
948 
950  pos.setElementID(_element.getID());
951 
952  auto const& b = _process_data.specific_body_force;
953 
956 
957  GlobalDimMatrixType const& I(
958  GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
959 
960  auto const& medium =
961  *_process_data.media_map->getMedium(_element.getID());
962  auto const& phase = medium.phase("AqueousLiquid");
963  // Hydraulic process id is 0 and thus transport process id starts
964  // from 1.
965  auto const component_id = transport_process_id - 1;
966  auto const& component = phase.component(
967  _transport_process_variables[component_id].get().getName());
968 
969  for (unsigned ip(0); ip < n_integration_points; ++ip)
970  {
971  pos.setIntegrationPoint(ip);
972 
973  auto& ip_data = _ip_data[ip];
974  auto const& N = ip_data.N;
975  auto const& dNdx = ip_data.dNdx;
976  auto const& w = ip_data.integration_weight;
977  auto& porosity = ip_data.porosity;
978  auto const& porosity_prev = ip_data.porosity_prev;
979 
980  double C_int_pt = 0.0;
981  double p_int_pt = 0.0;
982 
983  NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
984  NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
985 
986  vars.concentration = C_int_pt;
987  vars.phase_pressure = p_int_pt;
988 
990  {
991  vars.temperature = N.dot(local_T);
992  }
993 
994  // porosity
995  {
996  vars_prev.porosity = porosity_prev;
997 
998  porosity =
1000  ? porosity_prev
1002  .template value<double>(vars, vars_prev, pos, t,
1003  dt);
1004 
1005  vars.porosity = porosity;
1006  }
1007 
1008  auto const& retardation_factor =
1010  .template value<double>(vars, pos, t, dt);
1011 
1012  auto const& solute_dispersivity_transverse = medium.template value<
1013  double>(
1015  auto const& solute_dispersivity_longitudinal =
1016  medium.template value<double>(
1019 
1020  // Use the fluid density model to compute the density
1021  auto const density =
1023  .template value<double>(vars, pos, t, dt);
1024  auto const decay_rate =
1026  .template value<double>(vars, pos, t, dt);
1027 
1028  auto const& pore_diffusion_coefficient =
1029  MaterialPropertyLib::formEigenTensor<GlobalDim>(
1031  .value(vars, pos, t, dt));
1032 
1033  auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1035  vars, pos, t, dt));
1036  // Use the viscosity model to compute the viscosity
1037  auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1038  .template value<double>(vars, pos, t, dt);
1039 
1040  GlobalDimMatrixType const K_over_mu = K / mu;
1041  GlobalDimVectorType const velocity =
1043  ? GlobalDimVectorType(-K_over_mu *
1044  (dNdx * local_p - density * b))
1045  : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1046 
1047  GlobalDimMatrixType const hydrodynamic_dispersion =
1048  getHydrodynamicDispersion(I, pore_diffusion_coefficient,
1049  velocity, porosity,
1050  solute_dispersivity_transverse,
1051  solute_dispersivity_longitudinal);
1052 
1053  double const R_times_phi = retardation_factor * porosity;
1054  auto const N_t_N = (N.transpose() * N).eval();
1055 
1057  {
1058  const double drho_dC =
1060  .template dValue<double>(
1062  pos, t, dt);
1063  local_M.noalias() +=
1064  N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1065  }
1066 
1067  local_M.noalias() += N_t_N * (R_times_phi * density * w);
1068 
1069  // coupling term
1071  {
1072  double dot_p_int_pt = 0.0;
1073 
1074  NumLib::shapeFunctionInterpolate(local_pdot, N, dot_p_int_pt);
1075  const double drho_dp =
1077  .template dValue<double>(
1079  pos, t, dt);
1080 
1081  local_K.noalias() +=
1082  N_t_N * ((R_times_phi * drho_dp * dot_p_int_pt) * w) -
1083  dNdx.transpose() * velocity * N * (density * w);
1084  }
1085  else
1086  {
1087  local_K.noalias() +=
1088  N.transpose() * velocity.transpose() * dNdx * (density * w);
1089  }
1090  local_K.noalias() +=
1091  dNdx.transpose() * hydrodynamic_dispersion * dNdx *
1092  (density * w) +
1093  N_t_N * (decay_rate * R_times_phi * density * w);
1094  }
1095  }
1096 
1098  double const t, double const dt, Eigen::VectorXd const& local_x,
1099  Eigen::VectorXd const& local_xdot, int const process_id,
1100  std::vector<double>& /*local_M_data*/,
1101  std::vector<double>& /*local_K_data*/,
1102  std::vector<double>& local_b_data,
1103  std::vector<double>& local_Jac_data) override
1104  {
1105  if (process_id == _process_data.hydraulic_process_id)
1106  {
1107  assembleWithJacobianHydraulicEquation(t, dt, local_x, local_xdot,
1108  local_b_data, local_Jac_data);
1109  }
1110  else
1111  {
1112  int const component_id = process_id - 1;
1114  t, dt, local_x, local_xdot, local_b_data, local_Jac_data,
1115  component_id);
1116  }
1117  }
1118 
1120  double const t, double const dt, Eigen::VectorXd const& local_x,
1121  Eigen::VectorXd const& local_xdot, std::vector<double>& local_b_data,
1122  std::vector<double>& local_Jac_data)
1123  {
1124  auto const p = local_x.template segment<pressure_size>(pressure_index);
1125  auto const c = local_x.template segment<concentration_size>(
1127 
1128  auto const pdot = local_xdot.segment<pressure_size>(pressure_index);
1129  auto const cdot =
1130  local_xdot.segment<concentration_size>(first_concentration_index);
1131 
1132  auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1133  local_Jac_data, pressure_size, pressure_size);
1134  auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1135  local_b_data, pressure_size);
1136 
1137  unsigned const n_integration_points =
1139 
1141  pos.setElementID(_element.getID());
1142 
1143  auto const& b = _process_data.specific_body_force;
1144 
1145  auto const& medium =
1146  *_process_data.media_map->getMedium(_element.getID());
1147  auto const& phase = medium.phase("AqueousLiquid");
1148 
1151 
1152  for (unsigned ip(0); ip < n_integration_points; ++ip)
1153  {
1154  pos.setIntegrationPoint(ip);
1155 
1156  auto& ip_data = _ip_data[ip];
1157  auto const& N = ip_data.N;
1158  auto const& dNdx = ip_data.dNdx;
1159  auto const& w = ip_data.integration_weight;
1160  auto& phi = ip_data.porosity;
1161  auto const& phi_prev = ip_data.porosity_prev;
1162 
1163  double const p_ip = N.dot(p);
1164  double const c_ip = N.dot(c);
1165 
1166  double const cdot_ip = N.dot(cdot);
1167 
1168  vars.phase_pressure = p_ip;
1169  vars.concentration = c_ip;
1170 
1171  // porosity
1172  {
1173  vars_prev.porosity = phi_prev;
1174 
1176  ? phi_prev
1178  .template value<double>(vars, vars_prev, pos, t,
1179  dt);
1180 
1181  vars.porosity = phi;
1182  }
1183 
1184  auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1185  .template value<double>(vars, pos, t, dt);
1186 
1187  auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1189  vars, pos, t, dt));
1190 
1191  auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1192  .template value<double>(vars, pos, t, dt);
1193 
1194  auto const drho_dp =
1196  .template dValue<double>(
1198  pos, t, dt);
1199  auto const drho_dc =
1201  .template dValue<double>(
1203  t, dt);
1204 
1205  // matrix assembly
1206  local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1207  w * dNdx.transpose() * rho * k / mu * dNdx;
1208 
1209  local_rhs.noalias() -=
1210  w * N.transpose() * phi *
1211  (drho_dp * N * pdot + drho_dc * cdot_ip) +
1212  w * rho * dNdx.transpose() * k / mu * dNdx * p;
1213 
1215  {
1216  local_rhs.noalias() +=
1217  w * rho * dNdx.transpose() * k / mu * rho * b;
1218  }
1219  }
1220  }
1221 
1223  double const t, double const dt, Eigen::VectorXd const& local_x,
1224  Eigen::VectorXd const& local_xdot, std::vector<double>& local_b_data,
1225  std::vector<double>& local_Jac_data, int const component_id)
1226  {
1227  auto const concentration_index =
1229 
1230  auto const p = local_x.template segment<pressure_size>(pressure_index);
1231  auto const c =
1232  local_x.template segment<concentration_size>(concentration_index);
1233  auto const cdot =
1234  local_xdot.segment<concentration_size>(concentration_index);
1235 
1236  NodalVectorType T;
1238  {
1240  }
1241 
1242  auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1243  local_Jac_data, concentration_size, concentration_size);
1244  auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1245  local_b_data, concentration_size);
1246 
1247  unsigned const n_integration_points =
1249 
1251  pos.setElementID(_element.getID());
1252 
1253  auto const& b = _process_data.specific_body_force;
1254 
1257 
1258  GlobalDimMatrixType const& I(
1259  GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
1260 
1261  auto const& medium =
1262  *_process_data.media_map->getMedium(_element.getID());
1263  auto const& phase = medium.phase("AqueousLiquid");
1264  auto const& component = phase.component(
1265  _transport_process_variables[component_id].get().getName());
1266 
1267  for (unsigned ip(0); ip < n_integration_points; ++ip)
1268  {
1269  pos.setIntegrationPoint(ip);
1270 
1271  auto& ip_data = _ip_data[ip];
1272  auto const& N = ip_data.N;
1273  auto const& dNdx = ip_data.dNdx;
1274  auto const& w = ip_data.integration_weight;
1275  auto& phi = ip_data.porosity;
1276  auto const& phi_prev = ip_data.porosity_prev;
1277 
1278  double const p_ip = N.dot(p);
1279  double const c_ip = N.dot(c);
1280 
1281  vars.phase_pressure = p_ip;
1282  vars.concentration = c_ip;
1283 
1285  {
1286  vars.temperature = N.dot(T);
1287  }
1288 
1289  // porosity
1290  {
1291  vars_prev.porosity = phi_prev;
1292 
1294  ? phi_prev
1296  .template value<double>(vars, vars_prev, pos, t,
1297  dt);
1298 
1299  vars.porosity = phi;
1300  }
1301 
1302  auto const R =
1304  .template value<double>(vars, pos, t, dt);
1305 
1306  auto const alpha_T = medium.template value<double>(
1308  auto const alpha_L = medium.template value<double>(
1310 
1311  auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1312  .template value<double>(vars, pos, t, dt);
1313  // first-order decay constant
1314  auto const alpha =
1316  .template value<double>(vars, pos, t, dt);
1317 
1318  auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1320  .value(vars, pos, t, dt));
1321 
1322  auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1324  vars, pos, t, dt));
1325  auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1326  .template value<double>(vars, pos, t, dt);
1327  // Darcy flux
1328  GlobalDimVectorType const q =
1330  ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1331  : GlobalDimVectorType(-k / mu * dNdx * p);
1332 
1333  GlobalDimMatrixType const D =
1334  getHydrodynamicDispersion(I, Dp, q, phi, alpha_T, alpha_L);
1335 
1336  // matrix assembly
1337  local_Jac.noalias() +=
1338  w * rho *
1339  (N.transpose() * phi * R * (alpha + 1 / dt) * N +
1340  dNdx.transpose() * D * dNdx);
1341 
1342  local_rhs.noalias() -=
1343  w * rho *
1344  (N.transpose() * phi * R * N * (cdot + alpha * c) +
1345  dNdx.transpose() * D * dNdx * c);
1346 
1348  local_Jac.noalias() +=
1349  w * rho * N.transpose() * q.transpose() * dNdx;
1350 
1351  local_rhs.noalias() -=
1352  w * rho * N.transpose() * q.transpose() * dNdx * c;
1353  }
1354  }
1355 
1357  double const t, double const dt, Eigen::VectorXd const& local_x,
1358  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1359  std::vector<double>& local_b_data,
1360  int const transport_process_id) override
1361  {
1362  auto const local_C = local_x.template segment<concentration_size>(
1364  (transport_process_id - 1) * concentration_size);
1365 
1366  auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1367  local_M_data, concentration_size, concentration_size);
1368  auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1369  local_K_data, concentration_size, concentration_size);
1370  auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1371  local_b_data, concentration_size);
1372 
1373  unsigned const n_integration_points =
1375 
1377  pos.setElementID(_element.getID());
1378 
1381 
1382  auto const& medium =
1383  *_process_data.media_map->getMedium(_element.getID());
1384  auto const component_id = transport_process_id - 1;
1385  for (unsigned ip(0); ip < n_integration_points; ++ip)
1386  {
1387  pos.setIntegrationPoint(ip);
1388 
1389  auto& ip_data = _ip_data[ip];
1390  auto const& N = ip_data.N;
1391  auto const w = ip_data.integration_weight;
1392  auto& porosity = ip_data.porosity;
1393  auto const& porosity_prev = ip_data.porosity_prev;
1394  auto const chemical_system_id = ip_data.chemical_system_id;
1395 
1396  double C_int_pt = 0.0;
1397  NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
1398 
1399  vars.concentration = C_int_pt;
1400 
1401  auto const porosity_dot = (porosity - porosity_prev) / dt;
1402 
1403  // porosity
1404  {
1405  vars_prev.porosity = porosity_prev;
1406 
1407  porosity =
1409  ? porosity_prev
1411  .template value<double>(vars, vars_prev, pos, t,
1412  dt);
1413  }
1414 
1415  local_M.noalias() += w * N.transpose() * porosity * N;
1416 
1417  local_K.noalias() += w * N.transpose() * porosity_dot * N;
1418 
1419  if (chemical_system_id == -1)
1420  {
1421  continue;
1422  }
1423 
1424  auto const C_post_int_pt =
1426  component_id, chemical_system_id);
1427 
1428  local_b.noalias() +=
1429  w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1430  }
1431  }
1432 
1433  std::vector<double> const& getIntPtDarcyVelocity(
1434  const double t,
1435  std::vector<GlobalVector*> const& x,
1436  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
1437  std::vector<double>& cache) const override
1438  {
1439  assert(x.size() == dof_table.size());
1440 
1441  auto const n_processes = x.size();
1442  std::vector<std::vector<double>> local_x;
1443  local_x.reserve(n_processes);
1444 
1445  for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1446  {
1447  auto const indices =
1448  NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1449  assert(!indices.empty());
1450  local_x.push_back(x[process_id]->get(indices));
1451  }
1452 
1453  // only one process, must be monolithic.
1454  if (n_processes == 1)
1455  {
1456  auto const local_p = Eigen::Map<const NodalVectorType>(
1457  &local_x[0][pressure_index], pressure_size);
1458  auto const local_C = Eigen::Map<const NodalVectorType>(
1460  return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1461  }
1462 
1463  // multiple processes, must be staggered.
1464  {
1465  constexpr int pressure_process_id = 0;
1466  constexpr int concentration_process_id = 1;
1467  auto const local_p = Eigen::Map<const NodalVectorType>(
1468  &local_x[pressure_process_id][0], pressure_size);
1469  auto const local_C = Eigen::Map<const NodalVectorType>(
1470  &local_x[concentration_process_id][0], concentration_size);
1471  return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1472  }
1473  }
1474 
1475  std::vector<double> const& calculateIntPtDarcyVelocity(
1476  const double t,
1477  Eigen::Ref<const NodalVectorType> const& p_nodal_values,
1478  Eigen::Ref<const NodalVectorType> const& C_nodal_values,
1479  std::vector<double>& cache) const
1480  {
1481  auto const n_integration_points =
1483 
1484  cache.clear();
1485  auto cache_mat = MathLib::createZeroedMatrix<
1486  Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1487  cache, GlobalDim, n_integration_points);
1488 
1490  pos.setElementID(_element.getID());
1491 
1493 
1494  auto const& medium =
1495  *_process_data.media_map->getMedium(_element.getID());
1496  auto const& phase = medium.phase("AqueousLiquid");
1497 
1498  for (unsigned ip = 0; ip < n_integration_points; ++ip)
1499  {
1500  auto const& ip_data = _ip_data[ip];
1501  auto const& N = ip_data.N;
1502  auto const& dNdx = ip_data.dNdx;
1503  auto const& porosity = ip_data.porosity;
1504 
1505  pos.setIntegrationPoint(ip);
1506 
1507  double C_int_pt = 0.0;
1508  double p_int_pt = 0.0;
1509 
1510  NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1511  NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1512 
1513  vars.concentration = C_int_pt;
1514  vars.phase_pressure = p_int_pt;
1515  vars.porosity = porosity;
1516 
1517  // TODO (naumov) Temporary value not used by current material
1518  // models. Need extension of secondary variables interface.
1519  double const dt = std::numeric_limits<double>::quiet_NaN();
1520  auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1522  vars, pos, t, dt));
1523  auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1524  .template value<double>(vars, pos, t, dt);
1525  GlobalDimMatrixType const K_over_mu = K / mu;
1526 
1527  cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1529  {
1530  auto const rho_w =
1532  .template value<double>(vars, pos, t, dt);
1533  auto const b = _process_data.specific_body_force;
1534  // here it is assumed that the vector b is directed 'downwards'
1535  cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1536  }
1537  }
1538 
1539  return cache;
1540  }
1541 
1542  Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
1543  const unsigned integration_point) const override
1544  {
1545  auto const& N = _ip_data[integration_point].N;
1546 
1547  // assumes N is stored contiguously in memory
1548  return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1549  }
1550 
1551  Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
1552  double const t,
1553  std::vector<double> const& local_x) const override
1554  {
1555  auto const local_p = Eigen::Map<const NodalVectorType>(
1556  &local_x[pressure_index], pressure_size);
1557  auto const local_C = Eigen::Map<const NodalVectorType>(
1559 
1560  // Eval shape matrices at given point
1561  // Note: Axial symmetry is set to false here, because we only need dNdx
1562  // here, which is not affected by axial symmetry.
1563  auto const shape_matrices =
1565  GlobalDim>(
1566  _element, false /*is_axially_symmetric*/,
1567  std::array{pnt_local_coords})[0];
1568 
1570  pos.setElementID(_element.getID());
1571 
1573 
1574  auto const& medium =
1575  *_process_data.media_map->getMedium(_element.getID());
1576  auto const& phase = medium.phase("AqueousLiquid");
1577 
1578  // local_x contains the local concentration and pressure values
1579  double c_int_pt;
1580  NumLib::shapeFunctionInterpolate(local_C, shape_matrices.N, c_int_pt);
1581  vars.concentration = c_int_pt;
1582 
1583  double p_int_pt;
1584  NumLib::shapeFunctionInterpolate(local_p, shape_matrices.N, p_int_pt);
1585  vars.phase_pressure = p_int_pt;
1586 
1587  // TODO (naumov) Temporary value not used by current material models.
1588  // Need extension of secondary variables interface.
1589  double const dt = std::numeric_limits<double>::quiet_NaN();
1590  auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1592  vars, pos, t, dt));
1593 
1594  auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1595  .template value<double>(vars, pos, t, dt);
1596  GlobalDimMatrixType const K_over_mu = K / mu;
1597 
1598  GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
1599  auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
1600  .template value<double>(vars, pos, t, dt);
1602  {
1603  auto const b = _process_data.specific_body_force;
1604  q += K_over_mu * rho_w * b;
1605  }
1606  Eigen::Vector3d flux(0.0, 0.0, 0.0);
1607  flux.head<GlobalDim>() = rho_w * q;
1608  return flux;
1609  }
1610 
1612  double const t,
1613  double const /*dt*/,
1614  Eigen::VectorXd const& local_x,
1615  Eigen::VectorXd const& /*local_x_dot*/) override
1616  {
1617  auto const local_p =
1618  local_x.template segment<pressure_size>(pressure_index);
1619  auto const local_C = local_x.template segment<concentration_size>(
1621 
1622  std::vector<double> ele_velocity;
1623  calculateIntPtDarcyVelocity(t, local_p, local_C, ele_velocity);
1624 
1625  auto const n_integration_points =
1627  auto const ele_velocity_mat =
1628  MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
1629 
1630  auto const ele_id = _element.getID();
1631  Eigen::Map<LocalVectorType>(
1632  &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
1633  GlobalDim) =
1634  ele_velocity_mat.rowwise().sum() / n_integration_points;
1635  }
1636 
1638  std::size_t const ele_id) override
1639  {
1640  auto const n_integration_points =
1642 
1644  {
1645  auto const& medium = *_process_data.media_map->getMedium(ele_id);
1646 
1647  for (auto& ip_data : _ip_data)
1648  {
1649  ip_data.porosity = ip_data.porosity_prev;
1650 
1652  ->updatePorosityPostReaction(ip_data.chemical_system_id,
1653  medium, ip_data.porosity);
1654  }
1655 
1656  (*_process_data.mesh_prop_porosity)[ele_id] =
1657  std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
1658  [](double const s, auto const& ip)
1659  { return s + ip.porosity; }) /
1660  n_integration_points;
1661  }
1662 
1663  std::vector<GlobalIndexType> chemical_system_indices;
1664  chemical_system_indices.reserve(n_integration_points);
1665  std::transform(_ip_data.begin(), _ip_data.end(),
1666  std::back_inserter(chemical_system_indices),
1667  [](auto const& ip_data)
1668  { return ip_data.chemical_system_id; });
1669 
1671  ele_id, chemical_system_indices);
1672  }
1673 
1674  std::vector<double> const& getIntPtMolarFlux(
1675  const double t,
1676  std::vector<GlobalVector*> const& x,
1677  std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
1678  std::vector<double>& cache) const override
1679  {
1680  std::vector<double> local_x_vec;
1681 
1682  auto const n_processes = x.size();
1683  for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1684  {
1685  auto const indices =
1686  NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
1687  assert(!indices.empty());
1688  auto const local_solution = x[process_id]->get(indices);
1689  local_x_vec.insert(std::end(local_x_vec),
1690  std::begin(local_solution),
1691  std::end(local_solution));
1692  }
1693  auto const local_x = MathLib::toVector(local_x_vec);
1694 
1695  auto const p = local_x.template segment<pressure_size>(pressure_index);
1696  auto const c = local_x.template segment<concentration_size>(
1698 
1699  auto const n_integration_points =
1701 
1702  cache.clear();
1703  auto cache_mat = MathLib::createZeroedMatrix<
1704  Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1705  cache, GlobalDim, n_integration_points);
1706 
1708  pos.setElementID(_element.getID());
1709 
1710  auto const& b = _process_data.specific_body_force;
1711 
1713 
1714  GlobalDimMatrixType const& I(
1715  GlobalDimMatrixType::Identity(GlobalDim, GlobalDim));
1716 
1717  auto const& medium =
1718  *_process_data.media_map->getMedium(_element.getID());
1719  auto const& phase = medium.phase("AqueousLiquid");
1720 
1721  int const component_id = 0;
1722  auto const& component = phase.component(
1723  _transport_process_variables[component_id].get().getName());
1724 
1725  for (unsigned ip = 0; ip < n_integration_points; ++ip)
1726  {
1727  auto const& ip_data = _ip_data[ip];
1728  auto const& N = ip_data.N;
1729  auto const& dNdx = ip_data.dNdx;
1730  auto const& phi = ip_data.porosity;
1731 
1732  pos.setIntegrationPoint(ip);
1733 
1734  double const p_ip = N.dot(p);
1735  double const c_ip = N.dot(c);
1736 
1737  vars.concentration = c_ip;
1738  vars.phase_pressure = p_ip;
1739  vars.porosity = phi;
1740 
1741  double const dt = std::numeric_limits<double>::quiet_NaN();
1742 
1743  auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1745  vars, pos, t, dt));
1746  auto const mu = phase[MaterialPropertyLib::PropertyType::viscosity]
1747  .template value<double>(vars, pos, t, dt);
1748  auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1749  .template value<double>(vars, pos, t, dt);
1750 
1751  // Darcy flux
1752  GlobalDimVectorType const q =
1754  ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1755  : GlobalDimVectorType(-k / mu * dNdx * p);
1756 
1757  auto const alpha_T = medium.template value<double>(
1759  auto const alpha_L = medium.template value<double>(
1761  auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1763  .value(vars, pos, t, dt));
1764 
1765  double const q_magnitude = q.norm();
1766  // hydrodymanic dispersion
1767  GlobalDimMatrixType const D =
1768  q_magnitude != 0.0
1769  ? GlobalDimMatrixType(phi * Dp + alpha_T * q_magnitude * I +
1770  (alpha_L - alpha_T) / q_magnitude *
1771  q * q.transpose())
1772  : GlobalDimMatrixType(phi * Dp);
1773 
1774  cache_mat.col(ip).noalias() = q * c_ip - phi * D * dNdx * c;
1775  }
1776 
1777  return cache;
1778  }
1779 
1780  void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
1781  double const /*t*/, double const /*dt*/) override
1782  {
1783  unsigned const n_integration_points =
1785 
1786  for (unsigned ip = 0; ip < n_integration_points; ip++)
1787  {
1788  _ip_data[ip].pushBackState();
1789  }
1790  }
1791 
1792 private:
1795 
1797  std::vector<std::reference_wrapper<ProcessVariable>> const
1799 
1800  std::vector<
1802  Eigen::aligned_allocator<
1805 };
1806 
1807 } // namespace ComponentTransport
1808 } // namespace ProcessLib
GlobalMatrix::IndexType GlobalIndexType
std::string getName(std::string const &line)
Returns the name/title from the "Zone"-description.
virtual Eigen::SparseMatrix< double > const * getStoichiometricMatrix() const
virtual void updatePorosityPostReaction(GlobalIndexType const &, MaterialPropertyLib::Medium const &, double &)
virtual void computeSecondaryVariable(std::size_t const, std::vector< GlobalIndexType > const &)
virtual double getKineticPrefactor([[maybe_unused]] std::size_t reaction_id) const
virtual void setChemicalSystemConcrete(std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const *, MaterialPropertyLib::VariableArray const &, ParameterLib::SpatialPosition const &, double const, double const)
virtual void updateVolumeFractionPostReaction(GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const, double const, double const)
std::vector< GlobalIndexType > chemical_system_index_map
virtual double getConcentration(int const, GlobalIndexType const) const
virtual void initializeChemicalSystemConcrete(std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const)
int add(IndexType row, IndexType col, double val)
Definition: EigenMatrix.h:86
Global vector based on Eigen vector.
Definition: EigenVector.h:25
void add(IndexType rowId, double v)
add entry
Definition: EigenVector.h:76
double getWeight() const
Definition: WeightedPoint.h:80
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
void setElementID(std::size_t element_id)
void setIntegrationPoint(unsigned integration_point)
virtual void setChemicalSystemConcrete(Eigen::VectorXd const &, double const, double const)=0
virtual void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt)=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 computeReactionRelatedSecondaryVariable(std::size_t const ele_id)=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)
virtual void assembleReactionEquationConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id)=0
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 void initializeChemicalSystemConcrete(Eigen::VectorXd const &, double const)=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
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
GlobalDimMatrixType getHydrodynamicDispersion(GlobalDimMatrixType const &I, GlobalDimMatrixType const &pore_diffusion_coefficient, GlobalDimVectorType const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > LocalMatrixType
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt) override
NumLib::GenericIntegrationMethod const & _integration_method
void assembleWithJacobianComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, int const component_id)
void initializeChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t) override
void computeReactionRelatedSecondaryVariable(std::size_t const ele_id) override
typename ShapeMatricesType::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
std::vector< double > const & calculateIntPtDarcyVelocity(const double t, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< const NodalVectorType > const &C_nodal_values, std::vector< double > &cache) const
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
void assembleHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< double > &cache) const override
void postTimestepConcrete(Eigen::VectorXd const &, double const, double const) override
ComponentTransportProcessData const & _process_data
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 override
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
void assembleWithJacobianHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void computeSecondaryVariableConcrete(double const t, double const, Eigen::VectorXd const &local_x, Eigen::VectorXd const &) override
std::vector< std::reference_wrapper< ProcessVariable > > const _transport_process_variables
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > > > _ip_data
void assembleComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &, int const transport_process_id)
typename ShapeMatricesType::template MatrixType< pressure_size, pressure_size > LocalBlockMatrixType
LocalAssemblerData(MeshLib::Element const &element, std::size_t const local_matrix_size, NumLib::GenericIntegrationMethod const &integration_method, bool is_axially_symmetric, ComponentTransportProcessData const &process_data, std::vector< std::reference_wrapper< ProcessVariable >> const &transport_process_variables)
typename ShapeMatricesType::template VectorType< pressure_size > LocalSegmentVectorType
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, int const process_id, std::vector< double > &, std::vector< double > &, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
Eigen::Vector3d getFlux(MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
void assemble(double const t, double const dt, std::vector< double > const &local_x, std::vector< double > const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
Eigen::Matrix< double, Eigen::Dynamic, 1 > LocalVectorType
void assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_xdot, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
void assembleReactionEquationConcrete(double const t, double const dt, Eigen::VectorXd const &local_x, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data, int const transport_process_id) override
void assembleBlockMatrices(int const component_id, double const t, double const dt, Eigen::Ref< const NodalVectorType > const &C_nodal_values, Eigen::Ref< const NodalVectorType > const &p_nodal_values, Eigen::Ref< LocalBlockMatrixType > KCC, Eigen::Ref< LocalBlockMatrixType > MCC, Eigen::Ref< LocalBlockMatrixType > MCp, Eigen::Ref< LocalBlockMatrixType > MpC, Eigen::Ref< LocalBlockMatrixType > Kpp, Eigen::Ref< LocalBlockMatrixType > Mpp, Eigen::Ref< LocalSegmentVectorType > Bp)
void assembleKCmCn(int const component_id, double const t, double const dt, Eigen::Ref< LocalBlockMatrixType > KCmCn, double const stoichiometric_coefficient, double const kinetic_prefactor)
void setChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t, double dt) override
typename ShapeMatricesType::NodalVectorType NodalVectorType
@ longitudinal_dispersivity
used to compute the hydrodynamic dispersion tensor.
Definition: PropertyType.h:59
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
Definition: PropertyType.h:104
@ retardation_factor
specify retardation factor used in component transport process.
Definition: PropertyType.h:82
static const double q
Eigen::Map< const Vector > toVector(std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector.
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:32
static const double s
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:72
static const double t
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:80
std::vector< GlobalIndexType > getIndices(std::size_t const mesh_item_id, NumLib::LocalToGlobalIndexMap const &dof_table)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > initShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, IntegrationMethod const &integration_method)
NumLib::ShapeMatrices< NodalRowVectorType, DimNodalMatrixType, DimMatrixType, GlobalDimNodalMatrixType > ShapeMatrices
MatrixType< GlobalDim, ShapeFunction::NPOINTS > GlobalDimNodalMatrixType
MatrixType< GlobalDim, GlobalDim > GlobalDimMatrixType
VectorType< GlobalDim > GlobalDimVectorType
VectorType< ShapeFunction::NPOINTS > NodalVectorType
RowVectorType< ShapeFunction::NPOINTS > NodalRowVectorType
virtual Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > getNodalValuesOnElement(MeshLib::Element const &element, double const t) const
Returns a matrix of values for all nodes of the given element.
Definition: Parameter.h:163
ChemistryLib::ChemicalSolverInterface *const chemical_solver_interface
std::unique_ptr< NumLib::NumericalStabilization > stabilizer
std::unique_ptr< MaterialPropertyLib::MaterialSpatialDistributionMap > media_map
IntegrationPointData(NodalRowVectorType const &N_, GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)