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 () override
 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 () override
 This simply calls upwind. More...
 
virtual void computeOffDiagJacobian (MooseVariableFEBase &jvar) override
 this simply calls upwind More...
 
virtual void computeJacobian () override
 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 47 of file Q2PSaturationFlux.h.

Constructor & Destructor Documentation

◆ Q2PSaturationFlux()

Q2PSaturationFlux::Q2PSaturationFlux ( const InputParameters &  parameters)

Definition at line 44 of file Q2PSaturationFlux.C.

45  : Kernel(parameters),
46  _density(getUserObject<RichardsDensity>("fluid_density")),
47  _pp(coupledValue("porepressure_variable")),
48  _grad_pp(coupledGradient("porepressure_variable")),
49  _pp_nodal(coupledNodalValue("porepressure_variable")),
50  _pp_var(coupled("porepressure_variable")),
51  _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
52  _viscosity(getParam<Real>("fluid_viscosity")),
53  _gravity(getMaterialProperty<RealVectorValue>("gravity")),
54  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
55  _num_nodes(0),
56  _mobility(0),
57  _dmobility_dp(0),
58  _dmobility_ds(0)
59 {
60 }
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

◆ computeJacobian()

void Q2PSaturationFlux::computeJacobian ( )
overrideprotectedvirtual

this simply calls upwind

Definition at line 107 of file Q2PSaturationFlux.C.

108 {
109  upwind(false, true, _var.number());
110 }
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...

◆ computeOffDiagJacobian()

void Q2PSaturationFlux::computeOffDiagJacobian ( MooseVariableFEBase &  jvar)
overrideprotectedvirtual

this simply calls upwind

Definition at line 113 of file Q2PSaturationFlux.C.

114 {
115  upwind(false, true, jvar.number());
116 }
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...

◆ computeQpJac()

Real Q2PSaturationFlux::computeQpJac ( unsigned int  dvar)
protected

the derivative of the flux without the upstream mobility terms

Definition at line 119 of file Q2PSaturationFlux.C.

Referenced by upwind().

120 {
121  // this is just the derivative of the flux WITHOUT the upstream mobility terms
122  // Those terms get added in during computeJacobian()
123  if (dvar == _pp_var)
124  return _grad_test[_i][_qp] *
125  (_permeability[_qp] *
126  (_grad_phi[_j][_qp] - _density.ddensity(_pp[_qp]) * _gravity[_qp] * _phi[_j][_qp]));
127  else
128  return 0;
129 }
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
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 MaterialProperty< RealTensorValue > & _permeability
permeability

◆ computeQpResidual()

Real Q2PSaturationFlux::computeQpResidual ( )
overrideprotectedvirtual

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 92 of file Q2PSaturationFlux.C.

Referenced by upwind().

93 {
94  // note this is not the complete residual:
95  // the upwind mobility parts get added in computeResidual
96  return _grad_test[_i][_qp] *
97  (_permeability[_qp] * (_grad_pp[_qp] - _density.density(_pp[_qp]) * _gravity[_qp]));
98 }
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

◆ computeResidual()

void Q2PSaturationFlux::computeResidual ( )
overrideprotectedvirtual

This simply calls upwind.

Definition at line 101 of file Q2PSaturationFlux.C.

102 {
103  upwind(true, false, 0);
104 }
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...

◆ prepareNodalValues()

void Q2PSaturationFlux::prepareNodalValues ( )
protected

calculates the nodal values of mobility, and derivatives thereof

Definition at line 63 of file Q2PSaturationFlux.C.

Referenced by upwind().

