www.mooseframework.org
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
RichardsBorehole Class Reference

Approximates a borehole by a sequence of Dirac Points. More...

#include <RichardsBorehole.h>

Inheritance diagram for RichardsBorehole:
[legend]

Public Member Functions

 RichardsBorehole (const InputParameters &parameters)
 Creates a new RichardsBorehole This sets all the variables, but also reads the file containing the lines of the form radius x y z that defines the borehole geometry. More...
 
virtual void computeResidual ()
 Computes the residual. More...
 
virtual Real computeQpResidual ()
 Computes the Qp residual. More...
 
virtual void computeJacobian ()
 Computes the Jacobian. More...
 
virtual Real computeQpJacobian ()
 Computes the diagonal part of the jacobian. More...
 
virtual Real computeQpOffDiagJacobian (unsigned int jvar)
 Computes the off-diagonal part of the jacobian Note: at March2014 this is never called since moose does not have this functionality. More...
 

Protected Member Functions

void prepareNodalValues ()
 calculates the nodal values of pressure, mobility, and derivatives thereof More...
 
Real jac (unsigned int wrt_num)
 Calculates Jacobian. More...
 
virtual void addPoints ()
 Add Dirac Points to the borehole. More...
 
bool parseNextLineReals (std::ifstream &ifs, std::vector< Real > &myvec)
 reads a space-separated line of floats from ifs and puts in myvec More...
 
Real wellConstant (const RealTensorValue &perm, const RealTensorValue &rot, const Real &half_len, const Elem *ele, const Real &rad)
 Calculates Peaceman's form of the borehole well constant Z Chen, Y Zhang, Well flow models for various numerical methods, Int J Num Analysis and Modeling, 3 (2008) 375-388. More...
 

Protected Attributes

const bool _fully_upwind
 Whether to use full upwinding. More...
 
const RichardsVarNames_richards_name_UO
 Defines the richards variables in the simulation. More...
 
const unsigned int _num_p
 number of richards variables More...
 
const unsigned int _pvar
 The moose internal variable number of the richards variable of this Dirac Kernel. More...
 
const RichardsDensity_density_UO
 user object defining the density. Only used if _fully_upwind = true More...
 
const RichardsSeff_seff_UO
 user object defining the effective saturation. Only used if _fully_upwind = true More...
 
const RichardsRelPerm_relperm_UO
 user object defining the relative permeability. Only used if _fully_upwind = true More...
 
unsigned int _num_nodes
 number of nodes in this element. Only used if _fully_upwind = true More...
 
std::vector< Real > _mobility
 nodal values of mobility = density*relperm/viscosity These are used if _fully_upwind = true More...
 
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 calculations if _fully_upwind = true More...
 
const MaterialProperty< std::vector< Real > > & _pp
 fluid porepressure (or porepressures in case of multiphase) More...
 
const MaterialProperty< std::vector< std::vector< Real > > > & _dpp_dv
 d(porepressure_i)/d(variable_j) More...
 
const MaterialProperty< std::vector< Real > > & _viscosity
 fluid viscosity More...
 
const MaterialProperty< RealTensorValue > & _permeability
 material permeability More...
 
const MaterialProperty< std::vector< std::vector< Real > > > & _dseff_dv
 deriviatves of Seff wrt variables More...
 
const MaterialProperty< std::vector< Real > > & _rel_perm
 relative permeability More...
 
const MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
 d(relperm_i)/d(variable_j) More...
 
const MaterialProperty< std::vector< Real > > & _density
 fluid density More...
 
const MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
 d(density_i)/d(variable_j) More...
 
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: _ps_at_nodes[_pvar] is a pointer to this variable's nodal porepressure values So: (*_ps_at_nodes[_pvar])[i] = _var.dofValues()[i] = porepressure of pressure-variable _pvar at node i. More...
 
Function & _character
 If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure, and does nothing otherwise If negative then the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise The flow rate to/from the borehole is multiplied by |character|, so usually character = +/- 1. More...
 
const Real _p_bot
 bottomhole pressure of borehole More...
 
const RealVectorValue _unit_weight
 unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point) More...
 
RichardsSumQuantity_total_outflow_mass
 This is used to hold the total fluid flowing into the borehole Hence, it is positive for production wells where fluid is flowing from porespace into the borehole and removed from the model. More...
 
std::vector< Real > _rs
 radii of the borehole More...
 
std::vector< Real > _xs
 x points of the borehole More...
 
std::vector< Real > _ys
 y points of the borehole More...
 
std::vector< Real > _zs
 z points of borehole More...
 
Point _bottom_point
 the bottom point of the borehole (where bottom_pressure is defined) More...
 
std::vector< Real > _half_seg_len
 0.5*(length of polyline segments between points) More...
 
