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 11208 : VolumetricFlowRate::validParams()
21 : {
22 11208 : InputParameters params = SideIntegralPostprocessor::validParams();
23 11208 : params.addClassDescription(
24 : "Computes the volumetric flow rate of an advected quantity through a sideset.");
25 22416 : params.addRequiredCoupledVar("vel_x", "The x-axis velocity");
26 22416 : params.addCoupledVar("vel_y", 0, "The y-axis velocity");
27 22416 : params.addCoupledVar("vel_z", 0, "The z-axis velocity");
28 22416 : params.addCoupledVar("advected_variable",
29 : 0,
30 : "The advected variable quantity of which to study the flow; useful for "
31 : "finite element simulations");
32 22416 : params.addParam<MooseFunctorName>("advected_mat_prop",
33 22416 : 0,
34 : "The advected material property of which to study the flow; "
35 : "useful for finite element simulations");
36 22416 : 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 11208 : params += Moose::FV::interpolationParameters();
40 22416 : params.addParam<UserObjectName>("rhie_chow_user_object", "The rhie-chow user-object");
41 22416 : 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 11208 : return params;
45 0 : }
46 :
47 6023 : VolumetricFlowRate::VolumetricFlowRate(const InputParameters & parameters)
48 : : SideIntegralPostprocessor(parameters),
49 6023 : _vel_x(coupledValue("vel_x")),
50 6023 : _vel_y(coupledValue("vel_y")),
51 6023 : _vel_z(coupledValue("vel_z")),
52 6023 : _advected_variable_supplied(parameters.isParamSetByUser("advected_variable")),
53 6023 : _advected_variable(coupledValue("advected_variable")),
54 6023 : _advected_mat_prop_supplied(parameters.isParamSetByUser("advected_mat_prop")),
55 12046 : _advected_material_property(getFunctor<ADReal>("advected_mat_prop")),
56 20616 : _adv_quant(isParamValid("advected_quantity") ? &getFunctor<ADReal>("advected_quantity")
57 : : nullptr),
58 6023 : _rc_uo(isParamValid("rhie_chow_user_object")
59 10308 : ? &getUserObject<RhieChowFaceFluxProvider>("rhie_chow_user_object")
60 : : nullptr),
61 6023 : _subtract_mesh_velocity(isParamValid("subtract_mesh_velocity")
62 12366 : ? getParam<bool>("subtract_mesh_velocity")
63 11726 : : _fe_problem.haveDisplaced())
64 : {
65 : // Check that at most one advected quantity has been provided
66 6023 : 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 12046 : 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 6023 : _qp_integration = !getFieldVar("vel_x", 0)->isFV();
76 :
77 6023 : if (_advected_mat_prop_supplied)
78 176 : checkFunctorSupportsSideIntegration<ADReal>("advected_mat_prop", _qp_integration);
79 6023 : if (_adv_quant)
80 8570 : checkFunctorSupportsSideIntegration<ADReal>("advected_quantity", _qp_integration);
81 :
82 6023 : if (!_qp_integration)
83 : {
84 4285 : 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 4285 : if (!_adv_quant)
88 0 : mooseError("We were instructed to use finite volume, but no 'advected_quantity' parameter is "
89 : "supplied.");
90 :
91 4285 : Moose::FV::setInterpolationMethods(*this, _advected_interp_method, _velocity_interp_method);
92 : }
93 :
94 6023 : 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 6023 : 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 6023 : }
102 :
103 : void
104 6001 : VolumetricFlowRate::initialSetup()
105 : {
106 6001 : const auto * rc_base = dynamic_cast<const RhieChowInterpolatorBase *>(_rc_uo);
107 4263 : if (_rc_uo && rc_base &&
108 10098 : rc_base->velocityInterpolationMethod() == Moose::FV::InterpMethod::RhieChow &&
109 4097 : !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 8098 : for (const auto bid : boundaryIDs())
114 4049 : 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 4049 : if (getExecuteOnEnum().isValueSet(EXEC_INITIAL))
121 266 : for (const auto bid : boundaryIDs())
122 : {
123 134 : 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 5999 : }
135 :
136 : void
137 0 : VolumetricFlowRate::meshChanged()
138 : {
139 0 : initialSetup();
140 0 : }
141 :
142 : Real
143 26605 : 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 26605 : const auto state = determineState();
148 :
149 : // Get face value for velocity
150 26605 : const auto face_flux = MetaPhysicL::raw_value(_rc_uo->getVolumetricFaceFlux(
151 26605 : _velocity_interp_method, *fi, state, _tid, _subtract_mesh_velocity));
152 :
153 26605 : const bool correct_skewness =
154 26605 : _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 26605 : const auto * elem = _adv_quant->hasFaceSide(*fi, true) ? fi->elemPtr() : fi->neighborPtr();
160 :
161 : const auto adv_quant_face = MetaPhysicL::raw_value(
162 79815 : (*_adv_quant)(Moose::FaceArg({fi,
163 26605 : Moose::FV::limiterType(_advected_interp_method),
164 26605 : face_flux > 0,
165 : correct_skewness,
166 : elem,
167 : nullptr}),
168 : state));
169 26605 : return face_flux * adv_quant_face;
170 : }
171 :
172 : Real
173 11256 : VolumetricFlowRate::computeQpIntegral()
174 : {
175 11256 : if (_advected_variable_supplied)
176 2916 : return _advected_variable[_qp] * RealVectorValue(_vel_x[_qp], _vel_y[_qp], _vel_z[_qp]) *
177 2916 : _normals[_qp];
178 8340 : else if (_advected_mat_prop_supplied)
179 648 : return MetaPhysicL::raw_value(_advected_material_property(
180 648 : Moose::ElemQpArg{_current_elem, _qp, _qrule, _q_point[_qp]}, determineState())) *
181 648 : RealVectorValue(_vel_x[_qp], _vel_y[_qp], _vel_z[_qp]) * _normals[_qp];
182 : else
183 7692 : return RealVectorValue(_vel_x[_qp], _vel_y[_qp], _vel_z[_qp]) * _normals[_qp];
184 : }
|
