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 : #include "ExternalProblem.h"
12 : #include "PostprocessorInterface.h"
13 : #include "SubChannelApp.h"
14 : #include "QuadSubChannelMesh.h"
15 : #include "SolutionHandle.h"
16 : #include "SinglePhaseFluidProperties.h"
17 : #include <petscdm.h>
18 : #include <petscdmda.h>
19 : #include <petscksp.h>
20 : #include <petscsys.h>
21 : #include <petscvec.h>
22 : #include <petscsnes.h>
23 :
24 : /**
25 : * Base class for the 1-phase steady-state/transient subchannel solver.
26 : */
27 : class SubChannel1PhaseProblem : public ExternalProblem, public PostprocessorInterface
28 : {
29 : public:
30 : SubChannel1PhaseProblem(const InputParameters & params);
31 : virtual ~SubChannel1PhaseProblem();
32 :
33 : virtual void externalSolve() override;
34 : virtual void syncSolutions(Direction direction) override;
35 : virtual bool solverSystemConverged(const unsigned int) override;
36 : virtual void initialSetup() override;
37 :
38 : /// Function that computes the added heat coming from the fuel pins, for channel i_ch and cell iz
39 : virtual Real computeAddedHeatPin(unsigned int i_ch, unsigned int iz) = 0;
40 : /// Function that computes the heat added by the duct, for channel i_ch and cell iz
41 : Real computeAddedHeatDuct(unsigned int i_ch, unsigned int iz);
42 :
43 : protected:
44 : struct FrictionStruct
45 : {
46 : int i_ch;
47 : Real Re, S, w_perim;
48 : } _friction_args;
49 :
50 : /// Returns friction factor
51 : virtual Real computeFrictionFactor(FrictionStruct friction_args) = 0;
52 : /// Computes diversion crossflow per gap for block iblock
53 : void computeWijFromSolve(int iblock);
54 : /// Computes net diversion crossflow per channel for block iblock
55 : void computeSumWij(int iblock);
56 : /// Computes mass flow per channel for block iblock
57 : void computeMdot(int iblock);
58 : /// Computes turbulent crossflow per gap for block iblock
59 : void computeWijPrime(int iblock);
60 : /// Computes turbulent mixing coefficient
61 : virtual Real computeBeta(unsigned int i_gap, unsigned int iz, bool enthalpy) = 0;
62 : /// Computes Pressure Drop per channel for block iblock
63 : void computeDP(int iblock);
64 : /// Computes Pressure per channel for block iblock
65 : void computeP(int iblock);
66 : /// Computes Enthalpy per channel for block iblock
67 : virtual void computeh(int iblock) = 0;
68 : /// Computes Temperature per channel for block iblock
69 : void computeT(int iblock);
70 : /// Computes Density per channel for block iblock
71 : void computeRho(int iblock);
72 : /// Computes Viscosity per channel for block iblock
73 : void computeMu(int iblock);
74 : /// Computes Residual Matrix based on the lateral momentum conservation equation for block iblock
75 : void computeWijResidual(int iblock);
76 : /// Function that computes the width of the duct cell that the peripheral subchannel i_ch sees
77 : virtual Real getSubChannelPeripheralDuctWidth(unsigned int i_ch) = 0;
78 : /// Computes Residual Vector based on the lateral momentum conservation equation for block iblock & updates flow variables based on current crossflow solution
79 : libMesh::DenseVector<Real> residualFunction(int iblock, libMesh::DenseVector<Real> solution);
80 : /// Computes solution of nonlinear equation using snes and provided a residual in a formFunction
81 : PetscErrorCode petscSnesSolver(int iblock,
82 : const libMesh::DenseVector<Real> & solution,
83 : libMesh::DenseVector<Real> & root);
84 : /// This is the residual Vector function in a form compatible with the SNES PETC solvers
85 : friend PetscErrorCode formFunction(SNES snes, Vec x, Vec f, void * ctx);
86 :
87 : /// Computes implicit solve using PetSc
88 : PetscErrorCode implicitPetscSolve(int iblock);
89 :
90 : /// Function to initialize the solution & geometry fields
91 : virtual void initializeSolution() = 0;
92 :
93 : /// Functions that computes the interpolation scheme given the Peclet number
94 : PetscScalar computeInterpolationCoefficients(PetscScalar Peclet = 0.0);
95 : PetscScalar
96 : computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet = 0.0);
97 :
98 : /// inline function that is used to define the gravity direction
99 381 : Real computeGravityDir(const MooseEnum & dir) const
100 : {
101 : switch (dir)
102 : {
103 : case 0: // counter_flow
104 : return 1.0;
105 : case 1: // co_flow
106 : return -1.0;
107 : case 2: // none
108 : return 0.0;
109 0 : default:
110 0 : mooseError(name(), ": Invalid gravity direction: expected counter_flow, co_flow, or none");
111 : }
112 : }
113 :
114 : PetscErrorCode cleanUp();
115 : SubChannelMesh & _subchannel_mesh;
116 : /// number of axial blocks
117 : unsigned int _n_blocks;
118 : libMesh::DenseMatrix<Real> & _Wij;
119 : libMesh::DenseMatrix<Real> _Wij_old;
120 : libMesh::DenseMatrix<Real> _WijPrime;
121 : libMesh::DenseMatrix<Real> _Wij_residual_matrix;
122 : const Real _g_grav;
123 : const Real & _kij;
124 : unsigned int _n_cells;
125 : unsigned int _n_gaps;
126 : unsigned int _n_pins;
127 : unsigned int _n_channels;
128 : unsigned int _block_size;
129 : /// axial location of nodes
130 : std::vector<Real> _z_grid;
131 : Real _one;
132 : /// Flag that activates or deactivates the transient parts of the equations we solve by multiplication
133 : Real _TR;
134 : /// Flag that activates or deactivates the calculation of density
135 : const bool _compute_density;
136 : /// Flag that activates or deactivates the calculation of viscosity
137 : const bool _compute_viscosity;
138 : /// Flag that informs if we need to solve the Enthalpy/Temperature equations or not
139 : const bool _compute_power;
140 : /// Flag that informs if there is a pin mesh or not
141 : const bool _pin_mesh_exist;
142 : /// Flag that informs if there is a pin mesh or not
143 : const bool _duct_mesh_exist;
144 : /// Variable that informs whether we exited external solve with a converged solution or not
145 : bool _converged;
146 : /// Time step
147 : Real _dt;
148 : /// Outlet Pressure
149 : const PostprocessorValue & _P_out;
150 : /// Turbulent modeling parameter used in axial momentum equation
151 : const Real & _CT;
152 : /// Convergence tolerance for the pressure loop in external solve
153 : const Real & _P_tol;
154 : /// Convergence tolerance for the temperature loop in external solve
155 : const Real & _T_tol;
156 : /// Maximum iterations for the inner temperature loop
157 : const int & _T_maxit;
158 : /// The relative convergence tolerance, (relative decrease) for the ksp linear solver
159 : const PetscReal & _rtol;
160 : /// The absolute convergence tolerance for the ksp linear solver
161 : const PetscReal & _atol;
162 : /// The divergence tolerance for the ksp linear solver
163 : const PetscReal & _dtol;
164 : /// The maximum number of iterations to use for the ksp linear solver
165 : const PetscInt & _maxit;
166 : /// The interpolation method used in constructing the systems
167 : const MooseEnum _interpolation_scheme;
168 : /// The direction of gravity
169 : const MooseEnum _gravity_direction;
170 : const Real _dir_grav;
171 : /// Flag to define the usage of a implicit or explicit solution
172 : const bool _implicit_bool;
173 : /// Flag to define the usage of staggered or collocated pressure
174 : const bool _staggered_pressure_bool;
175 : /// Segregated solve
176 : const bool _segregated_bool;
177 : /// Thermal monolithic bool
178 : const bool _monolithic_thermal_bool;
179 : /// Boolean to printout information related to subchannel solve
180 : const bool _verbose_subchannel;
181 : /// Flag that activates the effect of deformation (pin/duct) based on the auxvalues for displacement, Dpin
182 : const bool _deformation;
183 :
184 : /// Solutions handles and link to TH tables properties
185 : const SinglePhaseFluidProperties * _fp;
186 : std::unique_ptr<SolutionHandle> _mdot_soln;
187 : std::unique_ptr<SolutionHandle> _SumWij_soln;
188 : std::unique_ptr<SolutionHandle> _P_soln;
189 : std::unique_ptr<SolutionHandle> _DP_soln;
190 : std::unique_ptr<SolutionHandle> _h_soln;
191 : std::unique_ptr<SolutionHandle> _T_soln;
192 : std::unique_ptr<SolutionHandle> _Tpin_soln;
193 : std::unique_ptr<SolutionHandle> _Dpin_soln;
194 : std::unique_ptr<SolutionHandle> _rho_soln;
195 : std::unique_ptr<SolutionHandle> _mu_soln;
196 : std::unique_ptr<SolutionHandle> _S_flow_soln;
197 : std::unique_ptr<SolutionHandle> _w_perim_soln;
198 : std::unique_ptr<SolutionHandle> _q_prime_soln;
199 : std::unique_ptr<SolutionHandle> _duct_heat_flux_soln; // Only used for ducted assemblies
200 : std::unique_ptr<SolutionHandle> _Tduct_soln; // Only used for ducted assemblies
201 : std::unique_ptr<SolutionHandle> _displacement_soln;
202 :
203 : /// Petsc Functions
204 69912 : inline PetscErrorCode createPetscVector(Vec & v, PetscInt n)
205 : {
206 : PetscFunctionBegin;
207 69912 : LibmeshPetscCall(VecCreate(PETSC_COMM_SELF, &v));
208 69912 : LibmeshPetscCall(PetscObjectSetName((PetscObject)v, "Solution"));
209 69912 : LibmeshPetscCall(VecSetSizes(v, PETSC_DECIDE, n));
210 69912 : LibmeshPetscCall(VecSetFromOptions(v));
211 69912 : LibmeshPetscCall(VecZeroEntries(v));
212 69912 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
213 : }
214 :
215 7332 : inline PetscErrorCode createPetscMatrix(Mat & M, PetscInt n, PetscInt m)
216 : {
217 : PetscFunctionBegin;
218 7332 : LibmeshPetscCall(MatCreate(PETSC_COMM_SELF, &M));
219 7332 : LibmeshPetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, n, m));
220 7332 : LibmeshPetscCall(MatSetFromOptions(M));
221 7332 : LibmeshPetscCall(MatSetUp(M));
222 7332 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
223 : }
224 :
225 : template <class T>
226 : PetscErrorCode populateVectorFromDense(Vec & x,
227 : const T & solution,
228 : const unsigned int first_axial_level,
229 : const unsigned int last_axial_level,
230 : const unsigned int cross_dimension);
231 : template <class T>
232 : PetscErrorCode populateDenseFromVector(const Vec & x,
233 : T & solution,
234 : const unsigned int first_axial_level,
235 : const unsigned int last_axial_level,
236 : const unsigned int cross_dimension);
237 : template <class T>
238 : PetscErrorCode populateVectorFromHandle(Vec & x,
239 : const T & solution,
240 : const unsigned int first_axial_level,
241 : const unsigned int last_axial_level,
242 : const unsigned int cross_dimension);
243 : template <class T>
244 : PetscErrorCode populateSolutionChan(const Vec & x,
245 : T & solution,
246 : const unsigned int first_axial_level,
247 : const unsigned int last_axial_level,
248 : const unsigned int cross_dimension);
249 :
250 : //// Matrices and vectors to be used in implicit assembly
251 : /// Mass conservation
252 : /// Mass conservation - sum of cross fluxes
253 : Mat _mc_sumWij_mat;
254 : Vec _Wij_vec;
255 : Vec _prod;
256 : Vec _prodp;
257 : /// Mass conservation - axial convection
258 : Mat _mc_axial_convection_mat;
259 : Vec _mc_axial_convection_rhs;
260 : /// Mass conservation - density time derivative
261 : /// No implicit matrix
262 :
263 : /// Axial momentum
264 : /// Axial momentum conservation - compute turbulent cross fluxes
265 : Mat _amc_turbulent_cross_flows_mat;
266 : Vec _amc_turbulent_cross_flows_rhs;
267 : /// Axial momentum conservation - time derivative
268 : Mat _amc_time_derivative_mat;
269 : Vec _amc_time_derivative_rhs;
270 : /// Axial momentum conservation - advective (Eulerian) derivative
271 : Mat _amc_advective_derivative_mat;
272 : Vec _amc_advective_derivative_rhs;
273 : /// Axial momentum conservation - cross flux derivative
274 : Mat _amc_cross_derivative_mat;
275 : Vec _amc_cross_derivative_rhs;
276 : /// Axial momentum conservation - friction force
277 : Mat _amc_friction_force_mat;
278 : Vec _amc_friction_force_rhs;
279 : /// Axial momentum conservation - buoyancy force
280 : /// No implicit matrix
281 : Vec _amc_gravity_rhs;
282 : /// Axial momentum conservation - pressure force
283 : Mat _amc_pressure_force_mat;
284 : Vec _amc_pressure_force_rhs;
285 : /// Axial momentum system matrix
286 : Mat _amc_sys_mdot_mat;
287 : Vec _amc_sys_mdot_rhs;
288 :
289 : /// Cross momentum
290 : /// Cross momentum conservation - time derivative
291 : Mat _cmc_time_derivative_mat;
292 : Vec _cmc_time_derivative_rhs;
293 : /// Cross momentum conservation - advective (Eulerian) derivative
294 : Mat _cmc_advective_derivative_mat;
295 : Vec _cmc_advective_derivative_rhs;
296 : /// Cross momentum conservation - friction force
297 : Mat _cmc_friction_force_mat;
298 : Vec _cmc_friction_force_rhs;
299 : /// Cross momentum conservation - pressure force
300 : Mat _cmc_pressure_force_mat;
301 : Vec _cmc_pressure_force_rhs;
302 : /// Lateral momentum system matrix
303 : Mat _cmc_sys_Wij_mat;
304 : Vec _cmc_sys_Wij_rhs;
305 :
306 : /// Enthalpy
307 : /// Enthalpy conservation - time derivative
308 : Mat _hc_time_derivative_mat;
309 : Vec _hc_time_derivative_rhs;
310 : /// Enthalpy conservation - advective (Eulerian) derivative;
311 : Mat _hc_advective_derivative_mat;
312 : Vec _hc_advective_derivative_rhs;
313 : /// Enthalpy conservation - cross flux derivative
314 : Mat _hc_cross_derivative_mat;
315 : Vec _hc_cross_derivative_rhs;
316 : /// Enthalpy conservation - source and sink
317 : Vec _hc_added_heat_rhs;
318 : /// System matrices
319 : Mat _hc_sys_h_mat;
320 : Vec _hc_sys_h_rhs;
321 :
322 : /// Added resistances for monolithic convergence
323 : PetscScalar _added_K = 0.0;
324 : PetscScalar _added_K_old = 1000.0;
325 : PetscScalar _max_sumWij;
326 : PetscScalar _max_sumWij_new;
327 : PetscScalar _correction_factor = 1.0;
328 :
329 : public:
330 : static InputParameters validParams();
331 : };
332 :
333 : template <class T>
334 : PetscErrorCode
335 22446 : SubChannel1PhaseProblem::populateDenseFromVector(const Vec & x,
336 : T & loc_solution,
337 : const unsigned int first_axial_level,
338 : const unsigned int last_axial_level,
339 : const unsigned int cross_dimension)
340 : {
341 : PetscScalar * xx;
342 :
343 : PetscFunctionBegin;
344 22446 : LibmeshPetscCall(VecGetArray(x, &xx));
345 409470 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
346 : {
347 387024 : unsigned int iz_ind = iz - first_axial_level;
348 25327032 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
349 : {
350 24940008 : loc_solution(i_l, iz) = xx[iz_ind * cross_dimension + i_l];
351 : }
352 : }
353 22446 : LibmeshPetscCall(VecRestoreArray(x, &xx));
354 :
355 22446 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
356 : }
357 :
358 : template <class T>
359 : PetscErrorCode
360 59256 : SubChannel1PhaseProblem::populateVectorFromHandle(Vec & x,
361 : const T & loc_solution,
362 : const unsigned int first_axial_level,
363 : const unsigned int last_axial_level,
364 : const unsigned int cross_dimension)
365 : {
366 : PetscScalar * xx;
367 :
368 : PetscFunctionBegin;
369 59256 : LibmeshPetscCall(VecGetArray(x, &xx));
370 956232 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
371 : {
372 896976 : unsigned int iz_ind = iz - first_axial_level;
373 39024000 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
374 : {
375 38127024 : auto * loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
376 38127024 : xx[iz_ind * cross_dimension + i_l] = loc_solution(loc_node);
377 : }
378 : }
379 59256 : LibmeshPetscCall(VecRestoreArray(x, &xx));
380 :
381 59256 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
382 : }
383 :
384 : template <class T>
385 : PetscErrorCode
386 22446 : SubChannel1PhaseProblem::populateVectorFromDense(Vec & x,
387 : const T & loc_solution,
388 : const unsigned int first_axial_level,
389 : const unsigned int last_axial_level,
390 : const unsigned int cross_dimension)
391 : {
392 : PetscScalar * xx;
393 : PetscFunctionBegin;
394 22446 : LibmeshPetscCall(VecGetArray(x, &xx));
395 412164 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
396 : {
397 389718 : unsigned int iz_ind = iz - first_axial_level;
398 25510734 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
399 : {
400 25121016 : xx[iz_ind * cross_dimension + i_l] = loc_solution(i_l, iz);
401 : }
402 : }
403 22446 : LibmeshPetscCall(VecRestoreArray(x, &xx));
404 22446 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
405 : }
406 :
407 : template <class T>
408 : PetscErrorCode
409 34116 : SubChannel1PhaseProblem::populateSolutionChan(const Vec & x,
410 : T & loc_solution,
411 : const unsigned int first_axial_level,
412 : const unsigned int last_axial_level,
413 : const unsigned int cross_dimension)
414 : {
415 : PetscScalar * xx;
416 : PetscFunctionBegin;
417 34116 : LibmeshPetscCall(VecGetArray(x, &xx));
418 : Node * loc_node;
419 450648 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
420 : {
421 416532 : unsigned int iz_ind = iz - first_axial_level;
422 17513580 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
423 : {
424 17097048 : loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
425 17097048 : loc_solution.set(loc_node, xx[iz_ind * cross_dimension + i_l]);
426 : }
427 : }
428 34116 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
429 : }
|
