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;
338 std::vector<double> C_int_pt(n_component);
339 for (
unsigned component_id = 0; component_id < n_component;
342 auto const concentration_index =
346 local_x.template segment<concentration_size>(
347 concentration_index);
350 C_int_pt[component_id]);
360 double const t,
double dt)
override
377 unsigned const n_integration_points =
380 for (
unsigned ip = 0; ip < n_integration_points; ip++)
383 auto const& N = Ns[ip];
384 auto& porosity = ip_data.porosity;
385 auto const& porosity_prev = ip_data.porosity_prev;
386 auto const& chemical_system_id = ip_data.chemical_system_id;
389 std::vector<double> C_int_pt(n_component);
390 for (
unsigned component_id = 0; component_id < n_component;
393 auto const concentration_index =
397 local_x.template segment<concentration_size>(
398 concentration_index);
401 C_int_pt[component_id]);
413 .template value<double>(vars, vars_prev, pos, t,
420 C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
425 double const dt)
override
439 ip_data.porosity = ip_data.porosity_prev;
444 ip_data.porosity, t, dt);
447 ip_data.chemical_system_id, medium, ip_data.porosity);
452 std::vector<double>
const& local_x,
453 std::vector<double>
const& ,
454 std::vector<double>& local_M_data,
455 std::vector<double>& local_K_data,
456 std::vector<double>& local_b_data)
override
458 auto const local_matrix_size = local_x.size();
463 assert(local_matrix_size == ShapeFunction::NPOINTS * num_nodal_dof);
466 local_M_data, local_matrix_size, local_matrix_size);
468 local_K_data, local_matrix_size, local_matrix_size);
470 local_b_data, local_matrix_size);
473 auto Kpp = local_K.template block<pressure_size, pressure_size>(
475 auto Mpp = local_M.template block<pressure_size, pressure_size>(
477 auto Bp = local_b.template segment<pressure_size>(
pressure_index);
479 auto local_p = Eigen::Map<const NodalVectorType>(
486 auto const number_of_components = num_nodal_dof - 1;
487 for (
int component_id = 0; component_id < number_of_components;
499 auto concentration_index =
503 local_K.template block<concentration_size, concentration_size>(
504 concentration_index, concentration_index);
506 local_M.template block<concentration_size, concentration_size>(
507 concentration_index, concentration_index);
509 local_M.template block<concentration_size, pressure_size>(
512 local_M.template block<pressure_size, concentration_size>(
515 auto local_C = Eigen::Map<const NodalVectorType>(
519 MCC, MCp, MpC, Kpp, Mpp, Bp);
523 auto const stoichiometric_matrix =
527 assert(stoichiometric_matrix);
529 for (Eigen::SparseMatrix<double>::InnerIterator it(
530 *stoichiometric_matrix, component_id);
534 auto const stoichiometric_coefficient = it.value();
535 auto const coupled_component_id = it.row();
536 auto const kinetic_prefactor =
540 auto const concentration_index =
542 auto const coupled_concentration_index =
547 concentration_index, coupled_concentration_index);
551 stoichiometric_coefficient,
561 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
562 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
563 Eigen::Ref<LocalBlockMatrixType> KCC,
564 Eigen::Ref<LocalBlockMatrixType> MCC,
565 Eigen::Ref<LocalBlockMatrixType> MCp,
566 Eigen::Ref<LocalBlockMatrixType> MpC,
567 Eigen::Ref<LocalBlockMatrixType> Kpp,
568 Eigen::Ref<LocalBlockMatrixType> Mpp,
569 Eigen::Ref<LocalSegmentVectorType> Bp)
571 unsigned const n_integration_points =
583 auto const& phase = medium.
phase(
"AqueousLiquid");
594 std::vector<GlobalDimVectorType> ip_flux_vector;
595 double average_velocity_norm = 0.0;
598 ip_flux_vector.reserve(n_integration_points);
605 for (
unsigned ip(0); ip < n_integration_points; ++ip)
608 auto const& dNdx = ip_data.dNdx;
609 auto const& N = Ns[ip];
610 auto const& w = ip_data.integration_weight;
611 auto& porosity = ip_data.porosity;
613 double C_int_pt = 0.0;
614 double p_int_pt = 0.0;
624 .template value<double>(vars, pos, t, dt);
627 auto const& retardation_factor =
629 .template value<double>(vars, pos, t, dt);
631 auto const& solute_dispersivity_transverse = medium.template value<
635 auto const& solute_dispersivity_longitudinal =
636 medium.template value<double>(
638 longitudinal_dispersivity);
645 .template value<double>(vars, pos, t, dt);
647 auto const decay_rate =
649 .template value<double>(vars, pos, t, dt);
651 auto const& pore_diffusion_coefficient =
654 .value(vars, pos, t, dt));
662 .template value<double>(vars, pos, t, dt);
668 (dNdx * p_nodal_values - density * b))
671 const double drho_dp =
673 .template dValue<double>(
678 const double drho_dC =
680 .template dValue<double>(
687 pore_diffusion_coefficient, velocity, porosity,
688 solute_dispersivity_transverse,
689 solute_dispersivity_longitudinal);
691 const double R_times_phi(retardation_factor * porosity);
693 auto const N_t_N = (N.transpose() * N).eval();
697 MCp.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dp * w);
698 MCC.noalias() += N_t_N * (C_int_pt * R_times_phi * drho_dC * w);
699 KCC.noalias() -= dNdx.transpose() * mass_density_flow * N * w;
703 ip_flux_vector.emplace_back(mass_density_flow);
704 average_velocity_norm += velocity.norm();
706 MCC.noalias() += N_t_N * (R_times_phi * density * w);
707 KCC.noalias() += N_t_N * (decay_rate * R_times_phi * density * w);
708 KCC_Laplacian.noalias() +=
709 dNdx.transpose() * hydrodynamic_dispersion * dNdx * density * w;
711 MpC.noalias() += N_t_N * (porosity * drho_dC * w);
714 if (component_id == 0)
716 Mpp.noalias() += N_t_N * (porosity * drho_dp * w);
718 dNdx.transpose() * K_over_mu * dNdx * (density * w);
722 Bp.noalias() += dNdx.transpose() * K_over_mu * b *
723 (density * density * w);
731 typename ShapeFunction::MeshElement>(
736 average_velocity_norm /
737 static_cast<double>(n_integration_points),
741 KCC.noalias() += KCC_Laplacian;
745 Eigen::Ref<LocalBlockMatrixType> KCmCn,
746 double const stoichiometric_coefficient,
747 double const kinetic_prefactor)
749 unsigned const n_integration_points =
759 auto const& phase = medium.
phase(
"AqueousLiquid");
767 for (
unsigned ip(0); ip < n_integration_points; ++ip)
770 auto const& w = ip_data.integration_weight;
771 auto const& N = Ns[ip];
772 auto& porosity = ip_data.porosity;
774 auto const retardation_factor =
776 .template value<double>(vars, pos, t, dt);
779 .template value<double>(vars, pos, t, dt);
783 .template value<double>(vars, pos, t, dt);
785 KCmCn.noalias() -= w * N.transpose() * stoichiometric_coefficient *
786 kinetic_prefactor * retardation_factor *
787 porosity * density * N;
792 Eigen::VectorXd
const& local_x,
793 Eigen::VectorXd
const& local_x_prev,
794 int const process_id,
795 std::vector<double>& local_M_data,
796 std::vector<double>& local_K_data,
797 std::vector<double>& local_b_data)
override
802 local_M_data, local_K_data, local_b_data);
807 local_M_data, local_K_data,
814 local_M_data, local_K_data,
815 local_b_data, process_id);
821 Eigen::VectorXd
const& local_x,
822 Eigen::VectorXd
const& local_x_prev,
823 std::vector<double>& local_M_data,
824 std::vector<double>& local_K_data,
825 std::vector<double>& local_b_data)
829 auto const local_C = local_x.template segment<concentration_size>(
831 auto const local_C_prev =
843 unsigned const n_integration_points =
855 auto const& phase = medium.
phase(
"AqueousLiquid");
864 for (
unsigned ip(0); ip < n_integration_points; ++ip)
867 auto const& dNdx = ip_data.dNdx;
868 auto const& w = ip_data.integration_weight;
869 auto const& N = Ns[ip];
870 auto& porosity = ip_data.porosity;
871 auto const& porosity_prev = ip_data.porosity_prev;
873 double const C_int_pt = N.dot(local_C);
874 double const p_int_pt = N.dot(local_p);
875 double const T_int_pt = N.dot(local_T);
889 .template value<double>(vars, vars_prev, pos, t,
900 .template value<double>(vars, pos, t, dt);
908 .template value<double>(vars, pos, t, dt);
912 const double drho_dp =
914 .template dValue<double>(
918 const double drho_dC =
920 .template dValue<double>(
925 local_M.noalias() += w * N.transpose() * porosity * drho_dp * N;
927 w * dNdx.transpose() * density * K_over_mu * dNdx;
932 w * density * density * dNdx.transpose() * K_over_mu * b;
937 double const C_dot = (C_int_pt - N.dot(local_C_prev)) / dt;
940 w * N.transpose() * porosity * drho_dC * C_dot;
946 Eigen::VectorXd
const& local_x,
947 Eigen::VectorXd
const& ,
948 std::vector<double>& local_M_data,
949 std::vector<double>& local_K_data,
950 std::vector<double>& )
952 assert(local_x.size() ==
971 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
979 unsigned const n_integration_points =
982 std::vector<GlobalDimVectorType> ip_flux_vector;
983 double average_velocity_norm = 0.0;
984 ip_flux_vector.reserve(n_integration_points);
990 for (
unsigned ip(0); ip < n_integration_points; ip++)
992 auto const& ip_data = this->
_ip_data[ip];
993 auto const& dNdx = ip_data.dNdx;
994 auto const& w = ip_data.integration_weight;
995 auto const& N = Ns[ip];
1007 auto const porosity =
1009 .template value<double>(vars, pos, t, dt);
1013 auto const fluid_density =
1016 .template value<double>(vars, pos, t, dt);
1018 auto const specific_heat_capacity_fluid =
1021 .template value<double>(vars, pos, t, dt);
1024 local_M.noalias() += w *
1026 vars, porosity, fluid_density,
1027 specific_heat_capacity_fluid, pos, t, dt) *
1031 auto const viscosity =
1034 .template value<double>(vars, pos, t, dt);
1036 auto const intrinsic_permeability =
1041 .value(vars, pos, t, dt));
1044 intrinsic_permeability / viscosity;
1046 process_data.has_gravity
1048 (dNdx * local_p - fluid_density * b))
1053 vars, fluid_density, specific_heat_capacity_fluid, velocity,
1056 local_K.noalias() +=
1057 w * dNdx.transpose() * thermal_conductivity_dispersivity * dNdx;
1059 ip_flux_vector.emplace_back(velocity * fluid_density *
1060 specific_heat_capacity_fluid);
1061 average_velocity_norm += velocity.norm();
1065 process_data.stabilizer, this->_ip_data,
1067 average_velocity_norm /
static_cast<double>(n_integration_points),
1072 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1073 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_M_data,
1074 std::vector<double>& local_K_data,
1075 std::vector<double>& ,
int const transport_process_id)
1077 assert(
static_cast<int>(local_x.size()) ==
1083 auto const local_p =
1088 auto const local_C = local_x.template segment<concentration_size>(
1092 auto const local_p_prev =
1103 unsigned const n_integration_points =
1106 std::vector<GlobalDimVectorType> ip_flux_vector;
1107 double average_velocity_norm = 0.0;
1110 ip_flux_vector.reserve(n_integration_points);
1123 auto const& medium =
1125 auto const& phase = medium.
phase(
"AqueousLiquid");
1126 auto const component_id =
1128 auto const& component = phase.
component(
1135 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1138 auto const& dNdx = ip_data.dNdx;
1139 auto const& w = ip_data.integration_weight;
1140 auto const& N = Ns[ip];
1141 auto& porosity = ip_data.porosity;
1142 auto const& porosity_prev = ip_data.porosity_prev;
1144 double const C_int_pt = N.dot(local_C);
1145 double const p_int_pt = N.dot(local_p);
1146 double const T_int_pt = N.dot(local_T);
1159 vars_prev.
porosity = porosity_prev;
1165 .template value<double>(vars, vars_prev, pos, t,
1171 auto const& retardation_factor =
1173 .template value<double>(vars, pos, t, dt);
1175 auto const& solute_dispersivity_transverse = medium.template value<
1178 auto const& solute_dispersivity_longitudinal =
1179 medium.template value<double>(
1181 longitudinal_dispersivity);
1184 auto const density =
1186 .template value<double>(vars, pos, t, dt);
1187 auto const decay_rate =
1189 .template value<double>(vars, pos, t, dt);
1191 auto const& pore_diffusion_coefficient =
1194 .value(vars, pos, t, dt));
1201 .template value<double>(vars, pos, t, dt);
1207 (dNdx * local_p - density * b))
1213 pore_diffusion_coefficient, velocity, porosity,
1214 solute_dispersivity_transverse,
1215 solute_dispersivity_longitudinal);
1217 double const R_times_phi = retardation_factor * porosity;
1218 auto const N_t_N = (N.transpose() * N).eval();
1222 const double drho_dC =
1224 .template dValue<double>(
1227 local_M.noalias() +=
1228 N_t_N * (R_times_phi * C_int_pt * drho_dC * w);
1231 local_M.noalias() += N_t_N * (R_times_phi * density * w);
1236 double const p_dot = (p_int_pt - N.dot(local_p_prev)) / dt;
1238 const double drho_dp =
1240 .template dValue<double>(vars,
1242 liquid_phase_pressure,
1245 local_K.noalias() +=
1246 N_t_N * ((R_times_phi * drho_dp * p_dot) * w) -
1247 dNdx.transpose() * velocity * N * (density * w);
1251 ip_flux_vector.emplace_back(velocity * density);
1252 average_velocity_norm += velocity.norm();
1254 local_K.noalias() +=
1255 N_t_N * (decay_rate * R_times_phi * density * w);
1257 KCC_Laplacian.noalias() += dNdx.transpose() *
1258 hydrodynamic_dispersion * dNdx *
1265 typename ShapeFunction::MeshElement>(
1268 average_velocity_norm /
1269 static_cast<double>(n_integration_points),
1272 local_K.noalias() += KCC_Laplacian;
1276 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1277 Eigen::VectorXd
const& local_x_prev,
int const process_id,
1278 std::vector<double>& local_b_data,
1279 std::vector<double>& local_Jac_data)
override
1284 local_b_data, local_Jac_data);
1288 int const component_id = process_id - 1;
1290 t, dt, local_x, local_x_prev, local_b_data, local_Jac_data,
1296 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1297 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
1298 std::vector<double>& local_Jac_data)
1300 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1301 auto const c = local_x.template segment<concentration_size>(
1313 unsigned const n_integration_points =
1322 auto const& medium =
1324 auto const& phase = medium.
phase(
"AqueousLiquid");
1333 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1336 auto const& dNdx = ip_data.dNdx;
1337 auto const& w = ip_data.integration_weight;
1338 auto const& N = Ns[ip];
1339 auto& phi = ip_data.porosity;
1340 auto const& phi_prev = ip_data.porosity_prev;
1342 double const p_ip = N.dot(p);
1343 double const c_ip = N.dot(c);
1345 double const cdot_ip = (c_ip - N.dot(c_prev)) / dt;
1357 .template value<double>(vars, vars_prev, pos, t,
1364 .template value<double>(vars, pos, t, dt);
1371 .template value<double>(vars, pos, t, dt);
1373 auto const drho_dp =
1375 .template dValue<double>(
1379 auto const drho_dc =
1381 .template dValue<double>(
1386 local_Jac.noalias() += w * N.transpose() * phi * drho_dp / dt * N +
1387 w * dNdx.transpose() * rho * k / mu * dNdx;
1389 local_rhs.noalias() -=
1390 w * N.transpose() * phi *
1391 (drho_dp * N * p_prev + drho_dc * cdot_ip) +
1392 w * rho * dNdx.transpose() * k / mu * dNdx * p;
1396 local_rhs.noalias() +=
1397 w * rho * dNdx.transpose() * k / mu * rho * b;
1403 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1404 Eigen::VectorXd
const& local_x_prev, std::vector<double>& local_b_data,
1405 std::vector<double>& local_Jac_data,
int const component_id)
1407 auto const concentration_index =
1410 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1412 local_x.template segment<concentration_size>(concentration_index);
1430 unsigned const n_integration_points =
1433 std::vector<GlobalDimVectorType> ip_flux_vector;
1434 double average_velocity_norm = 0.0;
1435 ip_flux_vector.reserve(n_integration_points);
1447 auto const& medium =
1449 auto const& phase = medium.
phase(
"AqueousLiquid");
1450 auto const& component = phase.
component(
1457 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1460 auto const& dNdx = ip_data.dNdx;
1461 auto const& w = ip_data.integration_weight;
1462 auto const& N = Ns[ip];
1463 auto& phi = ip_data.porosity;
1464 auto const& phi_prev = ip_data.porosity_prev;
1466 double const p_ip = N.dot(p);
1467 double const c_ip = N.dot(c);
1484 .template value<double>(vars, vars_prev, pos, t,
1492 .template value<double>(vars, pos, t, dt);
1494 auto const alpha_T = medium.template value<double>(
1496 auto const alpha_L = medium.template value<double>(
1500 .template value<double>(vars, pos, t, dt);
1504 .template value<double>(vars, pos, t, dt);
1508 .value(vars, pos, t, dt));
1514 .template value<double>(vars, pos, t, dt);
1526 local_Jac.noalias() +=
1527 w * rho * N.transpose() * phi * R * (alpha + 1 / dt) * N;
1529 KCC_Laplacian.noalias() += w * rho * dNdx.transpose() * D * dNdx;
1531 auto const cdot = (c - c_prev) / dt;
1532 local_rhs.noalias() -=
1533 w * rho * N.transpose() * phi * R * N * (cdot + alpha * c);
1535 ip_flux_vector.emplace_back(q * rho);
1536 average_velocity_norm += q.norm();
1542 average_velocity_norm /
static_cast<double>(n_integration_points),
1545 local_rhs.noalias() -= KCC_Laplacian * c;
1547 local_Jac.noalias() += KCC_Laplacian;
1551 double const t,
double const dt, Eigen::VectorXd
const& local_x,
1552 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
1553 std::vector<double>& local_b_data,
1554 int const transport_process_id)
override
1556 auto const local_C = local_x.template segment<concentration_size>(
1567 unsigned const n_integration_points =
1576 auto const& medium =
1578 auto const component_id = transport_process_id - 1;
1584 for (
unsigned ip(0); ip < n_integration_points; ++ip)
1587 auto const w = ip_data.integration_weight;
1588 auto const& N = Ns[ip];
1589 auto& porosity = ip_data.porosity;
1590 auto const& porosity_prev = ip_data.porosity_prev;
1591 auto const chemical_system_id = ip_data.chemical_system_id;
1593 double C_int_pt = 0.0;
1598 auto const porosity_dot = (porosity - porosity_prev) / dt;
1602 vars_prev.
porosity = porosity_prev;
1608 .template value<double>(vars, vars_prev, pos, t,
1612 local_M.noalias() += w * N.transpose() * porosity * N;
1614 local_K.noalias() += w * N.transpose() * porosity_dot * N;
1616 if (chemical_system_id == -1)
1621 auto const C_post_int_pt =
1623 component_id, chemical_system_id);
1625 local_b.noalias() +=
1626 w * N.transpose() * porosity * (C_post_int_pt - C_int_pt) / dt;
1632 std::vector<GlobalVector*>
const& x,
1633 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_table,
1634 std::vector<double>& cache)
const override
1636 assert(x.size() == dof_table.size());
1638 auto const n_processes = x.size();
1639 std::vector<std::vector<double>> local_x;
1640 local_x.reserve(n_processes);
1642 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1644 auto const indices =
1646 assert(!indices.empty());
1647 local_x.push_back(x[process_id]->get(indices));
1651 if (n_processes == 1)
1653 auto const local_p = Eigen::Map<const NodalVectorType>(
1655 auto const local_C = Eigen::Map<const NodalVectorType>(
1662 constexpr int pressure_process_id = 0;
1663 constexpr int concentration_process_id = 1;
1664 auto const local_p = Eigen::Map<const NodalVectorType>(
1666 auto const local_C = Eigen::Map<const NodalVectorType>(
1674 Eigen::Ref<const NodalVectorType>
const& p_nodal_values,
1675 Eigen::Ref<const NodalVectorType>
const& C_nodal_values,
1676 std::vector<double>& cache)
const
1678 auto const n_integration_points =
1683 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1684 cache, GlobalDim, n_integration_points);
1695 auto const& medium =
1697 auto const& phase = medium.
phase(
"AqueousLiquid");
1703 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1705 auto const& ip_data =
_ip_data[ip];
1706 auto const& dNdx = ip_data.dNdx;
1707 auto const& N = Ns[ip];
1708 auto const& porosity = ip_data.porosity;
1710 double C_int_pt = 0.0;
1711 double p_int_pt = 0.0;
1722 double const dt = std::numeric_limits<double>::quiet_NaN();
1727 .template value<double>(vars, pos, t, dt);
1730 cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
1735 .template value<double>(vars, pos, t, dt);
1737 cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
1745 const unsigned integration_point)
const override
1748 typename ShapeFunction::MeshElement>()[integration_point];
1751 return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
1756 std::vector<double>
const& local_x)
const override
1758 auto const local_p = Eigen::Map<const NodalVectorType>(
1760 auto const local_C = Eigen::Map<const NodalVectorType>(
1766 auto const shape_matrices =
1770 std::array{pnt_local_coords})[0];
1780 auto const& medium =
1782 auto const& phase = medium.
phase(
"AqueousLiquid");
1795 double const dt = std::numeric_limits<double>::quiet_NaN();
1801 .template value<double>(vars, pos, t, dt);
1806 .template value<double>(vars, pos, t, dt);
1809 q += K_over_mu * rho_w * b;
1811 Eigen::Vector3d flux(0.0, 0.0, 0.0);
1812 flux.head<GlobalDim>() = rho_w * q;
1819 Eigen::VectorXd
const& local_x,
1820 Eigen::VectorXd
const& )
override
1822 auto const local_p =
1824 auto const local_C = local_x.template segment<concentration_size>(
1827 std::vector<double> ele_velocity;
1830 auto const n_integration_points =
1832 auto const ele_velocity_mat =
1836 Eigen::Map<LocalVectorType>(
1839 ele_velocity_mat.rowwise().sum() / n_integration_points;
1843 std::size_t
const ele_id)
override
1845 auto const n_integration_points =
1854 ip_data.porosity = ip_data.porosity_prev;
1858 medium, ip_data.porosity);
1863 [](
double const s,
auto const& ip)
1864 { return s + ip.porosity; }) /
1865 n_integration_points;
1868 std::vector<GlobalIndexType> chemical_system_indices;
1869 chemical_system_indices.reserve(n_integration_points);
1871 std::back_inserter(chemical_system_indices),
1872 [](
auto const& ip_data)
1873 { return ip_data.chemical_system_id; });
1876 ele_id, chemical_system_indices);
1880 const double t, std::vector<GlobalVector*>
const& x,
1881 std::vector<NumLib::LocalToGlobalIndexMap const*>
const& dof_tables,
1882 std::vector<double>& cache,
int const component_id)
const override
1884 std::vector<double> local_x_vec;
1886 auto const n_processes = x.size();
1887 for (std::size_t process_id = 0; process_id < n_processes; ++process_id)
1889 auto const indices =
1891 assert(!indices.empty());
1892 auto const local_solution = x[process_id]->get(indices);
1893 local_x_vec.insert(std::end(local_x_vec),
1894 std::begin(local_solution),
1895 std::end(local_solution));
1899 auto const p = local_x.template segment<pressure_size>(
pressure_index);
1900 auto const c = local_x.template segment<concentration_size>(
1903 auto const n_integration_points =
1908 Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
1909 cache, GlobalDim, n_integration_points);
1920 auto const& medium =
1922 auto const& phase = medium.
phase(
"AqueousLiquid");
1924 auto const& component = phase.
component(
1931 for (
unsigned ip = 0; ip < n_integration_points; ++ip)
1933 auto const& ip_data =
_ip_data[ip];
1934 auto const& dNdx = ip_data.dNdx;
1935 auto const& N = Ns[ip];
1936 auto const& phi = ip_data.porosity;
1938 double const p_ip = N.dot(p);
1939 double const c_ip = N.dot(c);
1945 double const dt = std::numeric_limits<double>::quiet_NaN();
1951 .template value<double>(vars, pos, t, dt);
1953 .template value<double>(vars, pos, t, dt);
1961 auto const alpha_T = medium.template value<double>(
1963 auto const alpha_L = medium.template value<double>(
1967 .value(vars, pos, t, dt));
1974 cache_mat.col(ip).noalias() = q * c_ip - D * dNdx * c;
1981 Eigen::VectorXd
const& ,
1982 double const ,
double const ,
1983 int const )
override
1985 unsigned const n_integration_points =
1988 for (
unsigned ip = 0; ip < n_integration_points; ip++)
1999 std::vector<std::reference_wrapper<ProcessVariable>>
const
2002 std::vector<IntegrationPointData<GlobalDimNodalMatrixType>>
_ip_data;
2006 const double fluid_density,
const double specific_heat_capacity_fluid,
2010 auto const& medium =
2012 auto const& solid_phase = medium.
phase(
"Solid");
2014 auto const specific_heat_capacity_solid =
2018 .template value<double>(vars, pos, t, dt);
2020 auto const solid_density =
2022 .template value<double>(vars, pos, t, dt);
2024 return solid_density * specific_heat_capacity_solid * (1 - porosity) +
2025 fluid_density * specific_heat_capacity_fluid * porosity;
2030 const double fluid_density,
const double specific_heat_capacity_fluid,
2035 auto const& medium =
2038 auto thermal_conductivity =
2043 .value(vars, pos, t, dt));
2045 auto const thermal_dispersivity_transversal =
2048 thermal_transversal_dispersivity)
2049 .template value<double>();
2051 auto const thermal_dispersivity_longitudinal =
2054 thermal_longitudinal_dispersivity)
2055 .template value<double>();
2060 return thermal_conductivity +
2061 fluid_density * specific_heat_capacity_fluid *
2064 GlobalDimMatrixType::Zero(GlobalDim, GlobalDim),
2065 velocity, 0 , thermal_dispersivity_transversal,
2066 thermal_dispersivity_longitudinal);
2070 Eigen::VectorXd
const& local_x)
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 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)
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