OGS 6.2.0-97-g4a610c866
TwoPhaseFlowWithPrhoLocalAssembler-impl.h
Go to the documentation of this file.
1 
10 #pragma once
11 
13 
17 
18 namespace ProcessLib
19 {
20 namespace TwoPhaseFlowWithPrho
21 {
22 template <typename ShapeFunction, typename IntegrationMethod,
23  unsigned GlobalDim>
24 void TwoPhaseFlowWithPrhoLocalAssembler<
25  ShapeFunction, IntegrationMethod,
26  GlobalDim>::assemble(double const t, std::vector<double> const& local_x,
27  std::vector<double>& local_M_data,
28  std::vector<double>& local_K_data,
29  std::vector<double>& local_b_data)
30 {
31  auto const local_matrix_size = local_x.size();
32 
33  assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
34 
35  auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
36  local_M_data, local_matrix_size, local_matrix_size);
37  auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
38  local_K_data, local_matrix_size, local_matrix_size);
39  auto local_b = MathLib::createZeroedVector<LocalVectorType>(
40  local_b_data, local_matrix_size);
41 
42  auto Mgp =
43  local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
44  nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
45  auto Mgx = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
46  nonwet_pressure_matrix_index, cap_pressure_matrix_index);
47 
48  auto Mlp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
49  cap_pressure_matrix_index, nonwet_pressure_matrix_index);
50 
51  auto Mlx = local_M.template block<cap_pressure_size, cap_pressure_size>(
52  cap_pressure_matrix_index, cap_pressure_matrix_index);
53 
54  NodalMatrixType laplace_operator =
55  NodalMatrixType::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
56 
57  auto Kgp =
58  local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
59  nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
60 
61  auto Kgx = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
62  nonwet_pressure_matrix_index, cap_pressure_matrix_index);
63 
64  auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
65  cap_pressure_matrix_index, nonwet_pressure_matrix_index);
66 
67  auto Klx = local_K.template block<cap_pressure_size, cap_pressure_size>(
68  cap_pressure_matrix_index, cap_pressure_matrix_index);
69 
70  auto Bg = local_b.template segment<nonwet_pressure_size>(
71  nonwet_pressure_matrix_index);
72 
73  auto Bl =
74  local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
75 
76  unsigned const n_integration_points =
77  _integration_method.getNumberOfPoints();
78 
80  pos.setElementID(_element.getID());
81  const int material_id =
82  _process_data._material->getMaterialID(pos.getElementID().get());
83 
84  const Eigen::MatrixXd& perm = _process_data._material->getPermeability(
85  material_id, t, pos, _element.getDimension());
86  assert(perm.rows() == _element.getDimension() || perm.rows() == 1);
87  GlobalDimMatrixType permeability = GlobalDimMatrixType::Zero(
88  _element.getDimension(), _element.getDimension());
89  if (perm.rows() == _element.getDimension())
90  {
91  permeability = perm;
92  }
93  else if (perm.rows() == 1)
94  {
95  permeability.diagonal().setConstant(perm(0, 0));
96  }
97 
98  for (unsigned ip = 0; ip < n_integration_points; ip++)
99  {
100  auto const& sm = _shape_matrices[ip];
101 
102  double pl_int_pt = 0.;
103  double totalrho_int_pt =
104  0.; // total mass density of the light component
105  NumLib::shapeFunctionInterpolate(local_x, sm.N, pl_int_pt,
106  totalrho_int_pt);
107 
108  const double temperature = _process_data._temperature(t, pos)[0];
109 
110  double const rho_gas =
111  _process_data._material->getGasDensity(pl_int_pt, temperature);
112  double const rho_h2o =
113  _process_data._material->getLiquidDensity(pl_int_pt, temperature);
114 
115  double& Sw = _ip_data[ip].sw;
117  double const X_h2_nonwet = 1.0;
118  double& rho_h2_wet = _ip_data[ip].rho_m;
119  double& dSwdP = _ip_data[ip].dsw_dpg;
120  double& dSwdrho = _ip_data[ip].dsw_drho;
121  double& drhoh2wet = _ip_data[ip].drhom_dpg;
122  double& drhoh2wet_drho = _ip_data[ip].drhom_drho;
123  if (!_ip_data[ip].mat_property.computeConstitutiveRelation(
124  t,
125  pos,
126  material_id,
127  pl_int_pt,
128  totalrho_int_pt,
129  temperature,
130  Sw,
131  rho_h2_wet,
132  dSwdP,
133  dSwdrho,
134  drhoh2wet,
135  drhoh2wet_drho))
136  {
137  OGS_FATAL("Computation of local constitutive relation failed.");
138  }
139  double const pc = _process_data._material->getCapillaryPressure(
140  material_id, t, pos, pl_int_pt, temperature, Sw);
141 
142  double const rho_wet = rho_h2o + rho_h2_wet;
143  _saturation[ip] = Sw;
144  _pressure_nonwetting[ip] = pl_int_pt + pc;
145 
146  // Assemble M matrix
147  // nonwetting
148  double dPC_dSw =
149  _process_data._material->getCapillaryPressureDerivative(
150  material_id, t, pos, pl_int_pt, temperature, Sw);
151 
152  double const porosity = _process_data._material->getPorosity(
153  material_id, t, pos, pl_int_pt, temperature, 0);
154 
155  Mgx.noalias() += porosity * _ip_data[ip].massOperator;
156 
157  Mlp.noalias() += porosity * rho_h2o * dSwdP * _ip_data[ip].massOperator;
158 
159  Mlx.noalias() +=
160  porosity * (1 + dSwdrho * rho_h2o) * _ip_data[ip].massOperator;
161  double const k_rel_gas =
162  _process_data._material->getNonwetRelativePermeability(
163  t, pos, _pressure_nonwetting[ip], temperature, Sw);
164  double const mu_gas = _process_data._material->getGasViscosity(
165  _pressure_nonwetting[ip], temperature);
166  double const lambda_gas = k_rel_gas / mu_gas;
167  double const diffusion_coeff_component_h2 =
168  _process_data._diffusion_coeff_component_b(t, pos)[0];
169 
170  // wet
171  double const k_rel_wet =
172  _process_data._material->getWetRelativePermeability(
173  t, pos, pl_int_pt, temperature, Sw);
174  double const mu_wet =
175  _process_data._material->getLiquidViscosity(pl_int_pt, temperature);
176  double const lambda_wet = k_rel_wet / mu_wet;
177 
178  laplace_operator.noalias() = sm.dNdx.transpose() * permeability *
179  sm.dNdx * _ip_data[ip].integration_weight;
180 
181  Kgp.noalias() +=
182  (rho_gas * X_h2_nonwet * lambda_gas * (1 + dPC_dSw * dSwdP) +
183  rho_h2_wet * lambda_wet) *
184  laplace_operator +
185  (Sw * porosity * diffusion_coeff_component_h2 *
186  (rho_h2o / rho_wet) * drhoh2wet) *
187  _ip_data[ip].diffusionOperator;
188  Kgx.noalias() +=
189  (rho_gas * X_h2_nonwet * lambda_gas * dPC_dSw * dSwdrho) *
190  laplace_operator +
191  (Sw * porosity * diffusion_coeff_component_h2 *
192  (rho_h2o / rho_wet) * drhoh2wet_drho) *
193  _ip_data[ip].diffusionOperator;
194  Klp.noalias() += (rho_gas * lambda_gas * (1 + dPC_dSw * dSwdP) +
195  rho_wet * lambda_wet) *
196  laplace_operator;
197 
198  Klx.noalias() +=
199  (rho_gas * lambda_gas * dPC_dSw * dSwdrho) * laplace_operator;
200 
201  if (_process_data._has_gravity)
202  {
203  auto const& b = _process_data._specific_body_force;
204  Bg.noalias() += (rho_gas * rho_gas * lambda_gas +
205  rho_h2_wet * rho_wet * lambda_wet) *
206  sm.dNdx.transpose() * permeability * b *
207  _ip_data[ip].integration_weight;
208  Bl.noalias() += (rho_wet * lambda_wet * rho_wet +
209  rho_gas * rho_gas * lambda_gas) *
210  sm.dNdx.transpose() * permeability * b *
211  _ip_data[ip].integration_weight;
212 
213  } // end of has gravity
214  }
215  if (_process_data._has_mass_lumping)
216  {
217  for (unsigned row = 0; row < Mgp.cols(); row++)
218  {
219  for (unsigned column = 0; column < Mgp.cols(); column++)
220  {
221  if (row != column)
222  {
223  Mgx(row, row) += Mgx(row, column);
224  Mgx(row, column) = 0.0;
225  Mgp(row, row) += Mgp(row, column);
226  Mgp(row, column) = 0.0;
227  Mlx(row, row) += Mlx(row, column);
228  Mlx(row, column) = 0.0;
229  Mlp(row, row) += Mlp(row, column);
230  Mlp(row, column) = 0.0;
231  }
232  }
233  }
234  } // end of mass-lumping
235 }
236 
237 } // namespace TwoPhaseFlowWithPrho
238 } // namespace ProcessLib
void setElementID(std::size_t element_id)
Definition of the PiecewiseLinearInterpolation class.
void shapeFunctionInterpolate(const NodalValues &nodal_values, const ShapeMatrix &shape_matrix_N, double &interpolated_value, ScalarTypes &... interpolated_values)
Definition: Interpolation.h:73
boost::optional< std::size_t > getElementID() const
#define OGS_FATAL(fmt,...)
Definition: Error.h:63