www.mooseframework.org
RichardsPiecewiseLinearSink.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 
9 
10 // MOOSE includes
11 #include "MooseVariable.h"
12 
13 // C++ includes
14 #include <iostream>
15 
16 template <>
17 InputParameters
19 {
20  InputParameters params = validParams<IntegratedBC>();
21  params.addRequiredParam<bool>(
22  "use_mobility",
23  "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), "
24  "where the '_nn' indicates the component normal to the boundary. In this "
25  "case bare_flux is measured in Pa.s^-1. This can be used in conjunction "
26  "with use_relperm.");
27  params.addRequiredParam<bool>("use_relperm",
28  "If true, then fluxes are multiplied by relative "
29  "permeability. This can be used in conjunction "
30  "with use_mobility");
31  params.addParam<std::vector<UserObjectName>>("relperm_UO",
32  "List of names of user objects that "
33  "define relative permeability. Only "
34  "needed if fully_upwind is used");
35  params.addParam<std::vector<UserObjectName>>(
36  "seff_UO",
37  "List of name of user objects that define effective saturation as a function of "
38  "pressure list. Only needed if fully_upwind is used");
39  params.addParam<std::vector<UserObjectName>>("density_UO",
40  "List of names of user objects that "
41  "define the fluid density. Only "
42  "needed if fully_upwind is used");
43  params.addRequiredParam<std::vector<Real>>(
44  "pressures", "Tuple of pressure values. Must be monotonically increasing.");
45  params.addRequiredParam<std::vector<Real>>(
46  "bare_fluxes",
47  "Tuple of flux values (measured in kg.m^-2.s^-1 for use_mobility=false, and "
48  "in Pa.s^-1 if use_mobility=true). This flux is OUT of the medium: hence "
49  "positive values of flux means this will be a SINK, while negative values "
50  "indicate this flux will be a SOURCE. A piecewise-linear fit is performed to "
51  "the (pressure,bare_fluxes) pairs to obtain the flux at any arbitrary "
52  "pressure, and the first or last bare_flux values are used if the quad-point "
53  "pressure falls outside this range.");
54  params.addParam<FunctionName>("multiplying_fcn",
55  1.0,
56  "If this function is provided, the flux "
57  "will be multiplied by this function. "
58  "This is useful for spatially or "
59  "temporally varying sinks");
60  params.addRequiredParam<UserObjectName>(
61  "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
62  params.addParam<bool>("fully_upwind", false, "Use full upwinding");
63  params.addParam<PostprocessorName>(
64  "area_pp",
65  1,
66  "An area postprocessor. If given, the bare_fluxes will be divided by this "
67  "quantity. This means the bare fluxes are measured in kg.s^-1. This is "
68  "useful for the case when you wish to provide the *total* flux, and let MOOSE "
69  "proportion it uniformly across the boundary. In that case you would have "
70  "use_mobility=false=use_relperm, and only one bare flux should be specified");
71  return params;
72 }
73 
75  : IntegratedBC(parameters),
76  _use_mobility(getParam<bool>("use_mobility")),
77  _use_relperm(getParam<bool>("use_relperm")),
78  _fully_upwind(getParam<bool>("fully_upwind")),
79 
80  _sink_func(getParam<std::vector<Real>>("pressures"),
81  getParam<std::vector<Real>>("bare_fluxes")),
82 
83  _m_func(getFunction("multiplying_fcn")),
84 
85  _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
86  _num_p(_richards_name_UO.num_v()),
87  _pvar(_richards_name_UO.richards_var_num(_var.number())),
88 
89  // in the following, getUserObjectByName returns a reference (an alias) to a RichardsBLAH user
90  // object, and the & turns it into a pointer
91  _density_UO(_fully_upwind
92  ? &getUserObjectByName<RichardsDensity>(
93  getParam<std::vector<UserObjectName>>("density_UO")[_pvar])
94  : NULL),
95  _seff_UO(_fully_upwind
96  ? &getUserObjectByName<RichardsSeff>(
97  getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])
98  : NULL),
99  _relperm_UO(_fully_upwind
100  ? &getUserObjectByName<RichardsRelPerm>(
101  getParam<std::vector<UserObjectName>>("relperm_UO")[_pvar])
102  : NULL),
103 
104  _area_pp(getPostprocessorValue("area_pp")),
105 
106  _num_nodes(0),
107  _nodal_density(0),
108  _dnodal_density_dv(0),
109  _nodal_relperm(0),
110  _dnodal_relperm_dv(0),
111 
112  _pp(getMaterialProperty<std::vector<Real>>("porepressure")),
113  _dpp_dv(getMaterialProperty<std::vector<std::vector<Real>>>("dporepressure_dv")),
114 
115  _viscosity(getMaterialProperty<std::vector<Real>>("viscosity")),
116  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
117 
118  _dseff_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ds_eff_dv")),
119 
120  _rel_perm(getMaterialProperty<std::vector<Real>>("rel_perm")),
121  _drel_perm_dv(getMaterialProperty<std::vector<std::vector<Real>>>("drel_perm_dv")),
122 
123  _density(getMaterialProperty<std::vector<Real>>("density")),
124  _ddensity_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ddensity_dv"))
125 {
126  _ps_at_nodes.resize(_num_p);
127  for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
129 }
130 
131 void
133 {
134  // NOTE: i'm assuming that all the richards variables are pressure values
135 
136  _num_nodes = (*_ps_at_nodes[_pvar]).size();
137 
138  Real p;
139  Real seff;
140  std::vector<Real> dseff_dp;
141  Real drelperm_ds;
142 
143  _nodal_density.resize(_num_nodes);
145  _nodal_relperm.resize(_num_nodes);
147  dseff_dp.resize(_num_p);
148  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
149  {
150  // retrieve and calculate basic things at the node
151  p = (*_ps_at_nodes[_pvar])[nodenum]; // pressure of fluid _pvar at node nodenum
152 
153  _nodal_density[nodenum] = _density_UO->density(p); // density of fluid _pvar at node nodenum
154  _dnodal_density_dv[nodenum].resize(_num_p);
155  for (unsigned int ph = 0; ph < _num_p; ++ph)
156  _dnodal_density_dv[nodenum][ph] = 0;
157  _dnodal_density_dv[nodenum][_pvar] = _density_UO->ddensity(p); // d(density)/dP
158 
159  seff = _seff_UO->seff(_ps_at_nodes,
160  nodenum); // effective saturation of fluid _pvar at node nodenum
161  _seff_UO->dseff(
162  _ps_at_nodes, nodenum, dseff_dp); // d(seff)/d(P_ph), for ph = 0, ..., _num_p - 1
163 
164  _nodal_relperm[nodenum] =
165  _relperm_UO->relperm(seff); // relative permeability of fluid _pvar at node nodenum
166  drelperm_ds = _relperm_UO->drelperm(seff); // d(relperm)/dseff
167 
168  _dnodal_relperm_dv[nodenum].resize(_num_p);
169  for (unsigned int ph = 0; ph < _num_p; ++ph)
170  _dnodal_relperm_dv[nodenum][ph] = drelperm_ds * dseff_dp[ph];
171  }
172 }
173 
174 void
176 {
177  if (_fully_upwind)
179  IntegratedBC::computeResidual();
180 }
181 
182 Real
184 {
185  Real flux = 0;
186  Real k = 0;
187 
188  if (!_fully_upwind)
189  {
190  flux = _test[_i][_qp] * _sink_func.sample(_pp[_qp][_pvar]);
191  if (_use_mobility)
192  {
193  k = (_permeability[_qp] * _normals[_qp]) * _normals[_qp];
194  flux *= _density[_qp][_pvar] * k / _viscosity[_qp][_pvar];
195  }
196  if (_use_relperm)
197  flux *= _rel_perm[_qp][_pvar];
198  }
199  else
200  {
201  flux = _test[_i][_qp] * _sink_func.sample((*_ps_at_nodes[_pvar])[_i]);
202  if (_use_mobility)
203  {
204  k = (_permeability[0] * _normals[_qp]) * _normals[_qp]; // assume that _permeability is
205  // constant throughout element so
206  // doesn't need to be upwinded
207  flux *= _nodal_density[_i] * k /
208  _viscosity[0][_pvar]; // assume that viscosity is constant throughout element
209  }
210  if (_use_relperm)
211  flux *= _nodal_relperm[_i];
212  }
213 
214  flux *= _m_func.value(_t, _q_point[_qp]);
215 
216  if (_area_pp == 0.0)
217  {
218  if (flux != 0)
219  mooseError("RichardsPiecewiseLinearSink: flux is nonzero, but area is zero!\n");
220  // if flux == 0, then leave it as zero.
221  }
222  else
223  flux /= _area_pp;
224 
225  return flux;
226 }
227 
228 void
230 {
231  if (_fully_upwind)
233  IntegratedBC::computeJacobian();
234 }
235 
236 Real
238 {
239  return jac(_pvar);
240 }
241 
242 void
244 {
245  if (_fully_upwind)
247  IntegratedBC::computeJacobianBlock(jvar);
248 }
249 
250 Real
252 {
254  return 0.0;
255  unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
256  return jac(dvar);
257 }
258 
259 Real
260 RichardsPiecewiseLinearSink::jac(unsigned int wrt_num)
261 {
262  Real flux = 0;
263  Real deriv = 0;
264  Real k = 0;
265  Real mob = 0;
266  Real mobp = 0;
267  Real phi = 0;
268 
269  if (!_fully_upwind)
270  {
271  flux = _sink_func.sample(_pp[_qp][_pvar]);
272  deriv = _sink_func.sampleDerivative(_pp[_qp][_pvar]) * _dpp_dv[_qp][_pvar][wrt_num];
273  phi = _phi[_j][_qp];
274  if (_use_mobility)
275  {
276  k = (_permeability[_qp] * _normals[_qp]) * _normals[_qp];
277  mob = _density[_qp][_pvar] * k / _viscosity[_qp][_pvar];
278  mobp = _ddensity_dv[_qp][_pvar][wrt_num] * k / _viscosity[_qp][_pvar];
279  deriv = mob * deriv + mobp * flux;
280  flux *= mob;
281  }
282  if (_use_relperm)
283  deriv = _rel_perm[_qp][_pvar] * deriv + _drel_perm_dv[_qp][_pvar][wrt_num] * flux;
284  }
285  else
286  {
287  if (_i != _j)
288  return 0.0; // residual at node _i only depends on variables at that node
289  flux = _sink_func.sample((*_ps_at_nodes[_pvar])[_i]);
290  deriv = (_pvar == wrt_num ? _sink_func.sampleDerivative((*_ps_at_nodes[_pvar])[_i])
291  : 0); // NOTE: i'm assuming that the variables are pressure variables
292  phi = 1;
293  if (_use_mobility)
294  {
295  k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
296  mob = _nodal_density[_i] * k / _viscosity[0][_pvar];
297  mobp = _dnodal_density_dv[_i][wrt_num] * k / _viscosity[0][_pvar];
298  deriv = mob * deriv + mobp * flux;
299  flux *= mob;
300  }
301  if (_use_relperm)
302  deriv = _nodal_relperm[_i] * deriv + _dnodal_relperm_dv[_i][wrt_num] * flux;
303  }
304 
305  deriv *= _m_func.value(_t, _q_point[_qp]);
306 
307  if (_area_pp == 0.0)
308  {
309  if (deriv != 0)
310  mooseError("RichardsPiecewiseLinearSink: deriv is nonzero, but area is zero!\n");
311  // if deriv == 0, then leave it as zero.
312  }
313  else
314  deriv /= _area_pp;
315 
316  return _test[_i][_qp] * deriv * phi;
317 }
const MaterialProperty< std::vector< Real > > & _rel_perm
relative permeability (only the _pvar component is used)
virtual Real computeQpOffDiagJacobian(unsigned int jvar)
LinearInterpolation _sink_func
piecewise-linear function of porepressure (this defines the strength of the sink) ...
const MaterialProperty< std::vector< Real > > & _viscosity
viscosity (only the _pvar component is used)
RichardsPiecewiseLinearSink(const InputParameters &parameters)
virtual Real drelperm(Real seff) const =0
derivative of relative permeability wrt effective saturation This must be over-ridden in your derived...
bool _use_relperm
whether to multiply the sink flux by relative permeability
bool _use_mobility
whether to multiply the sink flux by permeability*density/viscosity
Base class for effective saturation as a function of porepressure(s) The functions seff...
Definition: RichardsSeff.h:22
Base class for Richards relative permeability classes that provide relative permeability as a functio...
std::vector< Real > _nodal_density
nodal values of fluid density These are used if _fully_upwind = true
const MaterialProperty< std::vector< Real > > & _density
fluid density (only the _pvar component is used)
const VariableValue * nodal_var(unsigned int richards_var_num) const
The nodal variable values for the given richards_var_num To extract a the value of pressure variable ...
bool not_richards_var(unsigned int moose_var_num) const
returns true if moose_var_num is not a richards var
std::vector< Real > _nodal_relperm
nodal values of relative permeability These are used if _fully_upwind = true
virtual void dseff(std::vector< const VariableValue * > p, unsigned int qp, std::vector< Real > &result) const =0
derivative(s) of effective saturation as a function of porepressure(s) at given quadpoint of the elem...
This holds maps between pressure_var or pressure_var, sat_var used in RichardsMaterial and kernels...
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...
std::vector< std::vector< Real > > _dnodal_density_dv
d(_nodal_density)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jac...
InputParameters validParams< RichardsPiecewiseLinearSink >()
const RichardsVarNames & _richards_name_UO
holds info about the names and values of richards variable in the simulation
const RichardsSeff * _seff_UO
user object defining the effective saturation. Only used if _fully_upwind = true
const PostprocessorValue & _area_pp
area postprocessor. if given then all bare_fluxes are divided by this quantity
std::vector< const VariableValue * > _ps_at_nodes
Holds the values of pressures at all the nodes of the element Only used if _fully_upwind = true Eg: _...
const MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
d(relperm_i)/d(variable_j)
virtual void computeJacobianBlock(unsigned int jvar)
unsigned int _pvar
the moose internal variable number corresponding to the porepressure of this sink flux ...
virtual Real ddensity(Real p) const =0
derivative of fluid density wrt porepressure This must be over-ridden in derived classes to provide a...
bool _fully_upwind
whether to use full upwinding
const MaterialProperty< std::vector< Real > > & _pp
porepressure values (only the _pvar component is used)
const MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
unsigned int _num_nodes
number of nodes in this element. Only used if _fully_upwind = true
const RichardsRelPerm * _relperm_UO
user object defining the relative permeability. Only used if _fully_upwind = true ...
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_p
number of richards variables
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
void prepareNodalValues()
calculates the nodal values of pressure, mobility, and derivatives thereof
std::vector< std::vector< Real > > _dnodal_relperm_dv
d(_nodal_relperm)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jac...
Real jac(unsigned int wrt_num)
derivative of residual wrt the wrt_num Richards variable
const MaterialProperty< RealTensorValue > & _permeability
permeability
const RichardsDensity * _density_UO
user object defining the density. Only used if _fully_upwind = true
const MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
d(density_i)/d(variable_j)
virtual Real seff(std::vector< const VariableValue * > p, unsigned int qp) const =0
effective saturation as a function of porepressure(s) at given quadpoint of the element ...
unsigned int richards_var_num(unsigned int moose_var_num) const
the richards variable number
Function & _m_func
sink flux gets multiplied by this function