64 {
65  _num_nodes = _pp_nodal.size();
66 
67  Real density;
68  Real ddensity_dp;
69  Real relperm;
70  Real drelperm_ds;
71 
72  _mobility.resize(_num_nodes);
73  _dmobility_dp.resize(_num_nodes);
74  _dmobility_ds.resize(_num_nodes);
75 
76  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
77  {
78  density = _density.density(_pp_nodal[nodenum]); // fluid density at the node
79  ddensity_dp = _density.ddensity(_pp_nodal[nodenum]); // d(fluid density)/dP at the node
80  relperm = _relperm.relperm(
81  _var.dofValues()[nodenum]); // relative permeability of the fluid at node nodenum
82  drelperm_ds = _relperm.drelperm(_var.dofValues()[nodenum]); // d(relperm)/dsat
83 
84  // calculate the mobility and its derivatives wrt P and S
85  _mobility[nodenum] = density * relperm / _viscosity;
86  _dmobility_dp[nodenum] = ddensity_dp * relperm / _viscosity;
87  _dmobility_ds[nodenum] = density * drelperm_ds / _viscosity;
88  }
89 }
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...
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 _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:17
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 relperm(Real seff) const =0
relative permeability as a function of effective saturation This must be over-ridden in your derived ...
const VariableValue & _pp_nodal
porepressure at the nodes

◆ upwind()

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 132 of file Q2PSaturationFlux.C.

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

133 {
134  if (compute_jac && !(jvar == _var.number() || jvar == _pp_var))
135  return;
136 
137  // calculate the mobility values and their derivatives
139 
140  // compute the residual without the mobility terms
141  // Even if we are computing the jacobian we still need this
142  // in order to see which nodes are upwind and which are downwind
143  DenseVector<Number> & re = _assembly.residualBlock(_var.number());
144  _local_re.resize(re.size());
145  _local_re.zero();
146 
147  for (_i = 0; _i < _test.size(); _i++)
148  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
149  _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
150 
151  DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar);
152  if (compute_jac)
153  {
154  _local_ke.resize(ke.m(), ke.n());
155  _local_ke.zero();
156 
157  for (_i = 0; _i < _test.size(); _i++)
158  for (_j = 0; _j < _phi.size(); _j++)
159  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
160  _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJac(jvar);
161  }
162 
163  // Now perform the upwinding by multiplying the residuals at the
164  // upstream nodes by their mobilities
165  //
166  // The residual for the kernel is the darcy flux.
167  // This is
168  // R_i = int{mobility*flux_no_mob} = int{mobility*grad(pot)*permeability*grad(test_i)}
169  // for node i. where int is the integral over the element.
170  // However, in fully-upwind, the first step is to take the mobility outside the
171  // integral, which was done in the _local_re calculation above.
172  //
173  // NOTE: Physically _local_re(_i) is a measure of fluid flowing out of node i
174  // If we had left in mobility, it would be exactly the mass flux flowing out of node i.
175  //
176  // This leads to the definition of upwinding:
177  // ***
178  // If _local_re(i) is positive then we use mobility_i. That is
179  // we use the upwind value of mobility.
180  // ***
181  //
182  // The final subtle thing is we must also conserve fluid mass: the total mass
183  // flowing out of node i must be the sum of the masses flowing
184  // into the other nodes.
185 
186  // FIRST:
187  // this is a dirty way of getting around precision loss problems
188  // and problems at steadystate where upwinding oscillates from
189  // node-to-node causing nonconvergence.
190  // I'm not sure if i actually need to do this in moose. Certainly
191  // in cosflow it is necessary.
192  // I will code a better algorithm if necessary
193  bool reached_steady = true;
194  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
195  {
196  if (_local_re(nodenum) >= 1E-20)
197  {
198  reached_steady = false;
199  break;
200  }
201  }
202  reached_steady = false;
203 
204  // DEFINE VARIABLES USED TO ENSURE MASS CONSERVATION
205  // total mass out - used for mass conservation
206  Real total_mass_out = 0;
207  // total flux in
208  Real total_in = 0;
209 
210  // the following holds derivatives of these
211  std::vector<Real> dtotal_mass_out;
212  std::vector<Real> dtotal_in;
213  if (compute_jac)
214  {
215  dtotal_mass_out.assign(_num_nodes, 0);
216  dtotal_in.assign(_num_nodes, 0);
217  }
218 
219  // PERFORM THE UPWINDING!
220  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
221  {
222  if (_local_re(nodenum) >= 0 || reached_steady) // upstream node
223  {
224  if (compute_jac)
225  {
226  for (_j = 0; _j < _phi.size(); _j++)
227  _local_ke(nodenum, _j) *= _mobility[nodenum];
228  if (jvar == _var.number())
229  // deriv wrt S
230  _local_ke(nodenum, nodenum) += _dmobility_ds[nodenum] * _local_re(nodenum);
231  else
232  // deriv wrt P
233  _local_ke(nodenum, nodenum) += _dmobility_dp[nodenum] * _local_re(nodenum);
234  for (_j = 0; _j < _phi.size(); _j++)
235  dtotal_mass_out[_j] += _local_ke(nodenum, _j);
236  }
237  _local_re(nodenum) *= _mobility[nodenum];
238  total_mass_out += _local_re(nodenum);
239  }
240  else
241  {
242  total_in -= _local_re(nodenum); // note the -= means the result is positive
243  if (compute_jac)
244  for (_j = 0; _j < _phi.size(); _j++)
245  dtotal_in[_j] -= _local_ke(nodenum, _j);
246  }
247  }
248 
249  // CONSERVE MASS
250  // proportion the total_mass_out mass to the inflow nodes, weighting by their _local_re values
251  if (!reached_steady)
252  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
253  if (_local_re(nodenum) < 0)
254  {
255  if (compute_jac)
256  for (_j = 0; _j < _phi.size(); _j++)
257  {
258  _local_ke(nodenum, _j) *= total_mass_out / total_in;
259  _local_ke(nodenum, _j) +=
260  _local_re(nodenum) * (dtotal_mass_out[_j] / total_in -
261  dtotal_in[_j] * total_mass_out / total_in / total_in);
262  }
263  _local_re(nodenum) *= total_mass_out / total_in;
264  }
265 
266  // ADD RESULTS TO RESIDUAL OR JACOBIAN
267  if (compute_res)
268  {
269  re += _local_re;
270 
271  if (_has_save_in)
272  {
273  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
274  for (unsigned int i = 0; i < _save_in.size(); i++)
275  _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
276  }
277  }
278 
279  if (compute_jac)
280  {
281  ke += _local_ke;
282 
283  if (_has_diag_save_in && jvar == _var.number())
284  {
285  const unsigned int rows = ke.m();
286  DenseVector<Number> diag(rows);
287  for (unsigned int i = 0; i < rows; i++)
288  diag(i) = _local_ke(i, i);
289 
290  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
291  for (unsigned int i = 0; i < _diag_save_in.size(); i++)
292  _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
293  }
294  }
295 }
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
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
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

