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 : /**
153 : * A map functor from faces to mass fluxes which are used in the advection terms.
154 : */
155 : FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> & _face_mass_flux;
156 :
157 : /// Pointer to the body force terms
158 : std::vector<std::vector<LinearFVElementalKernel *>> _body_force_kernels;
159 : /// Vector of body force term names
160 : std::vector<std::vector<std::string>> _body_force_kernel_names;
161 :
162 : /**
163 : * for a PISO iteration we need to hold on to the original pressure gradient field.
164 : * Should not be used in other conditions.
165 : */
166 : std::vector<std::unique_ptr<NumericVector<Number>>> _grad_p_current;
167 :
168 : /**
169 : * Functor describing the density of the fluid
170 : */
171 : const Moose::Functor<Real> & _rho;
172 :
173 : /// Pointers to the linear system(s) in moose corresponding to the momentum equation(s)
174 : std::vector<LinearSystem *> _momentum_systems;
175 :
176 : /// Numbers of the momentum system(s)
177 : std::vector<unsigned int> _momentum_system_numbers;
178 :
179 : /// Global numbers of the momentum system(s)
180 : std::vector<unsigned int> _global_momentum_system_numbers;
181 :
182 : /// Pointers to the momentum equation implicit system(s) from libmesh
183 : std::vector<libMesh::LinearImplicitSystem *> _momentum_implicit_systems;
184 :
185 : /// Pointer to the pressure system
186 : const LinearSystem * _pressure_system;
187 :
188 : /// Global number of the pressure system
189 : unsigned int _global_pressure_system_number;
190 :
191 : /// We will hold a vector of cell volumes to make sure we can do volume corrections rapidly
192 : std::unique_ptr<NumericVector<Number>> _cell_volumes;
193 :
194 : /// Enumerator for the method used for pressure projection
195 : const MooseEnum _pressure_projection_method;
196 :
197 : private:
198 : /// The subset of the FaceInfo objects that actually cover the subdomains which the
199 : /// flow field is defined on. Cached for performance optimization.
200 : std::vector<const FaceInfo *> _flow_face_info;
201 : };
202 :
203 : template <typename VarType>
204 : void
205 2637 : RhieChowMassFlux::checkBlocks(const VarType & var) const
206 : {
207 2637 : const auto & var_blocks = var.blockIDs();
208 2637 : const auto & uo_blocks = blockIDs();
209 :
210 : // Error if this UO has any blocks that the variable does not
211 : std::set<SubdomainID> uo_blocks_minus_var_blocks;
212 2637 : std::set_difference(uo_blocks.begin(),
213 : uo_blocks.end(),
214 : var_blocks.begin(),
215 : var_blocks.end(),
216 : std::inserter(uo_blocks_minus_var_blocks, uo_blocks_minus_var_blocks.end()));
217 2637 : if (uo_blocks_minus_var_blocks.size() > 0)
218 0 : mooseError("Block restriction of interpolator user object '",
219 : this->name(),
220 : "' (",
221 : Moose::stringify(blocks()),
222 : ") includes blocks not in the block restriction of variable '",
223 : var.name(),
224 : "' (",
225 : Moose::stringify(var.blocks()),
226 : ")");
227 :
228 : // Get the blocks in the variable but not this UO
229 : std::set<SubdomainID> var_blocks_minus_uo_blocks;
230 2637 : std::set_difference(var_blocks.begin(),
231 : var_blocks.end(),
232 : uo_blocks.begin(),
233 : uo_blocks.end(),
234 : std::inserter(var_blocks_minus_uo_blocks, var_blocks_minus_uo_blocks.end()));
235 :
236 : // For each block in the variable but not this UO, error if there is connection
237 : // to any blocks on the UO.
238 2637 : for (auto & block_id : var_blocks_minus_uo_blocks)
239 : {
240 0 : const auto connected_blocks = _moose_mesh.getBlockConnectedBlocks(block_id);
241 : std::set<SubdomainID> connected_blocks_on_uo;
242 0 : std::set_intersection(connected_blocks.begin(),
243 : connected_blocks.end(),
244 : uo_blocks.begin(),
245 : uo_blocks.end(),
246 : std::inserter(connected_blocks_on_uo, connected_blocks_on_uo.end()));
247 0 : if (connected_blocks_on_uo.size() > 0)
248 0 : mooseError("Block restriction of interpolator user object '",
249 : this->name(),
250 : "' (",
251 : Moose::stringify(uo_blocks),
252 : ") doesn't match the block restriction of variable '",
253 : var.name(),
254 : "' (",
255 : Moose::stringify(var_blocks),
256 : ")");
257 : }
258 2637 : }
|
