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