OGS 6.1.0-1699-ge946d4c5f
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 
79  SpatialPosition pos;
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  permeability = perm;
91  else if (perm.rows() == 1)
92  permeability.diagonal().setConstant(perm(0, 0));
93 
94  for (unsigned ip = 0; ip < n_integration_points; ip++)
95  {
96  auto const& sm = _shape_matrices[ip];
97 
98  double pl_int_pt = 0.;
99  double totalrho_int_pt =
100  0.; // total mass density of the light component
101  NumLib::shapeFunctionInterpolate(local_x, sm.N, pl_int_pt,
102  totalrho_int_pt);
103 
104  const double temperature = _process_data._temperature(t, pos)[0];
105 
106  double const rho_gas =
107  _process_data._material->getGasDensity(pl_int_pt, temperature);
108  double const rho_h2o =
109  _process_data._material->getLiquidDensity(pl_int_pt, temperature);
110 
111  double& Sw = _ip_data[ip].sw;
113  double const X_h2_nonwet = 1.0;
114  double& rho_h2_wet = _ip_data[ip].rho_m;
115  double& dSwdP = _ip_data[ip].dsw_dpg;
116  double& dSwdrho = _ip_data[ip].dsw_drho;
117  double& drhoh2wet = _ip_data[ip].drhom_dpg;
118  double& drhoh2wet_drho = _ip_data[ip].drhom_drho;
119  if (!_ip_data[ip].mat_property.computeConstitutiveRelation(
120  t,
121  pos,
122  material_id,
123  pl_int_pt,
124  totalrho_int_pt,
125  temperature,
126  Sw,
127  rho_h2_wet,
128  dSwdP,
129  dSwdrho,
130  drhoh2wet,
131  drhoh2wet_drho))
132  OGS_FATAL("Computation of local constitutive relation failed.");
133  double const pc = _process_data._material->getCapillaryPressure(
134  material_id, t, pos, pl_int_pt, temperature, Sw);
135 
136  double const rho_wet = rho_h2o + rho_h2_wet;
137  _saturation[ip] = Sw;
138  _pressure_nonwetting[ip] = pl_int_pt + pc;
139 
140  // Assemble M matrix
141  // nonwetting
142  double dPC_dSw =
143  _process_data._material->getCapillaryPressureDerivative(
144  material_id, t, pos, pl_int_pt, temperature, Sw);
145 
146  double const porosity = _process_data._material->getPorosity(
147  material_id, t, pos, pl_int_pt, temperature, 0);
148 
149  Mgx.noalias() += porosity * _ip_data[ip].massOperator;
150 
151  Mlp.noalias() += porosity * rho_h2o * dSwdP * _ip_data[ip].massOperator;
152 
153  Mlx.noalias() +=
154  porosity * (1 + dSwdrho * rho_h2o) * _ip_data[ip].massOperator;
155  double const k_rel_gas =
156  _process_data._material->getNonwetRelativePermeability(
157  t, pos, _pressure_nonwetting[ip], temperature, Sw);
158  double const mu_gas = _process_data._material->getGasViscosity(
159  _pressure_nonwetting[ip], temperature);
160  double const lambda_gas = k_rel_gas / mu_gas;
161  double const diffusion_coeff_component_h2 =
162  _process_data._diffusion_coeff_component_b(t, pos)[0];
163 
164  // wet
165  double const k_rel_wet =
166  _process_data._material->getWetRelativePermeability(
167  t, pos, pl_int_pt, temperature, Sw);
168  double const mu_wet =
169  _process_data._material->getLiquidViscosity(pl_int_pt, temperature);
170  double const lambda_wet = k_rel_wet / mu_wet;
171 
172  laplace_operator.noalias() = sm.dNdx.transpose() * permeability *
173  sm.dNdx * _ip_data[ip].integration_weight;
174 
175  Kgp.noalias() +=
176  (rho_gas * X_h2_nonwet * lambda_gas * (1 + dPC_dSw * dSwdP) +
177  rho_h2_wet * lambda_wet) *
178  laplace_operator +
179  (Sw * porosity * diffusion_coeff_component_h2 *
180  (rho_h2o / rho_wet) * drhoh2wet) *
181  _ip_data[ip].diffusionOperator;
182  Kgx.noalias() +=
183  (rho_gas * X_h2_nonwet * lambda_gas * dPC_dSw * dSwdrho) *
184  laplace_operator +
185  (Sw * porosity * diffusion_coeff_component_h2 *
186  (rho_h2o / rho_wet) * drhoh2wet_drho) *
187  _ip_data[ip].diffusionOperator;
188  Klp.noalias() += (rho_gas * lambda_gas * (1 + dPC_dSw * dSwdP) +
189  rho_wet * lambda_wet) *
190  laplace_operator;
191 
192  Klx.noalias() +=
193  (rho_gas * lambda_gas * dPC_dSw * dSwdrho) * laplace_operator;
194 
195  if (_process_data._has_gravity)
196  {
197  auto const& b = _process_data._specific_body_force;
198  Bg.noalias() += (rho_gas * rho_gas * lambda_gas +
199  rho_h2_wet * rho_wet * lambda_wet) *
200  sm.dNdx.transpose() * permeability * b *
201  _ip_data[ip].integration_weight;
202  Bl.noalias() += (rho_wet * lambda_wet * rho_wet +
203  rho_gas * rho_gas * lambda_gas) *
204  sm.dNdx.transpose() * permeability * b *
205  _ip_data[ip].integration_weight;
206 
207  } // end of has gravity
208  }
209  if (_process_data._has_mass_lumping)
210  {
211  for (unsigned row = 0; row < Mgp.cols(); row++)
212  {
213  for (unsigned column = 0; column < Mgp.cols(); column++)
214  {
215  if (row != column)
216  {
217  Mgx(row, row) += Mgx(row, column);
218  Mgx(row, column) = 0.0;
219  Mgp(row, row) += Mgp(row, column);
220  Mgp(row, column) = 0.0;
221  Mlx(row, row) += Mlx(row, column);
222  Mlx(row, column) = 0.0;
223  Mlp(row, row) += Mlp(row, column);
224  Mlp(row, column) = 0.0;
225  }
226  }
227  }
228  } // end of mass-lumping
229 }
230 
231 } // end of namespace
232 } // end of namespace
Definition of the PiecewiseLinearInterpolation class.
void setElementID(std::size_t element_id)
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:71