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 "PressureDrop.h"
11 : #include "MathFVUtils.h"
12 : #include "NSFVUtils.h"
13 : #include "NS.h"
14 :
15 : #include <cmath>
16 :
17 : registerMooseObject("NavierStokesApp", PressureDrop);
18 :
19 : InputParameters
20 1732 : PressureDrop::validParams()
21 : {
22 1732 : InputParameters params = SideIntegralPostprocessor::validParams();
23 1732 : params.addClassDescription(
24 : "Computes the pressure drop between an upstream and a downstream boundary.");
25 :
26 3464 : params.addRequiredParam<MooseFunctorName>("pressure", "The pressure functor");
27 3464 : params.addParam<MooseFunctorName>("weighting_functor",
28 : "A vector functor to compute a flux to weigh the pressure for "
29 : "the pressure average computations");
30 3464 : params.addRequiredParam<std::vector<BoundaryName>>(
31 : "upstream_boundary", "The upstream surface (must also be specified in 'boundary' parameter");
32 3464 : params.addRequiredParam<std::vector<BoundaryName>>(
33 : "downstream_boundary",
34 : "The downstream surface (must also be specified in 'boundary' parameter");
35 :
36 3464 : params.addParam<MooseEnum>("weighting_interp_method",
37 3464 : Moose::FV::interpolationMethods(),
38 : "The interpolation to use for the weighting functor.");
39 1732 : return params;
40 0 : }
41 :
42 915 : PressureDrop::PressureDrop(const InputParameters & parameters)
43 : : SideIntegralPostprocessor(parameters),
44 915 : _pressure(getFunctor<Real>("pressure")),
45 1830 : _weighting_functor(isParamValid("weighting_functor")
46 1669 : ? &getFunctor<RealVectorValue>("weighting_functor")
47 915 : : nullptr)
48 : {
49 : // Examine the variable type of pressure to determine the type of integration
50 915 : const auto pressure_name = getParam<MooseFunctorName>("pressure");
51 915 : if (!_subproblem.hasVariable(pressure_name))
52 0 : paramError("pressure", "Pressure must be a variable");
53 915 : const auto * const pressure_var = &_subproblem.getVariable(_tid, pressure_name);
54 915 : _qp_integration = dynamic_cast<const MooseVariableFE<Real> *>(pressure_var);
55 :
56 915 : checkFunctorSupportsSideIntegration<Real>("pressure", _qp_integration);
57 915 : if (_weighting_functor)
58 754 : checkFunctorSupportsSideIntegration<RealVectorValue>("weighting_functor", _qp_integration);
59 :
60 915 : if (!_qp_integration)
61 1258 : Moose::FV::setInterpolationMethod(*this, _weight_interp_method, "weighting_interp_method");
62 :
63 572 : else if (parameters.isParamSetByUser("weighting_interp_method"))
64 2 : paramError("weighting_interp_method", "Face interpolation only specified for finite volume");
65 :
66 : // Only keep the ids
67 1826 : auto upstream_bdies = getParam<std::vector<BoundaryName>>("upstream_boundary");
68 1826 : auto downstream_bdies = getParam<std::vector<BoundaryName>>("downstream_boundary");
69 913 : _upstream_boundaries.resize(upstream_bdies.size());
70 913 : _downstream_boundaries.resize(downstream_bdies.size());
71 1828 : for (auto i : make_range(upstream_bdies.size()))
72 915 : _upstream_boundaries[i] = _mesh.getBoundaryID(upstream_bdies[i]);
73 1826 : for (auto i : make_range(downstream_bdies.size()))
74 913 : _downstream_boundaries[i] = _mesh.getBoundaryID(downstream_bdies[i]);
75 :
76 : // Check that boundary restriction makes sense
77 1826 : for (auto bdy : _upstream_boundaries)
78 915 : if (!hasBoundary(bdy))
79 2 : paramError("boundary",
80 2 : "Upstream boundary '" + _mesh.getBoundaryName(bdy) +
81 : "' is not included in boundary restriction");
82 1820 : for (auto bdy : _downstream_boundaries)
83 911 : if (!hasBoundary(bdy))
84 2 : paramError("boundary",
85 2 : "Downstream boundary '" + _mesh.getBoundaryName(bdy) +
86 : "' is not included in boundary restriction");
87 2727 : for (auto bdy : boundaryIDs())
88 1820 : if (std::find(_upstream_boundaries.begin(), _upstream_boundaries.end(), bdy) ==
89 1820 : _upstream_boundaries.end() &&
90 1820 : std::find(_downstream_boundaries.begin(), _downstream_boundaries.end(), bdy) ==
91 : _downstream_boundaries.end())
92 2 : paramError("boundary",
93 2 : "Boundary restriction on boundary '" + _mesh.getBoundaryName(bdy) +
94 : "' is not part of upstream or downstream boundaries");
95 1812 : for (auto bdy : _upstream_boundaries)
96 907 : if (std::find(_downstream_boundaries.begin(), _downstream_boundaries.end(), bdy) !=
97 : _downstream_boundaries.end())
98 2 : paramError("upstream_boundary",
99 2 : "Upstream boundary '" + _mesh.getBoundaryName(bdy) +
100 : "' is also a downstream boundary");
101 1810 : }
102 :
103 : void
104 2698 : PressureDrop::initialize()
105 : {
106 2698 : _weighted_pressure_upstream = 0;
107 2698 : _weighted_pressure_downstream = 0;
108 2698 : _weight_upstream = 0;
109 2698 : _weight_downstream = 0;
110 :
111 : // Build the face infos in all cases, needed to detect upstream/downstream status
112 2698 : if (_mesh.isFiniteVolumeInfoDirty())
113 65 : _mesh.setupFiniteVolumeMeshData();
114 2698 : }
115 :
116 : void
117 0 : PressureDrop::meshChanged()
118 : {
119 0 : initialize();
120 0 : }
121 :
122 : void
123 15142 : PressureDrop::execute()
124 : {
125 : // Determine if upstream or downstream boundary
126 : bool upstream = false;
127 : bool status_known = false;
128 15142 : getFaceInfos();
129 :
130 15142 : for (auto & fi : _face_infos)
131 : {
132 24667 : for (const auto bdy : fi->boundaryIDs())
133 : {
134 24667 : if (std::find(_upstream_boundaries.begin(), _upstream_boundaries.end(), bdy) !=
135 : _upstream_boundaries.end())
136 : {
137 : upstream = true;
138 : status_known = true;
139 : #ifdef NDEBUG
140 : break;
141 : #else
142 : if (std::find(_downstream_boundaries.begin(), _downstream_boundaries.end(), bdy) !=
143 : _downstream_boundaries.end())
144 : mooseError("Boundary '", _mesh.getBoundaryName(bdy), "' is both upstream and downstream");
145 : #endif
146 : }
147 : #ifndef NDEBUG
148 : if (upstream &&
149 : std::find(_downstream_boundaries.begin(), _downstream_boundaries.end(), bdy) !=
150 : _downstream_boundaries.end())
151 : mooseError("Face info is part of both upstream and downstream boundaries");
152 : #endif
153 17096 : if (!upstream &&
154 17096 : std::find(_downstream_boundaries.begin(), _downstream_boundaries.end(), bdy) !=
155 : _downstream_boundaries.end())
156 : status_known = true;
157 :
158 : // in debug mode we will check all boundaries the face info is a part of
159 : // to make sure they are consistently upstream or downstream
160 : #ifdef NDEBUG
161 : if (status_known)
162 : break;
163 : #endif
164 : }
165 : // we ll assume all face infos for a given element and side have the same upstream status
166 : // it seems very unlikely that upstream and downstream boundaries would be touching
167 : // and the delimitation would be within a refined element's face
168 15142 : if (status_known)
169 : break;
170 : }
171 :
172 15142 : if (upstream)
173 : {
174 : // Integration loops are different for FE and FV at this point
175 7571 : if (_qp_integration)
176 : {
177 2960 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
178 : {
179 2220 : _weighted_pressure_upstream +=
180 2220 : _JxW[_qp] * _coord[_qp] * computeQpWeightedPressureIntegral();
181 2220 : _weight_upstream += _JxW[_qp] * _coord[_qp] * computeQpWeightIntegral();
182 : }
183 : }
184 : else
185 : {
186 13662 : for (auto & fi : _face_infos)
187 : {
188 6831 : _weighted_pressure_upstream +=
189 6831 : fi->faceArea() * fi->faceCoord() * computeFaceInfoWeightedPressureIntegral(fi);
190 6831 : _weight_upstream += fi->faceArea() * fi->faceCoord() * computeFaceInfoWeightIntegral(fi);
191 : }
192 : }
193 : }
194 : // Downstream contributions
195 : else
196 : {
197 7571 : if (_qp_integration)
198 : {
199 2960 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
200 : {
201 2220 : _weighted_pressure_downstream +=
202 2220 : _JxW[_qp] * _coord[_qp] * computeQpWeightedPressureIntegral();
203 2220 : _weight_downstream += _JxW[_qp] * _coord[_qp] * computeQpWeightIntegral();
204 : }
205 : }
206 : else
207 : {
208 13662 : for (auto & fi : _face_infos)
209 : {
210 6831 : _weighted_pressure_downstream +=
211 6831 : fi->faceArea() * fi->faceCoord() * computeFaceInfoWeightedPressureIntegral(fi);
212 6831 : _weight_downstream += fi->faceArea() * fi->faceCoord() * computeFaceInfoWeightIntegral(fi);
213 : }
214 : }
215 : }
216 15142 : }
217 :
218 : Real
219 13662 : PressureDrop::computeFaceInfoWeightedPressureIntegral(const FaceInfo * const fi) const
220 : {
221 : mooseAssert(fi, "We should have a face info in " + name());
222 :
223 13662 : const bool correct_skewness =
224 13662 : _weight_interp_method == Moose::FV::InterpMethod::SkewCorrectedAverage;
225 13662 : const Moose::ElemQpArg elem_arg = {_current_elem, _qp, _qrule, _q_point[_qp]};
226 13662 : const auto state = determineState();
227 : mooseAssert(_qp == 0, "Only one quadrature point");
228 : mooseAssert(_pressure.hasFaceSide(*fi, true) || _pressure.hasFaceSide(*fi, true),
229 : "Pressure must be defined at least on one side of the face!");
230 13662 : const auto * elem = _pressure.hasFaceSide(*fi, true) ? fi->elemPtr() : fi->neighborPtr();
231 :
232 13662 : if (_weighting_functor)
233 : {
234 : mooseAssert(_pressure.hasFaceSide(*fi, true) == _weighting_functor->hasFaceSide(*fi, true),
235 : "Pressure and weighting functor have to be defined on the same side of the face!");
236 11790 : auto ssf = Moose::FaceArg(
237 : {fi,
238 11790 : Moose::FV::limiterType(_weight_interp_method),
239 23580 : MetaPhysicL::raw_value((*_weighting_functor)(elem_arg, state)) * fi->normal() > 0,
240 : correct_skewness,
241 : elem,
242 11790 : nullptr});
243 11790 : const auto face_weighting = MetaPhysicL::raw_value((*_weighting_functor)(ssf, state));
244 : // Dont use upwinding or an advection limiter for pressure
245 11790 : ssf.limiter_type = Moose::FV::LimiterType::CentralDifference;
246 11790 : return fi->normal() * face_weighting * _pressure(ssf, state);
247 : }
248 : else
249 : {
250 1872 : const auto ssf = Moose::FaceArg(
251 1872 : {fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, elem, nullptr});
252 1872 : return _pressure(ssf, state);
253 : }
254 : }
255 :
256 : Real
257 13662 : PressureDrop::computeFaceInfoWeightIntegral(const FaceInfo * fi) const
258 : {
259 : mooseAssert(fi, "We should have a face info in " + name());
260 13662 : const Moose::ElemQpArg elem_arg = {_current_elem, _qp, _qrule, _q_point[_qp]};
261 : mooseAssert(_qp == 0, "Only one quadrature point");
262 13662 : const auto state = determineState();
263 :
264 13662 : const bool correct_skewness =
265 13662 : _weight_interp_method == Moose::FV::InterpMethod::SkewCorrectedAverage;
266 :
267 13662 : if (_weighting_functor)
268 : {
269 11790 : const auto ssf = Moose::FaceArg(
270 : {fi,
271 11790 : Moose::FV::limiterType(_weight_interp_method),
272 35370 : MetaPhysicL::raw_value((*_weighting_functor)(elem_arg, state)) * fi->normal() > 0,
273 : correct_skewness,
274 11790 : _weighting_functor->hasFaceSide(*fi, true) ? fi->elemPtr() : fi->neighborPtr(),
275 23580 : nullptr});
276 23580 : return fi->normal() * MetaPhysicL::raw_value((*_weighting_functor)(ssf, state));
277 : }
278 : else
279 : return 1;
280 : }
281 :
282 : Real
283 4440 : PressureDrop::computeQpWeightedPressureIntegral() const
284 : {
285 : const Moose::ElemSideQpArg elem_side_arg = {
286 4440 : _current_elem, _current_side, _qp, _qrule, _q_point[_qp]};
287 4440 : const auto state = determineState();
288 4440 : if (_weighting_functor)
289 1596 : return (*_weighting_functor)(elem_side_arg, state) * _normals[_qp] *
290 1596 : _pressure(elem_side_arg, state);
291 : else
292 2844 : return _pressure(elem_side_arg, state);
293 : }
294 :
295 : Real
296 4440 : PressureDrop::computeQpWeightIntegral() const
297 : {
298 : const Moose::ElemSideQpArg elem_side_arg = {
299 4440 : _current_elem, _current_side, _qp, _qrule, _q_point[_qp]};
300 4440 : const auto state = determineState();
301 4440 : if (_weighting_functor)
302 1596 : return (*_weighting_functor)(elem_side_arg, state) * _normals[_qp];
303 : else
304 : return 1;
305 : }
306 :
307 : void
308 469 : PressureDrop::threadJoin(const UserObject & y)
309 : {
310 : const auto & pps = static_cast<const PressureDrop &>(y);
311 469 : _weighted_pressure_upstream += pps._weighted_pressure_upstream;
312 469 : _weighted_pressure_downstream += pps._weighted_pressure_downstream;
313 469 : _weight_upstream += pps._weight_upstream;
314 469 : _weight_downstream += pps._weight_downstream;
315 469 : }
316 :
317 : void
318 2223 : PressureDrop::finalize()
319 : {
320 2223 : gatherSum(_weighted_pressure_upstream);
321 2223 : gatherSum(_weighted_pressure_downstream);
322 2223 : gatherSum(_weight_upstream);
323 2223 : gatherSum(_weight_downstream);
324 2223 : }
325 :
326 : Real
327 2223 : PressureDrop::getValue() const
328 : {
329 2223 : if (MooseUtils::absoluteFuzzyEqual(_weight_upstream, 0) ||
330 : MooseUtils::absoluteFuzzyEqual(_weight_downstream, 0))
331 : {
332 2 : mooseWarning("Weight integral is 0 (downstream or upstream), either :\n"
333 : "- pressure drop value is queried before being computed\n"
334 : "- the weight flux integral is simply 0, the weighting is not appropriate");
335 0 : return 0;
336 : }
337 2221 : return _weighted_pressure_upstream / _weight_upstream -
338 2221 : _weighted_pressure_downstream / _weight_downstream;
339 : }
|
