39namespace ComponentTransport
41template <
typename GlobalDimNodalMatrixType>
45 double const& integration_weight_)
51 GlobalDimNodalMatrixType
const dNdx;
58 double porosity = std::numeric_limits<double>::quiet_NaN();
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)
77 std::vector<double> local_x_vec;
79 auto const n_processes = x.size();
80 for (std::size_t process_id = 0; process_id < n_processes; ++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));
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)
100 std::vector<double> local_x_vec;
102 auto const n_processes = x.size();
103 for (std::size_t process_id = 0; process_id < n_processes; ++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));
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,
124 std::vector<double> local_x_vec;
126 auto const n_processes = x.size();
127 for (std::size_t pcs_id = 0; pcs_id < n_processes; ++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));
141 auto const num_r_c = indices.size();
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);
151 local_K_data, local_b_data,
154 auto const r_c_indices =
156 if (!local_M_data.empty())
160 M.
add(r_c_indices, local_M);
162 if (!local_K_data.empty())
166 K.
add(r_c_indices, local_K);
168 if (!local_b_data.empty())
170 b.
add(indices, local_b_data);
175 double const t,
double const dt) = 0;
178 std::size_t
const ele_id) = 0;
182 std::vector<GlobalVector*>
const& x,
183 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
184 std::vector<double>& cache)
const = 0;
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;
193 Eigen::VectorXd
const& ,
double const ) = 0;
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;
205template <
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)
250 : ShapeFunction::NPOINTS),
252 ? ShapeFunction::NPOINTS
253 : 2 * ShapeFunction::NPOINTS),
259 (void)local_matrix_size;
261 unsigned const n_integration_points =
263 _ip_data.reserve(n_integration_points);
270 auto const shape_matrices =
272 GlobalDim>(element, is_axially_symmetric,
276 for (
unsigned ip = 0; ip < n_integration_points; ip++)
279 shape_matrices[ip].dNdx,
281 shape_matrices[ip].integralMeasure *
282 shape_matrices[ip].detJ * aperture_size);
286 .template initialValue<double>(
287 pos, std::numeric_limits<double>::quiet_NaN() );
297 auto& chemical_system_index_map =
300 unsigned const n_integration_points =
302 for (
unsigned ip = 0; ip < n_integration_points; ip++)
305 chemical_system_index_map.empty()
307 : chemical_system_index_map.back() + 1;
308 chemical_system_index_map.push_back(
314 double const t)
override
328 unsigned const n_integration_points =
331 for (
unsigned ip = 0; ip < n_integration_points; ip++)
334 auto const& N = Ns[ip];
335 auto const& chemical_system_id = ip_data.chemical_system_id;
345 std::vector<double> C_int_pt(n_component);
346 for (
unsigned component_id = 0; component_id < n_component;
349 auto const concentration_index =
353 local_x.template segment<concentration_size>(
354 concentration_index);
357 C_int_pt[component_id]);
367 double const t,
double dt)
override
384 unsigned const n_integration_points =
387 for (
unsigned ip = 0; ip < n_integration_points; ip++)
390 auto const& N = Ns[ip];
391 auto& porosity = ip_data.porosity;
392 auto const& porosity_prev = ip_data.porosity_prev;
393 auto const& chemical_system_id = ip_data.chemical_system_id;
396 std::vector<double> C_int_pt(n_component);
397 for (
unsigned component_id = 0; component_id < n_component;
400 auto const concentration_index =
404 local_x.template segment<concentration_size>(
405 concentration_index);
408 C_int_pt[component_id]);
420 .template value<double>(vars, vars_prev, pos, t,
427 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
432 double const dt)
override
446 ip_data.porosity = ip_data.porosity_prev;
451 ip_data.porosity, t, dt);
454 ip_data.chemical_system_id, medium, ip_data.porosity);
459 std::vector<double>
const& local_x,
460 std::vector<double>
const& ,
461 std::vector<double>& local_M_data,
462 std::vector<double>& local_K_data,
463 std::vector<double>& local_b_data)
override
465 auto const local_matrix_size = local_x.size();
470 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
473 local_M_data, local_matrix_size, local_matrix_size);
475 local_K_data, local_matrix_size, local_matrix_size);
477 local_b_data, local_matrix_size);
480 auto Kpp = local_K.template block<pressure_size, pressure_size>(
482 auto Mpp = local_M.template block<pressure_size, pressure_size>(
484 auto Bp = local_b.template segment<pressure_size>(
pressure_index);
486 auto local_p = Eigen::Map<const NodalVectorType>(
493 auto const number_of_components = num_nodal_dof - 1;
494 for (
int component_id = 0; component_id < number_of_components;
506 auto concentration_index =
510 local_K.template block<concentration_size, concentration_size>(
511 concentration_index, concentration_index);
513 local_M.template block<concentration_size, concentration_size>(
514 concentration_index, concentration_index);
516 local_M.template block<concentration_size, pressure_size>(
519 local_M.template block<pressure_size, concentration_size>(
522 auto local_C = Eigen::Map<const NodalVectorType>(
526 MCC, MCp, MpC, Kpp, Mpp, Bp);
530 auto const stoichiometric_matrix =
534 assert(stoichiometric_matrix);
536 for (Eigen::SparseMatrix<double>::InnerIterator it(
537 *stoichiometric_matrix, component_id);
541 auto const stoichiometric_coefficient = it.value();
542 auto const coupled_component_id = it.row();
543 auto const kinetic_prefactor =
547 auto const concentration_index =
549 auto const coupled_concentration_index =
554 concentration_index, coupled_concentration_index);
558 stoichiometric_coefficient,
568 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
569 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
570 Eigen::Ref<LocalBlockMatrixType> KCC,
571 Eigen::Ref<LocalBlockMatrixType> MCC,
572 Eigen::Ref<LocalBlockMatrixType> MCp,
573 Eigen::Ref<LocalBlockMatrixType> MpC,
574 Eigen::Ref<LocalBlockMatrixType> Kpp,
575 Eigen::Ref<LocalBlockMatrixType> Mpp,
576 Eigen::Ref<LocalSegmentVectorType> Bp)
578 unsigned const n_integration_points =
590 auto const& phase = medium.
phase(
"AqueousLiquid");
601 std::vector<GlobalDimVectorType> ip_flux_vector;
602 double average_velocity_norm = 0.0;
605 ip_flux_vector.reserve(n_integration_points);
612 for (
unsigned ip(0); ip < n_integration_points; ++ip)
615 auto const& dNdx = ip_data.dNdx;
616 auto const& N = Ns[ip];
617 auto const& w = ip_data.integration_weight;
618 auto& porosity = ip_data.porosity;
620 double C_int_pt = 0.0;
621 double p_int_pt = 0.0;
638 .template value<double>(vars, pos, t, dt);
641 auto const& retardation_factor =
643 .template value<double>(vars, pos, t, dt);
645 auto const& solute_dispersivity_transverse = medium.template value<
649 auto const& solute_dispersivity_longitudinal =
650 medium.template value<double>(
652 longitudinal_dispersivity);
659 .template value<double>(vars, pos, t, dt);
661 auto const decay_rate =
663 .template value<double>(vars, pos, t, dt);
665 auto const& pore_diffusion_coefficient =
668 .value(vars, pos, t, dt));
676 .template value<double>(vars, pos, t, dt);
682 (dNdx * p_nodal_values - density * b))
685 const double drho_dp =
687 .template dValue<double>(
692 const double drho_dC =
694 .template dValue<double>(
701 pore_diffusion_coefficient, velocity, porosity,
702 solute_dispersivity_transverse,
703 solute_dispersivity_longitudinal);
705 const double R_times_phi(retardation_factor * porosity);
707 auto const N_t_N = (N.transpose() * N).eval();
711 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
712 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
713 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
717 ip_flux_vector.emplace_back(mass_density_flow);
718 average_velocity_norm += velocity.norm();
720 MCC.noalias() += N_t_N * (R_times_phi * density * w);
721 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
722 KCC_Laplacian.noalias() +=
723 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
725 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
728 if (component_id == 0)
730 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
732 dNdx.transpose() * K_over_mu * dNdx * (density * w);
736 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
737 (density * density * w);
745 typename ShapeFunction::MeshElement>(
750 average_velocity_norm /
751 static_cast<double>(n_integration_points),
755 KCC.noalias() += KCC_Laplacian;
759 Eigen::Ref<LocalBlockMatrixType> KCmCn,
760 double const stoichiometric_coefficient,
761 double const kinetic_prefactor)
763 unsigned const n_integration_points =
773 auto const& phase = medium.
phase(
"AqueousLiquid");
781 for (
unsigned ip(0); ip < n_integration_points; ++ip)
784 auto const& w = ip_data.integration_weight;
785 auto const& N = Ns[ip];
786 auto& porosity = ip_data.porosity;
795 auto const retardation_factor =
797 .template value<double>(vars, pos, t, dt);
800 .template value<double>(vars, pos, t, dt);
804 .template value<double>(vars, pos, t, dt);
806 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
807 kinetic_prefactor * retardation_factor *
808 porosity * density * N;
813 Eigen::VectorXd
const& local_x,
814 Eigen::VectorXd
const& local_x_prev,
815 int const process_id,
816 std::vector<double>& local_M_data,
817 std::vector<double>& local_K_data,
818 std::vector<double>& local_b_data)
override
823 local_M_data, local_K_data, local_b_data);
828 local_M_data, local_K_data,
835 local_M_data, local_K_data,
836 local_b_data, process_id);
842 Eigen::VectorXd
const& local_x,
843 Eigen::VectorXd
const& local_x_prev,
844 std::vector<double>& local_M_data,
845 std::vector<double>& local_K_data,
846 std::vector<double>& local_b_data)
850 auto const local_C = local_x.template segment<concentration_size>(
852 auto const local_C_prev =
864 unsigned const n_integration_points =
876 auto const& phase = medium.
phase(
"AqueousLiquid");
885 for (
unsigned ip(0); ip < n_integration_points; ++ip)
888 auto const& dNdx = ip_data.dNdx;
889 auto const& w = ip_data.integration_weight;
890 auto const& N = Ns[ip];
891 auto& porosity = ip_data.porosity;
892 auto const& porosity_prev = ip_data.porosity_prev;
894 double const C_int_pt = N.dot(local_C);
895 double const p_int_pt = N.dot(local_p);
896 double const T_int_pt = N.dot(local_T);
910 .template value<double>(vars, vars_prev, pos, t,
921 .template value<double>(vars, pos, t, dt);
929 .template value<double>(vars, pos, t, dt);
933 const double drho_dp =
935 .template dValue<double>(
939 const double drho_dC =
941 .template dValue<double>(
946 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
948 w * dNdx.transpose() * density * K_over_mu * dNdx;
953 w * density * density * dNdx.transpose() * K_over_mu * b;
958 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
961 w * N.transpose() * porosity * drho_dC * C_dot;
967 Eigen::VectorXd
const& local_x,
968 Eigen::VectorXd
const& ,
969 std::vector<double>& local_M_data,
970 std::vector<double>& local_K_data,
971 std::vector<double>& )
973 assert(local_x.size() ==
992 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
1000 unsigned const n_integration_points =
1003 std::vector<GlobalDimVectorType> ip_flux_vector;
1004 double average_velocity_norm = 0.0;
1005 ip_flux_vector.reserve(n_integration_points);
1011 for (
unsigned ip(0); ip < n_integration_points; ip++)
1013 auto const& ip_data = this->
_ip_data[ip];
1014 auto const& dNdx = ip_data.dNdx;
1015 auto const& w = ip_data.integration_weight;
1016 auto const& N = Ns[ip];
1018 double p_at_xi = 0.;
1020 double T_at_xi = 0.;
1028 auto const porosity =
1030 .template value<double>(vars, pos, t, dt);
1034 auto const fluid_density =
1037 .template value<double>(vars, pos, t, dt);
1039 auto const specific_heat_capacity_fluid =
1042 .template value<double>(vars, pos, t, dt);
1045 local_M.noalias() += w *
1047 vars, porosity, fluid_density,
1048 specific_heat_capacity_fluid, pos, t, dt) *
1052 auto const viscosity =
1055 .template value<double>(vars, pos, t, dt);
1057 auto const intrinsic_permeability =
1062 .value(vars, pos, t, dt));
1065 intrinsic_permeability / viscosity;
1067 process_data.has_gravity
1069 (dNdx * local_p - fluid_density * b))
1074 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1077 local_K.noalias() +=
1078 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1080 ip_flux_vector.emplace_back(velocity * fluid_density *
1081 specific_heat_capacity_fluid);
1082 average_velocity_norm += velocity.norm();
1086 process_data.stabilizer, this->_ip_data,
1088 average_velocity_norm /
static_cast<double>(n_integration_points),
1093 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1094 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_M_data,
1095 std::vector<double>& local_K_data,
1096 std::vector<double>& ,
int const transport_process_id)
1098 assert(
static_cast<int>(local_x.size()) ==
1104 auto const local_p =
1109 auto const local_C = local_x.template segment<concentration_size>(
1113 auto const local_p_prev =
1124 unsigned const n_integration_points =
1127 std::vector<GlobalDimVectorType> ip_flux_vector;
1128 double average_velocity_norm = 0.0;
1131 ip_flux_vector.reserve(n_integration_points);
1144 auto const& medium =
1146 auto const& phase = medium.
phase(
"AqueousLiquid");
1147 auto const component_id =
1149 auto const& component = phase.
component(
1156 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1159 auto const& dNdx = ip_data.dNdx;
1160 auto const& w = ip_data.integration_weight;
1161 auto const& N = Ns[ip];
1162 auto& porosity = ip_data.porosity;
1163 auto const& porosity_prev = ip_data.porosity_prev;
1165 double const C_int_pt = N.dot(local_C);
1166 double const p_int_pt = N.dot(local_p);
1167 double const T_int_pt = N.dot(local_T);
1180 vars_prev.
porosity = porosity_prev;
1186 .template value<double>(vars, vars_prev, pos, t,
1192 auto const& retardation_factor =
1194 .template value<double>(vars, pos, t, dt);
1196 auto const& solute_dispersivity_transverse = medium.template value<
1199 auto const& solute_dispersivity_longitudinal =
1200 medium.template value<double>(
1202 longitudinal_dispersivity);
1205 auto const density =
1207 .template value<double>(vars, pos, t, dt);
1208 auto const decay_rate =
1210 .template value<double>(vars, pos, t, dt);
1212 auto const& pore_diffusion_coefficient =
1215 .value(vars, pos, t, dt));
1222 .template value<double>(vars, pos, t, dt);
1228 (dNdx * local_p - density * b))
1234 pore_diffusion_coefficient, velocity, porosity,
1235 solute_dispersivity_transverse,
1236 solute_dispersivity_longitudinal);
1238 double const R_times_phi = retardation_factor * porosity;
1239 auto const N_t_N = (N.transpose() * N).eval();
1243 const double drho_dC =
1245 .template dValue<double>(
1248 local_M.noalias() +=
1249 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1252 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1257 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1259 const double drho_dp =
1261 .template dValue<double>(vars,
1263 liquid_phase_pressure,
1266 local_K.noalias() +=
1267 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1268 dNdx.transpose() * velocity * N * (density * w);
1272 ip_flux_vector.emplace_back(velocity * density);
1273 average_velocity_norm += velocity.norm();
1275 local_K.noalias() +=
1276 N_t_N * (decay_rate * R_times_phi * density * w);
1278 KCC_Laplacian.noalias() += dNdx.transpose() *
1279 hydrodynamic_dispersion * dNdx *
1286 typename ShapeFunction::MeshElement>(
1289 average_velocity_norm /
1290 static_cast<double>(n_integration_points),
1293 local_K.noalias() += KCC_Laplacian;
1297 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1298 Eigen::VectorXd
const& local_x_prev,
int const process_id,
1299 std::vector<double>& local_b_data,
1300 std::vector<double>& local_Jac_data)
override
1305 local_b_data, local_Jac_data);
1309 int const component_id = process_id - 1;
1311 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1317 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1318 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
1319 std::vector<double>& local_Jac_data)
1321 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1322 auto const c = local_x.template segment<concentration_size>(
1334 unsigned const n_integration_points =
1343 auto const& medium =
1345 auto const& phase = medium.
phase(
"AqueousLiquid");
1354 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1357 auto const& dNdx = ip_data.dNdx;
1358 auto const& w = ip_data.integration_weight;
1359 auto const& N = Ns[ip];
1360 auto& phi = ip_data.porosity;
1361 auto const& phi_prev = ip_data.porosity_prev;
1363 double const p_ip = N.dot(p);
1364 double const c_ip = N.dot(c);
1366 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1378 .template value<double>(vars, vars_prev, pos, t,
1385 .template value<double>(vars, pos, t, dt);
1392 .template value<double>(vars, pos, t, dt);
1394 auto const drho_dp =
1396 .template dValue<double>(
1400 auto const drho_dc =
1402 .template dValue<double>(
1407 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1408 w * dNdx.transpose() * rho * k / mu * dNdx;
1410 local_rhs.noalias() -=
1411 w * N.transpose() * phi *
1412 (drho_dp * N * p_prev + drho_dc * cdot_ip) +
1413 w * rho * dNdx.transpose() * k / mu * dNdx * p;
1417 local_rhs.noalias() +=
1418 w * rho * dNdx.transpose() * k / mu * rho * b;
1424 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1425 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
1426 std::vector<double>& local_Jac_data,
int const component_id)
1428 auto const concentration_index =
1431 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1433 local_x.template segment<concentration_size>(concentration_index);
1451 unsigned const n_integration_points =
1454 std::vector<GlobalDimVectorType> ip_flux_vector;
1455 double average_velocity_norm = 0.0;
1456 ip_flux_vector.reserve(n_integration_points);
1468 auto const& medium =
1470 auto const& phase = medium.
phase(
"AqueousLiquid");
1471 auto const& component = phase.
component(
1478 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1481 auto const& dNdx = ip_data.dNdx;
1482 auto const& w = ip_data.integration_weight;
1483 auto const& N = Ns[ip];
1484 auto& phi = ip_data.porosity;
1485 auto const& phi_prev = ip_data.porosity_prev;
1487 double const p_ip = N.dot(p);
1488 double const c_ip = N.dot(c);
1505 .template value<double>(vars, vars_prev, pos, t,
1513 .template value<double>(vars, pos, t, dt);
1515 auto const alpha_T = medium.template value<double>(
1517 auto const alpha_L = medium.template value<double>(
1521 .template value<double>(vars, pos, t, dt);
1525 .template value<double>(vars, pos, t, dt);
1529 .value(vars, pos, t, dt));
1535 .template value<double>(vars, pos, t, dt);
1547 local_Jac.noalias() +=
1548 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1550 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1552 auto const cdot = (c - c_prev) / dt;
1553 local_rhs.noalias() -=
1554 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1556 ip_flux_vector.emplace_back(q * rho);
1557 average_velocity_norm += q.norm();
1563 average_velocity_norm /
static_cast<double>(n_integration_points),
1566 local_rhs.noalias() -= KCC_Laplacian * c;
1568 local_Jac.noalias() += KCC_Laplacian;
1572 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1573 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1574 std::vector<double>& local_b_data,
1575 int const transport_process_id)
override
1577 auto const local_C = local_x.template segment<concentration_size>(
1588 unsigned const n_integration_points =
1597 auto const& medium =
1599 auto const component_id = transport_process_id - 1;
1605 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1608 auto const w = ip_data.integration_weight;
1609 auto const& N = Ns[ip];
1610 auto& porosity = ip_data.porosity;
1611 auto const& porosity_prev = ip_data.porosity_prev;
1612 auto const chemical_system_id = ip_data.chemical_system_id;
1614 double C_int_pt = 0.0;
1619 auto const porosity_dot = (porosity - porosity_prev) / dt;
1623 vars_prev.
porosity = porosity_prev;
1629 .template value<double>(vars, vars_prev, pos, t,
1633 local_M.noalias() += w * N.transpose() * porosity * N;
1635 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1637 if (chemical_system_id == -1)
1642 auto const C_post_int_pt =
1644 component_id, chemical_system_id);
1646 local_b.noalias() +=
1647 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1653 std::vector<GlobalVector*>
const& x,
1654 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
1655 std::vector<double>& cache)
const override
1657 assert(x.size() == dof_table.size());
1659 auto const n_processes = x.size();
1660 std::vector<std::vector<double>> local_x;
1661 local_x.reserve(n_processes);
1663 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1665 auto const indices =
1667 assert(!indices.empty());
1668 local_x.push_back(x[process_id]->get(indices));
1672 if (n_processes == 1)
1674 auto const local_p = Eigen::Map<const NodalVectorType>(
1676 auto const local_C = Eigen::Map<const NodalVectorType>(
1683 constexpr int pressure_process_id = 0;
1684 constexpr int concentration_process_id = 1;
1685 auto const local_p = Eigen::Map<const NodalVectorType>(
1687 auto const local_C = Eigen::Map<const NodalVectorType>(
1695 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
1696 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
1697 std::vector<double>& cache)
const
1699 auto const n_integration_points =
1704 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1705 cache, GlobalDim, n_integration_points);
1716 auto const& medium =
1718 auto const& phase = medium.
phase(
"AqueousLiquid");
1724 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1726 auto const& ip_data =
_ip_data[ip];
1727 auto const& dNdx = ip_data.dNdx;
1728 auto const& N = Ns[ip];
1729 auto const& porosity = ip_data.porosity;
1731 double C_int_pt = 0.0;
1732 double p_int_pt = 0.0;
1743 double const dt = std::numeric_limits<double>::quiet_NaN();
1748 .template value<double>(vars, pos, t, dt);
1751 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1756 .template value<double>(vars, pos, t, dt);
1758 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1766 const unsigned integration_point)
const override
1769 typename ShapeFunction::MeshElement>()[integration_point];
1772 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1777 std::vector<double>
const& local_x)
const override
1779 auto const local_p = Eigen::Map<const NodalVectorType>(
1781 auto const local_C = Eigen::Map<const NodalVectorType>(
1787 auto const shape_matrices =
1791 std::array{pnt_local_coords})[0];
1801 auto const& medium =
1803 auto const& phase = medium.
phase(
"AqueousLiquid");
1816 double const dt = std::numeric_limits<double>::quiet_NaN();
1822 .template value<double>(vars, pos, t, dt);
1827 .template value<double>(vars, pos, t, dt);
1830 q += K_over_mu * rho_w * b;
1832 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1833 flux.head<GlobalDim>() = rho_w * q;
1840 Eigen::VectorXd
const& local_x,
1841 Eigen::VectorXd
const& )
override
1843 auto const local_p =
1845 auto const local_C = local_x.template segment<concentration_size>(
1848 std::vector<double> ele_velocity;
1851 auto const n_integration_points =
1853 auto const ele_velocity_mat =
1857 Eigen::Map<LocalVectorType>(
1860 ele_velocity_mat.rowwise().sum() / n_integration_points;
1864 std::size_t
const ele_id)
override
1866 auto const n_integration_points =
1875 ip_data.porosity = ip_data.porosity_prev;
1879 medium, ip_data.porosity);
1884 [](
double const s,
auto const& ip)
1885 { return s + ip.porosity; }) /
1886 n_integration_points;
1889 std::vector<GlobalIndexType> chemical_system_indices;
1890 chemical_system_indices.reserve(n_integration_points);
1892 std::back_inserter(chemical_system_indices),
1893 [](
auto const& ip_data)
1894 { return ip_data.chemical_system_id; });
1897 ele_id, chemical_system_indices);
1901 const double t, std::vector<GlobalVector*>
const& x,
1902 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
1903 std::vector<double>& cache,
int const component_id)
const override
1905 std::vector<double> local_x_vec;
1907 auto const n_processes = x.size();
1908 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1910 auto const indices =
1912 assert(!indices.empty());
1913 auto const local_solution = x[process_id]->get(indices);
1914 local_x_vec.insert(std::end(local_x_vec),
1915 std::begin(local_solution),
1916 std::end(local_solution));
1920 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1921 auto const c = local_x.template segment<concentration_size>(
1924 auto const n_integration_points =
1929 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1930 cache, GlobalDim, n_integration_points);
1941 auto const& medium =
1943 auto const& phase = medium.
phase(
"AqueousLiquid");
1945 auto const& component = phase.
component(
1952 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1954 auto const& ip_data =
_ip_data[ip];
1955 auto const& dNdx = ip_data.dNdx;
1956 auto const& N = Ns[ip];
1957 auto const& phi = ip_data.porosity;
1959 double const p_ip = N.dot(p);
1960 double const c_ip = N.dot(c);
1966 double const dt = std::numeric_limits<double>::quiet_NaN();
1972 .template value<double>(vars, pos, t, dt);
1974 .template value<double>(vars, pos, t, dt);
1982 auto const alpha_T = medium.template value<double>(
1984 auto const alpha_L = medium.template value<double>(
1988 .value(vars, pos, t, dt));
1995 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
2002 Eigen::VectorXd
const& ,
2003 double const ,
double const ,
2004 int const )
override
2006 unsigned const n_integration_points =
2009 for (
unsigned ip = 0; ip < n_integration_points; ip++)
2020 std::vector<std::reference_wrapper<ProcessVariable>>
const
2023 std::vector<IntegrationPointData<GlobalDimNodalMatrixType>>
_ip_data;
2027 const double fluid_density,
const double specific_heat_capacity_fluid,
2031 auto const& medium =
2033 auto const& solid_phase = medium.
phase(
"Solid");
2035 auto const specific_heat_capacity_solid =
2039 .template value<double>(vars, pos, t, dt);
2041 auto const solid_density =
2043 .template value<double>(vars, pos, t, dt);
2045 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2046 fluid_density * specific_heat_capacity_fluid * porosity;
2051 const double fluid_density,
const double specific_heat_capacity_fluid,
2056 auto const& medium =
2059 auto thermal_conductivity =
2064 .value(vars, pos, t, dt));
2066 auto const thermal_dispersivity_transversal =
2069 thermal_transversal_dispersivity)
2070 .template value<double>();
2072 auto const thermal_dispersivity_longitudinal =
2075 thermal_longitudinal_dispersivity)
2076 .template value<double>();
2081 return thermal_conductivity +
2082 fluid_density * specific_heat_capacity_fluid *
2085 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2086 velocity, 0 , thermal_dispersivity_transversal,
2087 thermal_dispersivity_longitudinal);
2091 Eigen::VectorXd
const& local_x)
39namespace ComponentTransport {
…}
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
Property const & property(PropertyType const &p) const
Component const & component(std::size_t const &index) const
double liquid_phase_pressure
int add(IndexType row, IndexType col, double val)
Global vector based on Eigen vector.
void add(IndexType rowId, double v)
add entry
constexpr double getWeight() const
std::size_t getID() const
Returns the ID of the element.
MathLib::WeightedPoint const & getWeightedPoint(unsigned const igp) const
unsigned getNumberOfPoints() const
MathLib::RowColumnIndices< GlobalIndexType > RowColumnIndices
auto const & NsHigherOrder() const
void setCoordinates(MathLib::Point3d const &coordinates)
void setElementID(std::size_t element_id)
virtual std::vector< double > const & getIntPtDarcyVelocity(const double t, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_table, std::vector< double > &cache) const =0
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
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
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
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)
const int first_concentration_index
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
static const int temperature_size
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
const int temperature_index
static const int concentration_size
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
static const int pressure_size
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)
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
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
@ 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< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
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)
std::array< double, 3 > interpolateCoordinates(MeshLib::Element const &e, typename ShapeMatricesType::ShapeMatrices::ShapeType const &N)
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
NumLib::ShapeMatrixCache shape_matrix_cache
caches for each mesh element type the shape matrix
ParameterLib::Parameter< double > const *const temperature
const int thermal_process_id
bool const chemically_induced_porosity_change
MeshLib::PropertyVector< double > * mesh_prop_porosity
GlobalIndexType chemical_system_id
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
GlobalDimNodalMatrixType const dNdx
IntegrationPointData(GlobalDimNodalMatrixType const &dNdx_, double const &integration_weight_)
double const integration_weight