OGS
TwoPhaseFlowWithPrhoLocalAssembler-impl.h
Go to the documentation of this file.
1
11#pragma once
12
18
19namespace ProcessLib
20{
21namespace TwoPhaseFlowWithPrho
22{
23namespace MPL = MaterialPropertyLib;
24
25template <typename ShapeFunction, int GlobalDim>
27 double const t, double const dt, std::vector<double> const& local_x,
28 std::vector<double> const& /*local_x_prev*/,
29 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
30 std::vector<double>& local_b_data)
31{
32 auto const local_matrix_size = local_x.size();
33
34 assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
35
36 auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
37 local_M_data, local_matrix_size, local_matrix_size);
38 auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
39 local_K_data, local_matrix_size, local_matrix_size);
40 auto local_b = MathLib::createZeroedVector<LocalVectorType>(
41 local_b_data, local_matrix_size);
42
43 auto Mgp =
44 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
45 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
46 auto Mgx = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
47 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
48
49 auto Mlp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
50 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
51
52 auto Mlx = local_M.template block<cap_pressure_size, cap_pressure_size>(
53 cap_pressure_matrix_index, cap_pressure_matrix_index);
54
55 NodalMatrixType laplace_operator =
56 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
57
58 auto Kgp =
59 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
60 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
61
62 auto Kgx = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
63 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
64
65 auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
66 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
67
68 auto Klx = local_K.template block<cap_pressure_size, cap_pressure_size>(
69 cap_pressure_matrix_index, cap_pressure_matrix_index);
70
71 auto Bg = local_b.template segment<nonwet_pressure_size>(
72 nonwet_pressure_matrix_index);
73
74 auto Bl =
75 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
76
77 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
78 auto const& liquid_phase = medium.phase("AqueousLiquid");
79 auto const& gas_phase = medium.phase("Gas");
80
81 unsigned const n_integration_points =
82 _integration_method.getNumberOfPoints();
83
85 pos.setElementID(_element.getID());
86 const int material_id =
87 _process_data._material->getMaterialID(pos.getElementID().value());
88
89 for (unsigned ip = 0; ip < n_integration_points; ip++)
90 {
91 auto const& sm = _shape_matrices[ip];
92
93 double pl_int_pt = 0.;
94 double totalrho_int_pt =
95 0.; // total mass density of the light component
96 NumLib::shapeFunctionInterpolate(local_x, sm.N, pl_int_pt,
97 totalrho_int_pt);
98
99 const double temperature = _process_data._temperature(t, pos)[0];
100
101 MPL::VariableArray variables;
102
103 variables.liquid_phase_pressure = pl_int_pt;
104 variables.temperature = temperature;
105 variables.molar_mass =
106 gas_phase.property(MPL::PropertyType::molar_mass)
107 .template value<double>(variables, pos, t, dt);
108
109 auto const permeability = MPL::formEigenTensor<GlobalDim>(
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
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(
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 }
141 double const pc = _process_data._material->getCapillaryPressure(
142 material_id, t, pos, pl_int_pt, temperature, Sw);
143 variables.capillary_pressure = pc;
144 variables.liquid_saturation = Sw;
145
146 double const rho_wet = rho_h2o + rho_h2_wet;
147 _saturation[ip] = Sw;
148 variables.gas_phase_pressure = _pressure_nonwetting[ip] =
149 pl_int_pt + pc;
150
151 // Assemble M matrix
152 // nonwetting
153 double dPC_dSw =
154 _process_data._material->getCapillaryPressureDerivative(
155 material_id, t, pos, pl_int_pt, temperature, Sw);
156
157 double const porosity =
158 medium.property(MPL::PropertyType::porosity)
159 .template value<double>(variables, pos, t, dt);
160
161 Mgx.noalias() += porosity * _ip_data[ip].massOperator;
162
163 Mlp.noalias() += porosity * rho_h2o * dSwdP * _ip_data[ip].massOperator;
164
165 Mlx.noalias() +=
166 porosity * (1 + dSwdrho * rho_h2o) * _ip_data[ip].massOperator;
167 double const k_rel_gas =
168 _process_data._material->getNonwetRelativePermeability(
169 t, pos, _pressure_nonwetting[ip], temperature, Sw);
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 =
175 _process_data._diffusion_coeff_component_b(t, pos)[0];
176
177 // wet
178 double const k_rel_wet =
179 _process_data._material->getWetRelativePermeability(
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) *
194 _ip_data[ip].diffusionOperator;
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) *
200 _ip_data[ip].diffusionOperator;
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
208 if (_process_data._has_gravity)
209 {
210 auto const& b = _process_data._specific_body_force;
211 Bg.noalias() += (rho_gas * rho_gas * lambda_gas +
212 rho_h2_wet * rho_wet * lambda_wet) *
213 sm.dNdx.transpose() * permeability * b *
214 _ip_data[ip].integration_weight;
215 Bl.noalias() += (rho_wet * lambda_wet * rho_wet +
216 rho_gas * rho_gas * lambda_gas) *
217 sm.dNdx.transpose() * permeability * b *
218 _ip_data[ip].integration_weight;
219
220 } // end of has gravity
221 }
222 if (_process_data._has_mass_lumping)
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 } // end of mass-lumping
242}
243
244} // namespace TwoPhaseFlowWithPrho
245} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
Definition of the PiecewiseLinearInterpolation class.
std::optional< std::size_t > getElementID() const
void setElementID(std::size_t element_id)
void assemble(double const t, double const, std::vector< double > const &local_x, std::vector< double > const &local_x_prev, std::vector< double > &local_M_data, std::vector< double > &local_K_data, std::vector< double > &local_b_data) override
void shapeFunctionInterpolate(const NodalValues &, const ShapeMatrix &)