std::vector< RealTensorValue > _rot_matrix
 rotation matrix used in well_constant calculation More...
 

Detailed Description

Approximates a borehole by a sequence of Dirac Points.

Definition at line 28 of file RichardsBorehole.h.

Constructor & Destructor Documentation

◆ RichardsBorehole()

RichardsBorehole::RichardsBorehole ( const InputParameters &  parameters)

Creates a new RichardsBorehole This sets all the variables, but also reads the file containing the lines of the form radius x y z that defines the borehole geometry.

It also calculates segment-lengths and rotation matrices needed for computing the borehole well constant

Definition at line 41 of file RichardsBorehole.C.

42  : PeacemanBorehole(parameters),
43  _fully_upwind(getParam<bool>("fully_upwind")),
44  _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
47 
48  // in the following, getUserObjectByName returns a reference (an alias) to a RichardsBLAH user
49  // object, and the & turns it into a pointer
51  ? &getUserObjectByName<RichardsDensity>(
52  getParam<std::vector<UserObjectName>>("density_UO")[_pvar])
53  : NULL),
55  ? &getUserObjectByName<RichardsSeff>(
56  getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])
57  : NULL),
59  ? &getUserObjectByName<RichardsRelPerm>(
60  getParam<std::vector<UserObjectName>>("relperm_UO")[_pvar])
61  : NULL),
62 
63  _num_nodes(0),
64  _mobility(0),
65  _dmobility_dv(0),
66  _pp(getMaterialProperty<std::vector<Real>>("porepressure")),
67  _dpp_dv(getMaterialProperty<std::vector<std::vector<Real>>>("dporepressure_dv")),
68  _viscosity(getMaterialProperty<std::vector<Real>>("viscosity")),
69  _permeability(getMaterialProperty<RealTensorValue>("permeability")),
70  _dseff_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ds_eff_dv")),
71  _rel_perm(getMaterialProperty<std::vector<Real>>("rel_perm")),
72  _drel_perm_dv(getMaterialProperty<std::vector<std::vector<Real>>>("drel_perm_dv")),
73  _density(getMaterialProperty<std::vector<Real>>("density")),
74  _ddensity_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ddensity_dv"))
75 {
76  _ps_at_nodes.resize(_num_p);
77  for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
79 
80  // To correctly compute the Jacobian terms,
81  // tell MOOSE that this DiracKernel depends on all the Richards Vars
82  const std::vector<MooseVariableFEBase *> & coupled_vars = _richards_name_UO.getCoupledMooseVars();
83  for (unsigned int i = 0; i < coupled_vars.size(); i++)
84  addMooseVariableDependency(coupled_vars[i]);
85 }
const MaterialProperty< RealTensorValue > & _permeability
material permeability
const bool _fully_upwind
Whether to use full upwinding.
const RichardsDensity * _density_UO
user object defining the density. Only used if _fully_upwind = true
const MaterialProperty< std::vector< Real > > & _density
fluid density
PeacemanBorehole(const InputParameters &parameters)
Creates a new PeacemanBorehole This reads the file containing the lines of the form radius x y z that...
const MaterialProperty< std::vector< Real > > & _viscosity
fluid viscosity
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 ...
const RichardsSeff * _seff_UO
user object defining the effective saturation. Only used if _fully_upwind = true
const MaterialProperty< std::vector< std::vector< Real > > > & _dseff_dv
deriviatves of Seff wrt variables
unsigned int _num_nodes
number of nodes in this element. Only used if _fully_upwind = true
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: _...
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are used if _fully_upwind = true ...
unsigned int num_v() const
the number of porepressure variables
const MaterialProperty< std::vector< Real > > & _pp
fluid porepressure (or porepressures in case of multiphase)
const MaterialProperty< std::vector< Real > > & _rel_perm
relative permeability
const RichardsRelPerm * _relperm_UO
user object defining the relative permeability. Only used if _fully_upwind = true ...
const MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
d(relperm_i)/d(variable_j)
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. ...
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)
const RichardsVarNames & _richards_name_UO
Defines the richards variables in the simulation.
const unsigned int _num_p
number of richards variables
const MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
d(density_i)/d(variable_j)

Member Function Documentation

◆ addPoints()

void PeacemanBorehole::addPoints ( )
protectedvirtualinherited

Add Dirac Points to the borehole.

Definition at line 180 of file PeacemanBorehole.C.

181 {
182  // This function gets called just before the DiracKernel is evaluated
183  // so this is a handy place to zero this out.
185 
186  // Add point using the unique ID "i", let the DiracKernel take
187  // care of the caching. This should be fast after the first call,
188  // as long as the points don't move around.
189  for (unsigned int i = 0; i < _zs.size(); i++)
190  addPoint(Point(_xs[i], _ys[i], _zs[i]), i);
191 }
std::vector< Real > _xs
x points of the borehole
std::vector< Real > _ys
y points of the borehole
void zero()
sets _total = 0
RichardsSumQuantity & _total_outflow_mass
This is used to hold the total fluid flowing into the borehole Hence, it is positive for production w...
std::vector< Real > _zs
z points of borehole

◆ computeJacobian()

void RichardsBorehole::computeJacobian ( )
virtual

Computes the Jacobian.

This just calls prepareNodalValues, if _fully_upwind then calls DiracKernel::computeJacobian

Definition at line 204 of file RichardsBorehole.C.

205 {
206  if (_fully_upwind)
208  DiracKernel::computeJacobian();
209 }
const bool _fully_upwind
Whether to use full upwinding.
void prepareNodalValues()
calculates the nodal values of pressure, mobility, and derivatives thereof

◆ computeQpJacobian()

Real RichardsBorehole::computeQpJacobian ( )
virtual

Computes the diagonal part of the jacobian.

Definition at line 212 of file RichardsBorehole.C.

213 {
214  const Real character = _character.value(_t, _q_point[_qp]);
215  if (character == 0.0)
216  return 0.0;
217  return jac(_pvar);
218 }
Real jac(unsigned int wrt_num)
Calculates Jacobian.
const unsigned int _pvar
The moose internal variable number of the richards variable of this Dirac Kernel. ...
Function & _character
If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure...

◆ computeQpOffDiagJacobian()

Real RichardsBorehole::computeQpOffDiagJacobian ( unsigned int  jvar)
virtual

Computes the off-diagonal part of the jacobian Note: at March2014 this is never called since moose does not have this functionality.

Hence as of March2014 this has never been tested.

Definition at line 221 of file RichardsBorehole.C.

222 {
224  return 0.0;
225  const unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
226  return jac(dvar);
227 }
Real jac(unsigned int wrt_num)
Calculates Jacobian.
bool not_richards_var(unsigned int moose_var_num) const
returns true if moose_var_num is not a richards var
unsigned int richards_var_num(unsigned int moose_var_num) const
the richards variable number
const RichardsVarNames & _richards_name_UO
Defines the richards variables in the simulation.

◆ computeQpResidual()

Real RichardsBorehole::computeQpResidual ( )
virtual

Computes the Qp residual.

Definition at line 136 of file RichardsBorehole.C.

137 {
138  const Real character = _character.value(_t, _q_point[_qp]);
139  if (character == 0.0)
140  return 0.0;
141 
142  const Real bh_pressure =
143  _p_bot +
144  _unit_weight *
145  (_q_point[_qp] -
146  _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
147 
148  Real pp;
149  Real mob;
150  if (!_fully_upwind)
151  {
152  pp = _pp[_qp][_pvar];
153  mob = _rel_perm[_qp][_pvar] * _density[_qp][_pvar] / _viscosity[_qp][_pvar];
154  }
155  else
156  {
157  pp = (*_ps_at_nodes[_pvar])[_i];
158  mob = _mobility[_i];
159  }
160 
161  // Get the ID we initially assigned to this point
162  const unsigned current_dirac_ptid = currentPointCachedID();
163 
164  // If getting the ID failed, fall back to the old bodge!
165  // if (current_dirac_ptid == libMesh::invalid_uint)
166  // current_dirac_ptid = (_zs.size() > 2) ? 1 : 0;
167 
168  Real outflow(0.0); // this is the flow rate from porespace out of the system
169 
170  Real wc(0.0);
171  if (current_dirac_ptid > 0)
172  // contribution from half-segment "behind" this point (must have >1 point for
173  // current_dirac_ptid>0)
174  {
175  wc = wellConstant(_permeability[_qp],
176  _rot_matrix[current_dirac_ptid - 1],
177  _half_seg_len[current_dirac_ptid - 1],
178  _current_elem,
179  _rs[current_dirac_ptid]);
180  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
181  // injection, so outflow<0 || // production, so outflow>0
182  outflow += _test[_i][_qp] * std::abs(character) * wc * mob * (pp - bh_pressure);
183  }
184 
185  if (current_dirac_ptid + 1 < _zs.size() || _zs.size() == 1)
186  // contribution from half-segment "ahead of" this point, or we only have one point
187  {
188  wc = wellConstant(_permeability[_qp],
189  _rot_matrix[current_dirac_ptid],
190  _half_seg_len[current_dirac_ptid],
191  _current_elem,
192  _rs[current_dirac_ptid]);
193  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
194  // injection, so outflow<0 || // production, so outflow>0
195  outflow += _test[_i][_qp] * std::abs(character) * wc * mob * (pp - bh_pressure);
196  }
197 
199  outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
200  return outflow;
201 }
const MaterialProperty< RealTensorValue > & _permeability
material permeability
const bool _fully_upwind
Whether to use full upwinding.
const MaterialProperty< std::vector< Real > > & _density
fluid density
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
std::vector< Real > _rs
radii of the borehole
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: _...
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are used if _fully_upwind = true ...
const MaterialProperty< std::vector< Real > > & _pp
fluid porepressure (or porepressures in case of multiphase)
const MaterialProperty< std::vector< Real > > & _rel_perm
relative permeability
const RealVectorValue _unit_weight
unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point) ...
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 Real _p_bot
bottomhole pressure of borehole
const unsigned int _pvar
The moose internal variable number of the richards variable of this Dirac Kernel. ...
Point _bottom_point
the bottom point of the borehole (where bottom_pressure is defined)
std::vector< RealTensorValue > _rot_matrix
rotation matrix used in well_constant calculation
std::vector< Real > _half_seg_len
0.5*(length of polyline segments between points)
std::vector< Real > _zs
z points of borehole
Function & _character
If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure...

◆ computeResidual()

void RichardsBorehole::computeResidual ( )
virtual

Computes the residual.

This just calls prepareNodalValues, if _fully_upwind then calls DiracKernel::computeResidual

Definition at line 128 of file RichardsBorehole.C.

129 {
130  if (_fully_upwind)
132  DiracKernel::computeResidual();
133 }
const bool _fully_upwind
Whether to use full upwinding.
void prepareNodalValues()
calculates the nodal values of pressure, mobility, and derivatives thereof

◆ jac()

Real RichardsBorehole::jac ( unsigned int  wrt_num)
protected

Calculates Jacobian.

Parameters
wrt_numdifferentiate the residual wrt this Richards variable

Definition at line 230 of file RichardsBorehole.C.

Referenced by computeQpJacobian(), and computeQpOffDiagJacobian().

231 {
232  const Real character = _character.value(_t, _q_point[_qp]);
233  if (character == 0.0)
234  return 0.0;
235 
236  const Real bh_pressure =
237  _p_bot +
238  _unit_weight *
239  (_q_point[_qp] -
240  _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
241 
242  Real pp;
243  Real dpp_dv;
244  Real mob;
245  Real dmob_dv;
246  Real phi;
247  if (!_fully_upwind)
248  {
249  pp = _pp[_qp][_pvar];
250  dpp_dv = _dpp_dv[_qp][_pvar][wrt_num];
251  mob = _rel_perm[_qp][_pvar] * _density[_qp][_pvar] / _viscosity[_qp][_pvar];
252  dmob_dv = (_drel_perm_dv[_qp][_pvar][wrt_num] * _density[_qp][_pvar] +
253  _rel_perm[_qp][_pvar] * _ddensity_dv[_qp][_pvar][wrt_num]) /
254  _viscosity[_qp][_pvar];
255  phi = _phi[_j][_qp];
256  }
257  else
258  {
259  if (_i != _j)
260  return 0.0; // residual at node _i only depends on variables at that node
261  pp = (*_ps_at_nodes[_pvar])[_i];
262  dpp_dv =
263  (_pvar == wrt_num ? 1 : 0); // NOTE: i'm assuming that the variables are pressure variables
264  mob = _mobility[_i];
265  dmob_dv = _dmobility_dv[_i][wrt_num];
266  phi = 1;
267  }
268 
269  // Get the ID we initially assigned to this point
270  const unsigned current_dirac_ptid = currentPointCachedID();
271 
272  // If getting the ID failed, fall back to the old bodge!
273  // if (current_dirac_ptid == libMesh::invalid_uint)
274  // current_dirac_ptid = (_zs.size() > 2) ? 1 : 0;
275 
276  Real outflowp(0.0);
277 
278  Real wc(0.0);
279  if (current_dirac_ptid > 0)
280  // contribution from half-segment "behind" this point
281  {
282  wc = wellConstant(_permeability[_qp],
283  _rot_matrix[current_dirac_ptid - 1],
284  _half_seg_len[current_dirac_ptid - 1],
285  _current_elem,
286  _rs[current_dirac_ptid]);
287  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
288  // injection, so outflow<0 || // production, so outflow>0
289  outflowp += _test[_i][_qp] * std::abs(character) * wc *
290  (mob * phi * dpp_dv + dmob_dv * phi * (pp - bh_pressure));
291  }
292 
293  if (current_dirac_ptid < _zs.size() - 1 || _zs.size() == 1)
294  // contribution from half-segment "ahead of" this point
295  {
296  wc = wellConstant(_permeability[_qp],
297  _rot_matrix[current_dirac_ptid],
298  _half_seg_len[current_dirac_ptid],
299  _current_elem,
300  _rs[current_dirac_ptid]);
301  if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
302  // injection, so outflow<0 || // production, so outflow>0
303  outflowp += _test[_i][_qp] * std::abs(character) * wc *
304  (mob * phi * dpp_dv + dmob_dv * phi * (pp - bh_pressure));
305  }
306 
307  return outflowp;
308 }
const MaterialProperty< RealTensorValue > & _permeability
material permeability
const bool _fully_upwind
Whether to use full upwinding.
const MaterialProperty< std::vector< Real > > & _density
fluid density
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
std::vector< Real > _rs
radii of the borehole
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: _...
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are used if _fully_upwind = true ...
const MaterialProperty< std::vector< Real > > & _pp
fluid porepressure (or porepressures in case of multiphase)
const MaterialProperty< std::vector< Real > > & _rel_perm
relative permeability
const RealVectorValue _unit_weight
unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point) ...
const MaterialProperty< std::vector< std::vector< Real > > > & _drel_perm_dv
d(relperm_i)/d(variable_j)
const Real _p_bot
bottomhole pressure of borehole
const unsigned int _pvar
The moose internal variable number of the richards variable of this Dirac Kernel. ...
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)
Point _bottom_point
the bottom point of the borehole (where bottom_pressure is defined)
std::vector< RealTensorValue > _rot_matrix
rotation matrix used in well_constant calculation
std::vector< Real > _half_seg_len
0.5*(length of polyline segments between points)
std::vector< Real > _zs
z points of borehole
const MaterialProperty< std::vector< std::vector< Real > > > & _ddensity_dv
d(density_i)/d(variable_j)
Function & _character
If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure...

