21 namespace TwoPhaseFlowWithPrho
23 template <
typename ShapeFunction,
typename IntegrationMethod,
int GlobalDim>
24 void TwoPhaseFlowWithPrhoLocalAssembler<
25 ShapeFunction, IntegrationMethod,
26 GlobalDim>::assemble(
double const t,
double const ,
27 std::vector<double>
const& local_x,
28 std::vector<double>
const& ,
29 std::vector<double>& local_M_data,
30 std::vector<double>& local_K_data,
31 std::vector<double>& local_b_data)
33 auto const local_matrix_size = local_x.size();
35 assert(local_matrix_size == ShapeFunction::NPOINTS *
NUM_NODAL_DOF);
37 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
38 local_M_data, local_matrix_size, local_matrix_size);
39 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
40 local_K_data, local_matrix_size, local_matrix_size);
41 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
42 local_b_data, local_matrix_size);
45 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
46 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
47 auto Mgx = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
48 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
50 auto Mlp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
51 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
53 auto Mlx = local_M.template block<cap_pressure_size, cap_pressure_size>(
54 cap_pressure_matrix_index, cap_pressure_matrix_index);
57 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
60 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
61 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
63 auto Kgx = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
64 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
66 auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
67 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
69 auto Klx = local_K.template block<cap_pressure_size, cap_pressure_size>(
70 cap_pressure_matrix_index, cap_pressure_matrix_index);
72 auto Bg = local_b.template segment<nonwet_pressure_size>(
73 nonwet_pressure_matrix_index);
76 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
78 unsigned const n_integration_points =
79 _integration_method.getNumberOfPoints();
83 const int material_id =
84 _process_data._material->getMaterialID(pos.
getElementID().value());
86 const Eigen::MatrixXd& perm = _process_data._material->getPermeability(
87 material_id, t, pos, _element.getDimension());
88 assert(perm.rows() == _element.getDimension() || perm.rows() == 1);
90 _element.getDimension(), _element.getDimension());
91 if (perm.rows() == _element.getDimension())
95 else if (perm.rows() == 1)
100 for (
unsigned ip = 0; ip < n_integration_points; ip++)
102 auto const& sm = _shape_matrices[ip];
104 double pl_int_pt = 0.;
105 double totalrho_int_pt =
110 const double temperature = _process_data._temperature(t, pos)[0];
112 double const rho_gas =
113 _process_data._material->getGasDensity(pl_int_pt, temperature);
114 double const rho_h2o =
115 _process_data._material->getLiquidDensity(pl_int_pt, temperature);
117 double& Sw = _ip_data[ip].sw;
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(
139 OGS_FATAL(
"Computation of local constitutive relation failed.");
141 double const pc = _process_data._material->getCapillaryPressure(
142 material_id, t, pos, pl_int_pt, temperature, Sw);
144 double const rho_wet = rho_h2o + rho_h2_wet;
145 _saturation[ip] = Sw;
146 _pressure_nonwetting[ip] = pl_int_pt + pc;
151 _process_data._material->getCapillaryPressureDerivative(
152 material_id, t, pos, pl_int_pt, temperature, Sw);
154 double const porosity = _process_data._material->getPorosity(
155 material_id, t, pos, pl_int_pt, temperature, 0);
157 Mgx.noalias() += porosity * _ip_data[ip].massOperator;
159 Mlp.noalias() += porosity * rho_h2o * dSwdP * _ip_data[ip].massOperator;
162 porosity * (1 + dSwdrho * rho_h2o) * _ip_data[ip].massOperator;
163 double const k_rel_gas =
164 _process_data._material->getNonwetRelativePermeability(
165 t, pos, _pressure_nonwetting[ip], temperature, Sw);
166 double const mu_gas = _process_data._material->getGasViscosity(
167 _pressure_nonwetting[ip], temperature);
168 double const lambda_gas = k_rel_gas / mu_gas;
169 double const diffusion_coeff_component_h2 =
170 _process_data._diffusion_coeff_component_b(t, pos)[0];
173 double const k_rel_wet =
174 _process_data._material->getWetRelativePermeability(
175 t, pos, pl_int_pt, temperature, Sw);
176 double const mu_wet =
177 _process_data._material->getLiquidViscosity(pl_int_pt, temperature);
178 double const lambda_wet = k_rel_wet / mu_wet;
180 laplace_operator.noalias() = sm.dNdx.transpose() *
permeability *
181 sm.dNdx * _ip_data[ip].integration_weight;
184 (rho_gas * X_h2_nonwet * lambda_gas * (1 + dPC_dSw * dSwdP) +
185 rho_h2_wet * lambda_wet) *
187 (Sw * porosity * diffusion_coeff_component_h2 *
188 (rho_h2o / rho_wet) * drhoh2wet) *
189 _ip_data[ip].diffusionOperator;
191 (rho_gas * X_h2_nonwet * lambda_gas * dPC_dSw * dSwdrho) *
193 (Sw * porosity * diffusion_coeff_component_h2 *
194 (rho_h2o / rho_wet) * drhoh2wet_drho) *
195 _ip_data[ip].diffusionOperator;
196 Klp.noalias() += (rho_gas * lambda_gas * (1 + dPC_dSw * dSwdP) +
197 rho_wet * lambda_wet) *
201 (rho_gas * lambda_gas * dPC_dSw * dSwdrho) * laplace_operator;
203 if (_process_data._has_gravity)
205 auto const& b = _process_data._specific_body_force;
206 Bg.noalias() += (rho_gas * rho_gas * lambda_gas +
207 rho_h2_wet * rho_wet * lambda_wet) *
209 _ip_data[ip].integration_weight;
210 Bl.noalias() += (rho_wet * lambda_wet * rho_wet +
211 rho_gas * rho_gas * lambda_gas) *
213 _ip_data[ip].integration_weight;
217 if (_process_data._has_mass_lumping)
219 for (
unsigned row = 0; row < Mgp.cols(); row++)
221 for (
unsigned column = 0; column < Mgp.cols(); column++)
225 Mgx(row, row) += Mgx(row, column);
226 Mgx(row, column) = 0.0;
227 Mgp(row, row) += Mgp(row, column);
228 Mgp(row, column) = 0.0;
229 Mlx(row, row) += Mlx(row, column);
230 Mlx(row, column) = 0.0;
231 Mlp(row, row) += Mlp(row, column);
232 Mlp(row, column) = 0.0;
Definition of the PiecewiseLinearInterpolation class.
std::optional< std::size_t > getElementID() const
void setElementID(std::size_t element_id)
typename ShapeMatricesType::GlobalDimMatrixType GlobalDimMatrixType
typename ShapeMatricesType::NodalMatrixType NodalMatrixType
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
const unsigned NUM_NODAL_DOF