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 "Q2PPiecewiseLinearSink.h"
11 :
12 : // MOOSE includes
13 : #include "MooseVariable.h"
14 :
15 : // C++ includes
16 : #include <iostream>
17 :
18 : registerMooseObject("RichardsApp", Q2PPiecewiseLinearSink);
19 :
20 : InputParameters
21 55 : Q2PPiecewiseLinearSink::validParams()
22 : {
23 55 : InputParameters params = IntegratedBC::validParams();
24 110 : 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 110 : 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 110 : params.addRequiredParam<std::vector<Real>>(
35 : "pressures", "Tuple of pressure values. Must be monotonically increasing.");
36 110 : 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 110 : params.addParam<FunctionName>("multiplying_fcn",
46 110 : 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 110 : params.addRequiredParam<UserObjectName>(
52 : "fluid_density",
53 : "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
54 110 : 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 110 : 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 110 : 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 110 : params.addRequiredParam<Real>("fluid_viscosity", "The fluid dynamic viscosity");
69 55 : params.addClassDescription("Sink of fluid, controlled by (pressure, bare_fluxes) interpolation. "
70 : "This is for use in Q2P models");
71 55 : return params;
72 0 : }
73 :
74 20 : Q2PPiecewiseLinearSink::Q2PPiecewiseLinearSink(const InputParameters & parameters)
75 : : IntegratedBC(parameters),
76 20 : _use_mobility(getParam<bool>("use_mobility")),
77 40 : _use_relperm(getParam<bool>("use_relperm")),
78 80 : _sink_func(getParam<std::vector<Real>>("pressures"),
79 : getParam<std::vector<Real>>("bare_fluxes")),
80 20 : _m_func(getFunction("multiplying_fcn")),
81 20 : _density(getUserObject<RichardsDensity>("fluid_density")),
82 20 : _relperm(getUserObject<RichardsRelPerm>("fluid_relperm")),
83 20 : _other_var_nodal(coupledDofValues("other_var")),
84 20 : _other_var_num(coupled("other_var")),
85 40 : _var_is_pp(getParam<bool>("var_is_porepressure")),
86 40 : _viscosity(getParam<Real>("fluid_viscosity")),
87 40 : _permeability(getMaterialProperty<RealTensorValue>("permeability")),
88 20 : _num_nodes(0),
89 20 : _pp(0),
90 20 : _sat(0),
91 20 : _nodal_density(0),
92 20 : _dnodal_density_dp(0),
93 20 : _nodal_relperm(0),
94 40 : _dnodal_relperm_ds(0)
95 : {
96 20 : }
97 :
98 : void
99 1724 : Q2PPiecewiseLinearSink::prepareNodalValues()
100 : {
101 1724 : _num_nodes = _other_var_nodal.size();
102 :
103 : // set _pp and _sat variables
104 1724 : _pp.resize(_num_nodes);
105 1724 : _sat.resize(_num_nodes);
106 1724 : if (_var_is_pp)
107 : {
108 3134 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
109 : {
110 2544 : _pp[nodenum] = _var.dofValues()[nodenum];
111 2544 : _sat[nodenum] = _other_var_nodal[nodenum];
112 : }
113 : }
114 : else
115 : {
116 5854 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
117 : {
118 4720 : _pp[nodenum] = _other_var_nodal[nodenum];
119 4720 : _sat[nodenum] = _var.dofValues()[nodenum];
120 : }
121 : }
122 :
123 : // calculate derived things
124 1724 : if (_use_mobility)
125 : {
126 1180 : _nodal_density.resize(_num_nodes);
127 1180 : _dnodal_density_dp.resize(_num_nodes);
128 6268 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
129 : {
130 5088 : _nodal_density[nodenum] = _density.density(_pp[nodenum]);
131 5088 : _dnodal_density_dp[nodenum] = _density.ddensity(_pp[nodenum]);
132 : }
133 : }
134 :
135 1724 : if (_use_relperm)
136 : {
137 1180 : _nodal_relperm.resize(_num_nodes);
138 1180 : _dnodal_relperm_ds.resize(_num_nodes);
139 6268 : for (unsigned int nodenum = 0; nodenum < _num_nodes; ++nodenum)
140 : {
141 5088 : _nodal_relperm[nodenum] = _relperm.relperm(_sat[nodenum]);
142 5088 : _dnodal_relperm_ds[nodenum] = _relperm.drelperm(_sat[nodenum]);
143 : }
144 : }
145 1724 : }
146 :
147 : void
148 983 : Q2PPiecewiseLinearSink::computeResidual()
149 : {
150 983 : prepareNodalValues();
151 983 : IntegratedBC::computeResidual();
152 983 : }
153 :
154 : Real
155 9784 : Q2PPiecewiseLinearSink::computeQpResidual()
156 : {
157 : Real flux = 0;
158 : Real k = 0;
159 :
160 9784 : flux = _test[_i][_qp] * _sink_func.sample(_pp[_i]);
161 9784 : if (_use_mobility)
162 : {
163 7376 : k = (_permeability[0] * _normals[_qp]) * _normals[_qp]; // assume that _permeability is constant
164 : // throughout element so doesn't need to
165 : // be upwinded
166 7376 : flux *= _nodal_density[_i] * k / _viscosity;
167 : }
168 9784 : if (_use_relperm)
169 7376 : flux *= _nodal_relperm[_i];
170 :
171 9784 : flux *= _m_func.value(_t, _q_point[_qp]);
172 :
173 9784 : return flux;
174 : }
175 :
176 : void
177 313 : Q2PPiecewiseLinearSink::computeJacobian()
178 : {
179 313 : prepareNodalValues();
180 313 : IntegratedBC::computeJacobian();
181 313 : }
182 :
183 : Real
184 10912 : Q2PPiecewiseLinearSink::computeQpJacobian()
185 : {
186 10912 : return jac(_var.number());
187 : }
188 :
189 : void
190 428 : Q2PPiecewiseLinearSink::computeOffDiagJacobian(const unsigned int jvar)
191 : {
192 428 : prepareNodalValues();
193 428 : IntegratedBC::computeOffDiagJacobian(jvar);
194 428 : }
195 :
196 : Real
197 7744 : Q2PPiecewiseLinearSink::computeQpOffDiagJacobian(unsigned int jvar)
198 : {
199 7744 : if (jvar == _var.number() || jvar == _other_var_num)
200 7744 : return jac(jvar);
201 :
202 : return 0.0;
203 : }
204 :
205 : Real
206 18656 : 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 18656 : if (_i != _j)
216 : return 0.0; // residual at node _i only depends on variables at that node
217 :
218 4408 : flux = _sink_func.sample(_pp[_i]);
219 :
220 4408 : if (_var_is_pp)
221 : {
222 : // derivative of the _sink_func
223 1512 : if (wrt_num == _var.number())
224 888 : deriv = _sink_func.sampleDerivative(_pp[_i]);
225 : else
226 : deriv = 0;
227 :
228 : // add derivative of the mobility
229 1512 : if (_use_mobility)
230 : {
231 1512 : k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
232 1512 : mob = _nodal_density[_i] * k / _viscosity;
233 1512 : if (wrt_num == _var.number())
234 888 : mobp = _dnodal_density_dp[_i] * k / _viscosity; // else mobp = 0
235 1512 : deriv = mob * deriv + mobp * flux;
236 1512 : flux *= mob;
237 : }
238 :
239 : // add derivative of the relperm
240 1512 : if (_use_relperm)
241 : {
242 1512 : if (wrt_num == _other_var_num)
243 624 : deriv = _nodal_relperm[_i] * deriv + _dnodal_relperm_ds[_i] * flux;
244 : else
245 888 : deriv = _nodal_relperm[_i] * deriv;
246 : }
247 : }
248 : else
249 : {
250 : // derivative of the _sink_func
251 2896 : if (wrt_num == _other_var_num)
252 1184 : deriv = _sink_func.sampleDerivative(_pp[_i]);
253 : else
254 : deriv = 0;
255 :
256 : // add derivative of the mobility
257 2896 : if (_use_mobility)
258 : {
259 1512 : k = (_permeability[0] * _normals[_qp]) * _normals[_qp];
260 1512 : mob = _nodal_density[_i] * k / _viscosity;
261 1512 : if (wrt_num == _other_var_num)
262 624 : mobp = _dnodal_density_dp[_i] * k / _viscosity; // else mobp = 0
263 1512 : deriv = mob * deriv + mobp * flux;
264 1512 : flux *= mob;
265 : }
266 :
267 : // add derivative of the relperm
268 2896 : if (_use_relperm)
269 : {
270 1512 : if (wrt_num == _var.number())
271 888 : deriv = _nodal_relperm[_i] * deriv + _dnodal_relperm_ds[_i] * flux;
272 : else
273 624 : deriv = _nodal_relperm[_i] * deriv;
274 : }
275 : }
276 :
277 4408 : deriv *= _m_func.value(_t, _q_point[_qp]);
278 :
279 4408 : return _test[_i][_qp] * deriv * phi;
280 : }
|
