OGS
ComponentTransportFEM.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <numeric>
14#include <vector>
15
36
37namespace ProcessLib
38{
39namespace ComponentTransport
40{
41template <typename GlobalDimNodalMatrixType>
43{
44 IntegrationPointData(GlobalDimNodalMatrixType const& dNdx_,
45 double const& integration_weight_)
46 : dNdx(dNdx_), integration_weight(integration_weight_)
47 {
48 }
49
51 GlobalDimNodalMatrixType const dNdx;
52 double const integration_weight;
53
54 // -1 indicates that no chemical reaction takes place in the element to
55 // which the integration point belongs.
57
58 double porosity = std::numeric_limits<double>::quiet_NaN();
59 double porosity_prev = std::numeric_limits<double>::quiet_NaN();
61};
62
66{
67public:
69
70 virtual void setChemicalSystemID(std::size_t const /*mesh_item_id*/) = 0;
71
73 std::size_t const mesh_item_id,
74 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
75 std::vector<GlobalVector*> const& x, double const t)
76 {
77 std::vector<double> local_x_vec;
78
79 auto const n_processes = x.size();
80 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
81 {
82 auto const indices =
83 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
84 assert(!indices.empty());
85 auto const local_solution = x[process_id]->get(indices);
86 local_x_vec.insert(std::end(local_x_vec),
87 std::begin(local_solution),
88 std::end(local_solution));
89 }
90 auto const local_x = MathLib::toVector(local_x_vec);
91
93 }
94
96 std::size_t const mesh_item_id,
97 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
98 std::vector<GlobalVector*> const& x, double const t, double const dt)
99 {
100 std::vector<double> local_x_vec;
101
102 auto const n_processes = x.size();
103 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
104 {
105 auto const indices =
106 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
107 assert(!indices.empty());
108 auto const local_solution = x[process_id]->get(indices);
109 local_x_vec.insert(std::end(local_x_vec),
110 std::begin(local_solution),
111 std::end(local_solution));
112 }
113 auto const local_x = MathLib::toVector(local_x_vec);
114
115 setChemicalSystemConcrete(local_x, t, dt);
116 }
117
119 std::size_t const mesh_item_id,
120 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
121 std::vector<GlobalVector*> const& x, double const t, double const dt,
122 GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, int const process_id)
123 {
124 std::vector<double> local_x_vec;
125
126 auto const n_processes = x.size();
127 for (std::size_t pcs_id = 0; pcs_id < n_processes; ++pcs_id)
128 {
129 auto const indices =
130 NumLib::getIndices(mesh_item_id, *dof_tables[pcs_id]);
131 assert(!indices.empty());
132 auto const local_solution = x[pcs_id]->get(indices);
133 local_x_vec.insert(std::end(local_x_vec),
134 std::begin(local_solution),
135 std::end(local_solution));
136 }
137 auto const local_x = MathLib::toVector(local_x_vec);
138
139 auto const indices =
140 NumLib::getIndices(mesh_item_id, *dof_tables[process_id]);
141 auto const num_r_c = indices.size();
142
143 std::vector<double> local_M_data;
144 local_M_data.reserve(num_r_c * num_r_c);
145 std::vector<double> local_K_data;
146 local_K_data.reserve(num_r_c * num_r_c);
147 std::vector<double> local_b_data;
148 local_b_data.reserve(num_r_c);
149
150 assembleReactionEquationConcrete(t, dt, local_x, local_M_data,
151 local_K_data, local_b_data,
152 process_id);
153
154 auto const r_c_indices =
156 if (!local_M_data.empty())
157 {
158 auto const local_M =
159 MathLib::toMatrix(local_M_data, num_r_c, num_r_c);
160 M.add(r_c_indices, local_M);
161 }
162 if (!local_K_data.empty())
163 {
164 auto const local_K =
165 MathLib::toMatrix(local_K_data, num_r_c, num_r_c);
166 K.add(r_c_indices, local_K);
167 }
168 if (!local_b_data.empty())
169 {
170 b.add(indices, local_b_data);
171 }
172 }
173
174 virtual void postSpeciationCalculation(std::size_t const ele_id,
175 double const t, double const dt) = 0;
176
178 std::size_t const ele_id) = 0;
179
180 virtual std::vector<double> const& getIntPtDarcyVelocity(
181 const double t,
182 std::vector<GlobalVector*> const& x,
183 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
184 std::vector<double>& cache) const = 0;
185
186 virtual std::vector<double> const& getIntPtMolarFlux(
187 const double t, std::vector<GlobalVector*> const& x,
188 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
189 std::vector<double>& cache, int const component_id) const = 0;
190
191private:
193 Eigen::VectorXd const& /*local_x*/, double const /*t*/) = 0;
194
195 virtual void setChemicalSystemConcrete(Eigen::VectorXd const& /*local_x*/,
196 double const /*t*/,
197 double const /*dt*/) = 0;
198
200 double const t, double const dt, Eigen::VectorXd const& local_x,
201 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
202 std::vector<double>& local_b_data, int const transport_process_id) = 0;
203};
204
205template <typename ShapeFunction, int GlobalDim>
207{
208 // When monolithic scheme is adopted, nodal pressure and nodal concentration
209 // are accessed by vector index.
210 static const int pressure_index = 0;
211 const int temperature_index = -1;
213
214 static const int pressure_size = ShapeFunction::NPOINTS;
215 static const int temperature_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 : temperature_index(process_data.isothermal ? -1
250 : ShapeFunction::NPOINTS),
251 first_concentration_index(process_data.isothermal
252 ? ShapeFunction::NPOINTS
253 : 2 * ShapeFunction::NPOINTS),
254 _element(element),
255 _process_data(process_data),
256 _integration_method(integration_method),
257 _transport_process_variables(transport_process_variables)
258 {
259 (void)local_matrix_size;
260
261 unsigned const n_integration_points =
263 _ip_data.reserve(n_integration_points);
264
267
268 double const aperture_size = _process_data.aperture_size(0.0, pos)[0];
269
270 auto const shape_matrices =
272 GlobalDim>(element, is_axially_symmetric,
274 auto const& medium =
276 for (unsigned ip = 0; ip < n_integration_points; ip++)
277 {
278 _ip_data.emplace_back(
279 shape_matrices[ip].dNdx,
281 shape_matrices[ip].integralMeasure *
282 shape_matrices[ip].detJ * aperture_size);
283
284 pos.setIntegrationPoint(ip);
285
286 _ip_data[ip].porosity =
288 .template initialValue<double>(
289 pos, std::numeric_limits<double>::quiet_NaN() /*t*/);
290
291 _ip_data[ip].pushBackState();
292 }
293 }
294
295 void setChemicalSystemID(std::size_t const /*mesh_item_id*/) override
296 {
298 // chemical system index map
299 auto& chemical_system_index_map =
301
302 unsigned const n_integration_points =
304 for (unsigned ip = 0; ip < n_integration_points; ip++)
305 {
306 _ip_data[ip].chemical_system_id =
307 chemical_system_index_map.empty()
308 ? 0
309 : chemical_system_index_map.back() + 1;
310 chemical_system_index_map.push_back(
311 _ip_data[ip].chemical_system_id);
312 }
313 }
314
315 void initializeChemicalSystemConcrete(Eigen::VectorXd const& local_x,
316 double const t) override
317 {
319
320 auto const& medium =
322
325
326 auto const& Ns =
328 .NsHigherOrder<typename ShapeFunction::MeshElement>();
329
330 unsigned const n_integration_points =
332
333 for (unsigned ip = 0; ip < n_integration_points; ip++)
334 {
335 auto& ip_data = _ip_data[ip];
336 auto const& N = Ns[ip];
337 auto const& chemical_system_id = ip_data.chemical_system_id;
338
339 auto const n_component = _transport_process_variables.size();
340 std::vector<double> C_int_pt(n_component);
341 for (unsigned component_id = 0; component_id < n_component;
342 ++component_id)
343 {
344 auto const concentration_index =
346 component_id * concentration_size;
347 auto const local_C =
348 local_x.template segment<concentration_size>(
349 concentration_index);
350
352 C_int_pt[component_id]);
353 }
354
356 ->initializeChemicalSystemConcrete(C_int_pt, chemical_system_id,
357 medium, pos, t);
358 }
359 }
360
361 void setChemicalSystemConcrete(Eigen::VectorXd const& local_x,
362 double const t, double dt) override
363 {
365
366 auto const& medium =
368
371
374
375 auto const& Ns =
377 .NsHigherOrder<typename ShapeFunction::MeshElement>();
378
379 unsigned const n_integration_points =
381
382 for (unsigned ip = 0; ip < n_integration_points; ip++)
383 {
384 auto& ip_data = _ip_data[ip];
385 auto const& N = Ns[ip];
386 auto& porosity = ip_data.porosity;
387 auto const& porosity_prev = ip_data.porosity_prev;
388 auto const& chemical_system_id = ip_data.chemical_system_id;
389
390 auto const n_component = _transport_process_variables.size();
391 std::vector<double> C_int_pt(n_component);
392 for (unsigned component_id = 0; component_id < n_component;
393 ++component_id)
394 {
395 auto const concentration_index =
397 component_id * concentration_size;
398 auto const local_C =
399 local_x.template segment<concentration_size>(
400 concentration_index);
401
403 C_int_pt[component_id]);
404 }
405
406 {
407 vars_prev.porosity = porosity_prev;
408
409 porosity =
411 ? porosity_prev
412 : medium
413 ->property(
415 .template value<double>(vars, vars_prev, pos, t,
416 dt);
417
418 vars.porosity = porosity;
419 }
420
422 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
423 }
424 }
425
426 void postSpeciationCalculation(std::size_t const ele_id, double const t,
427 double const dt) override
428 {
430 {
431 return;
432 }
433
434 auto const& medium = *_process_data.media_map.getMedium(ele_id);
435
437 pos.setElementID(ele_id);
438
439 for (auto& ip_data : _ip_data)
440 {
441 ip_data.porosity = ip_data.porosity_prev;
442
444 ->updateVolumeFractionPostReaction(ip_data.chemical_system_id,
445 medium, pos,
446 ip_data.porosity, t, dt);
447
449 ip_data.chemical_system_id, medium, ip_data.porosity);
450 }
451 }
452
453 void assemble(double const t, double const dt,
454 std::vector<double> const& local_x,
455 std::vector<double> const& /*local_x_prev*/,
456 std::vector<double>& local_M_data,
457 std::vector<double>& local_K_data,
458 std::vector<double>& local_b_data) override
459 {
460 auto const local_matrix_size = local_x.size();
461 // Nodal DOFs include pressure
462 int const num_nodal_dof = 1 + _transport_process_variables.size();
463 // This assertion is valid only if all nodal d.o.f. use the same shape
464 // matrices.
465 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
466
467 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
468 local_M_data, local_matrix_size, local_matrix_size);
469 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
470 local_K_data, local_matrix_size, local_matrix_size);
471 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
472 local_b_data, local_matrix_size);
473
474 // Get block matrices
475 auto Kpp = local_K.template block<pressure_size, pressure_size>(
477 auto Mpp = local_M.template block<pressure_size, pressure_size>(
479 auto Bp = local_b.template segment<pressure_size>(pressure_index);
480
481 auto local_p = Eigen::Map<const NodalVectorType>(
482 &local_x[pressure_index], pressure_size);
483
484 auto const& b =
487
488 auto const number_of_components = num_nodal_dof - 1;
489 for (int component_id = 0; component_id < number_of_components;
490 ++component_id)
491 {
492 /* Partitioned assembler matrix
493 * | pp | pc1 | pc2 | pc3 |
494 * |-----|-----|-----|-----|
495 * | c1p | c1c1| 0 | 0 |
496 * |-----|-----|-----|-----|
497 * | c2p | 0 | c2c2| 0 |
498 * |-----|-----|-----|-----|
499 * | c3p | 0 | 0 | c3c3|
500 */
501 auto concentration_index =
502 pressure_size + component_id * concentration_size;
503
504 auto KCC =
505 local_K.template block<concentration_size, concentration_size>(
506 concentration_index, concentration_index);
507 auto MCC =
508 local_M.template block<concentration_size, concentration_size>(
509 concentration_index, concentration_index);
510 auto MCp =
511 local_M.template block<concentration_size, pressure_size>(
512 concentration_index, pressure_index);
513 auto MpC =
514 local_M.template block<pressure_size, concentration_size>(
515 pressure_index, concentration_index);
516
517 auto local_C = Eigen::Map<const NodalVectorType>(
518 &local_x[concentration_index], concentration_size);
519
520 assembleBlockMatrices(b, component_id, t, dt, local_C, local_p, KCC,
521 MCC, MCp, MpC, Kpp, Mpp, Bp);
522
524 {
525 auto const stoichiometric_matrix =
528
529 assert(stoichiometric_matrix);
530
531 for (Eigen::SparseMatrix<double>::InnerIterator it(
532 *stoichiometric_matrix, component_id);
533 it;
534 ++it)
535 {
536 auto const stoichiometric_coefficient = it.value();
537 auto const coupled_component_id = it.row();
538 auto const kinetic_prefactor =
540 ->getKineticPrefactor(coupled_component_id);
541
542 auto const concentration_index =
543 pressure_size + component_id * concentration_size;
544 auto const coupled_concentration_index =
546 coupled_component_id * concentration_size;
547 auto KCmCn = local_K.template block<concentration_size,
549 concentration_index, coupled_concentration_index);
550
551 // account for the coupling between components
552 assembleKCmCn(component_id, t, dt, KCmCn,
553 stoichiometric_coefficient,
554 kinetic_prefactor);
555 }
556 }
557 }
558 }
559
561 GlobalDimVectorType const& b, int const component_id, double const t,
562 double const dt,
563 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
564 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
565 Eigen::Ref<LocalBlockMatrixType> KCC,
566 Eigen::Ref<LocalBlockMatrixType> MCC,
567 Eigen::Ref<LocalBlockMatrixType> MCp,
568 Eigen::Ref<LocalBlockMatrixType> MpC,
569 Eigen::Ref<LocalBlockMatrixType> Kpp,
570 Eigen::Ref<LocalBlockMatrixType> Mpp,
571 Eigen::Ref<LocalSegmentVectorType> Bp)
572 {
573 unsigned const n_integration_points =
575
578
580
581 // Get material properties
582 auto const& medium =
584 // Select the only valid for component transport liquid phase.
585 auto const& phase = medium.phase("AqueousLiquid");
586
587 // Assume that the component name is the same as the process variable
588 // name. Components are shifted by one because the first one is always
589 // pressure.
590 auto const& component = phase.component(
591 _transport_process_variables[component_id].get().getName());
592
593 LocalBlockMatrixType KCC_Laplacian =
594 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
595
596 std::vector<GlobalDimVectorType> ip_flux_vector;
597 double average_velocity_norm = 0.0;
599 {
600 ip_flux_vector.reserve(n_integration_points);
601 }
602
603 auto const& Ns =
605 .NsHigherOrder<typename ShapeFunction::MeshElement>();
606
607 for (unsigned ip(0); ip < n_integration_points; ++ip)
608 {
609 pos.setIntegrationPoint(ip);
610
611 auto& ip_data = _ip_data[ip];
612 auto const& dNdx = ip_data.dNdx;
613 auto const& N = Ns[ip];
614 auto const& w = ip_data.integration_weight;
615 auto& porosity = ip_data.porosity;
616
617 double C_int_pt = 0.0;
618 double p_int_pt = 0.0;
619
620 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
621 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
622
623 vars.concentration = C_int_pt;
624 vars.liquid_phase_pressure = p_int_pt;
625
626 // update according to a particular porosity model
628 .template value<double>(vars, pos, t, dt);
629 vars.porosity = porosity;
630
631 auto const& retardation_factor =
633 .template value<double>(vars, pos, t, dt);
634
635 auto const& solute_dispersivity_transverse = medium.template value<
636 double>(
638
639 auto const& solute_dispersivity_longitudinal =
640 medium.template value<double>(
642 longitudinal_dispersivity);
643
644 // Use the fluid density model to compute the density
645 // TODO (renchao): concentration of which component as the argument
646 // for calculation of fluid density
647 auto const density =
649 .template value<double>(vars, pos, t, dt);
650
651 auto const decay_rate =
653 .template value<double>(vars, pos, t, dt);
654
655 auto const& pore_diffusion_coefficient =
656 MaterialPropertyLib::formEigenTensor<GlobalDim>(
658 .value(vars, pos, t, dt));
659
660 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
662 vars, pos, t, dt));
663
664 // Use the viscosity model to compute the viscosity
666 .template value<double>(vars, pos, t, dt);
667
668 GlobalDimMatrixType const K_over_mu = K / mu;
669 GlobalDimVectorType const velocity =
671 ? GlobalDimVectorType(-K_over_mu *
672 (dNdx * p_nodal_values - density * b))
673 : GlobalDimVectorType(-K_over_mu * dNdx * p_nodal_values);
674
675 const double drho_dp =
677 .template dValue<double>(
678 vars,
680 pos, t, dt);
681
682 const double drho_dC =
684 .template dValue<double>(
686 t, dt);
687
688 GlobalDimMatrixType const hydrodynamic_dispersion =
691 pore_diffusion_coefficient, velocity, porosity,
692 solute_dispersivity_transverse,
693 solute_dispersivity_longitudinal);
694
695 const double R_times_phi(retardation_factor * porosity);
696 GlobalDimVectorType const mass_density_flow = velocity * density;
697 auto const N_t_N = (N.transpose() * N).eval();
698
700 {
701 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
702 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
703 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
704 }
705 else
706 {
707 ip_flux_vector.emplace_back(mass_density_flow);
708 average_velocity_norm += velocity.norm();
709 }
710 MCC.noalias() += N_t_N * (R_times_phi * density * w);
711 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
712 KCC_Laplacian.noalias() +=
713 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
714
715 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
716
717 // Calculate Mpp, Kpp, and bp in the first loop over components
718 if (component_id == 0)
719 {
720 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
721 Kpp.noalias() +=
722 dNdx.transpose() * K_over_mu * dNdx * (density * w);
723
725 {
726 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
727 (density * density * w);
728 }
729 }
730 }
731
733 {
735 typename ShapeFunction::MeshElement>(
737 _ip_data,
739 ip_flux_vector,
740 average_velocity_norm /
741 static_cast<double>(n_integration_points),
742 KCC_Laplacian);
743 }
744
745 KCC.noalias() += KCC_Laplacian;
746 }
747
748 void assembleKCmCn(int const component_id, double const t, double const dt,
749 Eigen::Ref<LocalBlockMatrixType> KCmCn,
750 double const stoichiometric_coefficient,
751 double const kinetic_prefactor)
752 {
753 unsigned const n_integration_points =
755
758
760
761 auto const& medium =
763 auto const& phase = medium.phase("AqueousLiquid");
764 auto const& component = phase.component(
765 _transport_process_variables[component_id].get().getName());
766
767 auto const& Ns =
769 .NsHigherOrder<typename ShapeFunction::MeshElement>();
770
771 for (unsigned ip(0); ip < n_integration_points; ++ip)
772 {
773 auto& ip_data = _ip_data[ip];
774 auto const& w = ip_data.integration_weight;
775 auto const& N = Ns[ip];
776 auto& porosity = ip_data.porosity;
777
778 auto const retardation_factor =
780 .template value<double>(vars, pos, t, dt);
781
783 .template value<double>(vars, pos, t, dt);
784
785 auto const density =
787 .template value<double>(vars, pos, t, dt);
788
789 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
790 kinetic_prefactor * retardation_factor *
791 porosity * density * N;
792 }
793 }
794
795 void assembleForStaggeredScheme(double const t, double const dt,
796 Eigen::VectorXd const& local_x,
797 Eigen::VectorXd const& local_x_prev,
798 int const process_id,
799 std::vector<double>& local_M_data,
800 std::vector<double>& local_K_data,
801 std::vector<double>& local_b_data) override
802 {
803 if (process_id == _process_data.hydraulic_process_id)
804 {
805 assembleHydraulicEquation(t, dt, local_x, local_x_prev,
806 local_M_data, local_K_data, local_b_data);
807 }
808 else if (process_id == _process_data.thermal_process_id)
809 {
810 assembleHeatTransportEquation(t, dt, local_x, local_x_prev,
811 local_M_data, local_K_data,
812 local_b_data);
813 }
814 else
815 {
816 // Go for assembling in an order of transport process id.
817 assembleComponentTransportEquation(t, dt, local_x, local_x_prev,
818 local_M_data, local_K_data,
819 local_b_data, process_id);
820 }
821 }
822
823 void assembleHydraulicEquation(double const t,
824 double const dt,
825 Eigen::VectorXd const& local_x,
826 Eigen::VectorXd const& local_x_prev,
827 std::vector<double>& local_M_data,
828 std::vector<double>& local_K_data,
829 std::vector<double>& local_b_data)
830 {
831 auto const local_p =
832 local_x.template segment<pressure_size>(pressure_index);
833 auto const local_C = local_x.template segment<concentration_size>(
835 auto const local_C_prev =
836 local_x_prev.segment<concentration_size>(first_concentration_index);
837
838 NodalVectorType local_T = getLocalTemperature(t, local_x);
839
840 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
841 local_M_data, pressure_size, pressure_size);
842 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
843 local_K_data, pressure_size, pressure_size);
844 auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
845 local_b_data, pressure_size);
846
847 unsigned const n_integration_points =
849
852
853 auto const& b =
856
857 auto const& medium =
859 auto const& phase = medium.phase("AqueousLiquid");
860
863
864 auto const& Ns =
866 .NsHigherOrder<typename ShapeFunction::MeshElement>();
867
868 for (unsigned ip(0); ip < n_integration_points; ++ip)
869 {
870 pos.setIntegrationPoint(ip);
871
872 auto& ip_data = _ip_data[ip];
873 auto const& dNdx = ip_data.dNdx;
874 auto const& w = ip_data.integration_weight;
875 auto const& N = Ns[ip];
876 auto& porosity = ip_data.porosity;
877 auto const& porosity_prev = ip_data.porosity_prev;
878
879 double const C_int_pt = N.dot(local_C);
880 double const p_int_pt = N.dot(local_p);
881 double const T_int_pt = N.dot(local_T);
882
883 vars.concentration = C_int_pt;
884 vars.liquid_phase_pressure = p_int_pt;
885 vars.temperature = T_int_pt;
886
887 // porosity
888 {
889 vars_prev.porosity = porosity_prev;
890
891 porosity =
893 ? porosity_prev
895 .template value<double>(vars, vars_prev, pos, t,
896 dt);
897
898 vars.porosity = porosity;
899 }
900
901 // Use the fluid density model to compute the density
902 // TODO: Concentration of which component as one of arguments for
903 // calculation of fluid density
904 auto const density =
906 .template value<double>(vars, pos, t, dt);
907
908 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
910 vars, pos, t, dt));
911
912 // Use the viscosity model to compute the viscosity
914 .template value<double>(vars, pos, t, dt);
915
916 GlobalDimMatrixType const K_over_mu = K / mu;
917
918 const double drho_dp =
920 .template dValue<double>(
921 vars,
923 pos, t, dt);
924 const double drho_dC =
926 .template dValue<double>(
928 t, dt);
929
930 // matrix assembly
931 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
932 local_K.noalias() +=
933 w * dNdx.transpose() * density * K_over_mu * dNdx;
934
936 {
937 local_b.noalias() +=
938 w * density * density * dNdx.transpose() * K_over_mu * b;
939 }
940
941 // coupling term
942 {
943 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
944
945 local_b.noalias() -=
946 w * N.transpose() * porosity * drho_dC * C_dot;
947 }
948 }
949 }
950
951 void assembleHeatTransportEquation(double const t, double const dt,
952 Eigen::VectorXd const& local_x,
953 Eigen::VectorXd const& /*local_x_prev*/,
954 std::vector<double>& local_M_data,
955 std::vector<double>& local_K_data,
956 std::vector<double>& /*local_b_data*/)
957 {
958 assert(local_x.size() ==
960
961 auto const local_p =
962 local_x.template segment<pressure_size>(pressure_index);
963 auto const local_T =
964 local_x.template segment<temperature_size>(temperature_index);
965
966 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
967 local_M_data, temperature_size, temperature_size);
968 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
969 local_K_data, temperature_size, temperature_size);
970
972 pos.setElementID(this->_element.getID());
973
974 auto const& process_data = this->_process_data;
975 auto const& medium =
976 *process_data.media_map.getMedium(this->_element.getID());
977 auto const& liquid_phase = medium.phase("AqueousLiquid");
978
979 auto const& b =
982
984
985 unsigned const n_integration_points =
987
988 std::vector<GlobalDimVectorType> ip_flux_vector;
989 double average_velocity_norm = 0.0;
990 ip_flux_vector.reserve(n_integration_points);
991
992 auto const& Ns =
994 .NsHigherOrder<typename ShapeFunction::MeshElement>();
995
996 for (unsigned ip(0); ip < n_integration_points; ip++)
997 {
998 pos.setIntegrationPoint(ip);
999
1000 auto const& ip_data = this->_ip_data[ip];
1001 auto const& dNdx = ip_data.dNdx;
1002 auto const& w = ip_data.integration_weight;
1003 auto const& N = Ns[ip];
1004
1005 double p_at_xi = 0.;
1006 NumLib::shapeFunctionInterpolate(local_p, N, p_at_xi);
1007 double T_at_xi = 0.;
1008 NumLib::shapeFunctionInterpolate(local_T, N, T_at_xi);
1009
1010 vars.temperature = T_at_xi;
1011 vars.liquid_phase_pressure = p_at_xi;
1012
1013 vars.liquid_saturation = 1.0;
1014
1015 auto const porosity =
1017 .template value<double>(vars, pos, t, dt);
1018 vars.porosity = porosity;
1019
1020 // Use the fluid density model to compute the density
1021 auto const fluid_density =
1022 liquid_phase
1024 .template value<double>(vars, pos, t, dt);
1025 vars.density = fluid_density;
1026 auto const specific_heat_capacity_fluid =
1027 liquid_phase
1029 .template value<double>(vars, pos, t, dt);
1030
1031 // Assemble mass matrix
1032 local_M.noalias() += w *
1034 vars, porosity, fluid_density,
1035 specific_heat_capacity_fluid, pos, t, dt) *
1036 N.transpose() * N;
1037
1038 // Assemble Laplace matrix
1039 auto const viscosity =
1040 liquid_phase
1042 .template value<double>(vars, pos, t, dt);
1043
1044 auto const intrinsic_permeability =
1045 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1046 medium
1047 .property(
1049 .value(vars, pos, t, dt));
1050
1051 GlobalDimMatrixType const K_over_mu =
1052 intrinsic_permeability / viscosity;
1053 GlobalDimVectorType const velocity =
1054 process_data.has_gravity
1055 ? GlobalDimVectorType(-K_over_mu *
1056 (dNdx * local_p - fluid_density * b))
1057 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1058
1059 GlobalDimMatrixType const thermal_conductivity_dispersivity =
1061 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1062 pos, t, dt);
1063
1064 local_K.noalias() +=
1065 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1066
1067 ip_flux_vector.emplace_back(velocity * fluid_density *
1068 specific_heat_capacity_fluid);
1069 average_velocity_norm += velocity.norm();
1070 }
1071
1072 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1073 process_data.stabilizer, this->_ip_data,
1074 _process_data.shape_matrix_cache, ip_flux_vector,
1075 average_velocity_norm / static_cast<double>(n_integration_points),
1076 local_K);
1077 }
1078
1080 double const t, double const dt, Eigen::VectorXd const& local_x,
1081 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_M_data,
1082 std::vector<double>& local_K_data,
1083 std::vector<double>& /*local_b_data*/, int const transport_process_id)
1084 {
1085 assert(static_cast<int>(local_x.size()) ==
1088 static_cast<int>(_transport_process_variables.size()) +
1090
1091 auto const local_p =
1092 local_x.template segment<pressure_size>(pressure_index);
1093
1094 NodalVectorType local_T = getLocalTemperature(t, local_x);
1095
1096 auto const local_C = local_x.template segment<concentration_size>(
1098 (transport_process_id - (_process_data.isothermal ? 1 : 2)) *
1100 auto const local_p_prev =
1101 local_x_prev.segment<pressure_size>(pressure_index);
1102
1103 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1104 local_M_data, concentration_size, concentration_size);
1105 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1106 local_K_data, concentration_size, concentration_size);
1107
1108 LocalBlockMatrixType KCC_Laplacian =
1109 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1110
1111 unsigned const n_integration_points =
1113
1114 std::vector<GlobalDimVectorType> ip_flux_vector;
1115 double average_velocity_norm = 0.0;
1117 {
1118 ip_flux_vector.reserve(n_integration_points);
1119 }
1120
1123
1124 auto const& b =
1127
1130
1131 auto const& medium =
1133 auto const& phase = medium.phase("AqueousLiquid");
1134 auto const component_id =
1135 transport_process_id - (_process_data.isothermal ? 1 : 2);
1136 auto const& component = phase.component(
1137 _transport_process_variables[component_id].get().getName());
1138
1139 auto const& Ns =
1141 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1142
1143 for (unsigned ip(0); ip < n_integration_points; ++ip)
1144 {
1145 pos.setIntegrationPoint(ip);
1146
1147 auto& ip_data = _ip_data[ip];
1148 auto const& dNdx = ip_data.dNdx;
1149 auto const& w = ip_data.integration_weight;
1150 auto const& N = Ns[ip];
1151 auto& porosity = ip_data.porosity;
1152 auto const& porosity_prev = ip_data.porosity_prev;
1153
1154 double const C_int_pt = N.dot(local_C);
1155 double const p_int_pt = N.dot(local_p);
1156 double const T_int_pt = N.dot(local_T);
1157
1158 vars.concentration = C_int_pt;
1159 vars.liquid_phase_pressure = p_int_pt;
1160 vars.temperature = T_int_pt;
1161
1163 {
1164 vars.temperature = N.dot(local_T);
1165 }
1166
1167 // porosity
1168 {
1169 vars_prev.porosity = porosity_prev;
1170
1171 porosity =
1173 ? porosity_prev
1175 .template value<double>(vars, vars_prev, pos, t,
1176 dt);
1177
1178 vars.porosity = porosity;
1179 }
1180
1181 auto const& retardation_factor =
1183 .template value<double>(vars, pos, t, dt);
1184
1185 auto const& solute_dispersivity_transverse = medium.template value<
1186 double>(
1188 auto const& solute_dispersivity_longitudinal =
1189 medium.template value<double>(
1191 longitudinal_dispersivity);
1192
1193 // Use the fluid density model to compute the density
1194 auto const density =
1196 .template value<double>(vars, pos, t, dt);
1197 auto const decay_rate =
1199 .template value<double>(vars, pos, t, dt);
1200
1201 auto const& pore_diffusion_coefficient =
1202 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1204 .value(vars, pos, t, dt));
1205
1206 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1208 vars, pos, t, dt));
1209 // Use the viscosity model to compute the viscosity
1211 .template value<double>(vars, pos, t, dt);
1212
1213 GlobalDimMatrixType const K_over_mu = K / mu;
1214 GlobalDimVectorType const velocity =
1216 ? GlobalDimVectorType(-K_over_mu *
1217 (dNdx * local_p - density * b))
1218 : GlobalDimVectorType(-K_over_mu * dNdx * local_p);
1219
1220 GlobalDimMatrixType const hydrodynamic_dispersion =
1223 pore_diffusion_coefficient, velocity, porosity,
1224 solute_dispersivity_transverse,
1225 solute_dispersivity_longitudinal);
1226
1227 double const R_times_phi = retardation_factor * porosity;
1228 auto const N_t_N = (N.transpose() * N).eval();
1229
1231 {
1232 const double drho_dC =
1234 .template dValue<double>(
1236 pos, t, dt);
1237 local_M.noalias() +=
1238 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1239 }
1240
1241 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1242
1243 // coupling term
1245 {
1246 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1247
1248 const double drho_dp =
1250 .template dValue<double>(vars,
1252 liquid_phase_pressure,
1253 pos, t, dt);
1254
1255 local_K.noalias() +=
1256 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1257 dNdx.transpose() * velocity * N * (density * w);
1258 }
1259 else
1260 {
1261 ip_flux_vector.emplace_back(velocity * density);
1262 average_velocity_norm += velocity.norm();
1263 }
1264 local_K.noalias() +=
1265 N_t_N * (decay_rate * R_times_phi * density * w);
1266
1267 KCC_Laplacian.noalias() += dNdx.transpose() *
1268 hydrodynamic_dispersion * dNdx *
1269 (density * w);
1270 }
1271
1273 {
1275 typename ShapeFunction::MeshElement>(
1277 _process_data.shape_matrix_cache, ip_flux_vector,
1278 average_velocity_norm /
1279 static_cast<double>(n_integration_points),
1280 KCC_Laplacian);
1281 }
1282 local_K.noalias() += KCC_Laplacian;
1283 }
1284
1286 double const t, double const dt, Eigen::VectorXd const& local_x,
1287 Eigen::VectorXd const& local_x_prev, int const process_id,
1288 std::vector<double>& local_b_data,
1289 std::vector<double>& local_Jac_data) override
1290 {
1291 if (process_id == _process_data.hydraulic_process_id)
1292 {
1293 assembleWithJacobianHydraulicEquation(t, dt, local_x, local_x_prev,
1294 local_b_data, local_Jac_data);
1295 }
1296 else
1297 {
1298 int const component_id = process_id - 1;
1300 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1301 component_id);
1302 }
1303 }
1304
1306 double const t, double const dt, Eigen::VectorXd const& local_x,
1307 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1308 std::vector<double>& local_Jac_data)
1309 {
1310 auto const p = local_x.template segment<pressure_size>(pressure_index);
1311 auto const c = local_x.template segment<concentration_size>(
1313
1314 auto const p_prev = local_x_prev.segment<pressure_size>(pressure_index);
1315 auto const c_prev =
1316 local_x_prev.segment<concentration_size>(first_concentration_index);
1317
1318 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1319 local_Jac_data, pressure_size, pressure_size);
1320 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1321 local_b_data, pressure_size);
1322
1323 unsigned const n_integration_points =
1325
1328 auto const& b =
1331
1332 auto const& medium =
1334 auto const& phase = medium.phase("AqueousLiquid");
1335
1338
1339 auto const& Ns =
1341 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1342
1343 for (unsigned ip(0); ip < n_integration_points; ++ip)
1344 {
1345 pos.setIntegrationPoint(ip);
1346
1347 auto& ip_data = _ip_data[ip];
1348 auto const& dNdx = ip_data.dNdx;
1349 auto const& w = ip_data.integration_weight;
1350 auto const& N = Ns[ip];
1351 auto& phi = ip_data.porosity;
1352 auto const& phi_prev = ip_data.porosity_prev;
1353
1354 double const p_ip = N.dot(p);
1355 double const c_ip = N.dot(c);
1356
1357 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1358
1359 vars.liquid_phase_pressure = p_ip;
1360 vars.concentration = c_ip;
1361
1362 // porosity
1363 {
1364 vars_prev.porosity = phi_prev;
1365
1367 ? phi_prev
1369 .template value<double>(vars, vars_prev, pos, t,
1370 dt);
1371
1372 vars.porosity = phi;
1373 }
1374
1375 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1376 .template value<double>(vars, pos, t, dt);
1377
1378 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1380 vars, pos, t, dt));
1381
1383 .template value<double>(vars, pos, t, dt);
1384
1385 auto const drho_dp =
1387 .template dValue<double>(
1388 vars,
1390 pos, t, dt);
1391 auto const drho_dc =
1393 .template dValue<double>(
1395 t, dt);
1396
1397 // matrix assembly
1398 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1399 w * dNdx.transpose() * rho * k / mu * dNdx;
1400
1401 local_rhs.noalias() -=
1402 w * N.transpose() * phi *
1403 (drho_dp * N * p_prev + drho_dc * cdot_ip) +
1404 w * rho * dNdx.transpose() * k / mu * dNdx * p;
1405
1407 {
1408 local_rhs.noalias() +=
1409 w * rho * dNdx.transpose() * k / mu * rho * b;
1410 }
1411 }
1412 }
1413
1415 double const t, double const dt, Eigen::VectorXd const& local_x,
1416 Eigen::VectorXd const& local_x_prev, std::vector<double>& local_b_data,
1417 std::vector<double>& local_Jac_data, int const component_id)
1418 {
1419 auto const concentration_index =
1421
1422 auto const p = local_x.template segment<pressure_size>(pressure_index);
1423 auto const c =
1424 local_x.template segment<concentration_size>(concentration_index);
1425 auto const c_prev =
1426 local_x_prev.segment<concentration_size>(concentration_index);
1427
1430 {
1432 }
1433
1434 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1435 local_Jac_data, concentration_size, concentration_size);
1436 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1437 local_b_data, concentration_size);
1438
1439 LocalBlockMatrixType KCC_Laplacian =
1440 LocalBlockMatrixType::Zero(concentration_size, concentration_size);
1441
1442 unsigned const n_integration_points =
1444
1445 std::vector<GlobalDimVectorType> ip_flux_vector;
1446 double average_velocity_norm = 0.0;
1447 ip_flux_vector.reserve(n_integration_points);
1448
1451
1452 auto const& b =
1455
1458
1459 auto const& medium =
1461 auto const& phase = medium.phase("AqueousLiquid");
1462 auto const& component = phase.component(
1463 _transport_process_variables[component_id].get().getName());
1464
1465 auto const& Ns =
1467 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1468
1469 for (unsigned ip(0); ip < n_integration_points; ++ip)
1470 {
1471 pos.setIntegrationPoint(ip);
1472
1473 auto& ip_data = _ip_data[ip];
1474 auto const& dNdx = ip_data.dNdx;
1475 auto const& w = ip_data.integration_weight;
1476 auto const& N = Ns[ip];
1477 auto& phi = ip_data.porosity;
1478 auto const& phi_prev = ip_data.porosity_prev;
1479
1480 double const p_ip = N.dot(p);
1481 double const c_ip = N.dot(c);
1482
1483 vars.liquid_phase_pressure = p_ip;
1484 vars.concentration = c_ip;
1485
1487 {
1488 vars.temperature = N.dot(T);
1489 }
1490
1491 // porosity
1492 {
1493 vars_prev.porosity = phi_prev;
1494
1496 ? phi_prev
1498 .template value<double>(vars, vars_prev, pos, t,
1499 dt);
1500
1501 vars.porosity = phi;
1502 }
1503
1504 auto const R =
1506 .template value<double>(vars, pos, t, dt);
1507
1508 auto const alpha_T = medium.template value<double>(
1510 auto const alpha_L = medium.template value<double>(
1512
1513 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1514 .template value<double>(vars, pos, t, dt);
1515 // first-order decay constant
1516 auto const alpha =
1518 .template value<double>(vars, pos, t, dt);
1519
1520 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1522 .value(vars, pos, t, dt));
1523
1524 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1526 vars, pos, t, dt));
1528 .template value<double>(vars, pos, t, dt);
1529 // Darcy flux
1530 GlobalDimVectorType const q =
1532 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1533 : GlobalDimVectorType(-k / mu * dNdx * p);
1534
1536 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1537 alpha_L);
1538
1539 // matrix assembly
1540 local_Jac.noalias() +=
1541 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1542
1543 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1544
1545 auto const cdot = (c - c_prev) / dt;
1546 local_rhs.noalias() -=
1547 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1548
1549 ip_flux_vector.emplace_back(q * rho);
1550 average_velocity_norm += q.norm();
1551 }
1552
1553 NumLib::assembleAdvectionMatrix<typename ShapeFunction::MeshElement>(
1555 _process_data.shape_matrix_cache, ip_flux_vector,
1556 average_velocity_norm / static_cast<double>(n_integration_points),
1557 KCC_Laplacian);
1558
1559 local_rhs.noalias() -= KCC_Laplacian * c;
1560
1561 local_Jac.noalias() += KCC_Laplacian;
1562 }
1563
1565 double const t, double const dt, Eigen::VectorXd const& local_x,
1566 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1567 std::vector<double>& local_b_data,
1568 int const transport_process_id) override
1569 {
1570 auto const local_C = local_x.template segment<concentration_size>(
1572 (transport_process_id - 1) * concentration_size);
1573
1574 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1575 local_M_data, concentration_size, concentration_size);
1576 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1577 local_K_data, concentration_size, concentration_size);
1578 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1579 local_b_data, concentration_size);
1580
1581 unsigned const n_integration_points =
1583
1586
1589
1590 auto const& medium =
1592 auto const component_id = transport_process_id - 1;
1593
1594 auto const& Ns =
1596 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1597
1598 for (unsigned ip(0); ip < n_integration_points; ++ip)
1599 {
1600 pos.setIntegrationPoint(ip);
1601
1602 auto& ip_data = _ip_data[ip];
1603 auto const w = ip_data.integration_weight;
1604 auto const& N = Ns[ip];
1605 auto& porosity = ip_data.porosity;
1606 auto const& porosity_prev = ip_data.porosity_prev;
1607 auto const chemical_system_id = ip_data.chemical_system_id;
1608
1609 double C_int_pt = 0.0;
1610 NumLib::shapeFunctionInterpolate(local_C, N, C_int_pt);
1611
1612 vars.concentration = C_int_pt;
1613
1614 auto const porosity_dot = (porosity - porosity_prev) / dt;
1615
1616 // porosity
1617 {
1618 vars_prev.porosity = porosity_prev;
1619
1620 porosity =
1622 ? porosity_prev
1624 .template value<double>(vars, vars_prev, pos, t,
1625 dt);
1626 }
1627
1628 local_M.noalias() += w * N.transpose() * porosity * N;
1629
1630 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1631
1632 if (chemical_system_id == -1)
1633 {
1634 continue;
1635 }
1636
1637 auto const C_post_int_pt =
1639 component_id, chemical_system_id);
1640
1641 local_b.noalias() +=
1642 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1643 }
1644 }
1645
1646 std::vector<double> const& getIntPtDarcyVelocity(
1647 const double t,
1648 std::vector<GlobalVector*> const& x,
1649 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
1650 std::vector<double>& cache) const override
1651 {
1652 assert(x.size() == dof_table.size());
1653
1654 auto const n_processes = x.size();
1655 std::vector<std::vector<double>> local_x;
1656 local_x.reserve(n_processes);
1657
1658 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1659 {
1660 auto const indices =
1661 NumLib::getIndices(_element.getID(), *dof_table[process_id]);
1662 assert(!indices.empty());
1663 local_x.push_back(x[process_id]->get(indices));
1664 }
1665
1666 // only one process, must be monolithic.
1667 if (n_processes == 1)
1668 {
1669 auto const local_p = Eigen::Map<const NodalVectorType>(
1670 &local_x[0][pressure_index], pressure_size);
1671 auto const local_C = Eigen::Map<const NodalVectorType>(
1673 return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1674 }
1675
1676 // multiple processes, must be staggered.
1677 {
1678 constexpr int pressure_process_id = 0;
1679 constexpr int concentration_process_id = 1;
1680 auto const local_p = Eigen::Map<const NodalVectorType>(
1681 &local_x[pressure_process_id][0], pressure_size);
1682 auto const local_C = Eigen::Map<const NodalVectorType>(
1683 &local_x[concentration_process_id][0], concentration_size);
1684 return calculateIntPtDarcyVelocity(t, local_p, local_C, cache);
1685 }
1686 }
1687
1688 std::vector<double> const& calculateIntPtDarcyVelocity(
1689 const double t,
1690 Eigen::Ref<const NodalVectorType> const& p_nodal_values,
1691 Eigen::Ref<const NodalVectorType> const& C_nodal_values,
1692 std::vector<double>& cache) const
1693 {
1694 auto const n_integration_points =
1696
1697 cache.clear();
1698 auto cache_mat = MathLib::createZeroedMatrix<
1699 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1700 cache, GlobalDim, n_integration_points);
1701
1704
1705 auto const& b =
1708
1710
1711 auto const& medium =
1713 auto const& phase = medium.phase("AqueousLiquid");
1714
1715 auto const& Ns =
1717 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1718
1719 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1720 {
1721 auto const& ip_data = _ip_data[ip];
1722 auto const& dNdx = ip_data.dNdx;
1723 auto const& N = Ns[ip];
1724 auto const& porosity = ip_data.porosity;
1725
1726 pos.setIntegrationPoint(ip);
1727
1728 double C_int_pt = 0.0;
1729 double p_int_pt = 0.0;
1730
1731 NumLib::shapeFunctionInterpolate(C_nodal_values, N, C_int_pt);
1732 NumLib::shapeFunctionInterpolate(p_nodal_values, N, p_int_pt);
1733
1734 vars.concentration = C_int_pt;
1735 vars.liquid_phase_pressure = p_int_pt;
1736 vars.porosity = porosity;
1737
1738 // TODO (naumov) Temporary value not used by current material
1739 // models. Need extension of secondary variables interface.
1740 double const dt = std::numeric_limits<double>::quiet_NaN();
1741 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1743 vars, pos, t, dt));
1745 .template value<double>(vars, pos, t, dt);
1746 GlobalDimMatrixType const K_over_mu = K / mu;
1747
1748 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1750 {
1751 auto const rho_w =
1753 .template value<double>(vars, pos, t, dt);
1754 // here it is assumed that the vector b is directed 'downwards'
1755 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1756 }
1757 }
1758
1759 return cache;
1760 }
1761
1762 Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
1763 const unsigned integration_point) const override
1764 {
1766 typename ShapeFunction::MeshElement>()[integration_point];
1767
1768 // assumes N is stored contiguously in memory
1769 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1770 }
1771
1772 Eigen::Vector3d getFlux(MathLib::Point3d const& pnt_local_coords,
1773 double const t,
1774 std::vector<double> const& local_x) const override
1775 {
1776 auto const local_p = Eigen::Map<const NodalVectorType>(
1777 &local_x[pressure_index], pressure_size);
1778 auto const local_C = Eigen::Map<const NodalVectorType>(
1780
1781 // Eval shape matrices at given point
1782 // Note: Axial symmetry is set to false here, because we only need dNdx
1783 // here, which is not affected by axial symmetry.
1784 auto const shape_matrices =
1786 GlobalDim>(
1787 _element, false /*is_axially_symmetric*/,
1788 std::array{pnt_local_coords})[0];
1789
1792 auto const& b =
1795
1797
1798 auto const& medium =
1800 auto const& phase = medium.phase("AqueousLiquid");
1801
1802 // local_x contains the local concentration and pressure values
1803 double c_int_pt;
1804 NumLib::shapeFunctionInterpolate(local_C, shape_matrices.N, c_int_pt);
1805 vars.concentration = c_int_pt;
1806
1807 double p_int_pt;
1808 NumLib::shapeFunctionInterpolate(local_p, shape_matrices.N, p_int_pt);
1809 vars.liquid_phase_pressure = p_int_pt;
1810
1811 // TODO (naumov) Temporary value not used by current material models.
1812 // Need extension of secondary variables interface.
1813 double const dt = std::numeric_limits<double>::quiet_NaN();
1814 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1816 vars, pos, t, dt));
1817
1819 .template value<double>(vars, pos, t, dt);
1820 GlobalDimMatrixType const K_over_mu = K / mu;
1821
1822 GlobalDimVectorType q = -K_over_mu * shape_matrices.dNdx * local_p;
1823 auto const rho_w = phase[MaterialPropertyLib::PropertyType::density]
1824 .template value<double>(vars, pos, t, dt);
1826 {
1827 q += K_over_mu * rho_w * b;
1828 }
1829 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1830 flux.head<GlobalDim>() = rho_w * q;
1831 return flux;
1832 }
1833
1835 double const t,
1836 double const /*dt*/,
1837 Eigen::VectorXd const& local_x,
1838 Eigen::VectorXd const& /*local_x_prev*/) override
1839 {
1840 auto const local_p =
1841 local_x.template segment<pressure_size>(pressure_index);
1842 auto const local_C = local_x.template segment<concentration_size>(
1844
1845 std::vector<double> ele_velocity;
1846 calculateIntPtDarcyVelocity(t, local_p, local_C, ele_velocity);
1847
1848 auto const n_integration_points =
1850 auto const ele_velocity_mat =
1851 MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
1852
1853 auto const ele_id = _element.getID();
1854 Eigen::Map<LocalVectorType>(
1855 &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
1856 GlobalDim) =
1857 ele_velocity_mat.rowwise().sum() / n_integration_points;
1858 }
1859
1861 std::size_t const ele_id) override
1862 {
1863 auto const n_integration_points =
1865
1867 {
1868 auto const& medium = *_process_data.media_map.getMedium(ele_id);
1869
1870 for (auto& ip_data : _ip_data)
1871 {
1872 ip_data.porosity = ip_data.porosity_prev;
1873
1875 ->updatePorosityPostReaction(ip_data.chemical_system_id,
1876 medium, ip_data.porosity);
1877 }
1878
1880 std::accumulate(_ip_data.begin(), _ip_data.end(), 0.,
1881 [](double const s, auto const& ip)
1882 { return s + ip.porosity; }) /
1883 n_integration_points;
1884 }
1885
1886 std::vector<GlobalIndexType> chemical_system_indices;
1887 chemical_system_indices.reserve(n_integration_points);
1888 std::transform(_ip_data.begin(), _ip_data.end(),
1889 std::back_inserter(chemical_system_indices),
1890 [](auto const& ip_data)
1891 { return ip_data.chemical_system_id; });
1892
1894 ele_id, chemical_system_indices);
1895 }
1896
1897 std::vector<double> const& getIntPtMolarFlux(
1898 const double t, std::vector<GlobalVector*> const& x,
1899 std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
1900 std::vector<double>& cache, int const component_id) const override
1901 {
1902 std::vector<double> local_x_vec;
1903
1904 auto const n_processes = x.size();
1905 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1906 {
1907 auto const indices =
1908 NumLib::getIndices(_element.getID(), *dof_tables[process_id]);
1909 assert(!indices.empty());
1910 auto const local_solution = x[process_id]->get(indices);
1911 local_x_vec.insert(std::end(local_x_vec),
1912 std::begin(local_solution),
1913 std::end(local_solution));
1914 }
1915 auto const local_x = MathLib::toVector(local_x_vec);
1916
1917 auto const p = local_x.template segment<pressure_size>(pressure_index);
1918 auto const c = local_x.template segment<concentration_size>(
1920
1921 auto const n_integration_points =
1923
1924 cache.clear();
1925 auto cache_mat = MathLib::createZeroedMatrix<
1926 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1927 cache, GlobalDim, n_integration_points);
1928
1931
1932 auto const& b =
1935
1937
1938 auto const& medium =
1940 auto const& phase = medium.phase("AqueousLiquid");
1941
1942 auto const& component = phase.component(
1943 _transport_process_variables[component_id].get().getName());
1944
1945 auto const& Ns =
1947 .NsHigherOrder<typename ShapeFunction::MeshElement>();
1948
1949 for (unsigned ip = 0; ip < n_integration_points; ++ip)
1950 {
1951 auto const& ip_data = _ip_data[ip];
1952 auto const& dNdx = ip_data.dNdx;
1953 auto const& N = Ns[ip];
1954 auto const& phi = ip_data.porosity;
1955
1956 pos.setIntegrationPoint(ip);
1957
1958 double const p_ip = N.dot(p);
1959 double const c_ip = N.dot(c);
1960
1961 vars.concentration = c_ip;
1962 vars.liquid_phase_pressure = p_ip;
1963 vars.porosity = phi;
1964
1965 double const dt = std::numeric_limits<double>::quiet_NaN();
1966
1967 auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1969 vars, pos, t, dt));
1971 .template value<double>(vars, pos, t, dt);
1972 auto const rho = phase[MaterialPropertyLib::PropertyType::density]
1973 .template value<double>(vars, pos, t, dt);
1974
1975 // Darcy flux
1976 GlobalDimVectorType const q =
1978 ? GlobalDimVectorType(-k / mu * (dNdx * p - rho * b))
1979 : GlobalDimVectorType(-k / mu * dNdx * p);
1980
1981 auto const alpha_T = medium.template value<double>(
1983 auto const alpha_L = medium.template value<double>(
1985 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1987 .value(vars, pos, t, dt));
1988
1989 // Hydrodynamic dispersion
1991 _process_data.stabilizer, _element.getID(), Dp, q, phi, alpha_T,
1992 alpha_L);
1993
1994 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
1995 }
1996
1997 return cache;
1998 }
1999
2000 void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
2001 Eigen::VectorXd const& /*local_x_prev*/,
2002 double const /*t*/, double const /*dt*/,
2003 int const /*process_id*/) override
2004 {
2005 unsigned const n_integration_points =
2007
2008 for (unsigned ip = 0; ip < n_integration_points; ip++)
2009 {
2010 _ip_data[ip].pushBackState();
2011 }
2012 }
2013
2014private:
2017
2019 std::vector<std::reference_wrapper<ProcessVariable>> const
2021
2022 std::vector<IntegrationPointData<GlobalDimNodalMatrixType>> _ip_data;
2023
2025 MaterialPropertyLib::VariableArray const& vars, const double porosity,
2026 const double fluid_density, const double specific_heat_capacity_fluid,
2027 ParameterLib::SpatialPosition const& pos, double const t,
2028 double const dt)
2029 {
2030 auto const& medium =
2031 *_process_data.media_map.getMedium(this->_element.getID());
2032 auto const& solid_phase = medium.phase("Solid");
2033
2034 auto const specific_heat_capacity_solid =
2035 solid_phase
2036 .property(
2038 .template value<double>(vars, pos, t, dt);
2039
2040 auto const solid_density =
2041 solid_phase.property(MaterialPropertyLib::PropertyType::density)
2042 .template value<double>(vars, pos, t, dt);
2043
2044 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2045 fluid_density * specific_heat_capacity_fluid * porosity;
2046 }
2047
2050 const double fluid_density, const double specific_heat_capacity_fluid,
2051 const GlobalDimVectorType& velocity,
2052 ParameterLib::SpatialPosition const& pos, double const t,
2053 double const dt)
2054 {
2055 auto const& medium =
2057
2058 auto thermal_conductivity =
2059 MaterialPropertyLib::formEigenTensor<GlobalDim>(
2060 medium
2061 .property(
2063 .value(vars, pos, t, dt));
2064
2065 auto const thermal_dispersivity_transversal =
2066 medium
2068 thermal_transversal_dispersivity)
2069 .template value<double>();
2070
2071 auto const thermal_dispersivity_longitudinal =
2072 medium
2074 thermal_longitudinal_dispersivity)
2075 .template value<double>();
2076
2077 // Thermal conductivity is moved outside and zero matrix is passed
2078 // instead due to multiplication with fluid's density times specific
2079 // heat capacity.
2080 return thermal_conductivity +
2081 fluid_density * specific_heat_capacity_fluid *
2084 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2085 velocity, 0 /* phi */, thermal_dispersivity_transversal,
2086 thermal_dispersivity_longitudinal);
2087 }
2088
2090 Eigen::VectorXd const& local_x)
2091 {
2092 NodalVectorType local_T;
2094 {
2096 {
2098 _element, t);
2099 }
2100 else
2101 {
2102 local_T = NodalVectorType::Zero(temperature_size);
2103 }
2104 }
2105 else
2106 {
2107 local_T =
2108 local_x.template segment<temperature_size>(temperature_index);
2109 }
2110 return local_T;
2111 }
2112};
2113
2114} // namespace ComponentTransport
2115} // 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
Property const & property(PropertyType const &p) const
Definition Phase.cpp:53
Component const & component(std::size_t const &index) const
Definition Phase.cpp:33
int add(IndexType row, IndexType col, double val)
Definition EigenMatrix.h:87
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
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
auto const & NsHigherOrder() const
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 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
virtual std::vector< double > const & getIntPtMolarFlux(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache, int const component_id) const =0
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
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 assembleHeatTransportEquation(double const t, double const dt, Eigen::VectorXd const &local_x, Eigen::VectorXd const &, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > LocalMatrixType
void postTimestepConcrete(Eigen::VectorXd const &, Eigen::VectorXd const &, double const, double const, int const) override
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, int const component_id) 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
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::template MatrixType< pressure_size, pressure_size > LocalBlockMatrixType
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 > &local_b_data, std::vector< double > &local_Jac_data) override
std::vector< IntegrationPointData< GlobalDimNodalMatrixType > > _ip_data
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::GlobalDimNodalMatrixType GlobalDimNodalMatrixType
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
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)
NodalVectorType getLocalTemperature(double const t, Eigen::VectorXd const &local_x)
Eigen::Vector3d getFlux(MathLib::Point3d const &pnt_local_coords, double const t, std::vector< double > const &local_x) const override
GlobalDimMatrixType getThermalConductivityDispersivity(MaterialPropertyLib::VariableArray const &vars, const double fluid_density, const double specific_heat_capacity_fluid, const GlobalDimVectorType &velocity, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
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
double getHeatEnergyCoefficient(MaterialPropertyLib::VariableArray const &vars, const double porosity, const double fluid_density, const double specific_heat_capacity_fluid, ParameterLib::SpatialPosition const &pos, double const t, double const dt)
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)
typename ShapeMatricesType::template VectorType< pressure_size > LocalSegmentVectorType
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.
@ transversal_dispersivity
used to compute the hydrodynamic dispersion tensor.
@ retardation_factor
specify retardation factor used in component transport process.
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)
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
void assembleAdvectionMatrix(IPData const &ip_data_vector, NumLib::ShapeMatrixCache const &shape_matrix_cache, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
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
NumLib::ShapeMatrixCache shape_matrix_cache
caches for each mesh element type the shape matrix
IntegrationPointData(GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)