www.mooseframework.org
Q2PPorepressureFlux.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "Q2PPorepressureFlux.h"
11 
12 // MOOSE includes
13 #include "Assembly.h"
14 #include "MooseVariable.h"
15 #include "SystemBase.h"
16 
17 #include "libmesh/quadrature.h"
18 
19 // C++ includes
20 #include <iostream>
21 
23 
26 {
28  params.addRequiredParam<UserObjectName>(
29  "fluid_density",
30  "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
31  params.addRequiredCoupledVar("saturation_variable", "The variable representing fluid saturation");
32  params.addRequiredParam<UserObjectName>(
33  "fluid_relperm",
34  "A RichardsRelPerm UserObject (eg RichardsRelPermPower) that defines the "
35  "fluid relative permeability as a function of the saturation variable");
36  params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
37  params.addClassDescription("Flux according to Darcy-Richards flow. The Variable for this Kernel "
38  "should be the porepressure.");
39  return params;
40 }
41 
43  : Kernel(parameters),
44  _density(getUserObject<RichardsDensity>("fluid_density")),
45  _sat(coupledDofValues("saturation_variable")),
46  _sat_var(coupled("saturation_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 = _sat.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(_var.dofValues()[nodenum]); // fluid density at the node
75  ddensity_dp = _density.ddensity(_var.dofValues()[nodenum]); // d(fluid density)/dP at the node
76  relperm = _relperm.relperm(_sat[nodenum]); // relative permeability of the fluid at node nodenum
77  drelperm_ds = _relperm.drelperm(_sat[nodenum]); // d(relperm)/dsat
78 
79  // calculate the mobility and its derivatives wrt P and S
80  _mobility[nodenum] = density * relperm / _viscosity;
81  _dmobility_dp[nodenum] = ddensity_dp * relperm / _viscosity;
82  _dmobility_ds[nodenum] = density * drelperm_ds / _viscosity;
83  }
84 }
85 
86 Real
88 {
89  // note this is not the complete residual:
90  // the upwind mobility parts get added in computeResidual
91  return _grad_test[_i][_qp] *
93 }
94 
95 void
97 {
98  upwind(true, false, 0);
99 }
100 
101 void
103 {
104  upwind(false, true, _var.number());
105 }
106 
107 void
109 {
110  upwind(false, true, jvar);
111 }
112 
113 Real
115 {
116  // this is just the derivative of the flux WITHOUT the upstream mobility terms
117  // Those terms get added in during computeJacobian()
118  if (dvar == _var.number())
119  return _grad_test[_i][_qp] *
120  (_permeability[_qp] *
122  else
123  return 0;
124 }
125 
126 void
127 Q2PPorepressureFlux::upwind(bool compute_res, bool compute_jac, unsigned int jvar)
128 {
129  if (compute_jac && !(jvar == _var.number() || jvar == _sat_var))
130  return;
131 
132  // calculate the mobility values and their derivatives
134 
135  // compute the residual without the mobility terms
136  // Even if we are computing the jacobian we still need this
137  // in order to see which nodes are upwind and which are downwind
139 
140  for (_i = 0; _i < _test.size(); _i++)
141  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
143 
144  if (compute_jac)
145  {
147  for (_i = 0; _i < _test.size(); _i++)
148  for (_j = 0; _j < _phi.size(); _j++)
149  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
150  _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJac(jvar);
151  }
152 
153  // Now perform the upwinding by multiplying the residuals at the
154  // upstream nodes by their mobilities
155  //
156  // The residual for the kernel is the darcy flux.
157  // This is
158  // R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
159  // for node i. where int is the integral over the element.
160  // However, in fully-upwind, the first step is to take the mobility outside the
161  // integral, which was done in the _local_re calculation above.
162  //
163  // NOTE: Physically _local_re(_i) is a measure of fluid flowing out of node i
164  // If we had left in mobility, it would be exactly the mass flux flowing out of node i.
165  //
166  // This leads to the definition of upwinding:
167  // ***
168  // If _local_re(i) is positive then we use mobility_i. That is
169  // we use the upwind value of mobility.
170  // ***
171  //
172  // The final subtle thing is we must also conserve fluid mass: the total mass
173  // flowing out of node i must be the sum of the masses flowing
174  // into the other nodes.
175 
176  // FIRST:
177  // this is a dirty way of getting around precision loss problems
178  // and problems at steadystate where upwinding oscillates from
179  // node-to-node causing nonconvergence.
180  // I'm not sure if i actually need to do this in moose. Certainly
181  // in cosflow it is necessary.
182  // I will code a better algorithm if necessary
183  bool reached_steady = true;
184  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
185  {
186  if (_local_re(nodenum) >= 1E-20)
187  {
188  reached_steady = false;
189  break;
190  }
191  }
192  reached_steady = false;
193 
194  // DEFINE VARIABLES USED TO ENSURE MASS CONSERVATION
195  // total mass out - used for mass conservation
196  Real total_mass_out = 0;
197  // total flux in
198  Real total_in = 0;
199 
200  // the following holds derivatives of these
201  std::vector<Real> dtotal_mass_out;
202  std::vector<Real> dtotal_in;
203  if (compute_jac)
204  {
205  dtotal_mass_out.assign(_num_nodes, 0);
206  dtotal_in.assign(_num_nodes, 0);
207  }
208 
209  // PERFORM THE UPWINDING!
210  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
211  {
212  if (_local_re(nodenum) >= 0 || reached_steady) // upstream node
213  {
214  if (compute_jac)
215  {
216  for (_j = 0; _j < _phi.size(); _j++)
217  _local_ke(nodenum, _j) *= _mobility[nodenum];
218  if (jvar == _var.number())
219  // deriv wrt P
220  _local_ke(nodenum, nodenum) += _dmobility_dp[nodenum] * _local_re(nodenum);
221  else
222  // deriv wrt S
223  _local_ke(nodenum, nodenum) += _dmobility_ds[nodenum] * _local_re(nodenum);
224  for (_j = 0; _j < _phi.size(); _j++)
225  dtotal_mass_out[_j] += _local_ke(nodenum, _j);
226  }
227  _local_re(nodenum) *= _mobility[nodenum];
228  total_mass_out += _local_re(nodenum);
229  }
230  else
231  {
232  total_in -= _local_re(nodenum); // note the -= means the result is positive
233  if (compute_jac)
234  for (_j = 0; _j < _phi.size(); _j++)
235  dtotal_in[_j] -= _local_ke(nodenum, _j);
236  }
237  }
238 
239  // CONSERVE MASS
240  // proportion the total_mass_out mass to the inflow nodes, weighting by their _local_re values
241  if (!reached_steady)
242  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
243  if (_local_re(nodenum) < 0)
244  {
245  if (compute_jac)
246  for (_j = 0; _j < _phi.size(); _j++)
247  {
248  _local_ke(nodenum, _j) *= total_mass_out / total_in;
249  _local_ke(nodenum, _j) +=
250  _local_re(nodenum) * (dtotal_mass_out[_j] / total_in -
251  dtotal_in[_j] * total_mass_out / total_in / total_in);
252  }
253  _local_re(nodenum) *= total_mass_out / total_in;
254  }
255 
256  // ADD RESULTS TO RESIDUAL OR JACOBIAN
257  if (compute_res)
258  {
260  if (_has_save_in)
261  for (unsigned int i = 0; i < _save_in.size(); i++)
262  _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
263  }
264 
265  if (compute_jac)
266  {
268  if (_has_diag_save_in && jvar == _var.number())
269  {
270  const unsigned int rows = _local_ke.m();
271  DenseVector<Number> diag(rows);
272  for (unsigned int i = 0; i < rows; i++)
273  diag(i) = _local_ke(i, i);
274 
275  for (unsigned int i = 0; i < _diag_save_in.size(); i++)
276  _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
277  }
278  }
279 }
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...
const VariableGradient & _grad_u
static InputParameters validParams()
void accumulateTaggedLocalResidual()
virtual Real drelperm(Real seff) const =0
derivative of relative permeability wrt effective saturation This must be over-ridden in your derived...
MooseVariable & _var
std::vector< MooseVariableFEBase *> _save_in
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 MooseArray< Real > & _JxW
const RichardsDensity & _density
fluid density
unsigned int number() const
Real _viscosity
fluid viscosity
Base class for Richards relative permeability classes that provide relative permeability as a functio...
const VariablePhiGradient & _grad_phi
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
static const std::string density
Definition: NS.h:33
This is a fully upwinded flux Kernel The Variable of this Kernel should be the porepressure.
unsigned int m() const
Q2PPorepressureFlux(const InputParameters &parameters)
const MaterialProperty< RealTensorValue > & _permeability
permeability
virtual void computeOffDiagJacobian(unsigned int jvar) override
this simply calls upwind
const RichardsRelPerm & _relperm
fluid relative permeability
bool _has_diag_save_in
std::vector< Real > _dmobility_ds
d(_mobility)/d(saturation) These are used in the jacobian calculations
DenseMatrix< Number > _local_ke
void addRequiredParam(const std::string &name, const std::string &doc_string)
static InputParameters validParams()
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...
TensorValue< Real > RealTensorValue
const VariableTestValue & _test
std::vector< MooseVariableFEBase *> _diag_save_in
const QBase *const & _qrule
bool _has_save_in
unsigned int _i
virtual Real relperm(Real seff) const =0
relative permeability as a function of effective saturation This must be over-ridden in your derived ...
void accumulateTaggedLocalMatrix()
virtual void computeResidual() override
This simply calls upwind.
const MooseArray< Real > & _coord
Assembly & _assembly
unsigned int _num_nodes
number of nodes in the element
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
virtual void computeJacobian() override
this simply calls upwind
unsigned int _j
virtual Real computeQpResidual() override
Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quad...
void prepareNodalValues()
calculates the nodal values of mobility, and derivatives thereof
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const VariableTestGradient & _grad_test
const DoFValue & dofValues() const override
DenseVector< Number > _local_re
registerMooseObject("RichardsApp", Q2PPorepressureFlux)
std::vector< Real > _dmobility_dp
d(_mobility)/d(porepressure) These are used in the jacobian calculations
void addClassDescription(const std::string &doc_string)
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
void prepareVectorTag(Assembly &assembly, unsigned int ivar)
void prepareMatrixTag(Assembly &assembly, unsigned int ivar, unsigned int jvar)
const VariablePhiValue & _phi
const VariableValue & _u
unsigned int _qp
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