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 :
53 424240 : struct NusseltStruct
54 : {
55 : Real Re, Pr;
56 : unsigned int i_pin, iz, i_ch;
57 : MooseEnum htc_correlation;
58 : // parameterized constructor
59 : NusseltStruct(Real Re_,
60 : Real Pr_,
61 : unsigned int i_pin_,
62 : unsigned int iz_,
63 : unsigned int i_ch_,
64 : const MooseEnum & htc_corr)
65 424240 : : Re(Re_), Pr(Pr_), i_pin(i_pin_), iz(iz_), i_ch(i_ch_), htc_correlation(htc_corr)
66 : {
67 17220 : }
68 : };
69 :
70 : /// The correlation used for computing the heat transfer correlation near the pin
71 : const MooseEnum _pin_htc_correlation;
72 : /// The correlation used for computing the heat transfer correlation near the duct
73 : const MooseEnum _duct_htc_correlation;
74 : NusseltStruct _nusselt_args;
75 :
76 : /// Function that computes the Nusselt number given a heat exchange correlation
77 : Real computeNusseltNumber(const NusseltStruct & nusselt_args);
78 :
79 : /// Computes diversion crossflow per gap for block iblock
80 : void computeWijFromSolve(int iblock);
81 : /// Computes net diversion crossflow per channel for block iblock
82 : void computeSumWij(int iblock);
83 : /// Computes mass flow per channel for block iblock
84 : void computeMdot(int iblock);
85 : /// Computes turbulent crossflow per gap for block iblock
86 : void computeWijPrime(int iblock);
87 : /// Computes turbulent mixing coefficient
88 : virtual Real computeBeta(unsigned int i_gap, unsigned int iz, bool enthalpy) = 0;
89 : /// Computes Pressure Drop per channel for block iblock
90 : void computeDP(int iblock);
91 : /// Computes Pressure per channel for block iblock
92 : void computeP(int iblock);
93 : /// Computes Enthalpy per channel for block iblock
94 : virtual void computeh(int iblock) = 0;
95 : /// Computes Temperature per channel for block iblock
96 : void computeT(int iblock);
97 : /// Computes Density per channel for block iblock
98 : void computeRho(int iblock);
99 : /// Computes Viscosity per channel for block iblock
100 : void computeMu(int iblock);
101 : /// Computes Residual Matrix based on the lateral momentum conservation equation for block iblock
102 : void computeWijResidual(int iblock);
103 : /// Function that computes the width of the duct cell that the peripheral subchannel i_ch sees
104 : virtual Real getSubChannelPeripheralDuctWidth(unsigned int i_ch) = 0;
105 : /// Computes Residual Vector based on the lateral momentum conservation equation for block iblock & updates flow variables based on current crossflow solution
106 : libMesh::DenseVector<Real> residualFunction(int iblock, libMesh::DenseVector<Real> solution);
107 : /// Computes solution of nonlinear equation using snes and provided a residual in a formFunction
108 : PetscErrorCode petscSnesSolver(int iblock,
109 : const libMesh::DenseVector<Real> & solution,
110 : libMesh::DenseVector<Real> & root);
111 : /// This is the residual Vector function in a form compatible with the SNES PETC solvers
112 : friend PetscErrorCode formFunction(SNES snes, Vec x, Vec f, void * ctx);
113 :
114 : /// Computes implicit solve using PetSc
115 : PetscErrorCode implicitPetscSolve(int iblock);
116 :
117 : /// Function to initialize the solution & geometry fields
118 : virtual void initializeSolution() = 0;
119 :
120 : /// Functions that computes the interpolation scheme given the Peclet number
121 : PetscScalar computeInterpolationCoefficients(PetscScalar Peclet = 0.0);
122 : PetscScalar
123 : computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet = 0.0);
124 :
125 : /// inline function that is used to define the gravity direction
126 372 : Real computeGravityDir(const MooseEnum & dir) const
127 : {
128 : switch (dir)
129 : {
130 : case 0: // counter_flow
131 : return 1.0;
132 : case 1: // co_flow
133 : return -1.0;
134 : case 2: // none
135 : return 0.0;
136 0 : default:
137 0 : mooseError(name(), ": Invalid gravity direction: expected counter_flow, co_flow, or none");
138 : }
139 : }
140 :
141 : /**
142 : * Solve a linear system (A * x = rhs) with a simple PCJACOBI KSP and populate the
143 : * enthalpy solution into _h_soln for nodes [first_node, last_node].
144 : *
145 : * Uses member tolerances (_rtol, _atol, _dtol, _maxit), mesh (_subchannel_mesh),
146 : * channel count (_n_channels), and error/solution handles (mooseError, _h_soln).
147 : *
148 : * @param A PETSc matrix (operators)
149 : * @param rhs PETSc vector (right-hand side)
150 : * @param first_node inclusive start axial node index
151 : * @param last_node inclusive end axial node index
152 : * @param ksp_prefix options prefix for KSP (e.g. "h_sys_"), may be nullptr
153 : */
154 : PetscErrorCode solveAndPopulateEnthalpy(
155 : Mat A, Vec rhs, unsigned int first_node, unsigned int last_node, const char * ksp_prefix);
156 :
157 : PetscErrorCode cleanUp();
158 : SubChannelMesh & _subchannel_mesh;
159 : /// number of axial blocks
160 : unsigned int _n_blocks;
161 : libMesh::DenseMatrix<Real> & _Wij;
162 : libMesh::DenseMatrix<Real> _Wij_old;
163 : libMesh::DenseMatrix<Real> _WijPrime;
164 : libMesh::DenseMatrix<Real> _Wij_residual_matrix;
165 : const Real _g_grav;
166 : const Real & _kij;
167 : unsigned int _n_cells;
168 : unsigned int _n_gaps;
169 : unsigned int _n_pins;
170 : unsigned int _n_channels;
171 : unsigned int _block_size;
172 : /// axial location of nodes
173 : std::vector<Real> _z_grid;
174 : Real _one;
175 : /// Flag that activates or deactivates the transient parts of the equations we solve by multiplication
176 : Real _TR;
177 : /// Flag that activates or deactivates the calculation of density
178 : const bool _compute_density;
179 : /// Flag that activates or deactivates the calculation of viscosity
180 : const bool _compute_viscosity;
181 : /// Flag that informs if we need to solve the Enthalpy/Temperature equations or not
182 : const bool _compute_power;
183 : /// Flag that informs if there is a pin mesh or not
184 : const bool _pin_mesh_exist;
185 : /// Flag that informs if there is a pin mesh or not
186 : const bool _duct_mesh_exist;
187 : /// Variable that informs whether we exited external solve with a converged solution or not
188 : bool _converged;
189 : /// Time step
190 : Real _dt;
191 : /// Outlet Pressure
192 : const PostprocessorValue & _P_out;
193 : /// Turbulent modeling parameter used in axial momentum equation
194 : const Real & _CT;
195 : /// Convergence tolerance for the pressure loop in external solve
196 : const Real & _P_tol;
197 : /// Convergence tolerance for the temperature loop in external solve
198 : const Real & _T_tol;
199 : /// Maximum iterations for the inner temperature loop
200 : const int & _T_maxit;
201 : /// The relative convergence tolerance, (relative decrease) for the ksp linear solver
202 : const PetscReal & _rtol;
203 : /// The absolute convergence tolerance for the ksp linear solver
204 : const PetscReal & _atol;
205 : /// The divergence tolerance for the ksp linear solver
206 : const PetscReal & _dtol;
207 : /// The maximum number of iterations to use for the ksp linear solver
208 : const PetscInt & _maxit;
209 : /// The interpolation method used in constructing the systems
210 : const MooseEnum _interpolation_scheme;
211 : /// The direction of gravity
212 : const MooseEnum _gravity_direction;
213 : const Real _dir_grav;
214 : /// Flag to define the usage of a implicit or explicit solution
215 : const bool _implicit_bool;
216 : /// Flag to define the usage of staggered or collocated pressure
217 : const bool _staggered_pressure_bool;
218 : /// Segregated solve
219 : const bool _segregated_bool;
220 : /// Boolean to printout information related to subchannel solve
221 : const bool _verbose_subchannel;
222 : /// Flag that activates the effect of deformation (pin/duct) based on the auxvalues for displacement, Dpin
223 : const bool _deformation;
224 :
225 : /// Solutions handles and link to TH tables properties
226 : const SinglePhaseFluidProperties * _fp;
227 : std::unique_ptr<SolutionHandle> _mdot_soln;
228 : std::unique_ptr<SolutionHandle> _SumWij_soln;
229 : std::unique_ptr<SolutionHandle> _P_soln;
230 : std::unique_ptr<SolutionHandle> _DP_soln;
231 : std::unique_ptr<SolutionHandle> _h_soln;
232 : std::unique_ptr<SolutionHandle> _T_soln;
233 : std::unique_ptr<SolutionHandle> _Tpin_soln;
234 : std::unique_ptr<SolutionHandle> _Dpin_soln;
235 : std::unique_ptr<SolutionHandle> _rho_soln;
236 : std::unique_ptr<SolutionHandle> _mu_soln;
237 : std::unique_ptr<SolutionHandle> _S_flow_soln;
238 : std::unique_ptr<SolutionHandle> _w_perim_soln;
239 : std::unique_ptr<SolutionHandle> _q_prime_soln;
240 : std::unique_ptr<SolutionHandle> _duct_heat_flux_soln; // Only used for ducted assemblies
241 : std::unique_ptr<SolutionHandle> _Tduct_soln; // Only used for ducted assemblies
242 : std::unique_ptr<SolutionHandle> _displacement_soln;
243 :
244 : /// Petsc Functions
245 55719 : inline PetscErrorCode createPetscVector(Vec & v, PetscInt n)
246 : {
247 : PetscFunctionBegin;
248 55719 : LibmeshPetscCall(VecCreate(PETSC_COMM_SELF, &v));
249 55719 : LibmeshPetscCall(PetscObjectSetName((PetscObject)v, "Solution"));
250 55719 : LibmeshPetscCall(VecSetSizes(v, PETSC_DECIDE, n));
251 55719 : LibmeshPetscCall(VecSetFromOptions(v));
252 55719 : LibmeshPetscCall(VecZeroEntries(v));
253 55719 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
254 : }
255 :
256 7197 : inline PetscErrorCode createPetscMatrix(Mat & M, PetscInt n, PetscInt m)
257 : {
258 : PetscFunctionBegin;
259 7197 : LibmeshPetscCall(MatCreate(PETSC_COMM_SELF, &M));
260 7197 : LibmeshPetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, n, m));
261 7197 : LibmeshPetscCall(MatSetFromOptions(M));
262 7197 : LibmeshPetscCall(MatSetUp(M));
263 7197 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
264 : }
265 :
266 : template <class T>
267 : PetscErrorCode populateVectorFromDense(Vec & x,
268 : const T & solution,
269 : const unsigned int first_axial_level,
270 : const unsigned int last_axial_level,
271 : const unsigned int cross_dimension);
272 : template <class T>
273 : PetscErrorCode populateDenseFromVector(const Vec & x,
274 : T & solution,
275 : const unsigned int first_axial_level,
276 : const unsigned int last_axial_level,
277 : const unsigned int cross_dimension);
278 : template <class T>
279 : PetscErrorCode populateVectorFromHandle(Vec & x,
280 : const T & solution,
281 : const unsigned int first_axial_level,
282 : const unsigned int last_axial_level,
283 : const unsigned int cross_dimension);
284 : template <class T>
285 : PetscErrorCode populateSolutionChan(const Vec & x,
286 : T & solution,
287 : const unsigned int first_axial_level,
288 : const unsigned int last_axial_level,
289 : const unsigned int cross_dimension);
290 :
291 : //// Matrices and vectors to be used in implicit assembly
292 : /// Mass conservation
293 : /// Mass conservation - sum of cross fluxes
294 : Mat _mc_sumWij_mat;
295 : Vec _Wij_vec;
296 : Vec _prod;
297 : Vec _prodp;
298 : /// Mass conservation - axial convection
299 : Mat _mc_axial_convection_mat;
300 : Vec _mc_axial_convection_rhs;
301 : /// Mass conservation - density time derivative
302 : /// No implicit matrix
303 :
304 : /// Axial momentum
305 : /// Axial momentum conservation - compute turbulent cross fluxes
306 : Mat _amc_turbulent_cross_flows_mat;
307 : Vec _amc_turbulent_cross_flows_rhs;
308 : /// Axial momentum conservation - time derivative
309 : Mat _amc_time_derivative_mat;
310 : Vec _amc_time_derivative_rhs;
311 : /// Axial momentum conservation - advective (Eulerian) derivative
312 : Mat _amc_advective_derivative_mat;
313 : Vec _amc_advective_derivative_rhs;
314 : /// Axial momentum conservation - cross flux derivative
315 : Mat _amc_cross_derivative_mat;
316 : Vec _amc_cross_derivative_rhs;
317 : /// Axial momentum conservation - friction force
318 : Mat _amc_friction_force_mat;
319 : Vec _amc_friction_force_rhs;
320 : /// Axial momentum conservation - buoyancy force
321 : /// No implicit matrix
322 : Vec _amc_gravity_rhs;
323 : /// Axial momentum conservation - pressure force
324 : Mat _amc_pressure_force_mat;
325 : Vec _amc_pressure_force_rhs;
326 : /// Axial momentum system matrix
327 : Mat _amc_sys_mdot_mat;
328 : Vec _amc_sys_mdot_rhs;
329 :
330 : /// Cross momentum
331 : /// Cross momentum conservation - time derivative
332 : Mat _cmc_time_derivative_mat;
333 : Vec _cmc_time_derivative_rhs;
334 : /// Cross momentum conservation - advective (Eulerian) derivative
335 : Mat _cmc_advective_derivative_mat;
336 : Vec _cmc_advective_derivative_rhs;
337 : /// Cross momentum conservation - friction force
338 : Mat _cmc_friction_force_mat;
339 : Vec _cmc_friction_force_rhs;
340 : /// Cross momentum conservation - pressure force
341 : Mat _cmc_pressure_force_mat;
342 : Vec _cmc_pressure_force_rhs;
343 : /// Lateral momentum system matrix
344 : Mat _cmc_sys_Wij_mat;
345 : Vec _cmc_sys_Wij_rhs;
346 :
347 : /// Enthalpy
348 : /// Enthalpy conservation - time derivative
349 : Mat _hc_time_derivative_mat;
350 : Vec _hc_time_derivative_rhs;
351 : /// Enthalpy conservation - advective (Eulerian) derivative;
352 : Mat _hc_advective_derivative_mat;
353 : Vec _hc_advective_derivative_rhs;
354 : /// Enthalpy conservation - cross flux derivative
355 : Mat _hc_cross_derivative_mat;
356 : Vec _hc_cross_derivative_rhs;
357 : /// Enthalpy conservation - source and sink
358 : Vec _hc_added_heat_rhs;
359 : /// System matrices
360 : Mat _hc_sys_h_mat;
361 : Vec _hc_sys_h_rhs;
362 :
363 : /// Added resistances for monolithic convergence
364 : PetscScalar _added_K = 0.0;
365 : PetscScalar _added_K_old = 1000.0;
366 : PetscScalar _max_sumWij;
367 : PetscScalar _max_sumWij_new;
368 : PetscScalar _correction_factor = 1.0;
369 :
370 : public:
371 : static InputParameters validParams();
372 : };
373 :
374 : template <class T>
375 : PetscErrorCode
376 20466 : SubChannel1PhaseProblem::populateDenseFromVector(const Vec & x,
377 : T & loc_solution,
378 : const unsigned int first_axial_level,
379 : const unsigned int last_axial_level,
380 : const unsigned int cross_dimension)
381 : {
382 : PetscScalar * xx;
383 :
384 : PetscFunctionBegin;
385 20466 : LibmeshPetscCall(VecGetArray(x, &xx));
386 344826 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
387 : {
388 324360 : unsigned int iz_ind = iz - first_axial_level;
389 23383656 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
390 : {
391 23059296 : loc_solution(i_l, iz) = xx[iz_ind * cross_dimension + i_l];
392 : }
393 : }
394 20466 : LibmeshPetscCall(VecRestoreArray(x, &xx));
395 :
396 20466 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
397 : }
398 :
399 : template <class T>
400 : PetscErrorCode
401 53262 : SubChannel1PhaseProblem::populateVectorFromHandle(Vec & x,
402 : const T & loc_solution,
403 : const unsigned int first_axial_level,
404 : const unsigned int last_axial_level,
405 : const unsigned int cross_dimension)
406 : {
407 : PetscScalar * xx;
408 :
409 : PetscFunctionBegin;
410 53262 : LibmeshPetscCall(VecGetArray(x, &xx));
411 750582 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
412 : {
413 697320 : unsigned int iz_ind = iz - first_axial_level;
414 32191200 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
415 : {
416 31493880 : auto * loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
417 31493880 : xx[iz_ind * cross_dimension + i_l] = loc_solution(loc_node);
418 : }
419 : }
420 53262 : LibmeshPetscCall(VecRestoreArray(x, &xx));
421 :
422 53262 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
423 : }
424 :
425 : template <class T>
426 : PetscErrorCode
427 20466 : SubChannel1PhaseProblem::populateVectorFromDense(Vec & x,
428 : const T & loc_solution,
429 : const unsigned int first_axial_level,
430 : const unsigned int last_axial_level,
431 : const unsigned int cross_dimension)
432 : {
433 : PetscScalar * xx;
434 : PetscFunctionBegin;
435 20466 : LibmeshPetscCall(VecGetArray(x, &xx));
436 344826 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
437 : {
438 324360 : unsigned int iz_ind = iz - first_axial_level;
439 23383656 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
440 : {
441 23059296 : xx[iz_ind * cross_dimension + i_l] = loc_solution(i_l, iz);
442 : }
443 : }
444 20466 : LibmeshPetscCall(VecRestoreArray(x, &xx));
445 20466 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
446 : }
447 :
448 : template <class T>
449 : PetscErrorCode
450 32796 : SubChannel1PhaseProblem::populateSolutionChan(const Vec & x,
451 : T & loc_solution,
452 : const unsigned int first_axial_level,
453 : const unsigned int last_axial_level,
454 : const unsigned int cross_dimension)
455 : {
456 : PetscScalar * xx;
457 : PetscFunctionBegin;
458 32796 : LibmeshPetscCall(VecGetArray(x, &xx));
459 : Node * loc_node;
460 405756 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
461 : {
462 372960 : unsigned int iz_ind = iz - first_axial_level;
463 16480944 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
464 : {
465 16107984 : loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
466 16107984 : loc_solution.set(loc_node, xx[iz_ind * cross_dimension + i_l]);
467 : }
468 : }
469 32796 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
470 : }
|
