OGS
ComponentTransportFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <numeric>
14#include <vector>
15
35
36namespace ProcessLib
37{
38namespace ComponentTransport
39{
40template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
42{
43 IntegrationPointData(NodalRowVectorType const& N_,
44 GlobalDimNodalMatrixType const& dNdx_,
45 double const& integration_weight_)
46 : N(N_), dNdx(dNdx_), integration_weight(integration_weight_)
47 {
48 }
49
51 NodalRowVectorType const N;
52 GlobalDimNodalMatrixType const dNdx;
53 double const integration_weight;
54
55 // -1 indicates that no chemical reaction takes place in the element to
56 // which the integration point belongs.
58
59 double porosity = std::numeric_limits<double>::quiet_NaN();
60 double porosity_prev = std::numeric_limits<double>::quiet_NaN();
62};
63
67{
68public:
70
71 virtual void setChemicalSystemID(std::size_t const /*mesh_item_id*/) = 0;
72
74 std::size_t const mesh_item_id,
75 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
76 std::vector<GlobalVector*> const& x, double const t)
77 {
78 std::vector<double> local_x_vec;
79
80 auto const n_processes = x.size();
81 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
82 {
83 auto const indices =
84 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
85 assert(!indices.empty());
86 auto const local_solution = x[process_id]->get(indices);
87 local_x_vec.insert(std::end(local_x_vec),
88 std::begin(local_solution),
89 std::end(local_solution));
90 }
91 auto const local_x = MathLib::toVector(local_x_vec);
92
94 }
95
97 std::size_t const mesh_item_id,
98 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
99 std::vector<GlobalVector*> const& x, double const t, double const dt)
100 {
101 std::vector<double> local_x_vec;
102
103 auto const n_processes = x.size();
104 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
105 {
106 auto const indices =
107 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
108 assert(!indices.empty());
109 auto const local_solution = x[process_id]->get(indices);
110 local_x_vec.insert(std::end(local_x_vec),
111 std::begin(local_solution),
112 std::end(local_solution));
113 }
114 auto const local_x = MathLib::toVector(local_x_vec);
115
116 setChemicalSystemConcrete(local_x, t, dt);
117 }
118
120 std::size_t const mesh_item_id,
121 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
122 std::vector<GlobalVector*> const& x, double const t, double const dt,
123 GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, int const process_id)
124 {
125 std::vector<double> local_x_vec;
126
127 auto const n_processes = x.size();
128 for (std::size_t pcs_id = 0; pcs_id < n_processes; ++pcs_id)
129 {
130 auto const indices =
131 NumLib::getIndices(mesh_item_id, *dof_tables[pcs_id]);
132 assert(!indices.empty());
133 auto const local_solution = x[pcs_id]->get(indices);
134 local_x_vec.insert(std::end(local_x_vec),
135 std::begin(local_solution),
136 std::end(local_solution));
137 }
138 auto const local_x = MathLib::toVector(local_x_vec);
139
140 auto const indices =
141 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
142 auto const num_r_c = indices.size();
143
144 std::vector<double> local_M_data;
145 local_M_data.reserve(num_r_c * num_r_c);
146 std::vector<double> local_K_data;
147 local_K_data.reserve(num_r_c * num_r_c);
148 std::vector<double> local_b_data;
149 local_b_data.reserve(num_r_c);
150
151 assembleReactionEquationConcrete(t, dt, local_x, local_M_data,
152 local_K_data, local_b_data,
153 process_id);
154
155 auto const r_c_indices =
157 if (!local_M_data.empty())
158 {
159 auto const local_M =
160 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
161 M.add(r_c_indices, local_M);
162 }
163 if (!local_K_data.empty())
164 {
165 auto const local_K =
166 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
167 K.add(r_c_indices, local_K);
168 }
169 if (!local_b_data.empty())
170 {
171 b.add(indices, local_b_data);
172 }
173 }
174
175 virtual void postSpeciationCalculation(std::size_t const ele_id,
176 double const t, double const dt) = 0;
177
179 std::size_t const ele_id) = 0;
180
181 virtual std::vector<double> const& getIntPtDarcyVelocity(
182 const double t,
183 std::vector<GlobalVector*> const& x,
184 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
185 std::vector<double>& cache) const = 0;
186
187 virtual std::vector<double> const& getIntPtMolarFlux(
188 const double t,
189 std::vector<GlobalVector*> const& x,
190 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
191 std::vector<double>& cache) const = 0;
192
193private:
195 Eigen::VectorXd const& /*local_x*/, double const /*t*/) = 0;
196
197 virtual void setChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/,
198 double const /*t*/,
199 double const /*dt*/) = 0;
200
202 double const t, double const dt, Eigen::VectorXd const& local_x,
203 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
204 std::vector<double>& local_b_data, int const transport_process_id) = 0;
205};
206
207template <typename ShapeFunction, int GlobalDim>
209{
210 // When monolithic scheme is adopted, nodal pressure and nodal concentration
211 // are accessed by vector index.
212 static const int pressure_index = 0;
213 static const int first_concentration_index = ShapeFunction::NPOINTS;
214
215 static const int pressure_size = ShapeFunction::NPOINTS;
216 static const int concentration_size =
217 ShapeFunction::NPOINTS; // per component
218
221
223 typename ShapeMatricesType::template MatrixType<pressure_size,
226 typename ShapeMatricesType::template VectorType<pressure_size>;
227
229 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
230 using LocalVectorType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
231
234
239
240public:
242 MeshLib::Element const& element,
243 std::size_t const local_matrix_size,
244 NumLib::GenericIntegrationMethod const& integration_method,
245 bool is_axially_symmetric,
246 ComponentTransportProcessData const& process_data,
247 std::vector<std::reference_wrapper<ProcessVariable>> const&
248 transport_process_variables)
249 : _element(element),
250 _process_data(process_data),
251 _integration_method(integration_method),
252 _transport_process_variables(transport_process_variables)
253 {
254 (void)local_matrix_size;
255
256 unsigned const n_integration_points =
258 _ip_data.reserve(n_integration_points);
259
262
263 double const aperture_size = _process_data.aperture_size(0.0, pos)[0];
264
265 auto const shape_matrices =
267 GlobalDim>(element, is_axially_symmetric,
269 auto const& medium =
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 * aperture_size);
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 =
317
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 =
358
361
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_x_prev*/,
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& b =
472
473 auto const number_of_components = num_nodal_dof - 1;
474 for (int component_id = 0; component_id < number_of_components;
475 ++component_id)
476 {
477 /* Partitioned assembler matrix
478 * | pp | pc1 | pc2 | pc3 |
479 * |-----|-----|-----|-----|
480 * | c1p | c1c1| 0 | 0 |
481 * |-----|-----|-----|-----|
482 * | c2p | 0 | c2c2| 0 |
483 * |-----|-----|-----|-----|
484 * | c3p | 0 | 0 | c3c3|
485 */
486 auto concentration_index =
487 pressure_size + component_id * concentration_size;
488
489 auto KCC =
490 local_K.template block<concentration_size, concentration_size>(
491 concentration_index, concentration_index);
492 auto MCC =
493 local_M.template block<concentration_size, concentration_size>(
494 concentration_index, concentration_index);
495 auto MCp =
496 local_M.template block<concentration_size, pressure_size>(
497 concentration_index, pressure_index);
498 auto MpC =
499 local_M.template block<pressure_size, concentration_size>(
500 pressure_index, concentration_index);
501
502 auto local_C = Eigen::Map<const NodalVectorType>(
503 &local_x[concentration_index], concentration_size);
504
505 assembleBlockMatrices(b, component_id, t, dt, local_C, local_p, KCC,
506 MCC, MCp, MpC, Kpp, Mpp, Bp);
507
509 {
510 auto const stoichiometric_matrix =
513
514 assert(stoichiometric_matrix);
515
516 for (Eigen::SparseMatrix<double>::InnerIterator it(
517 *stoichiometric_matrix, component_id);
518 it;
519 ++it)
520 {
521 auto const stoichiometric_coefficient = it.value();
522 auto const coupled_component_id = it.row();
523 auto const kinetic_prefactor =
525 ->getKineticPrefactor(coupled_component_id);
526
527 auto const concentration_index =
528 pressure_size + component_id * concentration_size;
529 auto const coupled_concentration_index =
531 coupled_component_id * concentration_size;
532 auto KCmCn = local_K.template block<concentration_size,
534 concentration_index, coupled_concentration_index);
535
536 // account for the coupling between components
537 assembleKCmCn(component_id, t, dt, KCmCn,
538 stoichiometric_coefficient,
539 kinetic_prefactor);
540 }
541 }
542 }
543 }
544
546 GlobalDimVectorType const& b, int const component_id, double const t,
547 double const dt,
548 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
549 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
550 Eigen::Ref<LocalBlockMatrixType> KCC,
551 Eigen::Ref<LocalBlockMatrixType> MCC,
552 Eigen::Ref<LocalBlockMatrixType> MCp,
553 Eigen::Ref<LocalBlockMatrixType> MpC,
554 Eigen::Ref<LocalBlockMatrixType> Kpp,
555 Eigen::Ref<LocalBlockMatrixType> Mpp,
556 Eigen::Ref<LocalSegmentVectorType> Bp)
557 {
558 unsigned const n_integration_points =
560
563
565
566 // Get material properties
567 auto const& medium =
569 // Select the only valid for component transport liquid phase.
570 auto const& phase = medium.phase("AqueousLiquid");
571
572 // Assume that the component name is the same as the process variable
573 // name. Components are shifted by one because the first one is always
574 // pressure.
575 auto const& component = phase.component(
576 _transport_process_variables[component_id].get().getName());
577
578 LocalBlockMatrixType KCC_Laplacian =
579 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
580
581 std::vector<GlobalDimVectorType> ip_flux_vector;
582 double average_velocity_norm = 0.0;
584 {
585 ip_flux_vector.reserve(n_integration_points);
586 }
587
588 for (unsigned ip(0); ip < n_integration_points; ++ip)
589 {
590 pos.setIntegrationPoint(ip);
591
592 auto& ip_data = _ip_data[ip];
593 auto const& N = ip_data.N;
594 auto const& dNdx = ip_data.dNdx;
595 auto const& w = ip_data.integration_weight;
596 auto& porosity = ip_data.porosity;
597
598 double C_int_pt = 0.0;
599 double p_int_pt = 0.0;
600
601 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
602 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
603
604 vars.concentration = C_int_pt;
605 vars.phase_pressure = p_int_pt;
606
607 // update according to a particular porosity model
609 .template value<double>(vars, pos, t, dt);
610 vars.porosity = porosity;
611
612 auto const& retardation_factor =
614 .template value<double>(vars, pos, t, dt);
615
616 auto const& solute_dispersivity_transverse = medium.template value<
617 double>(
619
620 auto const& solute_dispersivity_longitudinal =
621 medium.template value<double>(
623 longitudinal_dispersivity);
624
625 // Use the fluid density model to compute the density
626 // TODO (renchao): concentration of which component as the argument
627 // for calculation of fluid density
628 auto const density =
630 .template value<double>(vars, pos, t, dt);
631
632 auto const decay_rate =
634 .template value<double>(vars, pos, t, dt);
635
636 auto const& pore_diffusion_coefficient =
637 MaterialPropertyLib::formEigenTensor<GlobalDim>(
639 .value(vars, pos, t, dt));
640
641 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
643 vars, pos, t, dt));
644
645 // Use the viscosity model to compute the viscosity
647 .template value<double>(vars, pos, t, dt);
648
649 GlobalDimMatrixType const K_over_mu = K / mu;
650 GlobalDimVectorType const velocity =
652 ? GlobalDimVectorType(-K_over_mu *
653 (dNdx * p_nodal_values - density * b))
654 : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
655
656 const double drho_dp =
658 .template dValue<double>(
660 pos, t, dt);
661
662 const double drho_dC =
664 .template dValue<double>(
666 t, dt);
667
668 GlobalDimMatrixType const hydrodynamic_dispersion =
671 pore_diffusion_coefficient, velocity, porosity,
672 solute_dispersivity_transverse,
673 solute_dispersivity_longitudinal);
674
675 const double R_times_phi(retardation_factor * porosity);
676 GlobalDimVectorType const mass_density_flow = velocity * density;
677 auto const N_t_N = (N.transpose() * N).eval();
678
680 {
681 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
682 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
683 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
684 }
685 else
686 {
687 ip_flux_vector.emplace_back(mass_density_flow);
688 average_velocity_norm += velocity.norm();
689 }
690 MCC.noalias() += N_t_N * (R_times_phi * density * w);
691 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
692 KCC_Laplacian.noalias() +=
693 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
694
695 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
696
697 // Calculate Mpp, Kpp, and bp in the first loop over components
698 if (component_id == 0)
699 {
700 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
701 Kpp.noalias() +=
702 dNdx.transpose() * K_over_mu * dNdx * (density * w);
703
705 {
706 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
707 (density * density * w);
708 }
709 }
710 }
711
713 {
716 _ip_data,
717 ip_flux_vector,
718 average_velocity_norm /
719 static_cast<double>(n_integration_points),
720 KCC_Laplacian);
721 }
722
723 KCC.noalias() += KCC_Laplacian;
724 }
725
726 void assembleKCmCn(int const component_id, double const t, double const dt,
727 Eigen::Ref<LocalBlockMatrixType> KCmCn,
728 double const stoichiometric_coefficient,
729 double const kinetic_prefactor)
730 {
731 unsigned const n_integration_points =
733
736
738
739 auto const& medium =
741 auto const& phase = medium.phase("AqueousLiquid");
742 auto const& component = phase.component(
743 _transport_process_variables[component_id].get().getName());
744
745 for (unsigned ip(0); ip < n_integration_points; ++ip)
746 {
747 auto& ip_data = _ip_data[ip];
748 auto const& N = ip_data.N;
749 auto const& w = ip_data.integration_weight;
750 auto& porosity = ip_data.porosity;
751
752 auto const retardation_factor =
754 .template value<double>(vars, pos, t, dt);
755
757 .template value<double>(vars, pos, t, dt);
758
759 auto const density =
761 .template value<double>(vars, pos, t, dt);
762
763 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
764 kinetic_prefactor * retardation_factor *
765 porosity * density * N;
766 }
767 }
768
769 void assembleForStaggeredScheme(double const t, double const dt,
770 Eigen::VectorXd const& local_x,
771 Eigen::VectorXd const& local_x_prev,
772 int const process_id,
773 std::vector<double>& local_M_data,
774 std::vector<double>& local_K_data,
775 std::vector<double>& local_b_data) override
776 {
777 if (process_id == _process_data.hydraulic_process_id)
778 {
779 assembleHydraulicEquation(t, dt, local_x, local_x_prev,
780 local_M_data, local_K_data, local_b_data);
781 }
782 else
783 {
784 // Go for assembling in an order of transport process id.
785 assembleComponentTransportEquation(t, dt, local_x, local_x_prev,
786 local_M_data, local_K_data,
787 local_b_data, process_id);
788 }
789 }
790
791 void assembleHydraulicEquation(double const t,
792 double const dt,
793 Eigen::VectorXd const& local_x,
794 Eigen::VectorXd const& local_x_prev,
795 std::vector<double>& local_M_data,
796 std::vector<double>& local_K_data,
797 std::vector<double>& local_b_data)
798 {
799 auto const local_p =
800 local_x.template segment<pressure_size>(pressure_index);
801 auto const local_C = local_x.template segment<concentration_size>(
803 auto const local_C_prev =
804 local_x_prev.segment<concentration_size>(first_concentration_index);
805
806 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
807 local_M_data, pressure_size, pressure_size);
808 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
809 local_K_data, pressure_size, pressure_size);
810 auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
811 local_b_data, pressure_size);
812
813 unsigned const n_integration_points =
815
818
819 auto const& b =
822
823 auto const& medium =
825 auto const& phase = medium.phase("AqueousLiquid");
826
829
830 for (unsigned ip(0); ip < n_integration_points; ++ip)
831 {
832 pos.setIntegrationPoint(ip);
833
834 auto& ip_data = _ip_data[ip];
835 auto const& N = ip_data.N;
836 auto const& dNdx = ip_data.dNdx;
837 auto const& w = ip_data.integration_weight;
838 auto& porosity = ip_data.porosity;
839 auto const& porosity_prev = ip_data.porosity_prev;
840
841 double C_int_pt = 0.0;
842 double p_int_pt = 0.0;
843
844 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
845 NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
846
847 vars.concentration = C_int_pt;
848 vars.phase_pressure = p_int_pt;
849
850 // porosity
851 {
852 vars_prev.porosity = porosity_prev;
853
854 porosity =
856 ? porosity_prev
858 .template value<double>(vars, vars_prev, pos, t,
859 dt);
860
861 vars.porosity = porosity;
862 }
863
864 // Use the fluid density model to compute the density
865 // TODO: Concentration of which component as one of arguments for
866 // calculation of fluid density
867 auto const density =
869 .template value<double>(vars, pos, t, dt);
870
871 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
873 vars, pos, t, dt));
874
875 // Use the viscosity model to compute the viscosity
877 .template value<double>(vars, pos, t, dt);
878
879 GlobalDimMatrixType const K_over_mu = K / mu;
880
881 const double drho_dp =
883 .template dValue<double>(
885 pos, t, dt);
886 const double drho_dC =
888 .template dValue<double>(
890 t, dt);
891
892 // matrix assembly
893 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
894 local_K.noalias() +=
895 w * dNdx.transpose() * density * K_over_mu * dNdx;
896
898 {
899 local_b.noalias() +=
900 w * density * density * dNdx.transpose() * K_over_mu * b;
901 }
902
903 // coupling term
904 {
905 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
906
907 local_b.noalias() -=
908 w * N.transpose() * porosity * drho_dC * C_dot;
909 }
910 }
911 }
912
914 double const t, double const dt, Eigen::VectorXd const& local_x,
915 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_M_data,
916 std::vector<double>& local_K_data,
917 std::vector<double>& /*local_b_data*/, int const transport_process_id)
918 {
919 auto const local_p =
920 local_x.template segment<pressure_size>(pressure_index);
921 auto const local_C = local_x.template segment<concentration_size>(
923 (transport_process_id - 1) * concentration_size);
924 auto const local_p_prev =
925 local_x_prev.segment<pressure_size>(pressure_index);
926
927 NodalVectorType local_T;
929 {
930 local_T =
932 }
933
934 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
936 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
938
939 LocalBlockMatrixType KCC_Laplacian =
940 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
941
942 unsigned const n_integration_points =
944
945 std::vector<GlobalDimVectorType> ip_flux_vector;
946 double average_velocity_norm = 0.0;
948 {
949 ip_flux_vector.reserve(n_integration_points);
950 }
951
954
955 auto const& b =
958
961
962 auto const& medium =
964 auto const& phase = medium.phase("AqueousLiquid");
965 // Hydraulic process id is 0 and thus transport process id starts
966 // from 1.
967 auto const component_id = transport_process_id - 1;
968 auto const& component = phase.component(
969 _transport_process_variables[component_id].get().getName());
970
971 for (unsigned ip(0); ip < n_integration_points; ++ip)
972 {
973 pos.setIntegrationPoint(ip);
974
975 auto& ip_data = _ip_data[ip];
976 auto const& N = ip_data.N;
977 auto const& dNdx = ip_data.dNdx;
978 auto const& w = ip_data.integration_weight;
979 auto& porosity = ip_data.porosity;
980 auto const& porosity_prev = ip_data.porosity_prev;
981
982 double C_int_pt = 0.0;
983 double p_int_pt = 0.0;
984
985 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
986 NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
987
988 vars.concentration = C_int_pt;
989 vars.phase_pressure = p_int_pt;
990
992 {
993 vars.temperature = N.dot(local_T);
994 }
995
996 // porosity
997 {
998 vars_prev.porosity = porosity_prev;
999
1000 porosity =
1002 ? porosity_prev
1004 .template value<double>(vars, vars_prev, pos, t,
1005 dt);
1006
1007 vars.porosity = porosity;
1008 }
1009
1010 auto const& retardation_factor =
1012 .template value<double>(vars, pos, t, dt);
1013
1014 auto const& solute_dispersivity_transverse = medium.template value<
1015 double>(
1017 auto const& solute_dispersivity_longitudinal =
1018 medium.template value<double>(
1020 longitudinal_dispersivity);
1021
1022 // Use the fluid density model to compute the density
1023 auto const density =
1025 .template value<double>(vars, pos, t, dt);
1026 auto const decay_rate =
1028 .template value<double>(vars, pos, t, dt);
1029
1030 auto const& pore_diffusion_coefficient =
1031 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1033 .value(vars, pos, t, dt));
1034
1035 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1037 vars, pos, t, dt));
1038 // Use the viscosity model to compute the viscosity
1040 .template value<double>(vars, pos, t, dt);
1041
1042 GlobalDimMatrixType const K_over_mu = K / mu;
1043 GlobalDimVectorType const velocity =
1045 ? GlobalDimVectorType(-K_over_mu *
1046 (dNdx * local_p - density * b))
1047 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1048
1049 GlobalDimMatrixType const hydrodynamic_dispersion =
1052 pore_diffusion_coefficient, velocity, porosity,
1053 solute_dispersivity_transverse,
1054 solute_dispersivity_longitudinal);
1055
1056 double const R_times_phi = retardation_factor * porosity;
1057 auto const N_t_N = (N.transpose() * N).eval();
1058
1060 {
1061 const double drho_dC =
1063 .template dValue<double>(
1065 pos, t, dt);
1066 local_M.noalias() +=
1067 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1068 }
1069
1070 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1071
1072 // coupling term
1074 {
1075 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1076
1077 const double drho_dp =
1079 .template dValue<double>(
1081 pos, t, dt);
1082
1083 local_K.noalias() +=
1084 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1085 dNdx.transpose() * velocity * N * (density * w);
1086 }
1087 else
1088 {
1089 ip_flux_vector.emplace_back(velocity * density);
1090 average_velocity_norm += velocity.norm();
1091 }
1092 local_K.noalias() +=
1093 N_t_N * (decay_rate * R_times_phi * density * w);
1094
1095 KCC_Laplacian.noalias() += dNdx.transpose() *
1096 hydrodynamic_dispersion * dNdx *
1097 (density * w);
1098 }
1099
1101 {
1104 _ip_data,
1105 ip_flux_vector,
1106 average_velocity_norm /
1107 static_cast<double>(n_integration_points),
1108 KCC_Laplacian);
1109 }
1110 local_K.noalias() += KCC_Laplacian;
1111 }
1112
1114 double const t, double const dt, Eigen::VectorXd const& local_x,
1115 Eigen::VectorXd const& local_x_prev, int const process_id,
1116 std::vector<double>& /*local_M_data*/,
1117 std::vector<double>& /*local_K_data*/,
1118 std::vector<double>& local_b_data,
1119 std::vector<double>& local_Jac_data) override
1120 {
1121 if (process_id == _process_data.hydraulic_process_id)
1122 {
1123 assembleWithJacobianHydraulicEquation(t, dt, local_x, local_x_prev,
1124 local_b_data, local_Jac_data);
1125 }
1126 else
1127 {
1128 int const component_id = process_id - 1;
1130 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1131 component_id);
1132 }
1133 }
1134
1136 double const t, double const dt, Eigen::VectorXd const& local_x,
1137 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1138 std::vector<double>& local_Jac_data)
1139 {
1140 auto const p = local_x.template segment<pressure_size>(pressure_index);
1141 auto const c = local_x.template segment<concentration_size>(
1143
1144 auto const p_prev = local_x_prev.segment<pressure_size>(pressure_index);
1145 auto const c_prev =
1146 local_x_prev.segment<concentration_size>(first_concentration_index);
1147
1148 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1149 local_Jac_data, pressure_size, pressure_size);
1150 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1151 local_b_data, pressure_size);
1152
1153 unsigned const n_integration_points =
1155
1158 auto const& b =
1161
1162 auto const& medium =
1164 auto const& phase = medium.phase("AqueousLiquid");
1165
1168
1169 for (unsigned ip(0); ip < n_integration_points; ++ip)
1170 {
1171 pos.setIntegrationPoint(ip);
1172
1173 auto& ip_data = _ip_data[ip];
1174 auto const& N = ip_data.N;
1175 auto const& dNdx = ip_data.dNdx;
1176 auto const& w = ip_data.integration_weight;
1177 auto& phi = ip_data.porosity;
1178 auto const& phi_prev = ip_data.porosity_prev;
1179
1180 double const p_ip = N.dot(p);
1181 double const c_ip = N.dot(c);
1182
1183 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1184
1185 vars.phase_pressure = p_ip;
1186 vars.concentration = c_ip;
1187
1188 // porosity
1189 {
1190 vars_prev.porosity = phi_prev;
1191
1193 ? phi_prev
1195 .template value<double>(vars, vars_prev, pos, t,
1196 dt);
1197
1198 vars.porosity = phi;
1199 }
1200
1201 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1202 .template value<double>(vars, pos, t, dt);
1203
1204 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1206 vars, pos, t, dt));
1207
1209 .template value<double>(vars, pos, t, dt);
1210
1211 auto const drho_dp =
1213 .template dValue<double>(
1215 pos, t, dt);
1216 auto const drho_dc =
1218 .template dValue<double>(
1220 t, dt);
1221
1222 // matrix assembly
1223 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1224 w * dNdx.transpose() * rho * k / mu * dNdx;
1225
1226 local_rhs.noalias() -=
1227 w * N.transpose() * phi *
1228 (drho_dp * N * p_prev + drho_dc * cdot_ip) +
1229 w * rho * dNdx.transpose() * k / mu * dNdx * p;
1230
1232 {
1233 local_rhs.noalias() +=
1234 w * rho * dNdx.transpose() * k / mu * rho * b;
1235 }
1236 }
1237 }
1238
1240 double const t, double const dt, Eigen::VectorXd const& local_x,
1241 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1242 std::vector<double>& local_Jac_data, int const component_id)
1243 {
1244 auto const concentration_index =
1246
1247 auto const p = local_x.template segment<pressure_size>(pressure_index);
1248 auto const c =
1249 local_x.template segment<concentration_size>(concentration_index);
1250 auto const c_prev =
1251 local_x_prev.segment<concentration_size>(concentration_index);
1252
1255 {
1257 }
1258
1259 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1260 local_Jac_data, concentration_size, concentration_size);
1261 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1262 local_b_data, concentration_size);
1263
1264 LocalBlockMatrixType KCC_Laplacian =
1265 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1266
1267 unsigned const n_integration_points =
1269
1270 std::vector<GlobalDimVectorType> ip_flux_vector;
1271 double average_velocity_norm = 0.0;
1272 ip_flux_vector.reserve(n_integration_points);
1273
1276
1277 auto const& b =
1280
1283
1284 auto const& medium =
1286 auto const& phase = medium.phase("AqueousLiquid");
1287 auto const& component = phase.component(
1288 _transport_process_variables[component_id].get().getName());
1289
1290 for (unsigned ip(0); ip < n_integration_points; ++ip)
1291 {
1292 pos.setIntegrationPoint(ip);
1293
1294 auto& ip_data = _ip_data[ip];
1295 auto const& N = ip_data.N;
1296 auto const& dNdx = ip_data.dNdx;
1297 auto const& w = ip_data.integration_weight;
1298 auto& phi = ip_data.porosity;
1299 auto const& phi_prev = ip_data.porosity_prev;
1300
1301 double const p_ip = N.dot(p);
1302 double const c_ip = N.dot(c);
1303
1304 vars.phase_pressure = p_ip;
1305 vars.concentration = c_ip;
1306
1308 {
1309 vars.temperature = N.dot(T);
1310 }
1311
1312 // porosity
1313 {
1314 vars_prev.porosity = phi_prev;
1315
1317 ? phi_prev
1319 .template value<double>(vars, vars_prev, pos, t,
1320 dt);
1321
1322 vars.porosity = phi;
1323 }
1324
1325 auto const R =
1327 .template value<double>(vars, pos, t, dt);
1328
1329 auto const alpha_T = medium.template value<double>(
1331 auto const alpha_L = medium.template value<double>(
1333
1334 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1335 .template value<double>(vars, pos, t, dt);
1336 // first-order decay constant
1337 auto const alpha =
1339 .template value<double>(vars, pos, t, dt);
1340
1341 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1343 .value(vars, pos, t, dt));
1344
1345 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1347 vars, pos, t, dt));
1349 .template value<double>(vars, pos, t, dt);
1350 // Darcy flux
1351 GlobalDimVectorType const q =
1353 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1354 : GlobalDimVectorType(-k / mu * dNdx * p);
1355
1357 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1358 alpha_L);
1359
1360 // matrix assembly
1361 local_Jac.noalias() +=
1362 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1363
1364 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1365
1366 auto const cdot = (c - c_prev) / dt;
1367 local_rhs.noalias() -=
1368 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1369
1370 ip_flux_vector.emplace_back(q * rho);
1371 average_velocity_norm += q.norm();
1372 }
1373
1375 _process_data.stabilizer, _ip_data, ip_flux_vector,
1376 average_velocity_norm / static_cast<double>(n_integration_points),
1377 KCC_Laplacian);
1378
1379 local_rhs.noalias() -= KCC_Laplacian * c;
1380
1381 local_Jac.noalias() += KCC_Laplacian;
1382 }
1383
1385 double const t, double const dt, Eigen::VectorXd const& local_x,
1386 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1387 std::vector<double>& local_b_data,
1388 int const transport_process_id) override
1389 {
1390 auto const local_C = local_x.template segment<concentration_size>(
1392 (transport_process_id - 1) * concentration_size);
1393
1394 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1395 local_M_data, concentration_size, concentration_size);
1396 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1397 local_K_data, concentration_size, concentration_size);
1398 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1399 local_b_data, concentration_size);
1400
1401 unsigned const n_integration_points =
1403
1406
1409
1410 auto const& medium =
1412 auto const component_id = transport_process_id - 1;
1413 for (unsigned ip(0); ip < n_integration_points; ++ip)
1414 {
1415 pos.setIntegrationPoint(ip);
1416
1417 auto& ip_data = _ip_data[ip];
1418 auto const& N = ip_data.N;
1419 auto const w = ip_data.integration_weight;
1420 auto& porosity = ip_data.porosity;
1421 auto const& porosity_prev = ip_data.porosity_prev;
1422 auto const chemical_system_id = ip_data.chemical_system_id;
1423
1424 double C_int_pt = 0.0;
1425 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
1426
1427 vars.concentration = C_int_pt;
1428
1429 auto const porosity_dot = (porosity - porosity_prev) / dt;
1430
1431 // porosity
1432 {
1433 vars_prev.porosity = porosity_prev;
1434
1435 porosity =
1437 ? porosity_prev
1439 .template value<double>(vars, vars_prev, pos, t,
1440 dt);
1441 }
1442
1443 local_M.noalias() += w * N.transpose() * porosity * N;
1444
1445 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1446
1447 if (chemical_system_id == -1)
1448 {
1449 continue;
1450 }
1451
1452 auto const C_post_int_pt =
1454 component_id, chemical_system_id);
1455
1456 local_b.noalias() +=
1457 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1458 }
1459 }
1460
1461 std::vector<double> const& getIntPtDarcyVelocity(
1462 const double t,
1463 std::vector<GlobalVector*> const& x,
1464 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
1465 std::vector<double>& cache) const override
1466 {
1467 assert(x.size() == dof_table.size());
1468
1469 auto const n_processes = x.size();
1470 std::vector<std::vector<double>> local_x;
1471 local_x.reserve(n_processes);
1472
1473 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1474 {
1475 auto const indices =
1476 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1477 assert(!indices.empty());
1478 local_x.push_back(x[process_id]->get(indices));
1479 }
1480
1481 // only one process, must be monolithic.
1482 if (n_processes == 1)
1483 {
1484 auto const local_p = Eigen::Map<const NodalVectorType>(
1485 &local_x[0][pressure_index], pressure_size);
1486 auto const local_C = Eigen::Map<const NodalVectorType>(
1488 return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1489 }
1490
1491 // multiple processes, must be staggered.
1492 {
1493 constexpr int pressure_process_id = 0;
1494 constexpr int concentration_process_id = 1;
1495 auto const local_p = Eigen::Map<const NodalVectorType>(
1496 &local_x[pressure_process_id][0], pressure_size);
1497 auto const local_C = Eigen::Map<const NodalVectorType>(
1498 &local_x[concentration_process_id][0], concentration_size);
1499 return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1500 }
1501 }
1502
1503 std::vector<double> const& calculateIntPtDarcyVelocity(
1504 const double t,
1505 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
1506 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
1507 std::vector<double>& cache) const
1508 {
1509 auto const n_integration_points =
1511
1512 cache.clear();
1513 auto cache_mat = MathLib::createZeroedMatrix<
1514 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1515 cache, GlobalDim, n_integration_points);
1516
1519
1520 auto const& b =
1523
1525
1526 auto const& medium =
1528 auto const& phase = medium.phase("AqueousLiquid");
1529
1530 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1531 {
1532 auto const& ip_data = _ip_data[ip];
1533 auto const& N = ip_data.N;
1534 auto const& dNdx = ip_data.dNdx;
1535 auto const& porosity = ip_data.porosity;
1536
1537 pos.setIntegrationPoint(ip);
1538
1539 double C_int_pt = 0.0;
1540 double p_int_pt = 0.0;
1541
1542 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1543 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1544
1545 vars.concentration = C_int_pt;
1546 vars.phase_pressure = p_int_pt;
1547 vars.porosity = porosity;
1548
1549 // TODO (naumov) Temporary value not used by current material
1550 // models. Need extension of secondary variables interface.
1551 double const dt = std::numeric_limits<double>::quiet_NaN();
1552 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1554 vars, pos, t, dt));
1556 .template value<double>(vars, pos, t, dt);
1557 GlobalDimMatrixType const K_over_mu = K / mu;
1558
1559 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1561 {
1562 auto const rho_w =
1564 .template value<double>(vars, pos, t, dt);
1565 // here it is assumed that the vector b is directed 'downwards'
1566 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1567 }
1568 }
1569
1570 return cache;
1571 }
1572
1573 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
1574 const unsigned integration_point) const override
1575 {
1576 auto const& N = _ip_data[integration_point].N;
1577
1578 // assumes N is stored contiguously in memory
1579 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1580 }
1581
1582 Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
1583 double const t,
1584 std::vector<double> const& local_x) const override
1585 {
1586 auto const local_p = Eigen::Map<const NodalVectorType>(
1587 &local_x[pressure_index], pressure_size);
1588 auto const local_C = Eigen::Map<const NodalVectorType>(
1590
1591 // Eval shape matrices at given point
1592 // Note: Axial symmetry is set to false here, because we only need dNdx
1593 // here, which is not affected by axial symmetry.
1594 auto const shape_matrices =
1596 GlobalDim>(
1597 _element, false /*is_axially_symmetric*/,
1598 std::array{pnt_local_coords})[0];
1599
1602 auto const& b =
1605
1607
1608 auto const& medium =
1610 auto const& phase = medium.phase("AqueousLiquid");
1611
1612 // local_x contains the local concentration and pressure values
1613 double c_int_pt;
1614 NumLib::shapeFunctionInterpolate(local_C, shape_matrices.N, c_int_pt);
1615 vars.concentration = c_int_pt;
1616
1617 double p_int_pt;
1618 NumLib::shapeFunctionInterpolate(local_p, shape_matrices.N, p_int_pt);
1619 vars.phase_pressure = p_int_pt;
1620
1621 // TODO (naumov) Temporary value not used by current material models.
1622 // Need extension of secondary variables interface.
1623 double const dt = std::numeric_limits<double>::quiet_NaN();
1624 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1626 vars, pos, t, dt));
1627
1629 .template value<double>(vars, pos, t, dt);
1630 GlobalDimMatrixType const K_over_mu = K / mu;
1631
1632 GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
1633 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
1634 .template value<double>(vars, pos, t, dt);
1636 {
1637 q += K_over_mu * rho_w * b;
1638 }
1639 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1640 flux.head<GlobalDim>() = rho_w * q;
1641 return flux;
1642 }
1643
1645 double const t,
1646 double const /*dt*/,
1647 Eigen::VectorXd const& local_x,
1648 Eigen::VectorXd const& /*local_x_prev*/) override
1649 {
1650 auto const local_p =
1651 local_x.template segment<pressure_size>(pressure_index);
1652 auto const local_C = local_x.template segment<concentration_size>(
1654
1655 std::vector<double> ele_velocity;
1656 calculateIntPtDarcyVelocity(t, local_p, local_C, ele_velocity);
1657
1658 auto const n_integration_points =
1660 auto const ele_velocity_mat =
1661 MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
1662
1663 auto const ele_id = _element.getID();
1664 Eigen::Map<LocalVectorType>(
1665 &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
1666 GlobalDim) =
1667 ele_velocity_mat.rowwise().sum() / n_integration_points;
1668 }
1669
1671 std::size_t const ele_id) override
1672 {
1673 auto const n_integration_points =
1675
1677 {
1678 auto const& medium = *_process_data.media_map.getMedium(ele_id);
1679
1680 for (auto& ip_data : _ip_data)
1681 {
1682 ip_data.porosity = ip_data.porosity_prev;
1683
1685 ->updatePorosityPostReaction(ip_data.chemical_system_id,
1686 medium, ip_data.porosity);
1687 }
1688
1690 std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
1691 [](double const s, auto const& ip)
1692 { return s + ip.porosity; }) /
1693 n_integration_points;
1694 }
1695
1696 std::vector<GlobalIndexType> chemical_system_indices;
1697 chemical_system_indices.reserve(n_integration_points);
1698 std::transform(_ip_data.begin(), _ip_data.end(),
1699 std::back_inserter(chemical_system_indices),
1700 [](auto const& ip_data)
1701 { return ip_data.chemical_system_id; });
1702
1704 ele_id, chemical_system_indices);
1705 }
1706
1707 std::vector<double> const& getIntPtMolarFlux(
1708 const double t,
1709 std::vector<GlobalVector*> const& x,
1710 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
1711 std::vector<double>& cache) const override
1712 {
1713 std::vector<double> local_x_vec;
1714
1715 auto const n_processes = x.size();
1716 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1717 {
1718 auto const indices =
1719 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
1720 assert(!indices.empty());
1721 auto const local_solution = x[process_id]->get(indices);
1722 local_x_vec.insert(std::end(local_x_vec),
1723 std::begin(local_solution),
1724 std::end(local_solution));
1725 }
1726 auto const local_x = MathLib::toVector(local_x_vec);
1727
1728 auto const p = local_x.template segment<pressure_size>(pressure_index);
1729 auto const c = local_x.template segment<concentration_size>(
1731
1732 auto const n_integration_points =
1734
1735 cache.clear();
1736 auto cache_mat = MathLib::createZeroedMatrix<
1737 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1738 cache, GlobalDim, n_integration_points);
1739
1742
1743 auto const& b =
1746
1748
1749 auto const& medium =
1751 auto const& phase = medium.phase("AqueousLiquid");
1752
1753 int const component_id = 0;
1754 auto const& component = phase.component(
1755 _transport_process_variables[component_id].get().getName());
1756
1757 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1758 {
1759 auto const& ip_data = _ip_data[ip];
1760 auto const& N = ip_data.N;
1761 auto const& dNdx = ip_data.dNdx;
1762 auto const& phi = ip_data.porosity;
1763
1764 pos.setIntegrationPoint(ip);
1765
1766 double const p_ip = N.dot(p);
1767 double const c_ip = N.dot(c);
1768
1769 vars.concentration = c_ip;
1770 vars.phase_pressure = p_ip;
1771 vars.porosity = phi;
1772
1773 double const dt = std::numeric_limits<double>::quiet_NaN();
1774
1775 auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1777 vars, pos, t, dt));
1779 .template value<double>(vars, pos, t, dt);
1780 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1781 .template value<double>(vars, pos, t, dt);
1782
1783 // Darcy flux
1784 GlobalDimVectorType const q =
1786 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1787 : GlobalDimVectorType(-k / mu * dNdx * p);
1788
1789 auto const alpha_T = medium.template value<double>(
1791 auto const alpha_L = medium.template value<double>(
1793 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1795 .value(vars, pos, t, dt));
1796
1797 // Hydrodynamic dispersion
1799 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1800 alpha_L);
1801
1802 cache_mat.col(ip).noalias() = q * c_ip - phi * D * dNdx * c;
1803 }
1804
1805 return cache;
1806 }
1807
1808 void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
1809 Eigen::VectorXd const& /*local_x_prev*/,
1810 double const /*t*/, double const /*dt*/,
1811 bool const /*use_monolithic_scheme*/,
1812 int const /*process_id*/) override
1813 {
1814 unsigned const n_integration_points =
1816
1817 for (unsigned ip = 0; ip < n_integration_points; ip++)
1818 {
1819 _ip_data[ip].pushBackState();
1820 }
1821 }
1822
1823private:
1826
1828 std::vector<std::reference_wrapper<ProcessVariable>> const
1830
1831 std::vector<
1833 Eigen::aligned_allocator<
1836};
1837
1838} // namespace ComponentTransport
1839} // namespace ProcessLib
GlobalMatrix::IndexType GlobalIndexType
std::string getName(std::string const &line)
Returns the name/title from the "Zone"-description.
virtual void updatePorosityPostReaction(GlobalIndexType const &, MaterialPropertyLib::Medium const &, double &)
virtual double getKineticPrefactor(std::size_t reaction_id) const
virtual void computeSecondaryVariable(std::size_t const, std::vector< GlobalIndexType > 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 Eigen::SparseMatrix< double > const * getStoichiometricMatrix() const
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)
Phase const & phase(std::size_t index) const
Definition: Medium.cpp:33
Component const & component(std::size_t const &index) const
Definition: Phase.cpp:33
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 std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
virtual void 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 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
typename ShapeMatricesType::GlobalDimVectorType GlobalDimVectorType
void assembleForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
typename ShapeMatricesType::ShapeMatrices ShapeMatrices
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > LocalMatrixType
void assembleWithJacobianHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data)
void postTimestepConcrete(Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, bool const, int const) override
void postSpeciationCalculation(std::size_t const ele_id, double const t, double const dt) override
NumLib::GenericIntegrationMethod const & _integration_method
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 & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< double > &cache) const override
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
void assembleBlockMatrices(GlobalDimVectorType const &b, 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)
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
void assembleWithJacobianForStaggeredScheme(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, int const process_id, std::vector< double > &, std::vector< double > &, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data) override
ComponentTransportProcessData const & _process_data
Eigen::Map< const Eigen::RowVectorXd > getShapeMatrix(const unsigned integration_point) const override
Provides the shape matrix at the given integration point.
ShapeMatrixPolicyType< ShapeFunction, GlobalDim > ShapeMatricesType
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
void assembleHydraulicEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data)
typename ShapeMatricesType::NodalRowVectorType NodalRowVectorType
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
std::vector< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType >, Eigen::aligned_allocator< IntegrationPointData< NodalRowVectorType, GlobalDimNodalMatrixType > > > _ip_data
typename ShapeMatricesType::template MatrixType< pressure_size, pressure_size > LocalBlockMatrixType
typename ShapeMatricesType::template VectorType< pressure_size > LocalSegmentVectorType
void assembleComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &, int const transport_process_id)
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
void assembleWithJacobianComponentTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &local_x_prev, std::vector< double > &local_b_data, std::vector< double > &local_Jac_data, int const component_id)
Eigen::Matrix< double, Eigen::Dynamic, 1 > LocalVectorType
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
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)
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:60
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
Definition: PropertyType.h:110
@ retardation_factor
specify retardation factor used in component transport process.
Definition: PropertyType.h:84
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
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
Definition: EigenMapTools.h:72
void assembleAdvectionMatrix(IPData const &ip_data_vector, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Definition: Interpolation.h:27
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)
std::vector< typename ShapeMatricesType::ShapeMatrices, Eigen::aligned_allocator< typename ShapeMatricesType::ShapeMatrices > > computeShapeMatrices(MeshLib::Element const &e, bool const is_axially_symmetric, PointContainer const &points)
Eigen::MatrixXd computeHydrodynamicDispersion(NumericalStabilization const &stabilizer, std::size_t const element_id, Eigen::MatrixXd const &pore_diffusion_coefficient, Eigen::VectorXd const &velocity, double const porosity, double const solute_dispersivity_transverse, double const solute_dispersivity_longitudinal)
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:164
std::vector< Eigen::VectorXd > const projected_specific_body_force_vectors
Projected specific body force vector: R * R^T * b.
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
ChemistryLib::ChemicalSolverInterface *const chemical_solver_interface
IntegrationPointData(NodalRowVectorType const &N_, GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)