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 "ExternalProblem.h"
13 : #include "PostprocessorInterface.h"
14 : #include "SubChannelApp.h"
15 : #include "QuadSubChannelMesh.h"
16 : #include "SolutionHandle.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 : #include <limits>
24 :
25 : class SinglePhaseFluidProperties;
26 : class SCMFrictionClosureBase;
27 : class SCMHTCClosureBase;
28 : class SCMMixingClosureBase;
29 :
30 : /**
31 : * Base class for the 1-phase steady-state/transient subchannel solver.
32 : */
33 : class SubChannel1PhaseProblem : public ExternalProblem, public PostprocessorInterface
34 : {
35 : public:
36 : SubChannel1PhaseProblem(const InputParameters & params);
37 : virtual ~SubChannel1PhaseProblem();
38 :
39 : virtual void externalSolve() override;
40 : virtual void syncSolutions(Direction direction) override;
41 : virtual bool solverSystemConverged(const unsigned int) override;
42 : virtual void initialSetup() override;
43 :
44 176 : const SCMHTCClosureBase * getDuctHTCClosure() const { return _duct_HTC_closure; }
45 : const SCMHTCClosureBase * getPinHTCClosure() const { return _pin_HTC_closure; } // optional
46 3791355 : const SCMFrictionClosureBase * getFrictionClosure() const { return _friction_closure; }
47 :
48 : /// structure with the needed information to compute the friction factor at a specific subchannel cell
49 : struct FrictionStruct
50 : {
51 : unsigned int i_ch = 0;
52 : Real Re = 1.0;
53 : Real S = 0.0;
54 : Real w_perim = 0.0;
55 :
56 : FrictionStruct() = delete;
57 : FrictionStruct(unsigned int i_ch_, Real Re_, Real S_, Real w_perim_)
58 3390301 : : i_ch(i_ch_), Re(Re_), S(S_), w_perim(w_perim_)
59 : {
60 : }
61 : } _friction_args;
62 :
63 : /// structure with the needed information to compute the Nusselt number at a specific subchannel cell and heated surface
64 : struct NusseltStruct
65 : {
66 : Real Re = 1.0;
67 : Real Pr = 1.0;
68 : unsigned int i_pin = std::numeric_limits<unsigned int>::max(); // sentinel (duct) default
69 : unsigned int iz = 0;
70 : unsigned int i_ch = 0;
71 :
72 : NusseltStruct() = delete;
73 : NusseltStruct(Real Re_, Real Pr_, unsigned int i_pin_, unsigned int iz_, unsigned int i_ch_)
74 301 : : Re(Re_), Pr(Pr_), i_pin(i_pin_), iz(iz_), i_ch(i_ch_)
75 : {
76 : }
77 : } _nusselt_args;
78 :
79 : /// Return the added heat coming from the fuel pins
80 : Real getAddedHeatPin(unsigned int i_ch, unsigned int iz) const
81 : {
82 32160 : return computeAddedHeatPin(i_ch, iz);
83 : }
84 :
85 : /// Return the added heat coming from the duct
86 : Real getAddedHeatDuct(unsigned int i_ch, unsigned int iz) const
87 : {
88 1050 : return computeAddedHeatDuct(i_ch, iz);
89 : }
90 :
91 : /// Outlet pressure postprocessor value
92 : const PostprocessorValue & _P_out;
93 :
94 : /// Get outlet pressure
95 3390000 : const PostprocessorValue & getOutletPressure() const { return _P_out; }
96 :
97 : /// Non-owning pointer to fluid properties user object
98 : const SinglePhaseFluidProperties * _fp;
99 :
100 : /// Get fluid properties object
101 3390000 : const SinglePhaseFluidProperties * getSinglePhaseFluidProperties() const { return _fp; }
102 :
103 : protected:
104 : /// Pure virtual: daughters provide different implementations
105 : virtual Real computeAddedHeatPin(unsigned int i_ch, unsigned int iz) const = 0;
106 :
107 : /// Non-pure: implemented in the base (or override in a child if needed)
108 : virtual Real computeAddedHeatDuct(unsigned int i_ch, unsigned int iz) const;
109 :
110 : /// Computes diversion crossflow per gap for block iblock
111 : void computeWijFromSolve(int iblock);
112 : /// Computes net diversion crossflow per channel for block iblock
113 : void computeSumWij(int iblock);
114 : /// Computes mass flow per channel for block iblock
115 : void computeMdot(int iblock);
116 : /// Computes turbulent crossflow per gap for block iblock
117 : void computeWijPrime(int iblock);
118 : /// Computes and validates the turbulent mixing parameter
119 : Real computeMixingParameter(unsigned int i_gap, unsigned int iz) const;
120 : /// Computes and validates the sweep-flow mixing parameter
121 : Real computeSweepFlowMixingParameter(unsigned int i_gap, unsigned int iz) const;
122 : /// Computes Pressure Drop per channel for block iblock
123 : void computeDP(int iblock);
124 : /// Computes Pressure per channel for block iblock
125 : void computeP(int iblock);
126 : /// Computes Enthalpy per channel for block iblock
127 : virtual void computeh(int iblock) = 0;
128 : /// Computes Temperature per channel for block iblock
129 : void computeT(int iblock);
130 : /// Computes Density per channel for block iblock
131 : void computeRho(int iblock);
132 : /// Computes Viscosity per channel for block iblock
133 : void computeMu(int iblock);
134 : /// Computes Residual Matrix based on the lateral momentum conservation equation for block iblock
135 : void computeWijResidual(int iblock);
136 : /// Function that computes the width of the duct cell that the peripheral subchannel i_ch sees
137 : virtual Real getSubChannelPeripheralDuctWidth(unsigned int i_ch) const = 0;
138 : /// Computes Residual Vector based on the lateral momentum conservation equation for block iblock & updates flow variables based on current crossflow solution
139 : libMesh::DenseVector<Real> residualFunction(int iblock, libMesh::DenseVector<Real> solution);
140 : /// Computes solution of nonlinear equation using snes and provided a residual in a formFunction
141 : PetscErrorCode petscSnesSolver(int iblock,
142 : const libMesh::DenseVector<Real> & solution,
143 : libMesh::DenseVector<Real> & root);
144 : /// This is the residual Vector function in a form compatible with the SNES PETC solvers
145 : friend PetscErrorCode formFunction(SNES snes, Vec x, Vec f, void * ctx);
146 :
147 : /// Computes implicit solve using PetSc
148 : PetscErrorCode implicitPetscSolve(int iblock);
149 :
150 : /// Function to initialize the solution & geometry fields
151 : virtual void initializeSolution() = 0;
152 : /// Detects whether pin diameter or duct displacement fields require geometry recalculation
153 : void detectDeformation();
154 :
155 : /// Functions that computes the interpolation scheme given the Peclet number
156 : PetscScalar computeInterpolationCoefficients(PetscScalar Peclet = 0.0);
157 : PetscScalar
158 : computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet = 0.0);
159 :
160 : /// inline function that is used to define the gravity direction
161 301 : Real computeGravityDir(const MooseEnum & dir) const
162 : {
163 : switch (dir)
164 : {
165 : case 0: // counter_flow
166 : return 1.0;
167 : case 1: // co_flow
168 : return -1.0;
169 : case 2: // none
170 : return 0.0;
171 0 : default:
172 0 : mooseError(name(), ": Invalid gravity direction: expected counter_flow, co_flow, or none");
173 : }
174 : }
175 :
176 : /**
177 : * Solve a linear system (A * x = rhs) with a simple PCJACOBI KSP and populate the
178 : * enthalpy solution into _h_soln for nodes [first_node, last_node].
179 : *
180 : * Uses member tolerances (_rtol, _atol, _dtol, _maxit), mesh (_subchannel_mesh),
181 : * channel count (_n_channels), and error/solution handles (mooseError, _h_soln).
182 : *
183 : * @param A PETSc matrix (operators)
184 : * @param rhs PETSc vector (right-hand side)
185 : * @param first_node inclusive start axial node index
186 : * @param last_node inclusive end axial node index
187 : * @param ksp_prefix options prefix for KSP (e.g. "h_sys_"), may be nullptr
188 : */
189 : PetscErrorCode solveAndPopulateEnthalpy(
190 : Mat A, Vec rhs, unsigned int first_node, unsigned int last_node, const char * ksp_prefix);
191 :
192 : PetscErrorCode cleanUp();
193 : SubChannelMesh & _subchannel_mesh;
194 : /// number of axial blocks
195 : unsigned int _n_blocks;
196 : libMesh::DenseMatrix<Real> _DP;
197 : libMesh::DenseMatrix<Real> & _Wij;
198 : libMesh::DenseMatrix<Real> _Wij_old;
199 : libMesh::DenseMatrix<Real> _WijPrime;
200 : libMesh::DenseMatrix<Real> _Wij_residual_matrix;
201 : const Real _g_grav;
202 : const Real & _kij;
203 : unsigned int _n_cells;
204 : unsigned int _n_gaps;
205 : unsigned int _n_pins;
206 : unsigned int _n_channels;
207 : unsigned int _block_size;
208 : /// axial location of nodes
209 : std::vector<Real> _z_grid;
210 : Real _one;
211 : /// Flag that activates or deactivates the transient parts of the equations we solve by multiplication
212 : Real _TR;
213 : /// Flag that activates or deactivates the calculation of density
214 : const bool _compute_density;
215 : /// Flag that activates or deactivates the calculation of viscosity
216 : const bool _compute_viscosity;
217 : /// Flag that informs if we need to solve the Enthalpy/Temperature equations or not
218 : const bool _compute_power;
219 : /// Flag that informs if there is a pin mesh or not
220 : const bool _pin_mesh_exist;
221 : /// Flag that informs if there is a duct mesh or not
222 : const bool _duct_mesh_exist;
223 : /// Variable that informs whether we exited external solve with a converged solution or not
224 : bool _converged;
225 : /// Time step
226 : Real _dt;
227 : /// Turbulent modeling parameter used in axial momentum equation
228 : Real _CT;
229 : /// Convergence tolerance for the pressure loop in external solve
230 : const Real & _P_tol;
231 : /// Convergence tolerance for the temperature loop in internal solve
232 : const Real & _T_tol;
233 : /// Maximum iterations for the inner temperature loop
234 : const int & _T_maxit;
235 : /// The relative convergence tolerance, (relative decrease) for the ksp linear solver
236 : const PetscReal & _rtol;
237 : /// The absolute convergence tolerance for the ksp linear solver
238 : const PetscReal & _atol;
239 : /// The divergence tolerance for the ksp linear solver
240 : const PetscReal & _dtol;
241 : /// The maximum number of iterations to use for the ksp linear solver
242 : const PetscInt & _maxit;
243 : /// The interpolation method used in constructing the systems
244 : const MooseEnum _interpolation_scheme;
245 : /// The direction of gravity
246 : const MooseEnum _gravity_direction;
247 : const Real _dir_grav;
248 : /// Flag to define the usage of a implicit or explicit solution
249 : const bool _implicit_bool;
250 : /// Flag to define the usage of staggered or collocated pressure
251 : const bool _staggered_pressure_bool;
252 : /// Segregated solve
253 : const bool _segregated_bool;
254 : /// Boolean to printout information related to subchannel solve
255 : const bool _verbose_subchannel;
256 : /// Flag that activates the effect of deformation (pin/duct) based on the auxvalues for displacement, Dpin
257 : bool _deformation = false;
258 : /// Friction closure object
259 : const SCMFrictionClosureBase * _friction_closure;
260 : /// Turbulent Mixing closure object
261 : const SCMMixingClosureBase * _mixing_closure;
262 : /// HTC closure objects
263 : const SCMHTCClosureBase * _pin_HTC_closure;
264 : const SCMHTCClosureBase * _duct_HTC_closure;
265 :
266 : /// Solutions handles and link to TH tables properties
267 : std::unique_ptr<SolutionHandle> _mdot_soln;
268 : std::unique_ptr<SolutionHandle> _SumWij_soln;
269 : std::unique_ptr<SolutionHandle> _P_soln;
270 : std::unique_ptr<SolutionHandle> _DP_soln;
271 : std::unique_ptr<SolutionHandle> _h_soln;
272 : std::unique_ptr<SolutionHandle> _T_soln;
273 : std::unique_ptr<SolutionHandle> _Tpin_soln;
274 : std::unique_ptr<SolutionHandle> _Dpin_soln;
275 : std::unique_ptr<SolutionHandle> _rho_soln;
276 : std::unique_ptr<SolutionHandle> _mu_soln;
277 : std::unique_ptr<SolutionHandle> _S_flow_soln;
278 : std::unique_ptr<SolutionHandle> _w_perim_soln;
279 : std::unique_ptr<SolutionHandle> _q_prime_soln;
280 : std::unique_ptr<SolutionHandle> _duct_heat_flux_soln; // Only used for ducted assemblies
281 : std::unique_ptr<SolutionHandle> _Tduct_soln; // Only used for ducted assemblies
282 : std::unique_ptr<SolutionHandle> _displacement_soln;
283 : std::unique_ptr<SolutionHandle> _ff_soln;
284 : std::unique_ptr<SolutionHandle> _HTC_soln;
285 :
286 : /// Petsc Functions
287 115394 : inline PetscErrorCode createPetscVector(Vec & v, PetscInt n)
288 : {
289 : PetscFunctionBegin;
290 115394 : LibmeshPetscCall(VecCreate(PETSC_COMM_SELF, &v));
291 115394 : LibmeshPetscCall(PetscObjectSetName((PetscObject)v, "Solution"));
292 115394 : LibmeshPetscCall(VecSetSizes(v, PETSC_DECIDE, n));
293 115394 : LibmeshPetscCall(VecSetFromOptions(v));
294 115394 : LibmeshPetscCall(VecZeroEntries(v));
295 115394 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
296 : }
297 :
298 5883 : inline PetscErrorCode createPetscMatrix(Mat & M, PetscInt n, PetscInt m)
299 : {
300 : PetscFunctionBegin;
301 5883 : LibmeshPetscCall(MatCreate(PETSC_COMM_SELF, &M));
302 5883 : LibmeshPetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, n, m));
303 5883 : LibmeshPetscCall(MatSetFromOptions(M));
304 5883 : LibmeshPetscCall(MatSetUp(M));
305 5883 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
306 : }
307 :
308 : template <class T>
309 : PetscErrorCode populateVectorFromDense(Vec & x,
310 : const T & solution,
311 : const unsigned int first_axial_level,
312 : const unsigned int last_axial_level,
313 : const unsigned int cross_dimension);
314 : template <class T>
315 : PetscErrorCode populateDenseFromVector(const Vec & x,
316 : T & solution,
317 : const unsigned int first_axial_level,
318 : const unsigned int last_axial_level,
319 : const unsigned int cross_dimension);
320 : template <class T>
321 : PetscErrorCode populateVectorFromHandle(Vec & x,
322 : const T & solution,
323 : const unsigned int first_axial_level,
324 : const unsigned int last_axial_level,
325 : const unsigned int cross_dimension);
326 : template <class T>
327 : PetscErrorCode populateSolutionChan(const Vec & x,
328 : T & solution,
329 : const unsigned int first_axial_level,
330 : const unsigned int last_axial_level,
331 : const unsigned int cross_dimension);
332 :
333 : //// Matrices and vectors to be used in implicit assembly
334 : /// Mass conservation
335 : /// Mass conservation - sum of cross fluxes
336 : Mat _mc_sumWij_mat;
337 : Vec _Wij_vec;
338 : Vec _prod;
339 : Vec _prodp;
340 : /// Mass conservation - axial convection
341 : Mat _mc_axial_convection_mat;
342 : Vec _mc_axial_convection_rhs;
343 : /// Mass conservation - density time derivative
344 : /// No implicit matrix
345 :
346 : /// Axial momentum
347 : /// Axial momentum conservation - compute turbulent cross fluxes
348 : Mat _amc_turbulent_cross_flows_mat;
349 : Vec _amc_turbulent_cross_flows_rhs;
350 : /// Axial momentum conservation - time derivative
351 : Mat _amc_time_derivative_mat;
352 : Vec _amc_time_derivative_rhs;
353 : /// Axial momentum conservation - advective (Eulerian) derivative
354 : Mat _amc_advective_derivative_mat;
355 : Vec _amc_advective_derivative_rhs;
356 : /// Axial momentum conservation - cross flux derivative
357 : Mat _amc_cross_derivative_mat;
358 : Vec _amc_cross_derivative_rhs;
359 : /// Axial momentum conservation - friction force
360 : Mat _amc_friction_force_mat;
361 : Vec _amc_friction_force_rhs;
362 : /// Axial momentum conservation - buoyancy force
363 : /// No implicit matrix
364 : Vec _amc_gravity_rhs;
365 : /// Axial momentum conservation - pressure force
366 : Mat _amc_pressure_force_mat;
367 : Vec _amc_pressure_force_rhs;
368 : /// Axial momentum system matrix
369 : Mat _amc_sys_mdot_mat;
370 : Vec _amc_sys_mdot_rhs;
371 :
372 : /// Cross momentum
373 : /// Cross momentum conservation - time derivative
374 : Mat _cmc_time_derivative_mat;
375 : Vec _cmc_time_derivative_rhs;
376 : /// Cross momentum conservation - advective (Eulerian) derivative
377 : Mat _cmc_advective_derivative_mat;
378 : Vec _cmc_advective_derivative_rhs;
379 : /// Cross momentum conservation - friction force
380 : Mat _cmc_friction_force_mat;
381 : Vec _cmc_friction_force_rhs;
382 : /// Cross momentum conservation - pressure force
383 : Mat _cmc_pressure_force_mat;
384 : Vec _cmc_pressure_force_rhs;
385 : /// Lateral momentum system matrix
386 : Mat _cmc_sys_Wij_mat;
387 : Vec _cmc_sys_Wij_rhs;
388 :
389 : /// Enthalpy
390 : /// Enthalpy conservation - time derivative
391 : Mat _hc_time_derivative_mat;
392 : Vec _hc_time_derivative_rhs;
393 : /// Enthalpy conservation - advective (Eulerian) derivative;
394 : Mat _hc_advective_derivative_mat;
395 : Vec _hc_advective_derivative_rhs;
396 : /// Enthalpy conservation - cross flux derivative
397 : Mat _hc_cross_derivative_mat;
398 : Vec _hc_cross_derivative_rhs;
399 : /// Enthalpy conservation - source and sink
400 : Vec _hc_added_heat_rhs;
401 : /// System matrices
402 : Mat _hc_sys_h_mat;
403 : Vec _hc_sys_h_rhs;
404 :
405 : /// Added resistances for monolithic convergence
406 : PetscScalar _added_K = 0.0;
407 : PetscScalar _added_K_old = 1000.0;
408 : PetscScalar _max_sumWij;
409 : PetscScalar _max_sumWij_new;
410 : PetscScalar _correction_factor = 1.0;
411 :
412 : public:
413 : static InputParameters validParams();
414 : };
415 :
416 : template <class T>
417 : PetscErrorCode
418 51323 : SubChannel1PhaseProblem::populateDenseFromVector(const Vec & x,
419 : T & loc_solution,
420 : const unsigned int first_axial_level,
421 : const unsigned int last_axial_level,
422 : const unsigned int cross_dimension)
423 : {
424 : PetscScalar * xx;
425 :
426 : PetscFunctionBegin;
427 51323 : LibmeshPetscCall(VecGetArray(x, &xx));
428 661398 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
429 : {
430 610075 : unsigned int iz_ind = iz - first_axial_level;
431 41714575 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
432 : {
433 41104500 : loc_solution(i_l, iz) = xx[iz_ind * cross_dimension + i_l];
434 : }
435 : }
436 51323 : LibmeshPetscCall(VecRestoreArray(x, &xx));
437 :
438 51323 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
439 : }
440 :
441 : template <class T>
442 : PetscErrorCode
443 146421 : SubChannel1PhaseProblem::populateVectorFromHandle(Vec & x,
444 : const T & loc_solution,
445 : const unsigned int first_axial_level,
446 : const unsigned int last_axial_level,
447 : const unsigned int cross_dimension)
448 : {
449 : PetscScalar * xx;
450 :
451 : PetscFunctionBegin;
452 146421 : LibmeshPetscCall(VecGetArray(x, &xx));
453 1744906 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
454 : {
455 1598485 : unsigned int iz_ind = iz - first_axial_level;
456 66150295 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
457 : {
458 64551810 : auto * loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
459 64551810 : xx[iz_ind * cross_dimension + i_l] = loc_solution(loc_node);
460 : }
461 : }
462 146421 : LibmeshPetscCall(VecRestoreArray(x, &xx));
463 :
464 146421 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
465 : }
466 :
467 : template <class T>
468 : PetscErrorCode
469 96985 : SubChannel1PhaseProblem::populateVectorFromDense(Vec & x,
470 : const T & loc_solution,
471 : const unsigned int first_axial_level,
472 : const unsigned int last_axial_level,
473 : const unsigned int cross_dimension)
474 : {
475 : PetscScalar * xx;
476 : PetscFunctionBegin;
477 96985 : LibmeshPetscCall(VecGetArray(x, &xx));
478 1143330 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
479 : {
480 1046345 : unsigned int iz_ind = iz - first_axial_level;
481 68327045 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
482 : {
483 67280700 : xx[iz_ind * cross_dimension + i_l] = loc_solution(i_l, iz);
484 : }
485 : }
486 96985 : LibmeshPetscCall(VecRestoreArray(x, &xx));
487 96985 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
488 : }
489 :
490 : template <class T>
491 : PetscErrorCode
492 95098 : SubChannel1PhaseProblem::populateSolutionChan(const Vec & x,
493 : T & loc_solution,
494 : const unsigned int first_axial_level,
495 : const unsigned int last_axial_level,
496 : const unsigned int cross_dimension)
497 : {
498 : PetscScalar * xx;
499 : PetscFunctionBegin;
500 95098 : LibmeshPetscCall(VecGetArray(x, &xx));
501 : Node * loc_node;
502 1083508 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
503 : {
504 988410 : unsigned int iz_ind = iz - first_axial_level;
505 39471390 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
506 : {
507 38482980 : loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
508 38482980 : loc_solution.set(loc_node, xx[iz_ind * cross_dimension + i_l]);
509 : }
510 : }
511 95098 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
512 : }
|
