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 "RichardsPiecewiseLinearSink.h"
11 :
12 : // MOOSE includes
13 : #include "MooseVariable.h"
14 :
15 : // C++ includes
16 : #include <iostream>
17 :
18 : registerMooseObject("RichardsApp", RichardsPiecewiseLinearSink);
19 :
20 : InputParameters
21 110 : RichardsPiecewiseLinearSink::validParams()
22 : {
23 110 : InputParameters params = IntegratedBC::validParams();
24 220 : 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 220 : 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 220 : params.addParam<std::vector<UserObjectName>>("relperm_UO",
35 : "List of names of user objects that "
36 : "define relative permeability. Only "
37 : "needed if fully_upwind is used");
38 220 : params.addParam<std::vector<UserObjectName>>(
39 : "seff_UO",
40 : "List of name of user objects that define effective saturation as a function of "
41 : "pressure list. Only needed if fully_upwind is used");
42 220 : params.addParam<std::vector<UserObjectName>>("density_UO",
43 : "List of names of user objects that "
44 : "define the fluid density. Only "
45 : "needed if fully_upwind is used");
46 220 : params.addRequiredParam<std::vector<Real>>(
47 : "pressures", "Tuple of pressure values. Must be monotonically increasing.");
48 220 : params.addRequiredParam<std::vector<Real>>(
49 : "bare_fluxes",
50 : "Tuple of flux values (measured in kg.m^-2.s^-1 for use_mobility=false, and "
51 : "in Pa.s^-1 if use_mobility=true). This flux is OUT of the medium: hence "
52 : "positive values of flux means this will be a SINK, while negative values "
53 : "indicate this flux will be a SOURCE. A piecewise-linear fit is performed to "
54 : "the (pressure,bare_fluxes) pairs to obtain the flux at any arbitrary "
55 : "pressure, and the first or last bare_flux values are used if the quad-point "
56 : "pressure falls outside this range.");
57 220 : params.addParam<FunctionName>("multiplying_fcn",
58 220 : 1.0,
59 : "If this function is provided, the flux "
60 : "will be multiplied by this function. "
61 : "This is useful for spatially or "
62 : "temporally varying sinks");
63 220 : params.addRequiredParam<UserObjectName>(
64 : "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
65 220 : params.addParam<bool>("fully_upwind", false, "Use full upwinding");
66 220 : params.addParam<PostprocessorName>(
67 : "area_pp",
68 220 : 1,
69 : "An area postprocessor. If given, the bare_fluxes will be divided by this "
70 : "quantity. This means the bare fluxes are measured in kg.s^-1. This is "
71 : "useful for the case when you wish to provide the *total* flux, and let MOOSE "
72 : "proportion it uniformly across the boundary. In that case you would have "
73 : "use_mobility=false=use_relperm, and only one bare flux should be specified");
74 110 : return params;
75 0 : }
76 :
77 52 : RichardsPiecewiseLinearSink::RichardsPiecewiseLinearSink(const InputParameters & parameters)
78 : : IntegratedBC(parameters),
79 52 : _use_mobility(getParam<bool>("use_mobility")),
80 104 : _use_relperm(getParam<bool>("use_relperm")),
81 104 : _fully_upwind(getParam<bool>("fully_upwind")),
82 :
83 208 : _sink_func(getParam<std::vector<Real>>("pressures"),
84 : getParam<std::vector<Real>>("bare_fluxes")),
85 :
86 52 : _m_func(getFunction("multiplying_fcn")),
87 :
88 52 : _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
89 52 : _num_p(_richards_name_UO.num_v()),
90 52 : _pvar(_richards_name_UO.richards_var_num(_var.number())),
91 :
92 : // in the following, getUserObjectByName returns a reference (an alias) to a RichardsBLAH user
93 : // object, and the & turns it into a pointer
94 75 : _density_UO(_fully_upwind ? &getUserObjectByName<RichardsDensity>(
95 75 : getParam<std::vector<UserObjectName>>("density_UO")[_pvar])
96 : : NULL),
97 75 : _seff_UO(_fully_upwind ? &getUserObjectByName<RichardsSeff>(
98 75 : getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])
99 : : NULL),
100 75 : _relperm_UO(_fully_upwind ? &getUserObjectByName<RichardsRelPerm>(
101 75 : getParam<std::vector<UserObjectName>>("relperm_UO")[_pvar])
102 : : NULL),
103 :
104 52 : _area_pp(getPostprocessorValue("area_pp")),
105 :
106 52 : _num_nodes(0),
107 52 : _nodal_density(0),
108 52 : _dnodal_density_dv(0),
109 52 : _nodal_relperm(0),
110 52 : _dnodal_relperm_dv(0),
111 :
112 104 : _pp(getMaterialProperty<std::vector<Real>>("porepressure")),
113 104 : _dpp_dv(getMaterialProperty<std::vector<std::vector<Real>>>("dporepressure_dv")),
114 :
115 104 : _viscosity(getMaterialProperty<std::vector<Real>>("viscosity")),
116 104 : _permeability(getMaterialProperty<RealTensorValue>("permeability")),
117 :
118 104 : _dseff_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ds_eff_dv")),
119 :
120 104 : _rel_perm(getMaterialProperty<std::vector<Real>>("rel_perm")),
121 104 : _drel_perm_dv(getMaterialProperty<std::vector<std::vector<Real>>>("drel_perm_dv")),
122 :
123 104 : _density(getMaterialProperty<std::vector<Real>>("density")),
124 156 : _ddensity_dv(getMaterialProperty<std::vector<std::vector<Real>>>("ddensity_dv"))
125 : {
126 52 : _ps_at_nodes.resize(_num_p);
127 110 : for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
128 58 : _ps_at_nodes[pnum] = _richards_name_UO.nodal_var(pnum);
129 52 : }
130 :
131 : void
132 16110 : RichardsPiecewiseLinearSink::prepareNodalValues()
133 : {
134 : // NOTE: i'm assuming that all the richards variables are pressure values
135 :
136 16110 : _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 16110 : _nodal_density.resize(_num_nodes);
144 16110 : _dnodal_density_dv.resize(_num_nodes);
145 16110 : _nodal_relperm.resize(_num_nodes);
146 16110 : _dnodal_relperm_dv.resize(_num_nodes);
147 16110 : dseff_dp.resize(_num_p);
148 80958 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
149 : {
150 : // retrieve and calculate basic things at the node
151 64848 : p = (*_ps_at_nodes[_pvar])[nodenum]; // pressure of fluid _pvar at node nodenum
152 :
153 64848 : _nodal_density[nodenum] = _density_UO->density(p); // density of fluid _pvar at node nodenum
154 64848 : _dnodal_density_dv[nodenum].resize(_num_p);
155 130064 : for (unsigned int ph = 0; ph < _num_p; ++ph)
156 65216 : _dnodal_density_dv[nodenum][ph] = 0;
157 64848 : _dnodal_density_dv[nodenum][_pvar] = _density_UO->ddensity(p); // d(density)/dP
158 :
159 64848 : seff = _seff_UO->seff(_ps_at_nodes,
160 : nodenum); // effective saturation of fluid _pvar at node nodenum
161 64848 : _seff_UO->dseff(
162 : _ps_at_nodes, nodenum, dseff_dp); // d(seff)/d(P_ph), for ph = 0, ..., _num_p - 1
163 :
164 64848 : _nodal_relperm[nodenum] =
165 64848 : _relperm_UO->relperm(seff); // relative permeability of fluid _pvar at node nodenum
166 64848 : drelperm_ds = _relperm_UO->drelperm(seff); // d(relperm)/dseff
167 :
168 64848 : _dnodal_relperm_dv[nodenum].resize(_num_p);
169 130064 : for (unsigned int ph = 0; ph < _num_p; ++ph)
170 65216 : _dnodal_relperm_dv[nodenum][ph] = drelperm_ds * dseff_dp[ph];
171 : }
172 16110 : }
173 :
174 : void
175 15350 : RichardsPiecewiseLinearSink::computeResidual()
176 : {
177 15350 : if (_fully_upwind)
178 8430 : prepareNodalValues();
179 15350 : IntegratedBC::computeResidual();
180 15349 : }
181 :
182 : Real
183 127017 : RichardsPiecewiseLinearSink::computeQpResidual()
184 : {
185 : Real flux = 0;
186 : Real k = 0;
187 :
188 127017 : if (!_fully_upwind)
189 : {
190 57472 : flux = _test[_i][_qp] * _sink_func.sample(_pp[_qp][_pvar]);
191 57472 : if (_use_mobility)
192 : {
193 15488 : k = (_permeability[_qp] * _normals[_qp]) * _normals[_qp];
194 15488 : flux *= _density[_qp][_pvar] * k / _viscosity[_qp][_pvar];
195 : }
196 57472 : if (_use_relperm)
197 15488 : flux *= _rel_perm[_qp][_pvar];
198 : }
199 : else
200 : {
201 69545 : flux = _test[_i][_qp] * _sink_func.sample((*_ps_at_nodes[_pvar])[_i]);
202 69545 : if (_use_mobility)
203 : {
204 15488 : k = (_permeability[0] * _normals[_qp]) * _normals[_qp]; // assume that _permeability is
205 : // constant throughout element so
206 : // doesn't need to be upwinded
207 15488 : flux *= _nodal_density[_i] * k /
208 15488 : _viscosity[0][_pvar]; // assume that viscosity is constant throughout element
209 : }
210 69545 : if (_use_relperm)
211 15488 : flux *= _nodal_relperm[_i];
212 : }
213 :
214 127017 : flux *= _m_func.value(_t, _q_point[_qp]);
215 :
216 127017 : if (_area_pp == 0.0)
217 : {
218 25 : if (flux != 0)
219 1 : mooseError("RichardsPiecewiseLinearSink: flux is nonzero, but area is zero!\n");
220 : // if flux == 0, then leave it as zero.
221 : }
222 : else
223 126992 : flux /= _area_pp;
224 :
225 127016 : return flux;
226 : }
227 :
228 : void
229 9620 : RichardsPiecewiseLinearSink::computeJacobian()
230 : {
231 9620 : if (_fully_upwind)
232 4466 : prepareNodalValues();
233 9620 : IntegratedBC::computeJacobian();
234 9619 : }
235 :
236 : Real
237 310497 : RichardsPiecewiseLinearSink::computeQpJacobian()
238 : {
239 310497 : return jac(_pvar);
240 : }
241 :
242 : void
243 8624 : RichardsPiecewiseLinearSink::computeOffDiagJacobian(const unsigned int jvar)
244 : {
245 8624 : if (_fully_upwind)
246 3214 : prepareNodalValues();
247 8624 : IntegratedBC::computeOffDiagJacobian(jvar);
248 8623 : }
249 :
250 : Real
251 9152 : RichardsPiecewiseLinearSink::computeQpOffDiagJacobian(unsigned int jvar)
252 : {
253 9152 : if (_richards_name_UO.not_richards_var(jvar))
254 : return 0.0;
255 9152 : unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
256 9152 : return jac(dvar);
257 : }
258 :
259 : Real
260 319649 : 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 319649 : if (!_fully_upwind)
270 : {
271 174912 : flux = _sink_func.sample(_pp[_qp][_pvar]);
272 174912 : deriv = _sink_func.sampleDerivative(_pp[_qp][_pvar]) * _dpp_dv[_qp][_pvar][wrt_num];
273 174912 : phi = _phi[_j][_qp];
274 174912 : if (_use_mobility)
275 : {
276 42432 : k = (_permeability[_qp] * _normals[_qp]) * _normals[_qp];
277 42432 : mob = _density[_qp][_pvar] * k / _viscosity[_qp][_pvar];
278 42432 : mobp = _ddensity_dv[_qp][_pvar][wrt_num] * k / _viscosity[_qp][_pvar];
279 42432 : deriv = mob * deriv + mobp * flux;
280 42432 : flux *= mob;
281 : }
282 174912 : if (_use_relperm)
283 42432 : deriv = _rel_perm[_qp][_pvar] * deriv + _drel_perm_dv[_qp][_pvar][wrt_num] * flux;
284 : }
285 : else
286 : {
287 144737 : if (_i != _j)
288 : return 0.0; // residual at node _i only depends on variables at that node
289 35929 : flux = _sink_func.sample((*_ps_at_nodes[_pvar])[_i]);
290 35929 : 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 35929 : if (_use_mobility)
294 : {
295 10416 : k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
296 10416 : mob = _nodal_density[_i] * k / _viscosity[0][_pvar];
297 10416 : mobp = _dnodal_density_dv[_i][wrt_num] * k / _viscosity[0][_pvar];
298 10416 : deriv = mob * deriv + mobp * flux;
299 10416 : flux *= mob;
300 : }
301 35929 : if (_use_relperm)
302 10416 : deriv = _nodal_relperm[_i] * deriv + _dnodal_relperm_dv[_i][wrt_num] * flux;
303 : }
304 :
305 210841 : deriv *= _m_func.value(_t, _q_point[_qp]);
306 :
307 210841 : if (_area_pp == 0.0)
308 : {
309 1 : if (deriv != 0)
310 1 : mooseError("RichardsPiecewiseLinearSink: deriv is nonzero, but area is zero!\n");
311 : // if deriv == 0, then leave it as zero.
312 : }
313 : else
314 210840 : deriv /= _area_pp;
315 :
316 210840 : return _test[_i][_qp] * deriv * phi;
317 : }
|
