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 : #pragma once
11 :
12 : #include "RhieChowFaceFluxProvider.h"
13 : #include "CellCenteredMapFunctor.h"
14 : #include "FaceCenteredMapFunctor.h"
15 : #include "VectorComponentFunctor.h"
16 : #include "LinearFVAnisotropicDiffusion.h"
17 : #include "LinearFVElementalKernel.h"
18 : #include <unordered_map>
19 : #include <set>
20 : #include <unordered_set>
21 :
22 : #include "libmesh/petsc_vector.h"
23 :
24 : class MooseMesh;
25 : class INSFVVelocityVariable;
26 : class INSFVPressureVariable;
27 : namespace libMesh
28 : {
29 : class Elem;
30 : class MeshBase;
31 : }
32 :
33 : /**
34 : * User object responsible for determining the face fluxes using the Rhie-Chow interpolation in a
35 : * segregated solver that uses the linear FV formulation.
36 : */
37 : class RhieChowMassFlux : public RhieChowFaceFluxProvider, public NonADFunctorInterface
38 : {
39 : public:
40 : static InputParameters validParams();
41 : RhieChowMassFlux(const InputParameters & params);
42 :
43 : /// Get the face velocity times density (used in advection terms)
44 : Real getMassFlux(const FaceInfo & fi) const;
45 :
46 : /// Get the volumetric face flux (used in advection terms)
47 : Real getVolumetricFaceFlux(const FaceInfo & fi) const;
48 :
49 : virtual Real getVolumetricFaceFlux(const Moose::FV::InterpMethod m,
50 : const FaceInfo & fi,
51 : const Moose::StateArg & time,
52 : const THREAD_ID tid,
53 : bool subtract_mesh_velocity) const override;
54 :
55 : /// Initialize the container for face velocities
56 : void initFaceMassFlux();
57 : /// Initialize the coupling fields (HbyA and Ainv)
58 : void initCouplingField();
59 : /// Update the values of the face velocities in the containers
60 : void computeFaceMassFlux();
61 : /// Update the cell values of the velocity variables
62 : void computeCellVelocity();
63 :
64 : virtual void meshChanged() override;
65 : virtual void initialize() override;
66 0 : virtual void execute() override {}
67 0 : virtual void finalize() override {}
68 : virtual void initialSetup() override;
69 :
70 : /**
71 : * Update the momentum system-related information
72 : * @param momentum_systems Pointers to the momentum systems which are solved for the momentum
73 : * vector components
74 : * @param pressure_system Reference to the pressure system
75 : * @param momentum_system_numbers The numbers of these systems
76 : */
77 : void linkMomentumPressureSystems(const std::vector<LinearSystem *> & momentum_systems,
78 : const LinearSystem & pressure_system,
79 : const std::vector<unsigned int> & momentum_system_numbers);
80 :
81 : /**
82 : * Computes the inverse of the diagonal (1/A) of the system matrix plus the H/A components for the
83 : * pressure equation plus Rhie-Chow interpolation.
84 : */
85 : void computeHbyA(const bool with_updated_pressure, const bool verbose);
86 :
87 : protected:
88 : /// Select the right pressure gradient field and return a reference to the container
89 : std::vector<std::unique_ptr<NumericVector<Number>>> &
90 : selectPressureGradient(const bool updated_pressure);
91 :
92 : /// Compute the cell volumes on the mesh
93 : void setupMeshInformation();
94 :
95 : /// Populate the face values of the H/A and 1/A fields
96 : void
97 : populateCouplingFunctors(const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_hbya,
98 : const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_Ainv);
99 :
100 : /**
101 : * Check the block consistency between the passed in \p var and us
102 : */
103 : template <typename VarType>
104 : void checkBlocks(const VarType & var) const;
105 :
106 0 : virtual bool supportMeshVelocity() const override { return false; }
107 :
108 : /// The \p MooseMesh that this user object operates on
109 : const MooseMesh & _moose_mesh;
110 :
111 : /// The \p libMesh mesh that this object acts on
112 : const libMesh::MeshBase & _mesh;
113 :
114 : /// The dimension of the mesh, e.g. 3 for hexes and tets, 2 for quads and tris
115 : const unsigned int _dim;
116 :
117 : /// The thread 0 copy of the pressure variable
118 : const MooseLinearVariableFVReal * const _p;
119 :
120 : /// The thread 0 copy of the x-velocity variable
121 : std::vector<const MooseLinearVariableFVReal *> _vel;
122 :
123 : /// Pointer to the pressure diffusion term in the pressure Poisson equation
124 : LinearFVAnisotropicDiffusion * _p_diffusion_kernel;
125 :
126 : /**
127 : * A map functor from faces to $HbyA_{ij} = (A_{offdiag}*\mathrm{(predicted~velocity)} -
128 : * \mathrm{Source})_{ij}/A_{ij}$. So this contains the off-diagonal part of the system matrix
129 : * multiplied by the predicted velocity minus the source terms from the right hand side of the
130 : * linearized momentum predictor step.
131 : */
132 : FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> _HbyA_flux;
133 :
134 : /**
135 : * We hold on to the cell-based HbyA vectors so that we can easily reconstruct the
136 : * cell velocities as well.
137 : */
138 : std::vector<std::unique_ptr<NumericVector<Number>>> _HbyA_raw;
139 :
140 : /**
141 : * A map functor from faces to $(1/A)_f$. Where $A_i$ is the diagonal of the system matrix
142 : * for the momentum equation.
143 : */
144 : FaceCenteredMapFunctor<RealVectorValue, std::unordered_map<dof_id_type, RealVectorValue>> _Ainv;
145 :
146 : /**
147 : * We hold on to the cell-based 1/A vectors so that we can easily reconstruct the
148 : * cell velocities as well.
149 : */
150 : std::vector<std::unique_ptr<NumericVector<Number>>> _Ainv_raw;
151 :
152 : std::unique_ptr<NumericVector<Number>> _A_avg;
153 :
154 : /**
155 : * A map functor from faces to mass fluxes which are used in the advection terms.
156 : */
157 : FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> & _face_mass_flux;
158 :
159 : /// Pointer to the body force terms
160 : std::vector<std::vector<LinearFVElementalKernel *>> _body_force_kernels;
161 : /// Vector of body force term names
162 : std::vector<std::vector<std::string>> _body_force_kernel_names;
163 :
164 : /**
165 : * for a PISO iteration we need to hold on to the original pressure gradient field.
166 : * Should not be used in other conditions.
167 : */
168 : std::vector<std::unique_ptr<NumericVector<Number>>> _grad_p_current;
169 :
170 : /**
171 : * Functor describing the density of the fluid
172 : */
173 : const Moose::Functor<Real> & _rho;
174 :
175 : /// Pointers to the linear system(s) in moose corresponding to the momentum equation(s)
176 : std::vector<LinearSystem *> _momentum_systems;
177 :
178 : /// Numbers of the momentum system(s)
179 : std::vector<unsigned int> _momentum_system_numbers;
180 :
181 : /// Global numbers of the momentum system(s)
182 : std::vector<unsigned int> _global_momentum_system_numbers;
183 :
184 : /// Pointers to the momentum equation implicit system(s) from libmesh
185 : std::vector<libMesh::LinearImplicitSystem *> _momentum_implicit_systems;
186 :
187 : /// Pointer to the pressure system
188 : const LinearSystem * _pressure_system;
189 :
190 : /// Global number of the pressure system
191 : unsigned int _global_pressure_system_number;
192 :
193 : /// We will hold a vector of cell volumes to make sure we can do volume corrections rapidly
194 : std::unique_ptr<NumericVector<Number>> _cell_volumes;
195 :
196 : /// Enumerator for the method used for pressure projection
197 : const MooseEnum _pressure_projection_method;
198 :
199 : private:
200 : /// The subset of the FaceInfo objects that actually cover the subdomains which the
201 : /// flow field is defined on. Cached for performance optimization.
202 : std::vector<const FaceInfo *> _flow_face_info;
203 : };
204 :
205 : template <typename VarType>
206 : void
207 1970 : RhieChowMassFlux::checkBlocks(const VarType & var) const
208 : {
209 1970 : const auto & var_blocks = var.blockIDs();
210 1970 : const auto & uo_blocks = blockIDs();
211 :
212 : // Error if this UO has any blocks that the variable does not
213 : std::set<SubdomainID> uo_blocks_minus_var_blocks;
214 1970 : std::set_difference(uo_blocks.begin(),
215 : uo_blocks.end(),
216 : var_blocks.begin(),
217 : var_blocks.end(),
218 : std::inserter(uo_blocks_minus_var_blocks, uo_blocks_minus_var_blocks.end()));
219 1970 : if (uo_blocks_minus_var_blocks.size() > 0)
220 0 : mooseError("Block restriction of interpolator user object '",
221 : this->name(),
222 : "' (",
223 : Moose::stringify(blocks()),
224 : ") includes blocks not in the block restriction of variable '",
225 : var.name(),
226 : "' (",
227 : Moose::stringify(var.blocks()),
228 : ")");
229 :
230 : // Get the blocks in the variable but not this UO
231 : std::set<SubdomainID> var_blocks_minus_uo_blocks;
232 1970 : std::set_difference(var_blocks.begin(),
233 : var_blocks.end(),
234 : uo_blocks.begin(),
235 : uo_blocks.end(),
236 : std::inserter(var_blocks_minus_uo_blocks, var_blocks_minus_uo_blocks.end()));
237 :
238 : // For each block in the variable but not this UO, error if there is connection
239 : // to any blocks on the UO.
240 1970 : for (auto & block_id : var_blocks_minus_uo_blocks)
241 : {
242 0 : const auto connected_blocks = _moose_mesh.getBlockConnectedBlocks(block_id);
243 : std::set<SubdomainID> connected_blocks_on_uo;
244 0 : std::set_intersection(connected_blocks.begin(),
245 : connected_blocks.end(),
246 : uo_blocks.begin(),
247 : uo_blocks.end(),
248 : std::inserter(connected_blocks_on_uo, connected_blocks_on_uo.end()));
249 0 : if (connected_blocks_on_uo.size() > 0)
250 0 : mooseError("Block restriction of interpolator user object '",
251 : this->name(),
252 : "' (",
253 : Moose::stringify(uo_blocks),
254 : ") doesn't match the block restriction of variable '",
255 : var.name(),
256 : "' (",
257 : Moose::stringify(var_blocks),
258 : ")");
259 : }
260 1970 : }
|