◆ _density

const RichardsDensity& Q2PSaturationFlux::_density
protected

fluid density

Definition at line 86 of file Q2PSaturationFlux.h.

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

◆ _dmobility_dp

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

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

Definition at line 125 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _dmobility_ds

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

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

Definition at line 131 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _grad_pp

const VariableGradient& Q2PSaturationFlux::_grad_pp
protected

grad(porepressure) at the quadpoints

Definition at line 92 of file Q2PSaturationFlux.h.

Referenced by computeQpResidual().

◆ _gravity

const MaterialProperty<RealVectorValue>& Q2PSaturationFlux::_gravity
protected

gravity

Definition at line 107 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

◆ _mobility

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 119 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _num_nodes

unsigned int Q2PSaturationFlux::_num_nodes
protected

number of nodes in the element

Definition at line 113 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _permeability

const MaterialProperty<RealTensorValue>& Q2PSaturationFlux::_permeability
protected

permeability

Definition at line 110 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

◆ _pp

const VariableValue& Q2PSaturationFlux::_pp
protected

porepressure at the quadpoints

Definition at line 89 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

◆ _pp_nodal

const VariableValue& Q2PSaturationFlux::_pp_nodal
protected

porepressure at the nodes

Definition at line 95 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().

◆ _pp_var

unsigned int Q2PSaturationFlux::_pp_var
protected

variable number of the porepressure variable

Definition at line 98 of file Q2PSaturationFlux.h.

Referenced by computeQpJac(), and upwind().

◆ _relperm

const RichardsRelPerm& Q2PSaturationFlux::_relperm
protected

fluid relative permeability

Definition at line 101 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().

◆ _viscosity

Real Q2PSaturationFlux::_viscosity
protected

fluid viscosity

Definition at line 104 of file Q2PSaturationFlux.h.

Referenced by prepareNodalValues().


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