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