www.mooseframework.org
RichardsFlux.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 
8 #include "RichardsFlux.h"
9 
10 // MOOSE includes
11 #include "Material.h"
12 #include "MooseVariable.h"
13 
14 // C++ includes
15 #include <iostream>
16 
17 template <>
18 InputParameters
20 {
21  InputParameters params = validParams<Kernel>();
22  params.addParam<bool>(
23  "linear_shape_fcns",
24  true,
25  "If you are using second-order Lagrange shape functions you need to set this to false.");
26  params.addRequiredParam<UserObjectName>(
27  "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
28  return params;
29 }
30 
31 RichardsFlux::RichardsFlux(const InputParameters & parameters)
32  : Kernel(parameters),
33  _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
34  _pvar(_richards_name_UO.richards_var_num(_var.number())),
35 
36  // This kernel gets lots of things from the material
37  _flux(getMaterialProperty<std::vector<RealVectorValue>>("flux")),
38  _dflux_dv(getMaterialProperty<std::vector<std::vector<RealVectorValue>>>("dflux_dv")),
39  _dflux_dgradv(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>("dflux_dgradv")),
40  _d2flux_dvdv(
41  getMaterialProperty<std::vector<std::vector<std::vector<RealVectorValue>>>>("d2flux_dvdv")),
42  _d2flux_dgradvdv(getMaterialProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>(
43  "d2flux_dgradvdv")),
44  _d2flux_dvdgradv(getMaterialProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>(
45  "d2flux_dvdgradv")),
46 
47  _second_u(getParam<bool>("linear_shape_fcns")
48  ? _second_zero
49  : (_is_implicit ? _var.secondSln() : _var.secondSlnOld())),
50  _second_phi(getParam<bool>("linear_shape_fcns") ? _second_phi_zero : secondPhi()),
51 
52  _tauvel_SUPG(getMaterialProperty<std::vector<RealVectorValue>>("tauvel_SUPG")),
53  _dtauvel_SUPG_dgradv(
54  getMaterialProperty<std::vector<std::vector<RealTensorValue>>>("dtauvel_SUPG_dgradv")),
55  _dtauvel_SUPG_dv(
56  getMaterialProperty<std::vector<std::vector<RealVectorValue>>>("dtauvel_SUPG_dv"))
57 {
58 }
59 
60 Real
62 {
63  Real flux_part = _grad_test[_i][_qp] * _flux[_qp][_pvar];
64 
65  Real supg_test = _tauvel_SUPG[_qp][_pvar] * _grad_test[_i][_qp];
66  Real supg_kernel = 0.0;
67  if (supg_test != 0)
68  // NOTE: Libmesh does not correctly calculate grad(_grad_u) correctly, so following might not be
69  // correct
70  // NOTE: The following is -divergence(flux)
71  // NOTE: The following must be generalised if a non PPPP formalism is used
72  // NOTE: The generalisation will look like
73  // supg_kernel = sum_j {-(_dflux_dgradv[_qp][_pvar][j]*_second_u[_qp][j]).tr() -
74  // _dflux_dv[_qp][_pvar][j]*_grad_u[_qp][j]}
75  // where _grad_u[_qp][j] is the gradient of the j^th variable at the quadpoint.
76  supg_kernel = -(_dflux_dgradv[_qp][_pvar][_pvar] * _second_u[_qp]).tr() -
77  _dflux_dv[_qp][_pvar][_pvar] * _grad_u[_qp];
78 
79  return flux_part + supg_test * supg_kernel;
80 }
81 
82 Real
83 RichardsFlux::computeQpJac(unsigned int wrt_num)
84 {
85  Real flux_prime = _grad_test[_i][_qp] * (_dflux_dgradv[_qp][_pvar][wrt_num] * _grad_phi[_j][_qp] +
86  _dflux_dv[_qp][_pvar][wrt_num] * _phi[_j][_qp]);
87 
88  Real supg_test = _tauvel_SUPG[_qp][_pvar] * _grad_test[_i][_qp];
89  Real supg_test_prime =
90  _grad_phi[_j][_qp] * (_dtauvel_SUPG_dgradv[_qp][_pvar][wrt_num] * _grad_test[_i][_qp]) +
91  _phi[_j][_qp] * _dtauvel_SUPG_dv[_qp][_pvar][wrt_num] * _grad_test[_i][_qp];
92  Real supg_kernel = 0.0;
93  Real supg_kernel_prime = 0.0;
94 
95  if (supg_test != 0)
96  {
97  // NOTE: since Libmesh does not correctly calculate grad(_grad_u) correctly, so following might
98  // not be correct
99  supg_kernel = -(_dflux_dgradv[_qp][_pvar][_pvar] * _second_u[_qp]).tr() -
100  _dflux_dv[_qp][_pvar][_pvar] * _grad_u[_qp];
101 
102  // NOTE: just like supg_kernel, this must be generalised for non-PPPP formulations
103  supg_kernel_prime =
104  -(_d2flux_dvdv[_qp][_pvar][_pvar][wrt_num] * _phi[_j][_qp] * _grad_u[_qp] +
105  _phi[_j][_qp] * (_d2flux_dgradvdv[_qp][_pvar][_pvar][wrt_num] * _second_u[_qp]).tr() +
106  (_d2flux_dvdgradv[_qp][_pvar][_pvar][wrt_num] * _grad_u[_qp]) * _grad_phi[_j][_qp]);
107  if (wrt_num == _pvar)
108  supg_kernel_prime -= _dflux_dv[_qp][_pvar][_pvar] * _grad_phi[_j][_qp];
109  supg_kernel_prime -= (_dflux_dgradv[_qp][_pvar][_pvar] * _second_phi[_j][_qp]).tr();
110  }
111 
112  return flux_prime + supg_test_prime * supg_kernel + supg_test * supg_kernel_prime;
113 }
114 
115 Real
117 {
118  return computeQpJac(_pvar);
119 }
120 
121 Real
123 {
125  return 0.0;
126  unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
127  return computeQpJac(dvar);
128 }
const MaterialProperty< std::vector< std::vector< std::vector< RealVectorValue > > > > & _d2flux_dvdv
d^2(Richards flux_i)/d(variable_j)/d(variable_k), here flux_i is the i_th flux, which is itself a Rea...
Definition: RichardsFlux.h:61
unsigned int _pvar
the index of this variable in the list of Richards variables held by _richards_name_UO.
Definition: RichardsFlux.h:49
const MaterialProperty< std::vector< RealVectorValue > > & _tauvel_SUPG
SUPGtau*SUPGvel (a vector of these if multiphase)
Definition: RichardsFlux.h:76
const RichardsVarNames & _richards_name_UO
holds info regarding the names of the Richards variables and methods for extracting values of these v...
Definition: RichardsFlux.h:40
virtual Real computeQpOffDiagJacobian(unsigned int jvar)
Definition: RichardsFlux.C:122
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dflux_dgradv
d(Richards flux_i)/d(grad(variable_j)), here flux_i is the i_th flux, which is itself a RealVectorVal...
Definition: RichardsFlux.h:58
bool not_richards_var(unsigned int moose_var_num) const
returns true if moose_var_num is not a richards var
RichardsFlux(const InputParameters &parameters)
Definition: RichardsFlux.C:31
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dtauvel_SUPG_dgradv
derivative of SUPGtau*SUPGvel_i wrt grad(variable_j)
Definition: RichardsFlux.h:79
This holds maps between pressure_var or pressure_var, sat_var used in RichardsMaterial and kernels...
const VariablePhiSecond & _second_phi
interpolation function for the _second_u
Definition: RichardsFlux.h:73
Real computeQpJac(unsigned int wrt_num)
Computes diagonal and off-diagonal jacobian entries.
Definition: RichardsFlux.C:83
const MaterialProperty< std::vector< RealVectorValue > > & _flux
Richards flux.
Definition: RichardsFlux.h:52
const MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dflux_dv
d(Richards flux_i)/d(variable_j), here flux_i is the i_th flux, which is itself a RealVectorValue ...
Definition: RichardsFlux.h:55
virtual Real computeQpJacobian()
Definition: RichardsFlux.C:116
InputParameters validParams< RichardsFlux >()
Definition: RichardsFlux.C:19
const VariableSecond & _second_u
grad_i grad_j porepressure. This is used in SUPG
Definition: RichardsFlux.h:70
const MaterialProperty< std::vector< std::vector< RealVectorValue > > > & _dtauvel_SUPG_dv
derivative of SUPGtau*SUPGvel_i wrt variable_j
Definition: RichardsFlux.h:82
const MaterialProperty< std::vector< std::vector< std::vector< RealTensorValue > > > > & _d2flux_dvdgradv
d^2(Richards flux_i)/d(variable_j)/d(grad(variable_k)), here flux_i is the i_th flux, which is itself a RealVectorValue
Definition: RichardsFlux.h:67
const MaterialProperty< std::vector< std::vector< std::vector< RealTensorValue > > > > & _d2flux_dgradvdv
d^2(Richards flux_i)/d(grad(variable_j))/d(variable_k), here flux_i is the i_th flux, which is itself a RealVectorValue
Definition: RichardsFlux.h:64
virtual Real computeQpResidual()
Definition: RichardsFlux.C:61
unsigned int richards_var_num(unsigned int moose_var_num) const
the richards variable number