www.mooseframework.org
RichardsBorehole.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 "RichardsBorehole.h"
11 #include "RotationMatrix.h"
12 
14 
17 {
19  params.addRequiredParam<UserObjectName>(
20  "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
21  params.addParam<std::vector<UserObjectName>>("relperm_UO",
22  "List of names of user objects that "
23  "define relative permeability. Only "
24  "needed if fully_upwind is used");
25  params.addParam<std::vector<UserObjectName>>(
26  "seff_UO",
27  "List of name of user objects that define effective saturation as a function of "
28  "pressure list. Only needed if fully_upwind is used");
29  params.addParam<std::vector<UserObjectName>>("density_UO",
30  "List of names of user objects that "
31  "define the fluid density. Only "
32  "needed if fully_upwind is used");
33  params.addParam<bool>("fully_upwind", false, "Fully upwind the flux");
34  params.addClassDescription("Approximates a borehole in the mesh with given bottomhole pressure, "
35  "and radii using a number of point sinks whose positions are read "
36  "from a file");
37  return params;
38 }
39 
41  : PeacemanBorehole(parameters),
42  _fully_upwind(getParam<bool>("fully_upwind")),
43  _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
44  _num_p(_richards_name_UO.num_v()),
45  _pvar(_richards_name_UO.richards_var_num(_var.number())),
46 
47  // in the following, getUserObjectByName returns a reference (an alias) to a RichardsBLAH user
48  // object, and the & turns it into a pointer
49  _density_UO(_fully_upwind ? &getUserObjectByName<RichardsDensity>(
50  getParam<std::vector<UserObjectName>>("density_UO")[_pvar])
51  : NULL),
52  _seff_UO(_fully_upwind ? &getUserObjectByName<RichardsSeff>(
53  getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])
54  : NULL),
55  _relperm_UO(_fully_upwind ? &getUserObjectByName<RichardsRelPerm>(
56  getParam<std::vector<UserObjectName>>("relperm_UO")[_pvar])
57  : NULL),
58 
59  _num_nodes(0),
60  _mobility(0),
61  _dmobility_dv(0),
62  _pp(getMaterialProperty<std::vector<Real>>("porepressure")),
63  _dpp_dv(getMaterialProperty<std::vector<std::vector<Real>>>("dporepressure_dv")),
64  _viscosity(getMaterialProperty<std::vector<Real>>("viscosity")),
65  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
66  _dseff_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ds_eff_dv")),
67  _rel_perm(getMaterialProperty<std::vector<Real>>("rel_perm")),
68  _drel_perm_dv(getMaterialProperty<std::vector<std::vector<Real>>>("drel_perm_dv")),
69  _density(getMaterialProperty<std::vector<Real>>("density")),
70  _ddensity_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ddensity_dv"))
71 {
72  _ps_at_nodes.resize(_num_p);
73  for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
75 
76  // To correctly compute the Jacobian terms,
77  // tell MOOSE that this DiracKernel depends on all the Richards Vars
78  const std::vector<MooseVariableFEBase *> & coupled_vars = _richards_name_UO.getCoupledMooseVars();
79  for (unsigned int i = 0; i < coupled_vars.size(); i++)
80  addMooseVariableDependency(coupled_vars[i]);
81 }
82 
83 void
85 {
86  // NOTE: i'm assuming that all the richards variables are pressure values
87 
88  _num_nodes = (*_ps_at_nodes[_pvar]).size();
89 
90  Real p;
91  Real density;
92  Real ddensity_dp;
93  Real seff;
94  std::vector<Real> dseff_dp;
95  Real relperm;
96  Real drelperm_ds;
97  _mobility.resize(_num_nodes);
98  _dmobility_dv.resize(_num_nodes);
99  dseff_dp.resize(_num_p);
100  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
101  {
102  // retrieve and calculate basic things at the node
103  p = (*_ps_at_nodes[_pvar])[nodenum]; // pressure of fluid _pvar at node nodenum
104  density = _density_UO->density(p); // density of fluid _pvar at node nodenum
105  ddensity_dp = _density_UO->ddensity(p); // d(density)/dP
106  seff = _seff_UO->seff(_ps_at_nodes,
107  nodenum); // effective saturation of fluid _pvar at node nodenum
108  _seff_UO->dseff(
109  _ps_at_nodes, nodenum, dseff_dp); // d(seff)/d(P_ph), for ph = 0, ..., _num_p - 1
110  relperm = _relperm_UO->relperm(seff); // relative permeability of fluid _pvar at node nodenum
111  drelperm_ds = _relperm_UO->drelperm(seff); // d(relperm)/dseff
112 
113  // calculate the mobility and its derivatives wrt (variable_ph = porepressure_ph)
114  _mobility[nodenum] =
115  density * relperm / _viscosity[0][_pvar]; // assume viscosity is constant throughout element
116  _dmobility_dv[nodenum].resize(_num_p);
117  for (unsigned int ph = 0; ph < _num_p; ++ph)
118  _dmobility_dv[nodenum][ph] = density * drelperm_ds * dseff_dp[ph] / _viscosity[0][_pvar];
119  _dmobility_dv[nodenum][_pvar] += ddensity_dp * relperm / _viscosity[0][_pvar];
120  }
121 }
122 
123 void
125 {
126  if (_fully_upwind)
129 }
130 
131 Real
133 {
134  const Real character = _character.value(_t, _q_point[_qp]);
135  if (character == 0.0)
136  return 0.0;
137 
138  const Real bh_pressure =
139  _p_bot +
140  _unit_weight *
141  (_q_point[_qp] -
142  _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
143 
144  Real pp;
145  Real mob;
146  if (!_fully_upwind)
147  {
148  pp = _pp[_qp][_pvar];
150  }
151  else
152  {
153  pp = (*_ps_at_nodes[_pvar])[_i];
154  mob = _mobility[_i];
155  }
156 
157  // Get the ID we initially assigned to this point
158  const unsigned current_dirac_ptid = currentPointCachedID();
159 
160  // If getting the ID failed, fall back to the old bodge!
161  // if (current_dirac_ptid == libMesh::invalid_uint)
162  // current_dirac_ptid = (_zs.size() > 2) ? 1 : 0;
163 
164  Real outflow(0.0); // this is the flow rate from porespace out of the system
165 
166  Real wc(0.0);
167  if (current_dirac_ptid > 0)
168  // contribution from half-segment "behind" this point (must have >1 point for
169  // current_dirac_ptid>0)
170  {
172  _rot_matrix[current_dirac_ptid - 1],
173  _half_seg_len[current_dirac_ptid - 1],
175  _rs[current_dirac_ptid]);
176  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
177  // injection, so outflow<0 || // production, so outflow>0
178  outflow += _test[_i][_qp] * std::abs(character) * wc * mob * (pp - bh_pressure);
179  }
180 
181  if (current_dirac_ptid + 1 < _zs.size() || _zs.size() == 1)
182  // contribution from half-segment "ahead of" this point, or we only have one point
183  {
185  _rot_matrix[current_dirac_ptid],
186  _half_seg_len[current_dirac_ptid],
188  _rs[current_dirac_ptid]);
189  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
190  // injection, so outflow<0 || // production, so outflow>0
191  outflow += _test[_i][_qp] * std::abs(character) * wc * mob * (pp - bh_pressure);
192  }
193 
195  outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
196  return outflow;
197 }
198 
199 void
201 {
202  if (_fully_upwind)
205 }
206 
207 Real
209 {
210  const Real character = _character.value(_t, _q_point[_qp]);
211  if (character == 0.0)
212  return 0.0;
213  return jac(_pvar);
214 }
215 
216 Real
218 {
220  return 0.0;
221  const unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
222  return jac(dvar);
223 }
224 
225 Real
226 RichardsBorehole::jac(unsigned int wrt_num)
227 {
228  const Real character = _character.value(_t, _q_point[_qp]);
229  if (character == 0.0)
230  return 0.0;
231 
232  const Real bh_pressure =
233  _p_bot +
234  _unit_weight *
235  (_q_point[_qp] -
236  _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
237 
238  Real pp;
239  Real dpp_dv;
240  Real mob;
241  Real dmob_dv;
242  Real phi;
243  if (!_fully_upwind)
244  {
245  pp = _pp[_qp][_pvar];
246  dpp_dv = _dpp_dv[_qp][_pvar][wrt_num];
248  dmob_dv = (_drel_perm_dv[_qp][_pvar][wrt_num] * _density[_qp][_pvar] +
249  _rel_perm[_qp][_pvar] * _ddensity_dv[_qp][_pvar][wrt_num]) /
250  _viscosity[_qp][_pvar];
251  phi = _phi[_j][_qp];
252  }
253  else
254  {
255  if (_i != _j)
256  return 0.0; // residual at node _i only depends on variables at that node
257  pp = (*_ps_at_nodes[_pvar])[_i];
258  dpp_dv =
259  (_pvar == wrt_num ? 1 : 0); // NOTE: i'm assuming that the variables are pressure variables
260  mob = _mobility[_i];
261  dmob_dv = _dmobility_dv[_i][wrt_num];
262  phi = 1;
263  }
264 
265  // Get the ID we initially assigned to this point
266  const unsigned current_dirac_ptid = currentPointCachedID();
267 
268  // If getting the ID failed, fall back to the old bodge!
269  // if (current_dirac_ptid == libMesh::invalid_uint)
270  // current_dirac_ptid = (_zs.size() > 2) ? 1 : 0;
271 
272  Real outflowp(0.0);
273 
274  Real wc(0.0);
275  if (current_dirac_ptid > 0)
276  // contribution from half-segment "behind" this point
277  {
279  _rot_matrix[current_dirac_ptid - 1],
280  _half_seg_len[current_dirac_ptid - 1],
282  _rs[current_dirac_ptid]);
283  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
284  // injection, so outflow<0 || // production, so outflow>0
285  outflowp += _test[_i][_qp] * std::abs(character) * wc *
286  (mob * phi * dpp_dv + dmob_dv * phi * (pp - bh_pressure));
287  }
288 
289  if (current_dirac_ptid < _zs.size() - 1 || _zs.size() == 1)
290  // contribution from half-segment "ahead of" this point
291  {
293  _rot_matrix[current_dirac_ptid],
294  _half_seg_len[current_dirac_ptid],
296  _rs[current_dirac_ptid]);
297  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
298  // injection, so outflow<0 || // production, so outflow>0
299  outflowp += _test[_i][_qp] * std::abs(character) * wc *
300  (mob * phi * dpp_dv + dmob_dv * phi * (pp - bh_pressure));
301  }
302 
303  return outflowp;
304 }
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...
const MaterialProperty< RealTensorValue > & _permeability
material permeability
const bool _fully_upwind
Whether to use full upwinding.
const MaterialProperty< std::vector< Real > > & _density
fluid density
Real jac(unsigned int wrt_num)
Calculates Jacobian.
const MooseArray< Point > & _q_point
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...
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
Real wellConstant(const RealTensorValue &perm, const RealTensorValue &rot, const Real &half_len, const Elem *ele, const Real &rad)
Calculates Peaceman&#39;s form of the borehole well constant Z Chen, Y Zhang, Well flow models for variou...
const MaterialProperty< std::vector< Real > > & _viscosity
fluid viscosity
Base class for effective saturation as a function of porepressure(s) The functions seff...
Definition: RichardsSeff.h:18
Base class for Richards relative permeability classes that provide relative permeability as a functio...
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 ...
std::vector< Real > _rs
radii of the borehole
const OutputTools< T >::VariableTestValue & _test
const OutputTools< T >::VariablePhiValue & _phi
virtual Real computeQpJacobian()
Computes the diagonal part of the jacobian.
static const std::string density
Definition: NS.h:33
unsigned int _num_nodes
number of nodes in this element. Only used if _fully_upwind = true
bool not_richards_var(unsigned int moose_var_num) const
returns true if moose_var_num is not a richards var
unsigned currentPointCachedID()
Approximates a borehole by a sequence of Dirac Points.
This holds maps between pressure_var or pressure_var, sat_var used in RichardsMaterial and kernels...
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: _...
static InputParameters validParams()
Creates a new PeacemanBorehole This reads the file containing the lines of the form radius x y z that...
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are used if _fully_upwind = true ...
void addRequiredParam(const std::string &name, const std::string &doc_string)
const MaterialProperty< std::vector< Real > > & _pp
fluid porepressure (or porepressures in case of multiphase)
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
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 RichardsRelPerm *const _relperm_UO
user object defining the relative permeability. Only used if _fully_upwind = true ...
TensorValue< Real > RealTensorValue
virtual void computeResidual()
Computes the residual.
const RichardsDensity *const _density_UO
user object defining the density. Only used if _fully_upwind = true
static InputParameters validParams()
Creates a new RichardsBorehole This sets all the variables, but also reads the file containing the li...
const MaterialProperty< std::vector< Real > > & _rel_perm
relative permeability
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 ...
const RealVectorValue _unit_weight
unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point) ...
const RichardsSeff *const _seff_UO
user object defining the effective saturation. 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 ...
virtual Real computeQpOffDiagJacobian(unsigned int jvar)
Computes the off-diagonal part of the jacobian Note: at March2014 this is never called since moose do...
void add(Real contrib)
adds contrib to _total
RichardsSumQuantity & _total_outflow_mass
This is used to hold the total fluid flowing into the borehole Hence, it is positive for production w...
const Elem *const & _current_elem
const std::vector< MooseVariableFieldBase *> & getCoupledMooseVars() const
const MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
d(relperm_i)/d(variable_j)
const Real _p_bot
bottomhole pressure of borehole
virtual void computeJacobian()
Computes the Jacobian.
void addMooseVariableDependency(MooseVariableFieldBase *var)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int _j
RichardsBorehole(const InputParameters &parameters)
unsigned int richards_var_num(unsigned int moose_var_num) const
the richards variable number
const unsigned int _pvar
The moose internal variable number of the richards variable of this Dirac Kernel. ...
unsigned int _qp
virtual void computeJacobian() override
std::vector< std::vector< Real > > _dmobility_dv
d(_mobility)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jacobian...
const MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
d(porepressure_i)/d(variable_j)
void addClassDescription(const std::string &doc_string)
Point _bottom_point
the bottom point of the borehole (where bottom_pressure is defined)
Base class for fluid density as a function of porepressure The functions density, ddensity and d2dens...
std::vector< RealTensorValue > _rot_matrix
rotation matrix used in well_constant calculation
const RichardsVarNames & _richards_name_UO
Defines the richards variables in the simulation.
std::vector< Real > _half_seg_len
0.5*(length of polyline segments between points)
virtual Real value(Real t, const Point &p) const
std::vector< Real > _zs
z points of borehole
unsigned int _i
virtual Real computeQpResidual()
Computes the Qp residual.
void prepareNodalValues()
calculates the nodal values of pressure, mobility, and derivatives thereof
registerMooseObject("RichardsApp", RichardsBorehole)
const unsigned int _num_p
number of richards variables
const MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
d(density_i)/d(variable_j)
virtual void computeResidual() override
const Function & _character
If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure...
Approximates a borehole by a sequence of Dirac Points.