38namespace ComponentTransport
40template <
typename NodalRowVectorType,
typename GlobalDimNodalMatrixType>
44 GlobalDimNodalMatrixType
const& dNdx_,
45 double const& integration_weight_)
51 NodalRowVectorType
const N;
52 GlobalDimNodalMatrixType
const dNdx;
59 double porosity = std::numeric_limits<double>::quiet_NaN();
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)
78 std::vector<double> local_x_vec;
80 auto const n_processes = x.size();
81 for (std::size_t process_id = 0; process_id < n_processes; ++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));
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)
101 std::vector<double> local_x_vec;
103 auto const n_processes = x.size();
104 for (std::size_t process_id = 0; process_id < n_processes; ++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));
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,
125 std::vector<double> local_x_vec;
127 auto const n_processes = x.size();
128 for (std::size_t pcs_id = 0; pcs_id < n_processes; ++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));
142 auto const num_r_c = indices.size();
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);
152 local_K_data, local_b_data,
155 auto const r_c_indices =
157 if (!local_M_data.empty())
161 M.
add(r_c_indices, local_M);
163 if (!local_K_data.empty())
167 K.
add(r_c_indices, local_K);
169 if (!local_b_data.empty())
171 b.
add(indices, local_b_data);
176 double const t,
double const dt) = 0;
179 std::size_t
const ele_id) = 0;
183 std::vector<GlobalVector*>
const& x,
184 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
185 std::vector<double>& cache)
const = 0;
189 std::vector<GlobalVector*>
const& x,
190 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
191 std::vector<double>& cache)
const = 0;
195 Eigen::VectorXd
const& ,
double const ) = 0;
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;
207template <
typename ShapeFunction,
int GlobalDim>
217 ShapeFunction::NPOINTS;
223 typename ShapeMatricesType::template MatrixType<
pressure_size,
226 typename ShapeMatricesType::template VectorType<pressure_size>;
229 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
243 std::size_t
const local_matrix_size,
245 bool is_axially_symmetric,
247 std::vector<std::reference_wrapper<ProcessVariable>>
const&
248 transport_process_variables)
254 (void)local_matrix_size;
256 unsigned const n_integration_points =
258 _ip_data.reserve(n_integration_points);
265 auto const shape_matrices =
267 GlobalDim>(element, is_axially_symmetric,
271 for (
unsigned ip = 0; ip < n_integration_points; ip++)
274 shape_matrices[ip].N, shape_matrices[ip].dNdx,
276 shape_matrices[ip].integralMeasure *
277 shape_matrices[ip].detJ * aperture_size);
283 .template initialValue<double>(
284 pos, std::numeric_limits<double>::quiet_NaN() );
294 auto& chemical_system_index_map =
297 unsigned const n_integration_points =
299 for (
unsigned ip = 0; ip < n_integration_points; ip++)
302 chemical_system_index_map.empty()
304 : chemical_system_index_map.back() + 1;
305 chemical_system_index_map.push_back(
311 double const t)
override
321 unsigned const n_integration_points =
323 for (
unsigned ip = 0; ip < n_integration_points; ip++)
326 auto const& N = ip_data.N;
327 auto const& chemical_system_id = ip_data.chemical_system_id;
330 std::vector<double> C_int_pt(n_component);
331 for (
unsigned component_id = 0; component_id < n_component;
334 auto const concentration_index =
338 local_x.template segment<concentration_size>(
339 concentration_index);
342 C_int_pt[component_id]);
352 double const t,
double dt)
override
365 unsigned const n_integration_points =
367 for (
unsigned ip = 0; ip < n_integration_points; 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;
376 std::vector<double> C_int_pt(n_component);
377 for (
unsigned component_id = 0; component_id < n_component;
380 auto const concentration_index =
384 local_x.template segment<concentration_size>(
385 concentration_index);
388 C_int_pt[component_id]);
400 .template value<double>(vars, vars_prev, pos, t,
407 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
412 double const dt)
override
426 ip_data.porosity = ip_data.porosity_prev;
431 ip_data.porosity, t, dt);
434 ip_data.chemical_system_id, medium, ip_data.porosity);
439 std::vector<double>
const& local_x,
440 std::vector<double>
const& ,
441 std::vector<double>& local_M_data,
442 std::vector<double>& local_K_data,
443 std::vector<double>& local_b_data)
override
445 auto const local_matrix_size = local_x.size();
450 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
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);
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);
466 auto local_p = Eigen::Map<const NodalVectorType>(
473 auto const number_of_components = num_nodal_dof - 1;
474 for (
int component_id = 0; component_id < number_of_components;
486 auto concentration_index =
490 local_K.template block<concentration_size, concentration_size>(
491 concentration_index, concentration_index);
493 local_M.template block<concentration_size, concentration_size>(
494 concentration_index, concentration_index);
496 local_M.template block<concentration_size, pressure_size>(
499 local_M.template block<pressure_size, concentration_size>(
502 auto local_C = Eigen::Map<const NodalVectorType>(
506 MCC, MCp, MpC, Kpp, Mpp, Bp);
510 auto const stoichiometric_matrix =
514 assert(stoichiometric_matrix);
516 for (Eigen::SparseMatrix<double>::InnerIterator it(
517 *stoichiometric_matrix, component_id);
521 auto const stoichiometric_coefficient = it.value();
522 auto const coupled_component_id = it.row();
523 auto const kinetic_prefactor =
527 auto const concentration_index =
529 auto const coupled_concentration_index =
534 concentration_index, coupled_concentration_index);
538 stoichiometric_coefficient,
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)
558 unsigned const n_integration_points =
570 auto const& phase = medium.
phase(
"AqueousLiquid");
581 std::vector<GlobalDimVectorType> ip_flux_vector;
582 double average_velocity_norm = 0.0;
585 ip_flux_vector.reserve(n_integration_points);
588 for (
unsigned ip(0); ip < n_integration_points; ++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;
598 double C_int_pt = 0.0;
599 double p_int_pt = 0.0;
609 .template value<double>(vars, pos, t, dt);
612 auto const& retardation_factor =
614 .template value<double>(vars, pos, t, dt);
616 auto const& solute_dispersivity_transverse = medium.template value<
620 auto const& solute_dispersivity_longitudinal =
621 medium.template value<double>(
623 longitudinal_dispersivity);
630 .template value<double>(vars, pos, t, dt);
632 auto const decay_rate =
634 .template value<double>(vars, pos, t, dt);
636 auto const& pore_diffusion_coefficient =
637 MaterialPropertyLib::formEigenTensor<GlobalDim>(
639 .value(vars, pos, t, dt));
641 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
647 .template value<double>(vars, pos, t, dt);
653 (dNdx * p_nodal_values - density * b))
656 const double drho_dp =
658 .template dValue<double>(
662 const double drho_dC =
664 .template dValue<double>(
671 pore_diffusion_coefficient, velocity, porosity,
672 solute_dispersivity_transverse,
673 solute_dispersivity_longitudinal);
675 const double R_times_phi(retardation_factor * porosity);
677 auto const N_t_N = (N.transpose() * N).eval();
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;
687 ip_flux_vector.emplace_back(mass_density_flow);
688 average_velocity_norm += velocity.norm();
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;
695 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
698 if (component_id == 0)
700 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
702 dNdx.transpose() * K_over_mu * dNdx * (density * w);
706 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
707 (density * density * w);
718 average_velocity_norm /
719 static_cast<double>(n_integration_points),
723 KCC.noalias() += KCC_Laplacian;
727 Eigen::Ref<LocalBlockMatrixType> KCmCn,
728 double const stoichiometric_coefficient,
729 double const kinetic_prefactor)
731 unsigned const n_integration_points =
741 auto const& phase = medium.
phase(
"AqueousLiquid");
745 for (
unsigned ip(0); ip < n_integration_points; ++ip)
748 auto const& N = ip_data.N;
749 auto const& w = ip_data.integration_weight;
750 auto& porosity = ip_data.porosity;
752 auto const retardation_factor =
754 .template value<double>(vars, pos, t, dt);
757 .template value<double>(vars, pos, t, dt);
761 .template value<double>(vars, pos, t, dt);
763 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
764 kinetic_prefactor * retardation_factor *
765 porosity * density * N;
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
780 local_M_data, local_K_data, local_b_data);
786 local_M_data, local_K_data,
787 local_b_data, process_id);
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)
801 auto const local_C = local_x.template segment<concentration_size>(
803 auto const local_C_prev =
806 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
808 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
810 auto local_b = MathLib::createZeroedVector<LocalSegmentVectorType>(
813 unsigned const n_integration_points =
825 auto const& phase = medium.
phase(
"AqueousLiquid");
830 for (
unsigned ip(0); ip < n_integration_points; ++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;
841 double C_int_pt = 0.0;
842 double p_int_pt = 0.0;
858 .template value<double>(vars, vars_prev, pos, t,
869 .template value<double>(vars, pos, t, dt);
871 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
877 .template value<double>(vars, pos, t, dt);
881 const double drho_dp =
883 .template dValue<double>(
886 const double drho_dC =
888 .template dValue<double>(
893 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
895 w * dNdx.transpose() * density * K_over_mu * dNdx;
900 w * density * density * dNdx.transpose() * K_over_mu * b;
905 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
908 w * N.transpose() * porosity * drho_dC * C_dot;
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>& ,
int const transport_process_id)
921 auto const local_C = local_x.template segment<concentration_size>(
924 auto const local_p_prev =
934 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
936 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
942 unsigned const n_integration_points =
945 std::vector<GlobalDimVectorType> ip_flux_vector;
946 double average_velocity_norm = 0.0;
949 ip_flux_vector.reserve(n_integration_points);
964 auto const& phase = medium.
phase(
"AqueousLiquid");
967 auto const component_id = transport_process_id - 1;
971 for (
unsigned ip(0); ip < n_integration_points; ++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;
982 double C_int_pt = 0.0;
983 double p_int_pt = 0.0;
1004 .template value<double>(vars, vars_prev, pos, t,
1010 auto const& retardation_factor =
1012 .template value<double>(vars, pos, t, dt);
1014 auto const& solute_dispersivity_transverse = medium.template value<
1017 auto const& solute_dispersivity_longitudinal =
1018 medium.template value<double>(
1020 longitudinal_dispersivity);
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);
1030 auto const& pore_diffusion_coefficient =
1031 MaterialPropertyLib::formEigenTensor<GlobalDim>(
1033 .value(vars, pos, t, dt));
1035 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1040 .template value<double>(vars, pos, t, dt);
1046 (dNdx * local_p - density * b))
1052 pore_diffusion_coefficient, velocity, porosity,
1053 solute_dispersivity_transverse,
1054 solute_dispersivity_longitudinal);
1056 double const R_times_phi = retardation_factor * porosity;
1057 auto const N_t_N = (N.transpose() * N).eval();
1061 const double drho_dC =
1063 .template dValue<double>(
1066 local_M.noalias() +=
1067 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1070 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1075 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1077 const double drho_dp =
1079 .template dValue<double>(
1083 local_K.noalias() +=
1084 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1085 dNdx.transpose() * velocity * N * (density * w);
1089 ip_flux_vector.emplace_back(velocity * density);
1090 average_velocity_norm += velocity.norm();
1092 local_K.noalias() +=
1093 N_t_N * (decay_rate * R_times_phi * density * w);
1095 KCC_Laplacian.noalias() += dNdx.transpose() *
1096 hydrodynamic_dispersion * dNdx *
1106 average_velocity_norm /
1107 static_cast<double>(n_integration_points),
1110 local_K.noalias() += KCC_Laplacian;
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>& ,
1117 std::vector<double>& ,
1118 std::vector<double>& local_b_data,
1119 std::vector<double>& local_Jac_data)
override
1124 local_b_data, local_Jac_data);
1128 int const component_id = process_id - 1;
1130 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
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)
1140 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1141 auto const c = local_x.template segment<concentration_size>(
1148 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1150 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1153 unsigned const n_integration_points =
1162 auto const& medium =
1164 auto const& phase = medium.
phase(
"AqueousLiquid");
1169 for (
unsigned ip(0); ip < n_integration_points; ++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;
1180 double const p_ip = N.dot(p);
1181 double const c_ip = N.dot(c);
1183 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1195 .template value<double>(vars, vars_prev, pos, t,
1202 .template value<double>(vars, pos, t, dt);
1204 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1209 .template value<double>(vars, pos, t, dt);
1211 auto const drho_dp =
1213 .template dValue<double>(
1216 auto const drho_dc =
1218 .template dValue<double>(
1223 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1224 w * dNdx.transpose() * rho * k / mu * dNdx;
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;
1233 local_rhs.noalias() +=
1234 w * rho * dNdx.transpose() * k / mu * rho * b;
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)
1244 auto const concentration_index =
1247 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1249 local_x.template segment<concentration_size>(concentration_index);
1259 auto local_Jac = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1261 auto local_rhs = MathLib::createZeroedVector<LocalSegmentVectorType>(
1267 unsigned const n_integration_points =
1270 std::vector<GlobalDimVectorType> ip_flux_vector;
1271 double average_velocity_norm = 0.0;
1272 ip_flux_vector.reserve(n_integration_points);
1284 auto const& medium =
1286 auto const& phase = medium.
phase(
"AqueousLiquid");
1287 auto const& component = phase.
component(
1290 for (
unsigned ip(0); ip < n_integration_points; ++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;
1301 double const p_ip = N.dot(p);
1302 double const c_ip = N.dot(c);
1319 .template value<double>(vars, vars_prev, pos, t,
1327 .template value<double>(vars, pos, t, dt);
1329 auto const alpha_T = medium.template value<double>(
1331 auto const alpha_L = medium.template value<double>(
1335 .template value<double>(vars, pos, t, dt);
1339 .template value<double>(vars, pos, t, dt);
1341 auto const Dp = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1343 .value(vars, pos, t, dt));
1345 auto const k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1349 .template value<double>(vars, pos, t, dt);
1361 local_Jac.noalias() +=
1362 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1364 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1366 auto const cdot = (c - c_prev) / dt;
1367 local_rhs.noalias() -=
1368 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1370 ip_flux_vector.emplace_back(q * rho);
1371 average_velocity_norm += q.norm();
1376 average_velocity_norm /
static_cast<double>(n_integration_points),
1379 local_rhs.noalias() -= KCC_Laplacian * c;
1381 local_Jac.noalias() += KCC_Laplacian;
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
1390 auto const local_C = local_x.template segment<concentration_size>(
1394 auto local_M = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1396 auto local_K = MathLib::createZeroedMatrix<LocalBlockMatrixType>(
1398 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
1401 unsigned const n_integration_points =
1410 auto const& medium =
1412 auto const component_id = transport_process_id - 1;
1413 for (
unsigned ip(0); ip < n_integration_points; ++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;
1424 double C_int_pt = 0.0;
1429 auto const porosity_dot = (porosity - porosity_prev) / dt;
1433 vars_prev.
porosity = porosity_prev;
1439 .template value<double>(vars, vars_prev, pos, t,
1443 local_M.noalias() += w * N.transpose() * porosity * N;
1445 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1447 if (chemical_system_id == -1)
1452 auto const C_post_int_pt =
1454 component_id, chemical_system_id);
1456 local_b.noalias() +=
1457 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1463 std::vector<GlobalVector*>
const& x,
1464 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
1465 std::vector<double>& cache)
const override
1467 assert(x.size() == dof_table.size());
1469 auto const n_processes = x.size();
1470 std::vector<std::vector<double>> local_x;
1471 local_x.reserve(n_processes);
1473 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1475 auto const indices =
1477 assert(!indices.empty());
1478 local_x.push_back(x[process_id]->get(indices));
1482 if (n_processes == 1)
1484 auto const local_p = Eigen::Map<const NodalVectorType>(
1486 auto const local_C = Eigen::Map<const NodalVectorType>(
1493 constexpr int pressure_process_id = 0;
1494 constexpr int concentration_process_id = 1;
1495 auto const local_p = Eigen::Map<const NodalVectorType>(
1497 auto const local_C = Eigen::Map<const NodalVectorType>(
1505 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
1506 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
1507 std::vector<double>& cache)
const
1509 auto const n_integration_points =
1514 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1515 cache, GlobalDim, n_integration_points);
1526 auto const& medium =
1528 auto const& phase = medium.
phase(
"AqueousLiquid");
1530 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
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;
1537 pos.setIntegrationPoint(ip);
1539 double C_int_pt = 0.0;
1540 double p_int_pt = 0.0;
1551 double const dt = std::numeric_limits<double>::quiet_NaN();
1552 auto const& K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1556 .template value<double>(vars, pos, t, dt);
1559 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1564 .template value<double>(vars, pos, t, dt);
1566 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1574 const unsigned integration_point)
const override
1576 auto const& N =
_ip_data[integration_point].N;
1579 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1584 std::vector<double>
const& local_x)
const override
1586 auto const local_p = Eigen::Map<const NodalVectorType>(
1588 auto const local_C = Eigen::Map<const NodalVectorType>(
1594 auto const shape_matrices =
1598 std::array{pnt_local_coords})[0];
1608 auto const& medium =
1610 auto const& phase = medium.
phase(
"AqueousLiquid");
1623 double const dt = std::numeric_limits<double>::quiet_NaN();
1624 auto const K = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1629 .template value<double>(vars, pos, t, dt);
1634 .template value<double>(vars, pos, t, dt);
1637 q += K_over_mu * rho_w * b;
1639 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1640 flux.head<GlobalDim>() = rho_w * q;
1647 Eigen::VectorXd
const& local_x,
1648 Eigen::VectorXd
const& )
override
1650 auto const local_p =
1652 auto const local_C = local_x.template segment<concentration_size>(
1655 std::vector<double> ele_velocity;
1658 auto const n_integration_points =
1660 auto const ele_velocity_mat =
1664 Eigen::Map<LocalVectorType>(
1667 ele_velocity_mat.rowwise().sum() / n_integration_points;
1671 std::size_t
const ele_id)
override
1673 auto const n_integration_points =
1682 ip_data.porosity = ip_data.porosity_prev;
1686 medium, ip_data.porosity);
1691 [](
double const s,
auto const& ip)
1692 {
return s + ip.porosity; }) /
1693 n_integration_points;
1696 std::vector<GlobalIndexType> chemical_system_indices;
1697 chemical_system_indices.reserve(n_integration_points);
1699 std::back_inserter(chemical_system_indices),
1700 [](
auto const& ip_data)
1701 { return ip_data.chemical_system_id; });
1704 ele_id, chemical_system_indices);
1709 std::vector<GlobalVector*>
const& x,
1710 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
1711 std::vector<double>& cache)
const override
1713 std::vector<double> local_x_vec;
1715 auto const n_processes = x.size();
1716 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1718 auto const indices =
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));
1728 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1729 auto const c = local_x.template segment<concentration_size>(
1732 auto const n_integration_points =
1737 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1738 cache, GlobalDim, n_integration_points);
1749 auto const& medium =
1751 auto const& phase = medium.
phase(
"AqueousLiquid");
1753 int const component_id = 0;
1754 auto const& component = phase.
component(
1757 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
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;
1764 pos.setIntegrationPoint(ip);
1766 double const p_ip = N.dot(p);
1767 double const c_ip = N.dot(c);
1773 double const dt = std::numeric_limits<double>::quiet_NaN();
1775 auto const& k = MaterialPropertyLib::formEigenTensor<GlobalDim>(
1779 .template value<double>(vars, pos, t, dt);
1781 .template value<double>(vars, pos, t, dt);
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));
1802 cache_mat.col(ip).noalias() = q * c_ip - phi * D * dNdx * c;
1809 Eigen::VectorXd
const& ,
1810 double const ,
double const ,
1812 int const )
override
1814 unsigned const n_integration_points =
1817 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1828 std::vector<std::reference_wrapper<ProcessVariable>>
const
1833 Eigen::aligned_allocator<
GlobalMatrix::IndexType GlobalIndexType
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)
Medium * getMedium(std::size_t element_id)
Phase const & phase(std::size_t index) const
Component const & component(std::size_t const &index) const
int add(IndexType row, IndexType col, double val)
Global vector based on Eigen vector.
void add(IndexType rowId, double v)
add entry
virtual std::size_t getID() const final
Returns the ID of the element.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() 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
virtual void setChemicalSystemID(std::size_t const)=0
void initializeChemicalSystem(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t)
void setChemicalSystem(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, std::vector< GlobalVector * > const &x, double const t, double const dt)
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
ComponentTransportLocalAssemblerInterface()=default
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
MeshLib::Element const & _element
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)
static const int first_concentration_index
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
static const int concentration_size
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
static const int pressure_size
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)
static const int pressure_index
void setChemicalSystemConcrete(Eigen::VectorXd const &local_x, double const t, double dt) override
void setChemicalSystemID(std::size_t const) 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 assembleAdvectionMatrix(IPData const &ip_data_vector, std::vector< FluxVectorType > const &ip_flux_vector, Eigen::MatrixBase< Derived > &laplacian_matrix)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
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.
NumLib::NumericalStabilization stabilizer
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
static const int hydraulic_process_id
ParameterLib::Parameter< double > const & aperture_size
MeshLib::PropertyVector< double > * mesh_prop_velocity
bool const non_advective_form
ParameterLib::Parameter< double > const *const temperature
bool const chemically_induced_porosity_change
MeshLib::PropertyVector< double > * mesh_prop_porosity
GlobalDimNodalMatrixType const dNdx
GlobalIndexType chemical_system_id
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
NodalRowVectorType const N
IntegrationPointData(NodalRowVectorType const &N_, GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)
double const integration_weight