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 3761355 : 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 3360278 : : 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 278 : : 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 16080 : 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 3360000 : 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 3360000 : 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 278 : 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 : /// Whether the time integrator has been checked for consistency with the implementation
226 : bool _time_integrator_checked = false;
227 : /// Time step
228 : Real _dt;
229 : /// Turbulent modeling parameter used in axial momentum equation
230 : Real _CT;
231 : /// Convergence tolerance for the pressure loop in external solve
232 : const Real & _P_tol;
233 : /// Convergence tolerance for the temperature loop in internal solve
234 : const Real & _T_tol;
235 : /// Maximum iterations for the inner temperature loop
236 : const int & _T_maxit;
237 : /// The relative convergence tolerance, (relative decrease) for the ksp linear solver
238 : const PetscReal & _rtol;
239 : /// The absolute convergence tolerance for the ksp linear solver
240 : const PetscReal & _atol;
241 : /// The divergence tolerance for the ksp linear solver
242 : const PetscReal & _dtol;
243 : /// The maximum number of iterations to use for the ksp linear solver
244 : const PetscInt & _maxit;
245 : /// The interpolation method used in constructing the systems
246 : const MooseEnum _interpolation_scheme;
247 : /// The direction of gravity
248 : const MooseEnum _gravity_direction;
249 : const Real _dir_grav;
250 : /// Flag to define the usage of a implicit or explicit solution
251 : const bool _implicit_bool;
252 : /// Flag to define the usage of staggered or collocated pressure
253 : const bool _staggered_pressure_bool;
254 : /// Segregated solve
255 : const bool _segregated_bool;
256 : /// Boolean to printout information related to subchannel solve
257 : const bool _verbose_subchannel;
258 : /// Flag that activates the effect of deformation (pin/duct) based on the auxvalues for displacement, Dpin
259 : bool _deformation = false;
260 : /// Friction closure object
261 : const SCMFrictionClosureBase * _friction_closure;
262 : /// Turbulent Mixing closure object
263 : const SCMMixingClosureBase * _mixing_closure;
264 : /// HTC closure objects
265 : const SCMHTCClosureBase * _pin_HTC_closure;
266 : const SCMHTCClosureBase * _duct_HTC_closure;
267 :
268 : /// Solutions handles and link to TH tables properties
269 : std::unique_ptr<SolutionHandle> _mdot_soln;
270 : std::unique_ptr<SolutionHandle> _SumWij_soln;
271 : std::unique_ptr<SolutionHandle> _P_soln;
272 : std::unique_ptr<SolutionHandle> _DP_soln;
273 : std::unique_ptr<SolutionHandle> _h_soln;
274 : std::unique_ptr<SolutionHandle> _T_soln;
275 : std::unique_ptr<SolutionHandle> _Tpin_soln;
276 : std::unique_ptr<SolutionHandle> _Dpin_soln;
277 : std::unique_ptr<SolutionHandle> _rho_soln;
278 : std::unique_ptr<SolutionHandle> _mu_soln;
279 : std::unique_ptr<SolutionHandle> _S_flow_soln;
280 : std::unique_ptr<SolutionHandle> _w_perim_soln;
281 : std::unique_ptr<SolutionHandle> _q_prime_soln;
282 : std::unique_ptr<SolutionHandle> _duct_heat_flux_soln; // Only used for ducted assemblies
283 : std::unique_ptr<SolutionHandle> _Tduct_soln; // Only used for ducted assemblies
284 : std::unique_ptr<SolutionHandle> _displacement_soln;
285 : std::unique_ptr<SolutionHandle> _ff_soln;
286 : std::unique_ptr<SolutionHandle> _HTC_soln;
287 :
288 : /// Petsc Functions
289 113271 : inline PetscErrorCode createPetscVector(Vec & v, PetscInt n)
290 : {
291 : PetscFunctionBegin;
292 113271 : LibmeshPetscCall(VecCreate(PETSC_COMM_SELF, &v));
293 113271 : LibmeshPetscCall(PetscObjectSetName((PetscObject)v, "Solution"));
294 113271 : LibmeshPetscCall(VecSetSizes(v, PETSC_DECIDE, n));
295 113271 : LibmeshPetscCall(VecSetFromOptions(v));
296 113271 : LibmeshPetscCall(VecZeroEntries(v));
297 113271 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
298 : }
299 :
300 5427 : inline PetscErrorCode createPetscMatrix(Mat & M, PetscInt n, PetscInt m)
301 : {
302 : PetscFunctionBegin;
303 5427 : LibmeshPetscCall(MatCreate(PETSC_COMM_SELF, &M));
304 5427 : LibmeshPetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, n, m));
305 5427 : LibmeshPetscCall(MatSetFromOptions(M));
306 5427 : LibmeshPetscCall(MatSetUp(M));
307 5427 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
308 : }
309 :
310 : template <class T>
311 : PetscErrorCode populateVectorFromDense(Vec & x,
312 : const T & solution,
313 : const unsigned int first_axial_level,
314 : const unsigned int last_axial_level,
315 : const unsigned int cross_dimension);
316 : template <class T>
317 : PetscErrorCode populateDenseFromVector(const Vec & x,
318 : T & solution,
319 : const unsigned int first_axial_level,
320 : const unsigned int last_axial_level,
321 : const unsigned int cross_dimension);
322 : template <class T>
323 : PetscErrorCode populateVectorFromHandle(Vec & x,
324 : const T & solution,
325 : const unsigned int first_axial_level,
326 : const unsigned int last_axial_level,
327 : const unsigned int cross_dimension);
328 : template <class T>
329 : PetscErrorCode populateSolutionChan(const Vec & x,
330 : T & solution,
331 : const unsigned int first_axial_level,
332 : const unsigned int last_axial_level,
333 : const unsigned int cross_dimension);
334 :
335 : //// Matrices and vectors to be used in implicit assembly
336 : /// Mass conservation
337 : /// Mass conservation - sum of cross fluxes
338 : Mat _mc_sumWij_mat;
339 : Vec _Wij_vec;
340 : Vec _prod;
341 : Vec _prodp;
342 : /// Mass conservation - axial convection
343 : Mat _mc_axial_convection_mat;
344 : Vec _mc_axial_convection_rhs;
345 : /// Mass conservation - density time derivative
346 : /// No implicit matrix
347 :
348 : /// Axial momentum
349 : /// Axial momentum conservation - compute turbulent cross fluxes
350 : Mat _amc_turbulent_cross_flows_mat;
351 : Vec _amc_turbulent_cross_flows_rhs;
352 : /// Axial momentum conservation - time derivative
353 : Mat _amc_time_derivative_mat;
354 : Vec _amc_time_derivative_rhs;
355 : /// Axial momentum conservation - advective (Eulerian) derivative
356 : Mat _amc_advective_derivative_mat;
357 : Vec _amc_advective_derivative_rhs;
358 : /// Axial momentum conservation - cross flux derivative
359 : Mat _amc_cross_derivative_mat;
360 : Vec _amc_cross_derivative_rhs;
361 : /// Axial momentum conservation - friction force
362 : Mat _amc_friction_force_mat;
363 : Vec _amc_friction_force_rhs;
364 : /// Axial momentum conservation - buoyancy force
365 : /// No implicit matrix
366 : Vec _amc_gravity_rhs;
367 : /// Axial momentum conservation - pressure force
368 : Mat _amc_pressure_force_mat;
369 : Vec _amc_pressure_force_rhs;
370 : /// Axial momentum system matrix
371 : Mat _amc_sys_mdot_mat;
372 : Vec _amc_sys_mdot_rhs;
373 :
374 : /// Cross momentum
375 : /// Cross momentum conservation - time derivative
376 : Mat _cmc_time_derivative_mat;
377 : Vec _cmc_time_derivative_rhs;
378 : /// Cross momentum conservation - advective (Eulerian) derivative
379 : Mat _cmc_advective_derivative_mat;
380 : Vec _cmc_advective_derivative_rhs;
381 : /// Cross momentum conservation - friction force
382 : Mat _cmc_friction_force_mat;
383 : Vec _cmc_friction_force_rhs;
384 : /// Cross momentum conservation - pressure force
385 : Mat _cmc_pressure_force_mat;
386 : Vec _cmc_pressure_force_rhs;
387 : /// Lateral momentum system matrix
388 : Mat _cmc_sys_Wij_mat;
389 : Vec _cmc_sys_Wij_rhs;
390 :
391 : /// Enthalpy
392 : /// Enthalpy conservation - time derivative
393 : Mat _hc_time_derivative_mat;
394 : Vec _hc_time_derivative_rhs;
395 : /// Enthalpy conservation - advective (Eulerian) derivative;
396 : Mat _hc_advective_derivative_mat;
397 : Vec _hc_advective_derivative_rhs;
398 : /// Enthalpy conservation - cross flux derivative
399 : Mat _hc_cross_derivative_mat;
400 : Vec _hc_cross_derivative_rhs;
401 : /// Enthalpy conservation - source and sink
402 : Vec _hc_added_heat_rhs;
403 : /// System matrices
404 : Mat _hc_sys_h_mat;
405 : Vec _hc_sys_h_rhs;
406 :
407 : /// Added resistances for monolithic convergence
408 : PetscScalar _added_K = 0.0;
409 : PetscScalar _added_K_old = 1000.0;
410 : PetscScalar _max_sumWij;
411 : PetscScalar _max_sumWij_new;
412 : PetscScalar _correction_factor = 1.0;
413 :
414 : public:
415 : static InputParameters validParams();
416 : };
417 :
418 : template <class T>
419 : PetscErrorCode
420 50798 : SubChannel1PhaseProblem::populateDenseFromVector(const Vec & x,
421 : T & loc_solution,
422 : const unsigned int first_axial_level,
423 : const unsigned int last_axial_level,
424 : const unsigned int cross_dimension)
425 : {
426 : PetscScalar * xx;
427 :
428 : PetscFunctionBegin;
429 50798 : LibmeshPetscCall(VecGetArray(x, &xx));
430 643908 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
431 : {
432 593110 : unsigned int iz_ind = iz - first_axial_level;
433 40735870 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
434 : {
435 40142760 : loc_solution(i_l, iz) = xx[iz_ind * cross_dimension + i_l];
436 : }
437 : }
438 50798 : LibmeshPetscCall(VecRestoreArray(x, &xx));
439 :
440 50798 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
441 : }
442 :
443 : template <class T>
444 : PetscErrorCode
445 145546 : SubChannel1PhaseProblem::populateVectorFromHandle(Vec & x,
446 : const T & loc_solution,
447 : const unsigned int first_axial_level,
448 : const unsigned int last_axial_level,
449 : const unsigned int cross_dimension)
450 : {
451 : PetscScalar * xx;
452 :
453 : PetscFunctionBegin;
454 145546 : LibmeshPetscCall(VecGetArray(x, &xx));
455 1715756 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
456 : {
457 1570210 : unsigned int iz_ind = iz - first_axial_level;
458 64998820 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
459 : {
460 63428610 : auto * loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
461 63428610 : xx[iz_ind * cross_dimension + i_l] = loc_solution(loc_node);
462 : }
463 : }
464 145546 : LibmeshPetscCall(VecRestoreArray(x, &xx));
465 :
466 145546 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
467 : }
468 :
469 : template <class T>
470 : PetscErrorCode
471 96460 : SubChannel1PhaseProblem::populateVectorFromDense(Vec & x,
472 : const T & loc_solution,
473 : const unsigned int first_axial_level,
474 : const unsigned int last_axial_level,
475 : const unsigned int cross_dimension)
476 : {
477 : PetscScalar * xx;
478 : PetscFunctionBegin;
479 96460 : LibmeshPetscCall(VecGetArray(x, &xx));
480 1125840 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
481 : {
482 1029380 : unsigned int iz_ind = iz - first_axial_level;
483 67348340 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
484 : {
485 66318960 : xx[iz_ind * cross_dimension + i_l] = loc_solution(i_l, iz);
486 : }
487 : }
488 96460 : LibmeshPetscCall(VecRestoreArray(x, &xx));
489 96460 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
490 : }
491 :
492 : template <class T>
493 : PetscErrorCode
494 94748 : SubChannel1PhaseProblem::populateSolutionChan(const Vec & x,
495 : T & loc_solution,
496 : const unsigned int first_axial_level,
497 : const unsigned int last_axial_level,
498 : const unsigned int cross_dimension)
499 : {
500 : PetscScalar * xx;
501 : PetscFunctionBegin;
502 94748 : LibmeshPetscCall(VecGetArray(x, &xx));
503 : Node * loc_node;
504 1071848 : for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
505 : {
506 977100 : unsigned int iz_ind = iz - first_axial_level;
507 39010800 : for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
508 : {
509 38033700 : loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
510 38033700 : loc_solution.set(loc_node, xx[iz_ind * cross_dimension + i_l]);
511 : }
512 : }
513 94748 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
514 : }
|
