www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
Q2PPorepressureFlux Class Reference

This is a fully upwinded flux Kernel The Variable of this Kernel should be the porepressure. More...

#include <Q2PPorepressureFlux.h>

Inheritance diagram for Q2PPorepressureFlux:
[legend]

Public Member Functions

 Q2PPorepressureFlux (const InputParameters &parameters)
 

Protected Member Functions

virtual Real computeQpResidual ()
 Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quadpoints and then add the upwind mobility parts. More...
 
virtual void computeResidual ()
 This simply calls upwind. More...
 
virtual void computeOffDiagJacobian (unsigned int jvar)
 this simply calls upwind More...
 
virtual void computeJacobian ()
 this simply calls upwind More...
 
Real computeQpJac (unsigned int dvar)
 the derivative of the flux without the upstream mobility terms More...
 
void upwind (bool compute_res, bool compute_jac, unsigned int jvar)
 Do the upwinding for both the residual and jacobian I've put both calculations in the same code to try to reduce code duplication. More...
 
void prepareNodalValues ()
 calculates the nodal values of mobility, and derivatives thereof More...
 

Protected Attributes

const RichardsDensity_density
 fluid density More...
 
const VariableValue & _sat
 saturation at the nodes More...
 
unsigned int _sat_var
 variable number of the saturation variable More...
 
const RichardsRelPerm_relperm
 fluid relative permeability More...
 
Real _viscosity
 fluid viscosity More...
 
const MaterialProperty< RealVectorValue > & _gravity
 gravity More...
 
const MaterialProperty< RealTensorValue > & _permeability
 permeability More...
 
unsigned int _num_nodes
 number of nodes in the element More...
 
std::vector< Real > _mobility
 nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the residual More...
 
std::vector< Real > _dmobility_dp
 d(_mobility)/d(porepressure) These are used in the jacobian calculations More...
 
std::vector< Real > _dmobility_ds
 d(_mobility)/d(saturation) These are used in the jacobian calculations More...
 

Detailed Description

This is a fully upwinded flux Kernel The Variable of this Kernel should be the porepressure.

The residual for the kernel is the darcy flux. This is R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)} for node i. where int is the integral over the element, and pot = Porepressure - density*gravity.x

However, in fully-upwind, the first step is to take the mobility outside the integral. R_i = mobility*int{flux_no_mob} = mobility*F_i NOTE: R_i is exactly the mass flux flowing out of node i. Similarly, F_i is a measure of fluid flowing out of node i.

This leads to the definition of upwinding:

If F_i is positive then R_i = mobility_i * F_i That is, we use the upwind value of mobility.

For the F_i<0 nodes we construct their R_i using mass conservation

Definition at line 45 of file Q2PPorepressureFlux.h.

Constructor & Destructor Documentation

Q2PPorepressureFlux::Q2PPorepressureFlux ( const InputParameters &  parameters)

Definition at line 39 of file Q2PPorepressureFlux.C.

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 }
const RichardsDensity & _density
fluid density
Real _viscosity
fluid viscosity
const VariableValue & _sat
saturation at the nodes
unsigned int _sat_var
variable number of the saturation variable
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
unsigned int _num_nodes
number of nodes in the element
std::vector< Real > _dmobility_dp
d(_mobility)/d(porepressure) These are used in the jacobian calculations
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

Member Function Documentation

void Q2PPorepressureFlux::computeJacobian ( )
protectedvirtual

this simply calls upwind

Definition at line 99 of file Q2PPorepressureFlux.C.

100 {
101  upwind(false, true, _var.number());
102 }
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...
void Q2PPorepressureFlux::computeOffDiagJacobian ( unsigned int  jvar)
protectedvirtual

this simply calls upwind

Definition at line 105 of file Q2PPorepressureFlux.C.

106 {
107  upwind(false, true, jvar);
108 }
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...
Real Q2PPorepressureFlux::computeQpJac ( unsigned int  dvar)
protected

the derivative of the flux without the upstream mobility terms

Definition at line 111 of file Q2PPorepressureFlux.C.

