www.mooseframework.org
RichardsFullyUpwindFlux.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 
9 
10 // MOOSE includes
11 #include "Assembly.h"
12 #include "MooseVariable.h"
13 #include "SystemBase.h"
14 
15 #include "libmesh/quadrature.h"
16 
17 template <>
18 InputParameters
20 {
21  InputParameters params = validParams<Kernel>();
22  params.addRequiredParam<std::vector<UserObjectName>>(
23  "relperm_UO", "List of names of user objects that define relative permeability");
24  params.addRequiredParam<std::vector<UserObjectName>>(
25  "seff_UO",
26  "List of name of user objects that define effective saturation as a function of "
27  "pressure list");
28  params.addRequiredParam<std::vector<UserObjectName>>(
29  "density_UO", "List of names of user objects that define the fluid density");
30  params.addRequiredParam<UserObjectName>(
31  "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
32  return params;
33 }
34 
35 RichardsFullyUpwindFlux::RichardsFullyUpwindFlux(const InputParameters & parameters)
36  : Kernel(parameters),
37  _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
38  _num_p(_richards_name_UO.num_v()),
39  _pvar(_richards_name_UO.richards_var_num(_var.number())),
40  _density_UO(getUserObjectByName<RichardsDensity>(
41  getParam<std::vector<UserObjectName>>("density_UO")[_pvar])),
42  _seff_UO(
43  getUserObjectByName<RichardsSeff>(getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])),
44  _relperm_UO(getUserObjectByName<RichardsRelPerm>(
45  getParam<std::vector<UserObjectName>>("relperm_UO")[_pvar])),
46  _viscosity(getMaterialProperty<std::vector<Real>>("viscosity")),
47  _flux_no_mob(getMaterialProperty<std::vector<RealVectorValue>>("flux_no_mob")),
48  _dflux_no_mob_dv(
49  getMaterialProperty<std::vector<std::vector<RealVectorValue>>>("dflux_no_mob_dv")),
50  _dflux_no_mob_dgradv(
51  getMaterialProperty<std::vector<std::vector<RealTensorValue>>>("dflux_no_mob_dgradv")),
52  _num_nodes(0),
53  _mobility(0),
54  _dmobility_dv(0)
55 {
56  _ps_at_nodes.resize(_num_p);
57  for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
59 }
60 
61 void
63 {
64  _num_nodes = (*_ps_at_nodes[_pvar]).size();
65 
66  Real p;
67  Real density;
68  Real ddensity_dp;
69  Real seff;
70  std::vector<Real> dseff_dp;
71  Real relperm;
72  Real drelperm_ds;
73  _mobility.resize(_num_nodes);
74  _dmobility_dv.resize(_num_nodes);
75  dseff_dp.resize(_num_p);
76  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
77  {
78  // retrieve and calculate basic things at the node
79  p = (*_ps_at_nodes[_pvar])[nodenum]; // pressure of fluid _pvar at node nodenum
80  density = _density_UO.density(p); // density of fluid _pvar at node nodenum
81  ddensity_dp = _density_UO.ddensity(p); // d(density)/dP
82  seff =
83  _seff_UO.seff(_ps_at_nodes, nodenum); // effective saturation of fluid _pvar at node nodenum
84  _seff_UO.dseff(_ps_at_nodes, nodenum, dseff_dp); // d(seff)/d(P_ph), for ph = 0, ..., _num_p - 1
85  relperm = _relperm_UO.relperm(seff); // relative permeability of fluid _pvar at node nodenum
86  drelperm_ds = _relperm_UO.drelperm(seff); // d(relperm)/dseff
87 
88  // calculate the mobility and its derivatives wrt (variable_ph = porepressure_ph)
89  _mobility[nodenum] =
90  density * relperm / _viscosity[0][_pvar]; // assume viscosity is constant throughout element
91  _dmobility_dv[nodenum].resize(_num_p);
92  for (unsigned int ph = 0; ph < _num_p; ++ph)
93  _dmobility_dv[nodenum][ph] = density * drelperm_ds * dseff_dp[ph] / _viscosity[0][_pvar];
94  _dmobility_dv[nodenum][_pvar] += ddensity_dp * relperm / _viscosity[0][_pvar];
95  }
96 }
97 
98 Real
100 {
101  // note this is not the complete residual:
102  // the upwind mobility parts get added in computeResidual
103  return _grad_test[_i][_qp] * _flux_no_mob[_qp][_pvar];
104 }
105 
106 void
108 {
109  upwind(true, false, 0);
110  return;
111 }
112 
113 void
115 {
116  upwind(false, true, jvar);
117  return;
118 }
119 
120 Real
122 {
123  // this is just the derivative of the flux WITHOUT the upstream mobility terms
124  // Those terms get added in during computeJacobian()
125  return _grad_test[_i][_qp] * (_dflux_no_mob_dgradv[_qp][_pvar][dvar] * _grad_phi[_j][_qp] +
126  _dflux_no_mob_dv[_qp][_pvar][dvar] * _phi[_j][_qp]);
127 }
128 
129 void
130 RichardsFullyUpwindFlux::upwind(bool compute_res, bool compute_jac, unsigned int jvar)
131 {
132  if (compute_jac && _richards_name_UO.not_richards_var(jvar))
133  return;
134 
135  // calculate the mobility values and their derivatives
137 
138  // compute the residual without the mobility terms
139  // Even if we are computing the jacobian we still need this
140  // in order to see which nodes are upwind and which are downwind
141  DenseVector<Number> & re = _assembly.residualBlock(_var.number());
142  _local_re.resize(re.size());
143  _local_re.zero();
144 
145  for (_i = 0; _i < _test.size(); _i++)
146  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
147  _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
148 
149  const unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
150  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
151  if (compute_jac)
152  {
153  _local_ke.resize(ke.m(), ke.n());
154  _local_ke.zero();
155 
156  for (_i = 0; _i < _test.size(); _i++)
157  for (_j = 0; _j < _phi.size(); _j++)
158  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
159  _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJac(dvar);
160  }
161 
162  // Now perform the upwinding by multiplying the residuals at the
163  // upstream nodes by their mobilities
164  //
165  // The residual for the kernel is the darcy flux.
166  // This is
167  // R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
168  // for node i. where int is the integral over the element.
169  // However, in fully-upwind, the first step is to take the mobility outside the
170  // integral, which was done in the _local_re calculation above.
171  //
172  // NOTE: Physically _local_re(_i) is a measure of fluid flowing out of node i
173  // If we had left in mobility, it would be exactly the mass flux flowing out of node i.
174  //
175  // This leads to the definition of upwinding:
176  // ***
177  // If _local_re(i) is positive then we use mobility_i. That is
178  // we use the upwind value of mobility.
179  // ***
180  //
181  // The final subtle thing is we must also conserve fluid mass: the total mass
182  // flowing out of node i must be the sum of the masses flowing
183  // into the other nodes.
184 
185  // FIRST:
186  // this is a dirty way of getting around precision loss problems
187  // and problems at steadystate where upwinding oscillates from
188  // node-to-node causing nonconvergence.
189  // I'm not sure if i actually need to do this in moose. Certainly
190  // in cosflow it is necessary.
191  // I will code a better algorithm if necessary
192  bool reached_steady = true;
193  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
194  {
195  if (_local_re(nodenum) >= 1E-20)
196  {
197  reached_steady = false;
198  break;
199  }
200  }
201 
202  // DEFINE VARIABLES USED TO ENSURE MASS CONSERVATION
203  // total mass out - used for mass conservation
204  Real total_mass_out = 0;
205  // total flux in
206  Real total_in = 0;
207 
208  // the following holds derivatives of these
209  std::vector<Real> dtotal_mass_out;
210  std::vector<Real> dtotal_in;
211  if (compute_jac)
212  {
213  dtotal_mass_out.resize(_num_nodes);
214  dtotal_in.resize(_num_nodes);
215  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
216  {
217  dtotal_mass_out[nodenum] = 0;
218  dtotal_in[nodenum] = 0;
219  }
220  }
221 
222  // PERFORM THE UPWINDING!
223  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
224  {
225  if (_local_re(nodenum) >= 0 || reached_steady) // upstream node
226  {
227  if (compute_jac)
228  {
229  for (_j = 0; _j < _phi.size(); _j++)
230  _local_ke(nodenum, _j) *= _mobility[nodenum];
231  _local_ke(nodenum, nodenum) += _dmobility_dv[nodenum][dvar] * _local_re(nodenum);
232  for (_j = 0; _j < _phi.size(); _j++)
233  dtotal_mass_out[_j] += _local_ke(nodenum, _j);
234  }
235  _local_re(nodenum) *= _mobility[nodenum];
236  total_mass_out += _local_re(nodenum);
237  }
238  else
239  {
240  total_in -= _local_re(nodenum); // note the -= means the result is positive
241  if (compute_jac)
242  for (_j = 0; _j < _phi.size(); _j++)
243  dtotal_in[_j] -= _local_ke(nodenum, _j);
244  }
245  }
246 
247  // CONSERVE MASS
248  // proportion the total_mass_out mass to the inflow nodes, weighting by their _local_re values
249  if (!reached_steady)
250  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
251  if (_local_re(nodenum) < 0)
252  {
253  if (compute_jac)
254  for (_j = 0; _j < _phi.size(); _j++)
255  {
256  _local_ke(nodenum, _j) *= total_mass_out / total_in;
257  _local_ke(nodenum, _j) +=
258  _local_re(nodenum) * (dtotal_mass_out[_j] / total_in -
259  dtotal_in[_j] * total_mass_out / total_in / total_in);
260  }
261  _local_re(nodenum) *= total_mass_out / total_in;
262  }
263 
264  // ADD RESULTS TO RESIDUAL OR JACOBIAN
265  if (compute_res)
266  {
267  re += _local_re;
268 
269  if (_has_save_in)
270  {
271  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
272  for (unsigned int i = 0; i < _save_in.size(); i++)
273  _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
274  }
275  }
276 
277  if (compute_jac)
278  {
279  ke += _local_ke;
280 
281  if (_has_diag_save_in && dvar == _pvar)
282  {
283  const unsigned int rows = ke.m();
284  DenseVector<Number> diag(rows);
285  for (unsigned int i = 0; i < rows; i++)
286  diag(i) = _local_ke(i, i);
287 
288  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
289  for (unsigned int i = 0; i < _diag_save_in.size(); i++)
290  _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
291  }
292  }
293 }
InputParameters validParams< RichardsFullyUpwindFlux >()
unsigned int _num_p
number of richards variables
virtual Real drelperm(Real seff) const =0
derivative of relative permeability wrt effective saturation This must be over-ridden in your derived...
const MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dflux_no_mob_dv
d(_flux_no_mob)/d(variable)
Base class for effective saturation as a function of porepressure(s) The functions seff...
Definition: RichardsSeff.h:22
Base class for Richards relative permeability classes that provide relative permeability as a functio...
std::vector< const VariableValue * > _ps_at_nodes
Holds the values of pressures at all the nodes of the element Eg: _ps_at_nodes[_pvar] is a pointer to...
std::vector< std::vector< Real > > _dmobility_dv
d(_mobility)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jacobian...
const std::string density
Definition: NS.h:15
const VariableValue * nodal_var(unsigned int richards_var_num) const
The nodal variable values for the given richards_var_num To extract a the value of pressure variable ...
bool not_richards_var(unsigned int moose_var_num) const
returns true if moose_var_num is not a richards var
virtual void dseff(std::vector< const VariableValue * > p, unsigned int qp, std::vector< Real > &result) const =0
derivative(s) of effective saturation as a function of porepressure(s) at given quadpoint of the elem...
This holds maps between pressure_var or pressure_var, sat_var used in RichardsMaterial and kernels...
virtual void computeResidual()
This simply calls upwind.
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 void computeOffDiagJacobian(unsigned int jvar)
this simply calls upwind
const RichardsRelPerm & _relperm_UO
user object defining the relative permeability
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 MaterialProperty< std::vector< RealVectorValue > > & _flux_no_mob
permeability*(grad(pressure) - density*gravity) (a vector of these in the multiphase case) ...
const RichardsSeff & _seff_UO
user object defining the effective saturation
RichardsFullyUpwindFlux(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...
unsigned int _pvar
the index of this variable in the list of Richards variables held by _richards_name_UO.
void prepareNodalValues()
calculates the nodal values of mobility, and derivatives thereof
virtual Real relperm(Real seff) const =0
relative permeability as a function of effective saturation This must be over-ridden in your derived ...
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dflux_no_mob_dgradv
d(_flux_no_mob)/d(grad(variable))
Real computeQpJac(unsigned int dvar)
the derivative of the flux without the upstream mobility terms
unsigned int _num_nodes
number of nodes in this element
const RichardsVarNames & _richards_name_UO
holds info regarding the names of the Richards variables and methods for extracting values of these v...
const MaterialProperty< std::vector< Real > > & _viscosity
viscosities
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
const RichardsDensity & _density_UO
user object defining the density
virtual Real computeQpResidual()
Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quad...
virtual Real seff(std::vector< const VariableValue * > p, unsigned int qp) const =0
effective saturation as a function of porepressure(s) at given quadpoint of the element ...
unsigned int richards_var_num(unsigned int moose_var_num) const
the richards variable number
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the...