◆ parseNextLineReals()

bool PeacemanBorehole::parseNextLineReals ( std::ifstream &  ifs,
std::vector< Real > &  myvec 
)
protectedinherited

reads a space-separated line of floats from ifs and puts in myvec

Definition at line 158 of file PeacemanBorehole.C.

Referenced by PeacemanBorehole::PeacemanBorehole().

160 {
161  std::string line;
162  myvec.clear();
163  bool gotline(false);
164  if (getline(ifs, line))
165  {
166  gotline = true;
167 
168  // Harvest floats separated by whitespace
169  std::istringstream iss(line);
170  Real f;
171  while (iss >> f)
172  {
173  myvec.push_back(f);
174  }
175  }
176  return gotline;
177 }

◆ prepareNodalValues()

void RichardsBorehole::prepareNodalValues ( )
protected

calculates the nodal values of pressure, mobility, and derivatives thereof

Definition at line 88 of file RichardsBorehole.C.

Referenced by computeJacobian(), and computeResidual().

89 {
90  // NOTE: i'm assuming that all the richards variables are pressure values
91 
92  _num_nodes = (*_ps_at_nodes[_pvar]).size();
93 
94  Real p;
95  Real density;
96  Real ddensity_dp;
97  Real seff;
98  std::vector<Real> dseff_dp;
99  Real relperm;
100  Real drelperm_ds;
101  _mobility.resize(_num_nodes);
102  _dmobility_dv.resize(_num_nodes);
103  dseff_dp.resize(_num_p);
104  for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
105  {
106  // retrieve and calculate basic things at the node
107  p = (*_ps_at_nodes[_pvar])[nodenum]; // pressure of fluid _pvar at node nodenum
108  density = _density_UO->density(p); // density of fluid _pvar at node nodenum
109  ddensity_dp = _density_UO->ddensity(p); // d(density)/dP
110  seff = _seff_UO->seff(_ps_at_nodes,
111  nodenum); // effective saturation of fluid _pvar at node nodenum
112  _seff_UO->dseff(
113  _ps_at_nodes, nodenum, dseff_dp); // d(seff)/d(P_ph), for ph = 0, ..., _num_p - 1
114  relperm = _relperm_UO->relperm(seff); // relative permeability of fluid _pvar at node nodenum
115  drelperm_ds = _relperm_UO->drelperm(seff); // d(relperm)/dseff
116 
117  // calculate the mobility and its derivatives wrt (variable_ph = porepressure_ph)
118  _mobility[nodenum] =
119  density * relperm / _viscosity[0][_pvar]; // assume viscosity is constant throughout element
120  _dmobility_dv[nodenum].resize(_num_p);
121  for (unsigned int ph = 0; ph < _num_p; ++ph)
122  _dmobility_dv[nodenum][ph] = density * drelperm_ds * dseff_dp[ph] / _viscosity[0][_pvar];
123  _dmobility_dv[nodenum][_pvar] += ddensity_dp * relperm / _viscosity[0][_pvar];
124  }
125 }
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 RichardsDensity * _density_UO
user object defining the density. Only used if _fully_upwind = true
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 MaterialProperty< std::vector< Real > > & _viscosity
fluid viscosity
const RichardsSeff * _seff_UO
user object defining the effective saturation. Only used if _fully_upwind = true
const std::string density
Definition: NS.h:17
unsigned int _num_nodes
number of nodes in this element. Only used if _fully_upwind = true
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: _...
std::vector< Real > _mobility
nodal values of mobility = density*relperm/viscosity These are used if _fully_upwind = true ...
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 * _relperm_UO
user object defining the relative permeability. Only used if _fully_upwind = true ...
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 ...
virtual Real relperm(Real seff) const =0
relative permeability as a function of effective saturation This must be over-ridden in your derived ...
const unsigned int _pvar
The moose internal variable number of the richards variable of this Dirac Kernel. ...
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 unsigned int _num_p
number of richards variables