Referenced by upwind().

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 }
const RichardsDensity & _density
fluid density
const MaterialProperty< RealTensorValue > & _permeability
permeability
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 MaterialProperty< RealVectorValue > & _gravity
gravity
Real Q2PPorepressureFlux::computeQpResidual ( )
protectedvirtual

Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quadpoints and then add the upwind mobility parts.

Definition at line 84 of file Q2PPorepressureFlux.C.

Referenced by upwind().

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 }
const RichardsDensity & _density
fluid density
const MaterialProperty< RealTensorValue > & _permeability
permeability
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 MaterialProperty< RealVectorValue > & _gravity
gravity
void Q2PPorepressureFlux::computeResidual ( )
protectedvirtual

This simply calls upwind.

Definition at line 93 of file Q2PPorepressureFlux.C.

94 {
95  upwind(true, false, 0);
96 }
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...
void Q2PPorepressureFlux::prepareNodalValues ( )
protected

calculates the nodal values of mobility, and derivatives thereof

Definition at line 56 of file Q2PPorepressureFlux.C.

Referenced by upwind().

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 }
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
Real _viscosity
fluid viscosity
const VariableValue & _sat
saturation at the nodes
const std::string density
Definition: NS.h:15
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...
unsigned int _num_nodes
number of nodes in the element
virtual Real ddensity(Real p) const =0
derivative of fluid density wrt porepressure This must be over-ridden in derived classes to provide a...
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
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the...
void Q2PPorepressureFlux::upwind ( bool  compute_res,
bool  compute_jac,
unsigned int  jvar 
)
protected

Do the upwinding for both the residual and jacobian I've put both calculations in the same code to try to reduce code duplication.

This is because when calculating the jacobian we need to calculate the residual to see which nodes are upwind and which are downwind

Definition at line 124 of file Q2PPorepressureFlux.C.

Referenced by computeJacobian(), computeOffDiagJacobian(), and computeResidual().

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 }
Real computeQpJac(unsigned int dvar)
the derivative of the flux without the upstream mobility terms
unsigned int _sat_var
variable number of the saturation variable
std::vector< Real > _dmobility_ds
d(_mobility)/d(saturation) These are used in the jacobian calculations
virtual Real computeQpResidual()
Note that this is not the complete residual for the quadpoint In computeResidual we sum over the quad...
unsigned int _num_nodes
number of nodes in the element
void prepareNodalValues()
calculates the nodal values of mobility, and derivatives thereof
std::vector< Real > _dmobility_dp
d(_mobility)/d(porepressure) These are used in the jacobian calculations
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the...

Member Data Documentation

const RichardsDensity& Q2PPorepressureFlux::_density
protected

fluid density

Definition at line 83 of file Q2PPorepressureFlux.h.

Referenced by computeQpJac(), computeQpResidual(), and prepareNodalValues().

std::vector<Real> Q2PPorepressureFlux::_dmobility_dp
protected

d(_mobility)/d(porepressure) These are used in the jacobian calculations

Definition at line 116 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

std::vector<Real> Q2PPorepressureFlux::_dmobility_ds
protected

d(_mobility)/d(saturation) These are used in the jacobian calculations

Definition at line 122 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

const MaterialProperty<RealVectorValue>& Q2PPorepressureFlux::_gravity
protected

gravity

Definition at line 98 of file Q2PPorepressureFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

std::vector<Real> Q2PPorepressureFlux::_mobility
protected

nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the residual

Definition at line 110 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

unsigned int Q2PPorepressureFlux::_num_nodes
protected

number of nodes in the element

Definition at line 104 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

const MaterialProperty<RealTensorValue>& Q2PPorepressureFlux::_permeability
protected

permeability

Definition at line 101 of file Q2PPorepressureFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

const RichardsRelPerm& Q2PPorepressureFlux::_relperm
protected

fluid relative permeability

Definition at line 92 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().

const VariableValue& Q2PPorepressureFlux::_sat
protected

saturation at the nodes

Definition at line 86 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().

unsigned int Q2PPorepressureFlux::_sat_var
protected

variable number of the saturation variable

Definition at line 89 of file Q2PPorepressureFlux.h.

Referenced by upwind().

Real Q2PPorepressureFlux::_viscosity
protected

fluid viscosity

Definition at line 95 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().


The documentation for this class was generated from the following files: