31{
32 auto const local_matrix_size = local_x.size();
33
34 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
35
37 local_M_data, local_matrix_size, local_matrix_size);
39 local_K_data, local_matrix_size, local_matrix_size);
41 local_b_data, local_matrix_size);
42
43 auto Mgp =
44 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
46 auto Mgx = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
48
49 auto Mlp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
51
52 auto Mlx = local_M.template block<cap_pressure_size, cap_pressure_size>(
54
56 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
57
58 auto Kgp =
59 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
61
62 auto Kgx = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
64
65 auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
67
68 auto Klx = local_K.template block<cap_pressure_size, cap_pressure_size>(
70
71 auto Bg = local_b.template segment<nonwet_pressure_size>(
73
74 auto Bl =
76
78 auto const& liquid_phase = medium.
phase(
"AqueousLiquid");
79 auto const& gas_phase = medium.phase("Gas");
80
81 unsigned const n_integration_points =
83
86 const int material_id =
88
89 for (unsigned ip = 0; ip < n_integration_points; ip++)
90 {
92
93 double pl_int_pt = 0.;
94 double totalrho_int_pt =
95 0.;
97 totalrho_int_pt);
98
100
102
106 gas_phase.property(MPL::PropertyType::molar_mass)
107 .template value<double>(variables, pos, t, dt);
108
110 medium.property(MPL::PropertyType::permeability)
111 .value(variables, pos, t, dt));
112 auto const rho_gas = gas_phase.property(MPL::PropertyType::density)
113 .template value<double>(variables, pos, t, dt);
114 auto const rho_h2o = liquid_phase.property(MPL::PropertyType::density)
115 .template value<double>(variables, pos, t, dt);
116
119 double const X_h2_nonwet = 1.0;
120 double& rho_h2_wet =
_ip_data[ip].rho_m;
121 double& dSwdP =
_ip_data[ip].dsw_dpg;
122 double& dSwdrho =
_ip_data[ip].dsw_drho;
123 double& drhoh2wet =
_ip_data[ip].drhom_dpg;
124 double& drhoh2wet_drho =
_ip_data[ip].drhom_drho;
125 if (!
_ip_data[ip].mat_property.computeConstitutiveRelation(
126 t,
127 pos,
128 material_id,
129 pl_int_pt,
130 totalrho_int_pt,
131 temperature,
132 Sw,
133 rho_h2_wet,
134 dSwdP,
135 dSwdrho,
136 drhoh2wet,
137 drhoh2wet_drho))
138 {
139 OGS_FATAL(
"Computation of local constitutive relation failed.");
140 }
142 material_id, t, pos, pl_int_pt, temperature, Sw);
145
146 double const rho_wet = rho_h2o + rho_h2_wet;
149 pl_int_pt + pc;
150
151
152
153 double dPC_dSw =
155 material_id, t, pos, pl_int_pt, temperature, Sw);
156
158 medium.property(MPL::PropertyType::porosity)
159 .template value<double>(variables, pos, t, dt);
160
162
164
165 Mlx.noalias() +=
167 double const k_rel_gas =
170 double const mu_gas =
171 gas_phase.property(MPL::PropertyType::viscosity)
172 .template value<double>(variables, pos, t, dt);
173 double const lambda_gas = k_rel_gas / mu_gas;
174 double const diffusion_coeff_component_h2 =
176
177
178 double const k_rel_wet =
180 t, pos, pl_int_pt, temperature, Sw);
181 auto const mu_wet = liquid_phase.property(MPL::PropertyType::viscosity)
182 .template value<double>(variables, pos, t, dt);
183 double const lambda_wet = k_rel_wet / mu_wet;
184
185 laplace_operator.noalias() = sm.dNdx.transpose() *
permeability *
186 sm.dNdx *
_ip_data[ip].integration_weight;
187
188 Kgp.noalias() +=
189 (rho_gas * X_h2_nonwet * lambda_gas * (1 + dPC_dSw * dSwdP) +
190 rho_h2_wet * lambda_wet) *
191 laplace_operator +
192 (Sw *
porosity * diffusion_coeff_component_h2 *
193 (rho_h2o / rho_wet) * drhoh2wet) *
195 Kgx.noalias() +=
196 (rho_gas * X_h2_nonwet * lambda_gas * dPC_dSw * dSwdrho) *
197 laplace_operator +
198 (Sw * porosity * diffusion_coeff_component_h2 *
199 (rho_h2o / rho_wet) * drhoh2wet_drho) *
201 Klp.noalias() += (rho_gas * lambda_gas * (1 + dPC_dSw * dSwdP) +
202 rho_wet * lambda_wet) *
203 laplace_operator;
204
205 Klx.noalias() +=
206 (rho_gas * lambda_gas * dPC_dSw * dSwdrho) * laplace_operator;
207
209 {
211 Bg.noalias() += (rho_gas * rho_gas * lambda_gas +
212 rho_h2_wet * rho_wet * lambda_wet) *
215 Bl.noalias() += (rho_wet * lambda_wet * rho_wet +
216 rho_gas * rho_gas * lambda_gas) *
219
220 }
221 }
223 {
224 for (unsigned row = 0; row < Mgp.cols(); row++)
225 {
226 for (unsigned column = 0; column < Mgp.cols(); column++)
227 {
228 if (row != column)
229 {
230 Mgx(row, row) += Mgx(row, column);
231 Mgx(row, column) = 0.0;
232 Mgp(row, row) += Mgp(row, column);
233 Mgp(row, column) = 0.0;
234 Mlx(row, row) += Mlx(row, column);
235 Mlx(row, column) = 0.0;
236 Mlp(row, row) += Mlp(row, column);
237 Mlp(row, column) = 0.0;
238 }
239 }
240 }
241 }
242}
Medium * getMedium(std::size_t element_id)
Phase const & phase(std::size_t index) const
double gas_phase_pressure
double capillary_pressure
double liquid_phase_pressure
std::size_t getID() const
Returns the ID of the element.
std::optional< std::size_t > getElementID() const
void setElementID(std::size_t element_id)
static const int cap_pressure_matrix_index
static const int nonwet_pressure_matrix_index
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
Eigen::Matrix< double, GlobalDim, GlobalDim > formEigenTensor(MaterialPropertyLib::PropertyDataType const &values)
Eigen::Map< Vector > createZeroedVector(std::vector< double > &data, Eigen::VectorXd::Index size)
Eigen::Map< Matrix > createZeroedMatrix(std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)
Eigen::VectorXd const _specific_body_force
ParameterLib::Parameter< double > const & _temperature
bool const _has_mass_lumping
ParameterLib::Parameter< double > const & _diffusion_coeff_component_b
MaterialPropertyLib::MaterialSpatialDistributionMap media_map
Public Member Functions inherited from ProcessLib::LocalAssemblerInterface
1.11.0