www.mooseframework.org
Q2PPorepressureFlux.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "Q2PPorepressureFlux.h"
9 
10 // MOOSE includes
11 #include "Assembly.h"
12 #include "MooseVariable.h"
13 #include "SystemBase.h"
14 
15 #include "libmesh/quadrature.h"
16 
17 // C++ includes
18 #include <iostream>
19 
20 template <>
21 InputParameters
23 {
24  InputParameters params = validParams<Kernel>();
25  params.addRequiredParam<UserObjectName>(
26  "fluid_density",
27  "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
28  params.addRequiredCoupledVar("saturation_variable", "The variable representing fluid saturation");
29  params.addRequiredParam<UserObjectName>(
30  "fluid_relperm",
31  "A RichardsRelPerm UserObject (eg RichardsRelPermPower) that defines the "
32  "fluid relative permeability as a function of the saturation variable");
33  params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
34  params.addClassDescription("Flux according to Darcy-Richards flow. The Variable for this Kernel "
35  "should be the porepressure.");
36  return params;
37 }
38 
39 Q2PPorepressureFlux::Q2PPorepressureFlux(const InputParameters & parameters)
40  : Kernel(parameters),
41  _density(getUserObject<RichardsDensity>("fluid_density")),
42  _sat(coupledNodalValue("saturation_variable")),
43  _sat_var(coupled("saturation_variable")),
44  _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
45  _viscosity(getParam<Real>("fluid_viscosity")),
46  _gravity(getMaterialProperty<RealVectorValue>("gravity")),
47  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
48  _num_nodes(0),
49  _mobility(0),
50  _dmobility_dp(0),
51  _dmobility_ds(0)
52 {
53 }
54 
55 void
57 {
58  _num_nodes = _sat.size();
59 
60  Real density;
61  Real ddensity_dp;
62  Real relperm;
63  Real drelperm_ds;
64 
65  _mobility.resize(_num_nodes);
66  _dmobility_dp.resize(_num_nodes);
67  _dmobility_ds.resize(_num_nodes);
68 
69  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
70  {
71  density = _density.density(_var.nodalSln()[nodenum]); // fluid density at the node
72  ddensity_dp = _density.ddensity(_var.nodalSln()[nodenum]); // d(fluid density)/dP at the node
73  relperm = _relperm.relperm(_sat[nodenum]); // relative permeability of the fluid at node nodenum
74  drelperm_ds = _relperm.drelperm(_sat[nodenum]); // d(relperm)/dsat
75 
76  // calculate the mobility and its derivatives wrt P and S
77  _mobility[nodenum] = density * relperm / _viscosity;
78  _dmobility_dp[nodenum] = ddensity_dp * relperm / _viscosity;
79  _dmobility_ds[nodenum] = density * drelperm_ds / _viscosity;
80  }
81 }
82 
83 Real
85 {
86  // note this is not the complete residual:
87  // the upwind mobility parts get added in computeResidual
88  return _grad_test[_i][_qp] *
89  (_permeability[_qp] * (_grad_u[_qp] - _density.density(_u[_qp]) * _gravity[_qp]));
90 }
91 
92 void
94 {
95  upwind(true, false, 0);
96 }
97 
98 void
100 {
101  upwind(false, true, _var.number());
102 }
103 
104 void
106 {
107  upwind(false, true, jvar);
108 }
109 
110 Real
112 {
113  // this is just the derivative of the flux WITHOUT the upstream mobility terms
114  // Those terms get added in during computeJacobian()
115  if (dvar == _var.number())
116  return _grad_test[_i][_qp] *
117  (_permeability[_qp] *
118  (_grad_phi[_j][_qp] - _density.ddensity(_u[_qp]) * _gravity[_qp] * _phi[_j][_qp]));
119  else
120  return 0;
121 }
122 
123 void
124 Q2PPorepressureFlux::upwind(bool compute_res, bool compute_jac, unsigned int jvar)
125 {
126  if (compute_jac && !(jvar == _var.number() || jvar == _sat_var))
127  return;
128 
129  // calculate the mobility values and their derivatives
131 
132  // compute the residual without the mobility terms
133  // Even if we are computing the jacobian we still need this
134  // in order to see which nodes are upwind and which are downwind
135  DenseVector<Number> & re = _assembly.residualBlock(_var.number());
136  _local_re.resize(re.size());
137  _local_re.zero();
138 
139  for (_i = 0; _i < _test.size(); _i++)
140  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
141  _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
142 
143  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
144  if (compute_jac)
145  {
146  _local_ke.resize(ke.m(), ke.n());
147  _local_ke.zero();
148 
149  for (_i = 0; _i < _test.size(); _i++)
150  for (_j = 0; _j < _phi.size(); _j++)
151  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
152  _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJac(jvar);
153  }
154 
155  // Now perform the upwinding by multiplying the residuals at the
156  // upstream nodes by their mobilities
157  //
158  // The residual for the kernel is the darcy flux.
159  // This is
160  // R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
161  // for node i. where int is the integral over the element.
162  // However, in fully-upwind, the first step is to take the mobility outside the
163  // integral, which was done in the _local_re calculation above.
164  //
165  // NOTE: Physically _local_re(_i) is a measure of fluid flowing out of node i
166  // If we had left in mobility, it would be exactly the mass flux flowing out of node i.
167  //
168  // This leads to the definition of upwinding:
169  // ***
170  // If _local_re(i) is positive then we use mobility_i. That is
171  // we use the upwind value of mobility.
172  // ***
173  //
174  // The final subtle thing is we must also conserve fluid mass: the total mass
175  // flowing out of node i must be the sum of the masses flowing
176  // into the other nodes.
177 
178  // FIRST:
179  // this is a dirty way of getting around precision loss problems
180  // and problems at steadystate where upwinding oscillates from
181  // node-to-node causing nonconvergence.
182  // I'm not sure if i actually need to do this in moose. Certainly
183  // in cosflow it is necessary.
184  // I will code a better algorithm if necessary
185  bool reached_steady = true;
186  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
187  {
188  if (_local_re(nodenum) >= 1E-20)
189  {
190  reached_steady = false;
191  break;
192  }
193  }
194  reached_steady = false;
195 
196  // DEFINE VARIABLES USED TO ENSURE MASS CONSERVATION
197  // total mass out - used for mass conservation
198  Real total_mass_out = 0;
199  // total flux in
200  Real total_in = 0;
201 
202  // the following holds derivatives of these
203  std::vector<Real> dtotal_mass_out;
204  std::vector<Real> dtotal_in;
205  if (compute_jac)
206  {
207  dtotal_mass_out.assign(_num_nodes, 0);
208  dtotal_in.assign(_num_nodes, 0);
209  }
210 
211  // PERFORM THE UPWINDING!
212  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
213  {
214  if (_local_re(nodenum) >= 0 || reached_steady) // upstream node
215  {
216  if (compute_jac)
217  {
218  for (_j = 0; _j < _phi.size(); _j++)
219  _local_ke(nodenum, _j) *= _mobility[nodenum];
220  if (jvar == _var.number())
221  // deriv wrt P
222  _local_ke(nodenum, nodenum) += _dmobility_dp[nodenum] * _local_re(nodenum);
223  else
224  // deriv wrt S
225  _local_ke(nodenum, nodenum) += _dmobility_ds[nodenum] * _local_re(nodenum);
226  for (_j = 0; _j < _phi.size(); _j++)
227  dtotal_mass_out[_j] += _local_ke(nodenum, _j);
228  }
229  _local_re(nodenum) *= _mobility[nodenum];
230  total_mass_out += _local_re(nodenum);
231  }
232  else
233  {
234  total_in -= _local_re(nodenum); // note the -= means the result is positive
235  if (compute_jac)
236  for (_j = 0; _j < _phi.size(); _j++)
237  dtotal_in[_j] -= _local_ke(nodenum, _j);
238  }
239  }
240 
241  // CONSERVE MASS
242  // proportion the total_mass_out mass to the inflow nodes, weighting by their _local_re values
243  if (!reached_steady)
244  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
245  if (_local_re(nodenum) < 0)
246  {
247  if (compute_jac)
248  for (_j = 0; _j < _phi.size(); _j++)
249  {
250  _local_ke(nodenum, _j) *= total_mass_out / total_in;
251  _local_ke(nodenum, _j) +=
252  _local_re(nodenum) * (dtotal_mass_out[_j] / total_in -
253  dtotal_in[_j] * total_mass_out / total_in / total_in);
254  }
255  _local_re(nodenum) *= total_mass_out / total_in;
256  }
257 
258  // ADD RESULTS TO RESIDUAL OR JACOBIAN
259  if (compute_res)
260  {
261  re += _local_re;
262 
263  if (_has_save_in)
264  {
265  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
266  for (unsigned int i = 0; i < _save_in.size(); i++)
267  _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
268  }
269  }
270 
271  if (compute_jac)
272  {
273  ke += _local_ke;
274 
275  if (_has_diag_save_in && jvar == _var.number())
276  {
277  const unsigned int rows = ke.m();
278  DenseVector<Number> diag(rows);
279  for (unsigned int i = 0; i < rows; i++)
280  diag(i) = _local_ke(i, i);
281 
282  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
283  for (unsigned int i = 0; i < _diag_save_in.size(); i++)
284  _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
285  }
286  }
287 }
void upwind(bool compute_res, bool compute_jac, unsigned int jvar)
Do the upwinding for both the residual and jacobian I&#39;ve put both calculations in the same code to tr...
virtual Real drelperm(Real seff) const =0
derivative of relative permeability wrt effective saturation This must be over-ridden in your derived...
const RichardsDensity & _density
fluid density
virtual void computeJacobian()
this simply calls upwind
Real _viscosity
fluid viscosity
Base class for Richards relative permeability classes that provide relative permeability as a functio...
Real computeQpJac(unsigned int dvar)
the derivative of the flux without the upstream mobility terms
const VariableValue & _sat
saturation at the nodes
unsigned int _sat_var
variable number of the saturation variable
const std::string density
Definition: NS.h:15
Q2PPorepressureFlux(const InputParameters &parameters)
const MaterialProperty< RealTensorValue > & _permeability
permeability
const RichardsRelPerm & _relperm
fluid relative permeability
std::vector< Real > _dmobility_ds
d(_mobility)/d(saturation) These are used in the jacobian calculations
virtual Real density(Real p) const =0
fluid density as a function of porepressure This must be over-ridden in derived classes to provide an...
virtual Real computeQpResidual()
Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quad...
InputParameters validParams< Q2PPorepressureFlux >()
unsigned int _num_nodes
number of nodes in the element
virtual void computeOffDiagJacobian(unsigned int jvar)
this simply calls upwind
virtual Real ddensity(Real p) const =0
derivative of fluid density wrt porepressure This must be over-ridden in derived classes to provide a...
void prepareNodalValues()
calculates the nodal values of mobility, and derivatives thereof
virtual void computeResidual()
This simply calls upwind.
virtual Real relperm(Real seff) const =0
relative permeability as a function of effective saturation This must be over-ridden in your derived ...
std::vector< Real > _dmobility_dp
d(_mobility)/d(porepressure) These are used in the jacobian calculations
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the...
const MaterialProperty< RealVectorValue > & _gravity
gravity