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 () 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 & _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 47 of file Q2PPorepressureFlux.h.

Constructor & Destructor Documentation

◆ Q2PPorepressureFlux()

Q2PPorepressureFlux::Q2PPorepressureFlux ( const InputParameters &  parameters)

Definition at line 43 of file Q2PPorepressureFlux.C.

44  : Kernel(parameters),
45  _density(getUserObject<RichardsDensity>("fluid_density")),
46  _sat(coupledNodalValue("saturation_variable")),
47  _sat_var(coupled("saturation_variable")),
48  _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
49  _viscosity(getParam<Real>("fluid_viscosity")),
50  _gravity(getMaterialProperty<RealVectorValue>("gravity")),
51  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
52  _num_nodes(0),
53  _mobility(0),
54  _dmobility_dp(0),
55  _dmobility_ds(0)
56 {
57 }
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

◆ computeJacobian()

void Q2PPorepressureFlux::computeJacobian ( )
overrideprotectedvirtual

this simply calls upwind

Definition at line 103 of file Q2PPorepressureFlux.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...

◆ computeOffDiagJacobian()

void Q2PPorepressureFlux::computeOffDiagJacobian ( MooseVariableFEBase &  jvar)
overrideprotectedvirtual

this simply calls upwind

Definition at line 109 of file Q2PPorepressureFlux.C.

110 {
111  upwind(false, true, jvar.number());
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...

◆ computeQpJac()

Real Q2PPorepressureFlux::computeQpJac ( unsigned int  dvar)
protected

the derivative of the flux without the upstream mobility terms

Definition at line 115 of file Q2PPorepressureFlux.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 == _var.number())
120  return _grad_test[_i][_qp] *
121  (_permeability[_qp] *
122  (_grad_phi[_j][_qp] - _density.ddensity(_u[_qp]) * _gravity[_qp] * _phi[_j][_qp]));
123  else
124  return 0;
125 }
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 RichardsDensity & _density
fluid density
const MaterialProperty< RealTensorValue > & _permeability
permeability
const MaterialProperty< RealVectorValue > & _gravity
gravity

◆ computeQpResidual()

Real Q2PPorepressureFlux::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 88 of file Q2PPorepressureFlux.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_u[_qp] - _density.density(_u[_qp]) * _gravity[_qp]));
94 }
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

◆ computeResidual()

void Q2PPorepressureFlux::computeResidual ( )
overrideprotectedvirtual

This simply calls upwind.

Definition at line 97 of file Q2PPorepressureFlux.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...

◆ prepareNodalValues()

void Q2PPorepressureFlux::prepareNodalValues ( )
protected

calculates the nodal values of mobility, and derivatives thereof

Definition at line 60 of file Q2PPorepressureFlux.C.

Referenced by upwind().

61 {
62  _num_nodes = _sat.size();
63 
64  Real density;
65  Real ddensity_dp;
66  Real relperm;
67  Real drelperm_ds;
68 
69  _mobility.resize(_num_nodes);
70  _dmobility_dp.resize(_num_nodes);
71  _dmobility_ds.resize(_num_nodes);
72 
73  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
74  {
75  density = _density.density(_var.dofValues()[nodenum]); // fluid density at the node
76  ddensity_dp = _density.ddensity(_var.dofValues()[nodenum]); // d(fluid density)/dP at the node
77  relperm = _relperm.relperm(_sat[nodenum]); // relative permeability of the fluid at node nodenum
78  drelperm_ds = _relperm.drelperm(_sat[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 }
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...
const RichardsDensity & _density
fluid density
Real _viscosity
fluid viscosity
const VariableValue & _sat
saturation at the nodes
const std::string density
Definition: NS.h:17
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...
virtual Real relperm(Real seff) const =0
relative permeability as a function of 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
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are multiplied by _flux_no_mob to give the...

◆ upwind()

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 128 of file Q2PPorepressureFlux.C.

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

129 {
130  if (compute_jac && !(jvar == _var.number() || jvar == _sat_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 P
226  _local_ke(nodenum, nodenum) += _dmobility_dp[nodenum] * _local_re(nodenum);
227  else
228  // deriv wrt S
229  _local_ke(nodenum, nodenum) += _dmobility_ds[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 }
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
unsigned int _num_nodes
number of nodes in the element
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 > _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

◆ _density

const RichardsDensity& Q2PPorepressureFlux::_density
protected

fluid density

Definition at line 86 of file Q2PPorepressureFlux.h.

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

◆ _dmobility_dp

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

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

Definition at line 119 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _dmobility_ds

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

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

Definition at line 125 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _gravity

const MaterialProperty<RealVectorValue>& Q2PPorepressureFlux::_gravity
protected

gravity

Definition at line 101 of file Q2PPorepressureFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

◆ _mobility

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 113 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _num_nodes

unsigned int Q2PPorepressureFlux::_num_nodes
protected

number of nodes in the element

Definition at line 107 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues(), and upwind().

◆ _permeability

const MaterialProperty<RealTensorValue>& Q2PPorepressureFlux::_permeability
protected

permeability

Definition at line 104 of file Q2PPorepressureFlux.h.

Referenced by computeQpJac(), and computeQpResidual().

◆ _relperm

const RichardsRelPerm& Q2PPorepressureFlux::_relperm
protected

fluid relative permeability

Definition at line 95 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().

◆ _sat

const VariableValue& Q2PPorepressureFlux::_sat
protected

saturation at the nodes

Definition at line 89 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().

◆ _sat_var

unsigned int Q2PPorepressureFlux::_sat_var
protected

variable number of the saturation variable

Definition at line 92 of file Q2PPorepressureFlux.h.

Referenced by upwind().

◆ _viscosity

Real Q2PPorepressureFlux::_viscosity
protected

fluid viscosity

Definition at line 98 of file Q2PPorepressureFlux.h.

Referenced by prepareNodalValues().


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