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 "VolumetricFlowRate.h"
11 : #include "MathFVUtils.h"
12 : #include "RhieChowInterpolatorBase.h"
13 : #include "NSFVUtils.h"
14 :
15 : #include <cmath>
16 :
17 : registerMooseObject("NavierStokesApp", VolumetricFlowRate);
18 :
19 : InputParameters
20 5532 : VolumetricFlowRate::validParams()
21 : {
22 5532 : InputParameters params = SideIntegralPostprocessor::validParams();
23 5532 : params.addClassDescription(
24 : "Computes the volumetric flow rate of an advected quantity through a sideset.");
25 11064 : params.addRequiredCoupledVar("vel_x", "The x-axis velocity");
26 11064 : params.addCoupledVar("vel_y", 0, "The y-axis velocity");
27 11064 : params.addCoupledVar("vel_z", 0, "The z-axis velocity");
28 11064 : params.addCoupledVar("advected_variable",
29 : 0,
30 : "The advected variable quantity of which to study the flow; useful for "
31 : "finite element simulations");
32 11064 : params.addParam<MooseFunctorName>("advected_mat_prop",
33 11064 : 0,
34 : "The advected material property of which to study the flow; "
35 : "useful for finite element simulations");
36 11064 : params.addParam<MooseFunctorName>("advected_quantity",
37 : "The quantity to advect. This is the canonical parameter to "
38 : "set the advected quantity when finite volume is being used.");
39 5532 : params += Moose::FV::interpolationParameters();
40 11064 : params.addParam<UserObjectName>("rhie_chow_user_object", "The rhie-chow user-object");
41 11064 : params.addParam<bool>("subtract_mesh_velocity",
42 : "To subtract the velocity of the potentially moving mesh. Defaults to true "
43 : "if a displaced problem exists, else false.");
44 5532 : return params;
45 0 : }
46 :
47 2911 : VolumetricFlowRate::VolumetricFlowRate(const InputParameters & parameters)
48 : : SideIntegralPostprocessor(parameters),
49 2911 : _vel_x(coupledValue("vel_x")),
50 2911 : _vel_y(coupledValue("vel_y")),
51 2911 : _vel_z(coupledValue("vel_z")),
52 2911 : _advected_variable_supplied(parameters.isParamSetByUser("advected_variable")),
53 2911 : _advected_variable(coupledValue("advected_variable")),
54 2911 : _advected_mat_prop_supplied(parameters.isParamSetByUser("advected_mat_prop")),
55 5822 : _advected_material_property(getFunctor<ADReal>("advected_mat_prop")),
56 10056 : _adv_quant(isParamValid("advected_quantity") ? &getFunctor<ADReal>("advected_quantity")
57 : : nullptr),
58 2911 : _rc_uo(isParamValid("rhie_chow_user_object")
59 5028 : ? &getUserObject<RhieChowFaceFluxProvider>("rhie_chow_user_object")
60 : : nullptr),
61 2911 : _subtract_mesh_velocity(isParamValid("subtract_mesh_velocity")
62 6032 : ? getParam<bool>("subtract_mesh_velocity")
63 5612 : : _fe_problem.haveDisplaced())
64 : {
65 : // Check that at most one advected quantity has been provided
66 2911 : if (_advected_variable_supplied && _advected_mat_prop_supplied)
67 0 : mooseError("VolumetricFlowRatePostprocessor should be provided either an advected variable "
68 : "or an advected material property");
69 :
70 : // Check that the user isn't trying to get face values for material properties
71 5822 : if (parameters.isParamSetByUser("advected_interp_method") && _advected_mat_prop_supplied)
72 0 : mooseWarning("Advected quantity interpolation methods are currently unavailable for "
73 : "advected material properties.");
74 :
75 2911 : _qp_integration = !getFieldVar("vel_x", 0)->isFV();
76 :
77 2911 : if (_advected_mat_prop_supplied)
78 80 : checkFunctorSupportsSideIntegration<ADReal>("advected_mat_prop", _qp_integration);
79 2911 : if (_adv_quant)
80 4234 : checkFunctorSupportsSideIntegration<ADReal>("advected_quantity", _qp_integration);
81 :
82 2911 : if (!_qp_integration)
83 : {
84 2117 : if (!_rc_uo)
85 0 : mooseError("We were instructed to use finite volume, but no Rhie-Chow user object is "
86 : "supplied. Please make sure to set the 'rhie_chow_user_object' parameter");
87 2117 : if (!_adv_quant)
88 0 : mooseError("We were instructed to use finite volume, but no 'advected_quantity' parameter is "
89 : "supplied.");
90 :
91 2117 : Moose::FV::setInterpolationMethods(*this, _advected_interp_method, _velocity_interp_method);
92 : }
93 :
94 2911 : if (_subtract_mesh_velocity && _rc_uo && !_rc_uo->supportMeshVelocity())
95 0 : paramError("subtract_mesh_velocity",
96 : "Rhie Chow user object does not support subtracting the mesh velocity");
97 2911 : if (_subtract_mesh_velocity && !_fe_problem.haveDisplaced())
98 0 : paramError(
99 : "subtract_mesh_velocity",
100 : "No displaced problem, thus the mesh velocity is 0 and does not need to be subtracted");
101 2911 : }
102 :
103 : void
104 2889 : VolumetricFlowRate::initialSetup()
105 : {
106 2889 : const auto * rc_base = dynamic_cast<const RhieChowInterpolatorBase *>(_rc_uo);
107 2095 : if (_rc_uo && rc_base &&
108 4826 : rc_base->velocityInterpolationMethod() == Moose::FV::InterpMethod::RhieChow &&
109 1937 : !rc_base->segregated())
110 : {
111 : // We must make sure the A coefficients in the Rhie Chow interpolator are present on
112 : // both sides of the boundaries so that interpolation coefficients may be computed
113 3810 : for (const auto bid : boundaryIDs())
114 1905 : const_cast<RhieChowInterpolatorBase *>(rc_base)->ghostADataOnBoundary(bid);
115 :
116 : // On INITIAL, we cannot compute Rhie Chow coefficients on internal surfaces because
117 : // - the time integrator is not ready to compute time derivatives
118 : // - the setup routine is called too early for porosity functions to be initialized
119 : // We must check that the boundaries requested are all external
120 1905 : if (getExecuteOnEnum().isValueSet(EXEC_INITIAL))
121 122 : for (const auto bid : boundaryIDs())
122 : {
123 62 : if (!_mesh.isBoundaryFullyExternalToSubdomains(bid, rc_base->blockIDs()))
124 2 : paramError(
125 : "execute_on",
126 : "Boundary '",
127 2 : _mesh.getBoundaryName(bid),
128 : "' (id=",
129 : bid,
130 : ") has been detected to be internal to the flow domain.\n"
131 : "Volumetric flow rates cannot be computed on internal flow boundaries on INITIAL");
132 : }
133 : }
134 2887 : }
135 :
136 : void
137 0 : VolumetricFlowRate::meshChanged()
138 : {
139 0 : initialSetup();
140 0 : }
141 :
142 : Real
143 20190 : VolumetricFlowRate::computeFaceInfoIntegral(const FaceInfo * fi)
144 : {
145 : mooseAssert(fi, "We should have a face info in " + name());
146 : mooseAssert(_adv_quant, "We should have an advected quantity in " + name());
147 20190 : const auto state = determineState();
148 :
149 : // Get face value for velocity
150 20190 : const auto face_flux = MetaPhysicL::raw_value(_rc_uo->getVolumetricFaceFlux(
151 20190 : _velocity_interp_method, *fi, state, _tid, _subtract_mesh_velocity));
152 :
153 20190 : const bool correct_skewness =
154 20190 : _advected_interp_method == Moose::FV::InterpMethod::SkewCorrectedAverage;
155 :
156 : mooseAssert(_adv_quant->hasFaceSide(*fi, true) || _adv_quant->hasFaceSide(*fi, false),
157 : "Advected quantity should be defined on one side of the face!");
158 :
159 20190 : const auto * elem = _adv_quant->hasFaceSide(*fi, true) ? fi->elemPtr() : fi->neighborPtr();
160 :
161 : const auto adv_quant_face = MetaPhysicL::raw_value(
162 60570 : (*_adv_quant)(Moose::FaceArg({fi,
163 20190 : Moose::FV::limiterType(_advected_interp_method),
164 20190 : face_flux > 0,
165 : correct_skewness,
166 : elem,
167 : nullptr}),
168 : state));
169 20190 : return face_flux * adv_quant_face;
170 : }
171 :
172 : Real
173 7536 : VolumetricFlowRate::computeQpIntegral()
174 : {
175 7536 : if (_advected_variable_supplied)
176 1944 : return _advected_variable[_qp] * RealVectorValue(_vel_x[_qp], _vel_y[_qp], _vel_z[_qp]) *
177 1944 : _normals[_qp];
178 5592 : else if (_advected_mat_prop_supplied)
179 432 : return MetaPhysicL::raw_value(_advected_material_property(
180 432 : Moose::ElemQpArg{_current_elem, _qp, _qrule, _q_point[_qp]}, determineState())) *
181 432 : RealVectorValue(_vel_x[_qp], _vel_y[_qp], _vel_z[_qp]) * _normals[_qp];
182 : else
183 5160 : return RealVectorValue(_vel_x[_qp], _vel_y[_qp], _vel_z[_qp]) * _normals[_qp];
184 : }
|
