OGS
TwoPhaseFlowWithPPLocalAssembler-impl.h
Go to the documentation of this file.
1
31#pragma once
32
38
39namespace ProcessLib
40{
41namespace TwoPhaseFlowWithPP
42{
43namespace MPL = MaterialPropertyLib;
44
45template <typename ShapeFunction, int GlobalDim>
47 setInitialConditionsConcrete(Eigen::VectorXd const local_x,
48 double const t,
49 int const /*process_id*/)
50{
51 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
52 unsigned const n_integration_points =
53 _integration_method.getNumberOfPoints();
54
55 auto const pc =
56 local_x.segment(cap_pressure_matrix_index, cap_pressure_size);
57 constexpr double dt = std::numeric_limits<double>::quiet_NaN();
59 pos.setElementID(_element.getID());
60 for (unsigned ip = 0; ip < n_integration_points; ip++)
61 {
62 auto const& N = _ip_data[ip].N;
63 MPL::VariableArray variables;
64 variables.capillary_pressure = N.dot(pc);
65 _saturation[ip] = medium.property(MPL::PropertyType::saturation)
66 .template value<double>(variables, pos, t, dt);
67 }
68}
69
70template <typename ShapeFunction, int GlobalDim>
72 double const t, double const dt, std::vector<double> const& local_x,
73 std::vector<double> const& /*local_x_prev*/,
74 std::vector<double>& local_M_data, std::vector<double>& local_K_data,
75 std::vector<double>& local_b_data)
76{
77 auto const local_matrix_size = local_x.size();
78
79 assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
80
82 local_M_data, local_matrix_size, local_matrix_size);
84 local_K_data, local_matrix_size, local_matrix_size);
86 local_b_data, local_matrix_size);
87
88 auto Mgp =
89 local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
90 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
91 auto Mgpc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
92 nonwet_pressure_matrix_index, cap_pressure_matrix_index);
93
94 auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
95 cap_pressure_matrix_index, cap_pressure_matrix_index);
96
97 NodalMatrixType laplace_operator =
98 NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
99
100 auto Kgp =
101 local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
102 nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
103
104 auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
105 cap_pressure_matrix_index, nonwet_pressure_matrix_index);
106
107 auto Klpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
108 cap_pressure_matrix_index, cap_pressure_matrix_index);
109
110 auto Bg = local_b.template segment<nonwet_pressure_size>(
111 nonwet_pressure_matrix_index);
112
113 auto Bl =
114 local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
115
116 auto const& medium = *_process_data.media_map.getMedium(_element.getID());
117 auto const& liquid_phase = medium.phase("AqueousLiquid");
118 auto const& gas_phase = medium.phase("Gas");
119
120 unsigned const n_integration_points =
121 _integration_method.getNumberOfPoints();
122
124 pos.setElementID(_element.getID());
125
126 for (unsigned ip = 0; ip < n_integration_points; ip++)
127 {
128 auto const& w = _ip_data[ip].integration_weight;
129 auto const& N = _ip_data[ip].N;
130 auto const& dNdx = _ip_data[ip].dNdx;
131
132 NodalMatrixType M = N.transpose() * N * w;
133
134 double pc_int_pt = 0.;
135 double pn_int_pt = 0.;
136 NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, pn_int_pt,
137 pc_int_pt);
138
139 _pressure_wet[ip] = pn_int_pt - pc_int_pt;
140 MPL::VariableArray variables;
141
142 variables.gas_phase_pressure = pn_int_pt;
143 variables.capillary_pressure = pc_int_pt;
144 variables.temperature = _process_data.temperature(t, pos)[0];
145 variables.molar_mass =
146 gas_phase.property(MPL::PropertyType::molar_mass)
147 .template value<double>(variables, pos, t, dt);
148
149 auto const rho_nonwet =
150 gas_phase.property(MPL::PropertyType::density)
151 .template value<double>(variables, pos, t, dt);
152
153 auto const rho_wet = liquid_phase.property(MPL::PropertyType::density)
154 .template value<double>(variables, pos, t, dt);
155 auto& Sw = _saturation[ip];
156 Sw = medium.property(MPL::PropertyType::saturation)
157 .template value<double>(variables, pos, t, dt);
158
159 variables.liquid_saturation = Sw;
160 auto const dSw_dpc =
161 medium.property(MPL::PropertyType::saturation)
162 .template dValue<double>(
163 variables, MPL::Variable::capillary_pressure, pos, t, dt);
164
165 auto const porosity =
166 medium.property(MPL::PropertyType::porosity)
167 .template value<double>(variables, pos, t, dt);
168
169 auto const drhononwet_dpn =
170 gas_phase.property(MPL::PropertyType::density)
171 .template dValue<double>(
172 variables, MPL::Variable::gas_phase_pressure, pos, t, dt);
173
174 auto const k_rel_wet =
175 medium.property(MPL::PropertyType::relative_permeability)
176 .template value<double>(variables, pos, t, dt);
177 auto const k_rel_nonwet =
178 medium
179 .property(
180 MPL::PropertyType::relative_permeability_nonwetting_phase)
181 .template value<double>(variables, pos, t, dt);
182
183 auto const mu_nonwet =
184 gas_phase.property(MPL::PropertyType::viscosity)
185 .template value<double>(variables, pos, t, dt);
186
187 auto const lambda_nonwet = k_rel_nonwet / mu_nonwet;
188
189 auto const mu_wet = liquid_phase.property(MPL::PropertyType::viscosity)
190 .template value<double>(variables, pos, t, dt);
191
192 auto const lambda_wet = k_rel_wet / mu_wet;
193
194 auto const K = MPL::formEigenTensor<GlobalDim>(
195 medium.property(MPL::PropertyType::permeability)
196 .template value<double>(variables, pos, t, dt));
197
198 Mgp.noalias() += porosity * (1 - Sw) * drhononwet_dpn * M;
199 Mgpc.noalias() += -porosity * rho_nonwet * dSw_dpc * M;
200 Mlpc.noalias() += porosity * dSw_dpc * rho_wet * M;
201
202 laplace_operator.noalias() = dNdx.transpose() * K * dNdx * w;
203
204 Kgp.noalias() += rho_nonwet * lambda_nonwet * laplace_operator;
205 Klp.noalias() += rho_wet * lambda_wet * laplace_operator;
206 Klpc.noalias() += -rho_wet * lambda_wet * laplace_operator;
207
208 if (_process_data.has_gravity)
209 {
210 auto const& b = _process_data.specific_body_force;
211
212 NodalVectorType gravity_operator = dNdx.transpose() * K * b * w;
213 Bg.noalias() +=
214 rho_nonwet * rho_nonwet * lambda_nonwet * gravity_operator;
215 Bl.noalias() += rho_wet * rho_wet * lambda_wet * gravity_operator;
216 } // end of has gravity
217 }
218 if (_process_data.has_mass_lumping)
219 {
220 for (unsigned row = 0; row < Mgpc.cols(); row++)
221 {
222 for (unsigned column = 0; column < Mgpc.cols(); column++)
223 {
224 if (row != column)
225 {
226 Mgpc(row, row) += Mgpc(row, column);
227 Mgpc(row, column) = 0.0;
228 Mgp(row, row) += Mgp(row, column);
229 Mgp(row, column) = 0.0;
230 Mlpc(row, row) += Mlpc(row, column);
231 Mlpc(row, column) = 0.0;
232 }
233 }
234 }
235 } // end of mass-lumping
236}
237
238} // namespace TwoPhaseFlowWithPP
239} // namespace ProcessLib
Definition of the PiecewiseLinearInterpolation class.
void setElementID(std::size_t element_id)
void setInitialConditionsConcrete(Eigen::VectorXd const local_x, double const t, int const process_id) override
void assemble(double const t, double const dt, 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
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 &)