◆ wellConstant()

Real PeacemanBorehole::wellConstant ( const RealTensorValue &  perm,
const RealTensorValue &  rot,
const Real &  half_len,
const Elem *  ele,
const Real &  rad 
)
protectedinherited

Calculates Peaceman's form of the borehole well constant Z Chen, Y Zhang, Well flow models for various numerical methods, Int J Num Analysis and Modeling, 3 (2008) 375-388.

Definition at line 194 of file PeacemanBorehole.C.

Referenced by Q2PBorehole::computeQpResidual(), computeQpResidual(), Q2PBorehole::jac(), and jac().

200 {
201  if (_well_constant > 0)
202  return _well_constant;
203 
204  // rot_perm has its "2" component lying along the half segment
205  // we want to determine the eigenvectors of rot(0:1, 0:1), since, when
206  // rotated back to the original frame we will determine the element
207  // lengths along these directions
208  const RealTensorValue rot_perm = (rot * perm) * rot.transpose();
209  const Real trace2D = rot_perm(0, 0) + rot_perm(1, 1);
210  const Real det2D = rot_perm(0, 0) * rot_perm(1, 1) - rot_perm(0, 1) * rot_perm(1, 0);
211  const Real sq = std::sqrt(std::max(0.25 * trace2D * trace2D - det2D,
212  0.0)); // the std::max accounts for wierdo precision loss
213  const Real eig_val1 = 0.5 * trace2D + sq;
214  const Real eig_val2 = 0.5 * trace2D - sq;
215  RealVectorValue eig_vec1, eig_vec2;
216  if (sq > std::abs(trace2D) * 1E-7) // matrix is not a multiple of the identity (1E-7 accounts for
217  // precision in a crude way)
218  {
219  if (rot_perm(1, 0) != 0)
220  {
221  eig_vec1(0) = eig_val1 - rot_perm(1, 1);
222  eig_vec1(1) = rot_perm(1, 0);
223  eig_vec2(0) = eig_val2 - rot_perm(1, 1);
224  eig_vec2(1) = rot_perm(1, 0);
225  }
226  else if (rot_perm(0, 1) != 0)
227  {
228  eig_vec1(0) = rot_perm(0, 1);
229  eig_vec1(1) = eig_val1 - rot_perm(0, 0);
230  eig_vec2(0) = rot_perm(0, 1);
231  eig_vec2(1) = eig_val2 - rot_perm(0, 0);
232  }
233  else // off diagonal terms are both zero
234  {
235  eig_vec1(0) = 1;
236  eig_vec2(1) = 1;
237  }
238  }
239  else // matrix is basically a multiple of the identity
240  {
241  eig_vec1(0) = 1;
242  eig_vec2(1) = 1;
243  }
244 
245  // finally, rotate these to original frame and normalise
246  eig_vec1 = rot.transpose() * eig_vec1;
247  eig_vec1 /= std::sqrt(eig_vec1 * eig_vec1);
248  eig_vec2 = rot.transpose() * eig_vec2;
249  eig_vec2 /= std::sqrt(eig_vec2 * eig_vec2);
250 
251  // find the "length" of the element in these directions
252  // TODO - probably better to use variance than max&min
253  Real max1 = eig_vec1 * ele->point(0);
254  Real max2 = eig_vec2 * ele->point(0);
255  Real min1 = max1;
256  Real min2 = max2;
257  Real proj;
258  for (unsigned int i = 1; i < ele->n_nodes(); i++)
259  {
260  proj = eig_vec1 * ele->point(i);
261  max1 = (max1 < proj) ? proj : max1;
262  min1 = (min1 < proj) ? min1 : proj;
263 
264  proj = eig_vec2 * ele->point(i);
265  max2 = (max2 < proj) ? proj : max2;
266  min2 = (min2 < proj) ? min2 : proj;
267  }
268  const Real ll1 = max1 - min1;
269  const Real ll2 = max2 - min2;
270 
271  const Real r0 = _re_constant * std::sqrt(std::sqrt(eig_val1 / eig_val2) * std::pow(ll2, 2) +
272  std::sqrt(eig_val2 / eig_val1) * std::pow(ll1, 2)) /
273  (std::pow(eig_val1 / eig_val2, 0.25) + std::pow(eig_val2 / eig_val1, 0.25));
274 
275  const Real effective_perm = std::sqrt(det2D);
276 
277  const Real halfPi = acos(0.0);
278 
279  if (r0 <= rad)
280  mooseError("The effective element size (about 0.2-times-true-ele-size) for an element "
281  "containing a Peaceman-type borehole must be (much) larger than the borehole radius "
282  "for the Peaceman formulation to be correct. Your element has effective size ",
283  r0,
284  " and the borehole radius is ",
285  rad,
286  "\n");
287 
288  return 4 * halfPi * effective_perm * half_len / std::log(r0 / rad);
289 }
const Real _well_constant
well constant
const Real _re_constant
borehole constant
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)

