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

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

#include <Q2PSaturationFlux.h>

Inheritance diagram for Q2PSaturationFlux:
[legend]

Public Member Functions

 Q2PSaturationFlux (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 & _pp
 porepressure at the quadpoints More...
 
const VariableGradient & _grad_pp
 grad(porepressure) at the quadpoints More...
 
const VariableValue & _pp_nodal
 porepressure at the nodes More...
 
unsigned int _pp_var
 variable number of the porepressure 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 saturation.

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 Q2PSaturationFlux.h.

Constructor & Destructor Documentation

Q2PSaturationFlux::Q2PSaturationFlux ( const InputParameters &  parameters)

Definition at line 40 of file Q2PSaturationFlux.C.

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 }
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...
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
Real _viscosity
fluid viscosity
std::vector< Real > _dmobility_ds
d(_mobility)/d(saturation) These are used in the jacobian calculations
const VariableGradient & _grad_pp
grad(porepressure) at the quadpoints
const VariableValue & _pp
porepressure at the quadpoints
unsigned int _pp_var
variable number of the porepressure variable
const RichardsDensity & _density
fluid density
const MaterialProperty< RealVectorValue > & _gravity
gravity
const VariableValue & _pp_nodal
porepressure at the nodes
const MaterialProperty< RealTensorValue > & _permeability
permeability

Member Function Documentation

void Q2PSaturationFlux::computeJacobian ( )
protectedvirtual

this simply calls upwind

Definition at line 103 of file Q2PSaturationFlux.C.

104 {
105  upwind(false, true, _var.number());
106 }
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 Q2PSaturationFlux::computeOffDiagJacobian ( unsigned int  jvar)
protectedvirtual

this simply calls upwind

Definition at line 109 of file Q2PSaturationFlux.C.

110 {
111  upwind(false, true, jvar);
112 }
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 Q2PSaturationFlux::computeQpJac ( unsigned int  dvar)
protected

the derivative of the flux without the upstream mobility terms

Definition at line 115 of file Q2PSaturationFlux.C.

Referenced by upwind().

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 }
const VariableValue & _pp
porepressure at the quadpoints
unsigned int _pp_var
variable number of the porepressure variable
const RichardsDensity & _density
fluid density
const MaterialProperty< RealVectorValue > & _gravity
gravity
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< RealTensorValue > & _permeability
permeability
Real Q2PSaturationFlux::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 88 of file Q2PSaturationFlux.C.

Referenced by upwind().

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 }
const VariableGradient & _grad_pp
grad(porepressure) at the quadpoints
const VariableValue & _pp
porepressure at the quadpoints
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
const MaterialProperty< RealTensorValue > & _permeability
permeability
void Q2PSaturationFlux::computeResidual ( )
protectedvirtual

This simply calls upwind.

Definition at line 97 of file Q2PSaturationFlux.C.

98 {
99  upwind(true, false, 0);
100 }
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 Q2PSaturationFlux::prepareNodalValues ( )
protected

calculates the nodal values of mobility, and derivatives thereof

Definition at line 59 of file Q2PSaturationFlux.C.

Referenced by upwind().

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 }
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...
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
Real _viscosity
fluid viscosity
std::vector< Real > _dmobility_ds
d(_mobility)/d(saturation) These are used in the jacobian calculations
const std::string density
Definition: NS.h:15
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
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 ...
void Q2PSaturationFlux::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 128 of file Q2PSaturationFlux.C.

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

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
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the...
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 > _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 _pp_var
variable number of the porepressure variable
Real computeQpJac(unsigned int dvar)
the derivative of the flux without the upstream mobility terms

Member Data Documentation

const RichardsDensity& Q2PSaturationFlux::_density
protected

fluid density

Definition at line 83 of file Q2PSaturationFlux.h.

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

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

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

Definition at line 122 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

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

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

Definition at line 128 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

const VariableGradient& Q2PSaturationFlux::_grad_pp
protected

grad(porepressure) at the quadpoints

Definition at line 89 of file Q2PSaturationFlux.h.

Referenced by computeQpResidual().

const MaterialProperty<RealVectorValue>& Q2PSaturationFlux::_gravity
protected

gravity

Definition at line 104 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

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

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

Definition at line 116 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

unsigned int Q2PSaturationFlux::_num_nodes
protected

number of nodes in the element

Definition at line 110 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

const MaterialProperty<RealTensorValue>& Q2PSaturationFlux::_permeability
protected

permeability

Definition at line 107 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

const VariableValue& Q2PSaturationFlux::_pp
protected

porepressure at the quadpoints

Definition at line 86 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

const VariableValue& Q2PSaturationFlux::_pp_nodal
protected

porepressure at the nodes

Definition at line 92 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().

unsigned int Q2PSaturationFlux::_pp_var
protected

variable number of the porepressure variable

Definition at line 95 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and upwind().

const RichardsRelPerm& Q2PSaturationFlux::_relperm
protected

fluid relative permeability

Definition at line 98 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().

Real Q2PSaturationFlux::_viscosity
protected

fluid viscosity

Definition at line 101 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().


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