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