Member Data Documentation

◆ _bottom_point

Point PeacemanBorehole::_bottom_point
protectedinherited

the bottom point of the borehole (where bottom_pressure is defined)

Definition at line 95 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), computeQpResidual(), Q2PBorehole::jac(), jac(), and PeacemanBorehole::PeacemanBorehole().

◆ _character

Function& PeacemanBorehole::_character
protectedinherited

If positive then the borehole acts as a sink (producion well) for porepressure > borehole pressure, and does nothing otherwise If negative then the borehole acts as a source (injection well) for porepressure < borehole pressure, and does nothing otherwise The flow rate to/from the borehole is multiplied by |character|, so usually character = +/- 1.

Definition at line 67 of file PeacemanBorehole.h.

Referenced by computeQpJacobian(), Q2PBorehole::computeQpResidual(), computeQpResidual(), Q2PBorehole::jac(), and jac().

◆ _ddensity_dv

const MaterialProperty<std::vector<std::vector<Real> > >& RichardsBorehole::_ddensity_dv
protected

d(density_i)/d(variable_j)

Definition at line 136 of file RichardsBorehole.h.

Referenced by jac().

◆ _density

const MaterialProperty<std::vector<Real> >& RichardsBorehole::_density
protected

fluid density

Definition at line 133 of file RichardsBorehole.h.

Referenced by computeQpResidual(), and jac().

◆ _density_UO

const RichardsDensity* RichardsBorehole::_density_UO
protected

user object defining the density. Only used if _fully_upwind = true

Definition at line 88 of file RichardsBorehole.h.

Referenced by prepareNodalValues().

◆ _dmobility_dv

std::vector<std::vector<Real> > RichardsBorehole::_dmobility_dv
protected

d(_mobility)/d(variable_ph) (variable_ph is the variable for phase=ph) These are used in the jacobian calculations if _fully_upwind = true

Definition at line 109 of file RichardsBorehole.h.

Referenced by jac(), and prepareNodalValues().

◆ _dpp_dv

const MaterialProperty<std::vector<std::vector<Real> > >& RichardsBorehole::_dpp_dv
protected

d(porepressure_i)/d(variable_j)

Definition at line 115 of file RichardsBorehole.h.

Referenced by jac().

◆ _drel_perm_dv

const MaterialProperty<std::vector<std::vector<Real> > >& RichardsBorehole::_drel_perm_dv
protected

d(relperm_i)/d(variable_j)

Definition at line 130 of file RichardsBorehole.h.

Referenced by jac().

◆ _dseff_dv

const MaterialProperty<std::vector<std::vector<Real> > >& RichardsBorehole::_dseff_dv
protected

deriviatves of Seff wrt variables

Definition at line 124 of file RichardsBorehole.h.

◆ _fully_upwind

const bool RichardsBorehole::_fully_upwind
protected

Whether to use full upwinding.

Definition at line 76 of file RichardsBorehole.h.

Referenced by computeJacobian(), computeQpResidual(), computeResidual(), and jac().

◆ _half_seg_len

std::vector<Real> PeacemanBorehole::_half_seg_len
protectedinherited

0.5*(length of polyline segments between points)

Definition at line 98 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), computeQpResidual(), Q2PBorehole::jac(), jac(), and PeacemanBorehole::PeacemanBorehole().

◆ _mobility

std::vector<Real> RichardsBorehole::_mobility
protected

nodal values of mobility = density*relperm/viscosity These are used if _fully_upwind = true

Definition at line 103 of file RichardsBorehole.h.

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

◆ _num_nodes

unsigned int RichardsBorehole::_num_nodes
protected

number of nodes in this element. Only used if _fully_upwind = true

Definition at line 97 of file RichardsBorehole.h.

Referenced by prepareNodalValues().

◆ _num_p

const unsigned int RichardsBorehole::_num_p
protected

number of richards variables

Definition at line 82 of file RichardsBorehole.h.

