www.mooseframework.org
Q2PPiecewiseLinearSink.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "Q2PPiecewiseLinearSink.h"
11 
12 // MOOSE includes
13 #include "MooseVariable.h"
14 
15 // C++ includes
16 #include <iostream>
17 
19 
22 {
24  params.addRequiredParam<bool>(
25  "use_mobility",
26  "If true, then fluxes are multiplied by (density*permeability_nn/viscosity), "
27  "where the '_nn' indicates the component normal to the boundary. In this "
28  "case bare_flux is measured in Pa.s^-1. This can be used in conjunction "
29  "with use_relperm.");
30  params.addRequiredParam<bool>("use_relperm",
31  "If true, then fluxes are multiplied by relative "
32  "permeability. This can be used in conjunction "
33  "with use_mobility");
34  params.addRequiredParam<std::vector<Real>>(
35  "pressures", "Tuple of pressure values. Must be monotonically increasing.");
36  params.addRequiredParam<std::vector<Real>>(
37  "bare_fluxes",
38  "Tuple of flux values (measured in kg.m^-2.s^-1 for use_mobility=false, and "
39  "in Pa.s^-1 if use_mobility=true). This flux is OUT of the medium: hence "
40  "positive values of flux means this will be a SINK, while negative values "
41  "indicate this flux will be a SOURCE. A piecewise-linear fit is performed to "
42  "the (pressure,bare_fluxes) pairs to obtain the flux at any arbitrary "
43  "pressure, and the first or last bare_flux values are used if the quad-point "
44  "pressure falls outside this range.");
45  params.addParam<FunctionName>("multiplying_fcn",
46  1.0,
47  "If this function is provided, the flux "
48  "will be multiplied by this function. "
49  "This is useful for spatially or "
50  "temporally varying sinks");
51  params.addRequiredParam<UserObjectName>(
52  "fluid_density",
53  "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
54  params.addRequiredParam<UserObjectName>(
55  "fluid_relperm",
56  "A RichardsRelPerm UserObject (eg RichardsRelPermPower) that defines the "
57  "fluid relative permeability as a function of the saturation Variable.");
58  params.addRequiredCoupledVar("other_var",
59  "The other variable in the 2-phase system. If "
60  "Variable=porepressure, the other_var=saturation, and "
61  "vice-versa.");
62  params.addRequiredParam<bool>("var_is_porepressure",
63  "This flag is needed to correctly calculate the Jacobian entries. "
64  "If set to true, this Sink will extract fluid from the phase with "
65  "porepressure as its Variable (usually the liquid phase). If set "
66  "to false, this Sink will extract fluid from the phase with "
67  "saturation as its variable (usually the gas phase)");
68  params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
69  params.addClassDescription("Sink of fluid, controlled by (pressure, bare_fluxes) interpolation. "
70  "This is for use in Q2P models");
71  return params;
72 }
73 
75  : IntegratedBC(parameters),
76  _use_mobility(getParam<bool>("use_mobility")),
77  _use_relperm(getParam<bool>("use_relperm")),
78  _sink_func(getParam<std::vector<Real>>("pressures"),
79  getParam<std::vector<Real>>("bare_fluxes")),
80  _m_func(getFunction("multiplying_fcn")),
81  _density(getUserObject<RichardsDensity>("fluid_density")),
82  _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
83  _other_var_nodal(coupledDofValues("other_var")),
84  _other_var_num(coupled("other_var")),
85  _var_is_pp(getParam<bool>("var_is_porepressure")),
86  _viscosity(getParam<Real>("fluid_viscosity")),
87  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
88  _num_nodes(0),
89  _pp(0),
90  _sat(0),
91  _nodal_density(0),
92  _dnodal_density_dp(0),
93  _nodal_relperm(0),
94  _dnodal_relperm_ds(0)
95 {
96 }
97 
98 void
100 {
101  _num_nodes = _other_var_nodal.size();
102 
103  // set _pp and _sat variables
104  _pp.resize(_num_nodes);
105  _sat.resize(_num_nodes);
106  if (_var_is_pp)
107  {
108  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
109  {
110  _pp[nodenum] = _var.dofValues()[nodenum];
111  _sat[nodenum] = _other_var_nodal[nodenum];
112  }
113  }
114  else
115  {
116  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
117  {
118  _pp[nodenum] = _other_var_nodal[nodenum];
119  _sat[nodenum] = _var.dofValues()[nodenum];
120  }
121  }
122 
123  // calculate derived things
124  if (_use_mobility)
125  {
126  _nodal_density.resize(_num_nodes);
128  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
129  {
130  _nodal_density[nodenum] = _density.density(_pp[nodenum]);
131  _dnodal_density_dp[nodenum] = _density.ddensity(_pp[nodenum]);
132  }
133  }
134 
135  if (_use_relperm)
136  {
137  _nodal_relperm.resize(_num_nodes);
139  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
140  {
141  _nodal_relperm[nodenum] = _relperm.relperm(_sat[nodenum]);
142  _dnodal_relperm_ds[nodenum] = _relperm.drelperm(_sat[nodenum]);
143  }
144  }
145 }
146 
147 void
149 {
152 }
153 
154 Real
156 {
157  Real flux = 0;
158  Real k = 0;
159 
160  flux = _test[_i][_qp] * _sink_func.sample(_pp[_i]);
161  if (_use_mobility)
162  {
163  k = (_permeability[0] * _normals[_qp]) * _normals[_qp]; // assume that _permeability is constant
164  // throughout element so doesn't need to
165  // be upwinded
166  flux *= _nodal_density[_i] * k / _viscosity;
167  }
168  if (_use_relperm)
169  flux *= _nodal_relperm[_i];
170 
171  flux *= _m_func.value(_t, _q_point[_qp]);
172 
173  return flux;
174 }
175 
176 void
178 {
181 }
182 
183 Real
185 {
186  return jac(_var.number());
187 }
188 
189 void
191 {
194 }
195 
196 Real
198 {
199  if (jvar == _var.number() || jvar == _other_var_num)
200  return jac(jvar);
201 
202  return 0.0;
203 }
204 
205 Real
206 Q2PPiecewiseLinearSink::jac(unsigned int wrt_num)
207 {
208  Real flux = 0;
209  Real deriv = 0;
210  Real k = 0;
211  Real mob = 0;
212  Real mobp = 0;
213  Real phi = 1;
214 
215  if (_i != _j)
216  return 0.0; // residual at node _i only depends on variables at that node
217 
218  flux = _sink_func.sample(_pp[_i]);
219 
220  if (_var_is_pp)
221  {
222  // derivative of the _sink_func
223  if (wrt_num == _var.number())
225  else
226  deriv = 0;
227 
228  // add derivative of the mobility
229  if (_use_mobility)
230  {
231  k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
232  mob = _nodal_density[_i] * k / _viscosity;
233  if (wrt_num == _var.number())
234  mobp = _dnodal_density_dp[_i] * k / _viscosity; // else mobp = 0
235  deriv = mob * deriv + mobp * flux;
236  flux *= mob;
237  }
238 
239  // add derivative of the relperm
240  if (_use_relperm)
241  {
242  if (wrt_num == _other_var_num)
244  else
246  }
247  }
248  else
249  {
250  // derivative of the _sink_func
251  if (wrt_num == _other_var_num)
253  else
254  deriv = 0;
255 
256  // add derivative of the mobility
257  if (_use_mobility)
258  {
259  k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
260  mob = _nodal_density[_i] * k / _viscosity;
261  if (wrt_num == _other_var_num)
262  mobp = _dnodal_density_dp[_i] * k / _viscosity; // else mobp = 0
263  deriv = mob * deriv + mobp * flux;
264  flux *= mob;
265  }
266 
267  // add derivative of the relperm
268  if (_use_relperm)
269  {
270  if (wrt_num == _var.number())
272  else
274  }
275  }
276 
278 
279  return _test[_i][_qp] * deriv * phi;
280 }
const VariableTestValue & _test
virtual void computeOffDiagJacobian(unsigned int jvar) override
void prepareNodalValues()
calculates the nodal values of pressure, mobility, and derivatives thereof
virtual void computeJacobian() override
virtual void computeResidual() override
unsigned int _j
virtual void computeOffDiagJacobian(unsigned int jvar) override
virtual Real drelperm(Real seff) const =0
derivative of relative permeability wrt effective saturation This must be over-ridden in your derived...
const MooseArray< Point > & _normals
std::vector< Real > _nodal_density
nodal values of fluid density
LinearInterpolation _sink_func
piecewise-linear function of porepressure (this defines the strength of the sink) ...
virtual Real ddensity(Real p) const =0
derivative of fluid density wrt porepressure This must be over-ridden in derived classes to provide a...
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
bool _use_mobility
whether to multiply the sink flux by permeability*density/viscosity
Q2PPiecewiseLinearSink(const InputParameters &parameters)
unsigned int number() const
Base class for Richards relative permeability classes that provide relative permeability as a functio...
static InputParameters validParams()
T sample(const T &x) const
unsigned int _i
Real jac(unsigned int wrt_num)
derivative of residual wrt the wrt_num variable
unsigned int _other_var_num
the variable number of the other variable
static InputParameters validParams()
Applies a fully-upwinded flux sink to a boundary The sink is a piecewise linear function of porepress...
void addRequiredParam(const std::string &name, const std::string &doc_string)
virtual Real computeQpResidual() override
unsigned int _qp
bool _var_is_pp
whether the Variable for this BC is porepressure or not
virtual void computeJacobian() override
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 MooseArray< Point > & _q_point
TensorValue< Real > RealTensorValue
T sampleDerivative(const T &x) const
Real deriv(unsigned n, unsigned alpha, unsigned beta, Real x)
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
registerMooseObject("RichardsApp", Q2PPiecewiseLinearSink)
virtual Real relperm(Real seff) const =0
relative permeability as a function of effective saturation This must be over-ridden in your derived ...
std::vector< Real > _dnodal_density_dp
d(_nodal_density)/d(porepressure)
std::vector< Real > _pp
nodal values of porepressure
const RichardsRelPerm & _relperm
fluid relative permeability
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
MooseVariable & _var
virtual Real computeQpJacobian() override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const DoFValue & dofValues() const override
const MaterialProperty< RealTensorValue > & _permeability
permeability
bool _use_relperm
whether to multiply the sink flux by relative permeability
void addClassDescription(const std::string &doc_string)
std::vector< Real > _dnodal_relperm_ds
d(_nodal_relperm)/d(saturation)
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
virtual void computeResidual() override
virtual Real value(Real t, const Point &p) const
unsigned int _num_nodes
number of nodes in this element.
std::vector< Real > _nodal_relperm
nodal values of relative permeability
std::vector< Real > _sat
nodal values of saturation
const RichardsDensity & _density
fluid density
const Function & _m_func
sink flux gets multiplied by this function
static const std::string k
Definition: NS.h:124
const VariableValue & _other_var_nodal
the other variable in the 2-phase system (this is saturation if Variable=porepressure, and viceversa)