Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 :
13 : registerMooseObject("RichardsApp", RichardsBorehole);
14 :
15 : InputParameters
16 84 : RichardsBorehole::validParams()
17 : {
18 84 : InputParameters params = PeacemanBorehole::validParams();
19 168 : params.addRequiredParam<UserObjectName>(
20 : "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
21 168 : 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 168 : 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 168 : 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 168 : params.addParam<bool>("fully_upwind", false, "Fully upwind the flux");
34 84 : 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 84 : return params;
38 0 : }
39 :
40 43 : RichardsBorehole::RichardsBorehole(const InputParameters & parameters)
41 : : PeacemanBorehole(parameters),
42 41 : _fully_upwind(getParam<bool>("fully_upwind")),
43 41 : _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
44 41 : _num_p(_richards_name_UO.num_v()),
45 41 : _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 56 : _density_UO(_fully_upwind ? &getUserObjectByName<RichardsDensity>(
50 56 : getParam<std::vector<UserObjectName>>("density_UO")[_pvar])
51 : : NULL),
52 56 : _seff_UO(_fully_upwind ? &getUserObjectByName<RichardsSeff>(
53 56 : getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])
54 : : NULL),
55 56 : _relperm_UO(_fully_upwind ? &getUserObjectByName<RichardsRelPerm>(
56 56 : getParam<std::vector<UserObjectName>>("relperm_UO")[_pvar])
57 : : NULL),
58 :
59 41 : _num_nodes(0),
60 41 : _mobility(0),
61 41 : _dmobility_dv(0),
62 82 : _pp(getMaterialProperty<std::vector<Real>>("porepressure")),
63 82 : _dpp_dv(getMaterialProperty<std::vector<std::vector<Real>>>("dporepressure_dv")),
64 82 : _viscosity(getMaterialProperty<std::vector<Real>>("viscosity")),
65 82 : _permeability(getMaterialProperty<RealTensorValue>("permeability")),
66 82 : _dseff_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ds_eff_dv")),
67 82 : _rel_perm(getMaterialProperty<std::vector<Real>>("rel_perm")),
68 82 : _drel_perm_dv(getMaterialProperty<std::vector<std::vector<Real>>>("drel_perm_dv")),
69 82 : _density(getMaterialProperty<std::vector<Real>>("density")),
70 125 : _ddensity_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ddensity_dv"))
71 : {
72 41 : _ps_at_nodes.resize(_num_p);
73 86 : for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
74 45 : _ps_at_nodes[pnum] = _richards_name_UO.nodal_var(pnum);
75 :
76 : // To correctly compute the Jacobian terms,
77 : // tell MOOSE that this DiracKernel depends on all the Richards Vars
78 41 : const std::vector<MooseVariableFEBase *> & coupled_vars = _richards_name_UO.getCoupledMooseVars();
79 86 : for (unsigned int i = 0; i < coupled_vars.size(); i++)
80 90 : addMooseVariableDependency(coupled_vars[i]);
81 41 : }
82 :
83 : void
84 2399 : RichardsBorehole::prepareNodalValues()
85 : {
86 : // NOTE: i'm assuming that all the richards variables are pressure values
87 :
88 2399 : _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 2399 : _mobility.resize(_num_nodes);
98 2399 : _dmobility_dv.resize(_num_nodes);
99 2399 : dseff_dp.resize(_num_p);
100 21591 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
101 : {
102 : // retrieve and calculate basic things at the node
103 19192 : p = (*_ps_at_nodes[_pvar])[nodenum]; // pressure of fluid _pvar at node nodenum
104 19192 : density = _density_UO->density(p); // density of fluid _pvar at node nodenum
105 19192 : ddensity_dp = _density_UO->ddensity(p); // d(density)/dP
106 19192 : seff = _seff_UO->seff(_ps_at_nodes,
107 : nodenum); // effective saturation of fluid _pvar at node nodenum
108 19192 : _seff_UO->dseff(
109 : _ps_at_nodes, nodenum, dseff_dp); // d(seff)/d(P_ph), for ph = 0, ..., _num_p - 1
110 19192 : relperm = _relperm_UO->relperm(seff); // relative permeability of fluid _pvar at node nodenum
111 19192 : drelperm_ds = _relperm_UO->drelperm(seff); // d(relperm)/dseff
112 :
113 : // calculate the mobility and its derivatives wrt (variable_ph = porepressure_ph)
114 19192 : _mobility[nodenum] =
115 19192 : density * relperm / _viscosity[0][_pvar]; // assume viscosity is constant throughout element
116 19192 : _dmobility_dv[nodenum].resize(_num_p);
117 42752 : for (unsigned int ph = 0; ph < _num_p; ++ph)
118 23560 : _dmobility_dv[nodenum][ph] = density * drelperm_ds * dseff_dp[ph] / _viscosity[0][_pvar];
119 19192 : _dmobility_dv[nodenum][_pvar] += ddensity_dp * relperm / _viscosity[0][_pvar];
120 : }
121 2399 : }
122 :
123 : void
124 4714 : RichardsBorehole::computeResidual()
125 : {
126 4714 : if (_fully_upwind)
127 1603 : prepareNodalValues();
128 4714 : DiracKernel::computeResidual();
129 4713 : }
130 :
131 : Real
132 75249 : RichardsBorehole::computeQpResidual()
133 : {
134 75249 : const Real character = _character.value(_t, _q_point[_qp]);
135 75249 : if (character == 0.0)
136 : return 0.0;
137 :
138 : const Real bh_pressure =
139 75249 : _p_bot +
140 : _unit_weight *
141 75249 : (_q_point[_qp] -
142 75249 : _bottom_point); // really want to use _q_point instaed of _current_point, i think?!
143 :
144 : Real pp;
145 : Real mob;
146 75249 : if (!_fully_upwind)
147 : {
148 49761 : pp = _pp[_qp][_pvar];
149 49761 : mob = _rel_perm[_qp][_pvar] * _density[_qp][_pvar] / _viscosity[_qp][_pvar];
150 : }
151 : else
152 : {
153 25488 : pp = (*_ps_at_nodes[_pvar])[_i];
154 25488 : mob = _mobility[_i];
155 : }
156 :
157 : // Get the ID we initially assigned to this point
158 75249 : 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 75249 : if (current_dirac_ptid > 0)
168 : // contribution from half-segment "behind" this point (must have >1 point for
169 : // current_dirac_ptid>0)
170 : {
171 37544 : wc = wellConstant(_permeability[_qp],
172 : _rot_matrix[current_dirac_ptid - 1],
173 37544 : _half_seg_len[current_dirac_ptid - 1],
174 37544 : _current_elem,
175 37544 : _rs[current_dirac_ptid]);
176 37544 : if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
177 : // injection, so outflow<0 || // production, so outflow>0
178 37114 : outflow += _test[_i][_qp] * std::abs(character) * wc * mob * (pp - bh_pressure);
179 : }
180 :
181 75249 : 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 : {
184 37705 : wc = wellConstant(_permeability[_qp],
185 : _rot_matrix[current_dirac_ptid],
186 : _half_seg_len[current_dirac_ptid],
187 37705 : _current_elem,
188 37705 : _rs[current_dirac_ptid]);
189 37704 : if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
190 : // injection, so outflow<0 || // production, so outflow>0
191 37274 : outflow += _test[_i][_qp] * std::abs(character) * wc * mob * (pp - bh_pressure);
192 : }
193 :
194 75248 : _total_outflow_mass.add(
195 75248 : outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
196 75248 : return outflow;
197 : }
198 :
199 : void
200 2594 : RichardsBorehole::computeJacobian()
201 : {
202 2594 : if (_fully_upwind)
203 796 : prepareNodalValues();
204 2594 : DiracKernel::computeJacobian();
205 2594 : }
206 :
207 : Real
208 331072 : RichardsBorehole::computeQpJacobian()
209 : {
210 331072 : const Real character = _character.value(_t, _q_point[_qp]);
211 331072 : if (character == 0.0)
212 : return 0.0;
213 331072 : return jac(_pvar);
214 : }
215 :
216 : Real
217 6656 : RichardsBorehole::computeQpOffDiagJacobian(unsigned int jvar)
218 : {
219 6656 : if (_richards_name_UO.not_richards_var(jvar))
220 : return 0.0;
221 6656 : const unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
222 6656 : return jac(dvar);
223 : }
224 :
225 : Real
226 337728 : RichardsBorehole::jac(unsigned int wrt_num)
227 : {
228 337728 : const Real character = _character.value(_t, _q_point[_qp]);
229 337728 : if (character == 0.0)
230 : return 0.0;
231 :
232 : const Real bh_pressure =
233 337728 : _p_bot +
234 : _unit_weight *
235 337728 : (_q_point[_qp] -
236 337728 : _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 337728 : if (!_fully_upwind)
244 : {
245 233472 : pp = _pp[_qp][_pvar];
246 233472 : dpp_dv = _dpp_dv[_qp][_pvar][wrt_num];
247 233472 : mob = _rel_perm[_qp][_pvar] * _density[_qp][_pvar] / _viscosity[_qp][_pvar];
248 233472 : dmob_dv = (_drel_perm_dv[_qp][_pvar][wrt_num] * _density[_qp][_pvar] +
249 233472 : _rel_perm[_qp][_pvar] * _ddensity_dv[_qp][_pvar][wrt_num]) /
250 : _viscosity[_qp][_pvar];
251 233472 : phi = _phi[_j][_qp];
252 : }
253 : else
254 : {
255 104256 : if (_i != _j)
256 : return 0.0; // residual at node _i only depends on variables at that node
257 13032 : pp = (*_ps_at_nodes[_pvar])[_i];
258 13032 : dpp_dv =
259 : (_pvar == wrt_num ? 1 : 0); // NOTE: i'm assuming that the variables are pressure variables
260 13032 : mob = _mobility[_i];
261 13032 : dmob_dv = _dmobility_dv[_i][wrt_num];
262 : phi = 1;
263 : }
264 :
265 : // Get the ID we initially assigned to this point
266 246504 : 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 246504 : if (current_dirac_ptid > 0)
276 : // contribution from half-segment "behind" this point
277 : {
278 123192 : wc = wellConstant(_permeability[_qp],
279 : _rot_matrix[current_dirac_ptid - 1],
280 123192 : _half_seg_len[current_dirac_ptid - 1],
281 123192 : _current_elem,
282 123192 : _rs[current_dirac_ptid]);
283 123192 : if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
284 : // injection, so outflow<0 || // production, so outflow>0
285 123192 : outflowp += _test[_i][_qp] * std::abs(character) * wc *
286 123192 : (mob * phi * dpp_dv + dmob_dv * phi * (pp - bh_pressure));
287 : }
288 :
289 246504 : if (current_dirac_ptid < _zs.size() - 1 || _zs.size() == 1)
290 : // contribution from half-segment "ahead of" this point
291 : {
292 123312 : wc = wellConstant(_permeability[_qp],
293 : _rot_matrix[current_dirac_ptid],
294 : _half_seg_len[current_dirac_ptid],
295 123312 : _current_elem,
296 : _rs[current_dirac_ptid]);
297 123312 : if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
298 : // injection, so outflow<0 || // production, so outflow>0
299 123312 : outflowp += _test[_i][_qp] * std::abs(character) * wc *
300 123312 : (mob * phi * dpp_dv + dmob_dv * phi * (pp - bh_pressure));
301 : }
302 :
303 : return outflowp;
304 : }
|