Referenced by prepareNodalValues(), and RichardsBorehole().

◆ _p_bot

const Real PeacemanBorehole::_p_bot
protectedinherited

bottomhole pressure of borehole

Definition at line 70 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), computeQpResidual(), Q2PBorehole::jac(), and jac().

◆ _permeability

const MaterialProperty<RealTensorValue>& RichardsBorehole::_permeability
protected

material permeability

Definition at line 121 of file RichardsBorehole.h.

Referenced by computeQpResidual(), and jac().

◆ _pp

const MaterialProperty<std::vector<Real> >& RichardsBorehole::_pp
protected

fluid porepressure (or porepressures in case of multiphase)

Definition at line 112 of file RichardsBorehole.h.

Referenced by computeQpResidual(), and jac().

◆ _ps_at_nodes

std::vector<const VariableValue *> RichardsBorehole::_ps_at_nodes
protected

Holds the values of pressures at all the nodes of the element Only used if _fully_upwind = true Eg: _ps_at_nodes[_pvar] is a pointer to this variable's nodal porepressure values So: (*_ps_at_nodes[_pvar])[i] = _var.dofValues()[i] = porepressure of pressure-variable _pvar at node i.

Definition at line 146 of file RichardsBorehole.h.

Referenced by computeQpResidual(), jac(), prepareNodalValues(), and RichardsBorehole().

◆ _pvar

const unsigned int RichardsBorehole::_pvar
protected

The moose internal variable number of the richards variable of this Dirac Kernel.

Definition at line 85 of file RichardsBorehole.h.

Referenced by computeQpJacobian(), computeQpResidual(), jac(), and prepareNodalValues().

◆ _rel_perm

const MaterialProperty<std::vector<Real> >& RichardsBorehole::_rel_perm
protected

relative permeability

Definition at line 127 of file RichardsBorehole.h.

Referenced by computeQpResidual(), and jac().

◆ _relperm_UO

const RichardsRelPerm* RichardsBorehole::_relperm_UO
protected

user object defining the relative permeability. Only used if _fully_upwind = true

Definition at line 94 of file RichardsBorehole.h.

Referenced by prepareNodalValues().

◆ _richards_name_UO

const RichardsVarNames& RichardsBorehole::_richards_name_UO
protected

Defines the richards variables in the simulation.

Definition at line 79 of file RichardsBorehole.h.

Referenced by computeQpOffDiagJacobian(), and RichardsBorehole().

◆ _rot_matrix

std::vector<RealTensorValue> PeacemanBorehole::_rot_matrix
protectedinherited

rotation matrix used in well_constant calculation

Definition at line 101 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), computeQpResidual(), Q2PBorehole::jac(), jac(), and PeacemanBorehole::PeacemanBorehole().

◆ _rs

std::vector<Real> PeacemanBorehole::_rs
protectedinherited

◆ _seff_UO

const RichardsSeff* RichardsBorehole::_seff_UO
protected

user object defining the effective saturation. Only used if _fully_upwind = true

Definition at line 91 of file RichardsBorehole.h.

Referenced by prepareNodalValues().

◆ _total_outflow_mass

RichardsSumQuantity& PeacemanBorehole::_total_outflow_mass
protectedinherited

This is used to hold the total fluid flowing into the borehole Hence, it is positive for production wells where fluid is flowing from porespace into the borehole and removed from the model.

Definition at line 80 of file PeacemanBorehole.h.

Referenced by PeacemanBorehole::addPoints(), Q2PBorehole::computeQpResidual(), computeQpResidual(), and PeacemanBorehole::PeacemanBorehole().

◆ _unit_weight

const RealVectorValue PeacemanBorehole::_unit_weight
protectedinherited

unit weight of fluid in borehole (for calculating bottomhole pressure at each Dirac Point)

Definition at line 73 of file PeacemanBorehole.h.

Referenced by Q2PBorehole::computeQpResidual(), computeQpResidual(), Q2PBorehole::jac(), and jac().

◆ _viscosity

const MaterialProperty<std::vector<Real> >& RichardsBorehole::_viscosity
protected

fluid viscosity

Definition at line 118 of file RichardsBorehole.h.

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

◆ _xs

std::vector<Real> PeacemanBorehole::_xs
protectedinherited

x points of the borehole

Definition at line 86 of file PeacemanBorehole.h.

Referenced by PeacemanBorehole::addPoints(), and PeacemanBorehole::PeacemanBorehole().

◆ _ys

std::vector<Real> PeacemanBorehole::_ys
protectedinherited

y points of the borehole

Definition at line 89 of file PeacemanBorehole.h.

Referenced by PeacemanBorehole::addPoints(), and PeacemanBorehole::PeacemanBorehole().

◆ _zs

std::vector<Real> PeacemanBorehole::_zs
protectedinherited

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