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 : #include "SubChannel1PhaseProblem.h"
11 : #include "SystemBase.h"
12 : #include "libmesh/petsc_vector.h"
13 : #include "libmesh/dense_matrix.h"
14 : #include "libmesh/dense_vector.h"
15 : #include <iostream>
16 : #include <cmath>
17 : #include "AuxiliarySystem.h"
18 : #include "SCM.h"
19 :
20 : struct Ctx
21 : {
22 : int iblock;
23 : SubChannel1PhaseProblem * schp;
24 : };
25 :
26 : PetscErrorCode
27 40634 : formFunction(SNES, Vec x, Vec f, void * ctx)
28 : {
29 : const PetscScalar * xx;
30 : PetscScalar * ff;
31 : PetscInt size;
32 :
33 : PetscFunctionBegin;
34 : Ctx * cc = static_cast<Ctx *>(ctx);
35 40634 : LibmeshPetscCallQ(VecGetSize(x, &size));
36 :
37 40634 : libMesh::DenseVector<Real> solution_seed(size, 0.0);
38 40634 : LibmeshPetscCallQ(VecGetArrayRead(x, &xx));
39 22994954 : for (PetscInt i = 0; i < size; i++)
40 22954320 : solution_seed(i) = xx[i];
41 :
42 40634 : LibmeshPetscCallQ(VecRestoreArrayRead(x, &xx));
43 :
44 : libMesh::DenseVector<Real> Wij_residual_vector =
45 40634 : cc->schp->residualFunction(cc->iblock, solution_seed);
46 :
47 40634 : LibmeshPetscCallQ(VecGetArray(f, &ff));
48 22994954 : for (int i = 0; i < size; i++)
49 22954320 : ff[i] = Wij_residual_vector(i);
50 :
51 40634 : LibmeshPetscCallQ(VecRestoreArray(f, &ff));
52 40634 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
53 : }
54 :
55 : InputParameters
56 488 : SubChannel1PhaseProblem::validParams()
57 : {
58 976 : MooseEnum schemes("upwind downwind central_difference exponential", "central_difference");
59 488 : InputParameters params = ExternalProblem::validParams();
60 488 : params += PostprocessorInterface::validParams();
61 488 : params.addClassDescription("Base class of the subchannel solvers");
62 976 : params.addRequiredParam<unsigned int>("n_blocks", "The number of blocks in the axial direction");
63 976 : params.addRequiredParam<Real>("CT", "Turbulent modeling parameter");
64 976 : params.addParam<Real>("P_tol", 1e-6, "Pressure tolerance");
65 976 : params.addParam<Real>("T_tol", 1e-6, "Temperature tolerance");
66 976 : params.addParam<int>("T_maxit", 100, "Maximum number of iterations for inner temperature loop");
67 976 : params.addParam<PetscReal>("rtol", 1e-6, "Relative tolerance for ksp solver");
68 976 : params.addParam<PetscReal>("atol", 1e-6, "Absolute tolerance for ksp solver");
69 976 : params.addParam<PetscReal>("dtol", 1e5, "Divergence tolerance or ksp solver");
70 976 : params.addParam<PetscInt>("maxit", 1e4, "Maximum number of iterations for ksp solver");
71 976 : params.addParam<MooseEnum>("interpolation_scheme",
72 : schemes,
73 : "Interpolation scheme used for the method. Default is exponential");
74 976 : params.addParam<bool>(
75 976 : "implicit", false, "Boolean to define the use of explicit or implicit solution.");
76 976 : params.addParam<bool>(
77 976 : "staggered_pressure", false, "Boolean to define the use of explicit or implicit solution.");
78 976 : params.addParam<bool>(
79 976 : "segregated", true, "Boolean to define whether to use a segregated solution.");
80 976 : params.addParam<bool>(
81 976 : "monolithic_thermal", false, "Boolean to define whether to use thermal monolithic solve.");
82 976 : params.addParam<bool>(
83 976 : "verbose_subchannel", false, "Boolean to print out information related to subchannel solve.");
84 976 : params.addParam<bool>(
85 : "deformation",
86 976 : false,
87 : "Boolean that activates the deformation effect based on values for: displacement, Dpin");
88 976 : params.addRequiredParam<bool>("compute_density", "Flag that enables the calculation of density");
89 976 : params.addRequiredParam<bool>("compute_viscosity",
90 : "Flag that enables the calculation of viscosity");
91 976 : params.addRequiredParam<bool>(
92 : "compute_power",
93 : "Flag that informs whether we solve the Enthalpy/Temperature equations or not");
94 976 : params.addRequiredParam<PostprocessorName>(
95 : "P_out", "The postprocessor (or scalar) that provides the value of outlet pressure [Pa]");
96 976 : params.addRequiredParam<UserObjectName>("fp", "Fluid properties user object name");
97 488 : return params;
98 488 : }
99 :
100 244 : SubChannel1PhaseProblem::SubChannel1PhaseProblem(const InputParameters & params)
101 : : ExternalProblem(params),
102 : PostprocessorInterface(this),
103 488 : _subchannel_mesh(SCM::getMesh<SubChannelMesh>(_mesh)),
104 488 : _n_blocks(getParam<unsigned int>("n_blocks")),
105 488 : _Wij(declareRestartableData<libMesh::DenseMatrix<Real>>("Wij")),
106 244 : _g_grav(9.81),
107 244 : _kij(_subchannel_mesh.getKij()),
108 244 : _one(1.0),
109 488 : _compute_density(getParam<bool>("compute_density")),
110 488 : _compute_viscosity(getParam<bool>("compute_viscosity")),
111 488 : _compute_power(getParam<bool>("compute_power")),
112 244 : _pin_mesh_exist(_subchannel_mesh.pinMeshExist()),
113 244 : _duct_mesh_exist(_subchannel_mesh.ductMeshExist()),
114 244 : _P_out(getPostprocessorValue("P_out")),
115 488 : _CT(getParam<Real>("CT")),
116 488 : _P_tol(getParam<Real>("P_tol")),
117 488 : _T_tol(getParam<Real>("T_tol")),
118 488 : _T_maxit(getParam<int>("T_maxit")),
119 488 : _rtol(getParam<PetscReal>("rtol")),
120 488 : _atol(getParam<PetscReal>("atol")),
121 488 : _dtol(getParam<PetscReal>("dtol")),
122 488 : _maxit(getParam<PetscInt>("maxit")),
123 488 : _interpolation_scheme(getParam<MooseEnum>("interpolation_scheme")),
124 488 : _implicit_bool(getParam<bool>("implicit")),
125 488 : _staggered_pressure_bool(getParam<bool>("staggered_pressure")),
126 488 : _segregated_bool(getParam<bool>("segregated")),
127 488 : _monolithic_thermal_bool(getParam<bool>("monolithic_thermal")),
128 488 : _verbose_subchannel(getParam<bool>("verbose_subchannel")),
129 488 : _deformation(getParam<bool>("deformation")),
130 244 : _fp(nullptr),
131 : _Tpin_soln(nullptr),
132 : _q_prime_duct_soln(nullptr),
133 732 : _Tduct_soln(nullptr)
134 : {
135 : // Require only one process
136 244 : if (n_processors() > 1)
137 0 : mooseError("Cannot use more than one MPI process.");
138 244 : _n_cells = _subchannel_mesh.getNumOfAxialCells();
139 244 : _n_gaps = _subchannel_mesh.getNumOfGapsPerLayer();
140 244 : _n_pins = _subchannel_mesh.getNumOfPins();
141 244 : _n_channels = _subchannel_mesh.getNumOfChannels();
142 244 : _z_grid = _subchannel_mesh.getZGrid();
143 244 : _block_size = _n_cells / _n_blocks;
144 : // Turbulent crossflow (stuff that live on the gaps)
145 244 : if (!_app.isRestarting() && !_app.isRecovering())
146 : {
147 238 : _Wij.resize(_n_gaps, _n_cells + 1);
148 238 : _Wij.zero();
149 : }
150 244 : _Wij_old.resize(_n_gaps, _n_cells + 1);
151 : _Wij_old.zero();
152 244 : _WijPrime.resize(_n_gaps, _n_cells + 1);
153 : _WijPrime.zero();
154 244 : _Wij_residual_matrix.resize(_n_gaps, _block_size);
155 : _Wij_residual_matrix.zero();
156 244 : _converged = true;
157 :
158 244 : _outer_channels = 0.0;
159 10106 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
160 : {
161 9862 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
162 9862 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
163 3988 : _outer_channels += 1.0;
164 : }
165 :
166 : // Mass conservation components
167 244 : LibmeshPetscCall(
168 : createPetscMatrix(_mc_sumWij_mat, _block_size * _n_channels, _block_size * _n_gaps));
169 244 : LibmeshPetscCall(createPetscVector(_Wij_vec, _block_size * _n_gaps));
170 244 : LibmeshPetscCall(createPetscVector(_prod, _block_size * _n_channels));
171 244 : LibmeshPetscCall(createPetscVector(_prodp, _block_size * _n_channels));
172 244 : LibmeshPetscCall(createPetscMatrix(
173 : _mc_axial_convection_mat, _block_size * _n_channels, _block_size * _n_channels));
174 244 : LibmeshPetscCall(createPetscVector(_mc_axial_convection_rhs, _block_size * _n_channels));
175 :
176 : // Axial momentum conservation components
177 244 : LibmeshPetscCall(createPetscMatrix(
178 : _amc_turbulent_cross_flows_mat, _block_size * _n_gaps, _block_size * _n_channels));
179 244 : LibmeshPetscCall(createPetscVector(_amc_turbulent_cross_flows_rhs, _block_size * _n_gaps));
180 244 : LibmeshPetscCall(createPetscMatrix(
181 : _amc_time_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
182 244 : LibmeshPetscCall(createPetscVector(_amc_time_derivative_rhs, _block_size * _n_channels));
183 244 : LibmeshPetscCall(createPetscMatrix(
184 : _amc_advective_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
185 244 : LibmeshPetscCall(createPetscVector(_amc_advective_derivative_rhs, _block_size * _n_channels));
186 244 : LibmeshPetscCall(createPetscMatrix(
187 : _amc_cross_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
188 244 : LibmeshPetscCall(createPetscVector(_amc_cross_derivative_rhs, _block_size * _n_channels));
189 244 : LibmeshPetscCall(createPetscMatrix(
190 : _amc_friction_force_mat, _block_size * _n_channels, _block_size * _n_channels));
191 244 : LibmeshPetscCall(createPetscVector(_amc_friction_force_rhs, _block_size * _n_channels));
192 244 : LibmeshPetscCall(createPetscVector(_amc_gravity_rhs, _block_size * _n_channels));
193 244 : LibmeshPetscCall(createPetscMatrix(
194 : _amc_pressure_force_mat, _block_size * _n_channels, _block_size * _n_channels));
195 244 : LibmeshPetscCall(createPetscVector(_amc_pressure_force_rhs, _block_size * _n_channels));
196 244 : LibmeshPetscCall(
197 : createPetscMatrix(_amc_sys_mdot_mat, _block_size * _n_channels, _block_size * _n_channels));
198 244 : LibmeshPetscCall(createPetscVector(_amc_sys_mdot_rhs, _block_size * _n_channels));
199 :
200 : // Lateral momentum conservation components
201 244 : LibmeshPetscCall(
202 : createPetscMatrix(_cmc_time_derivative_mat, _block_size * _n_gaps, _block_size * _n_gaps));
203 244 : LibmeshPetscCall(createPetscVector(_cmc_time_derivative_rhs, _block_size * _n_gaps));
204 244 : LibmeshPetscCall(createPetscMatrix(
205 : _cmc_advective_derivative_mat, _block_size * _n_gaps, _block_size * _n_gaps));
206 244 : LibmeshPetscCall(createPetscVector(_cmc_advective_derivative_rhs, _block_size * _n_gaps));
207 244 : LibmeshPetscCall(
208 : createPetscMatrix(_cmc_friction_force_mat, _block_size * _n_gaps, _block_size * _n_gaps));
209 244 : LibmeshPetscCall(createPetscVector(_cmc_friction_force_rhs, _block_size * _n_gaps));
210 244 : LibmeshPetscCall(
211 : createPetscMatrix(_cmc_pressure_force_mat, _block_size * _n_gaps, _block_size * _n_channels));
212 244 : LibmeshPetscCall(createPetscVector(_cmc_pressure_force_rhs, _block_size * _n_gaps));
213 244 : LibmeshPetscCall(
214 : createPetscMatrix(_cmc_sys_Wij_mat, _block_size * _n_gaps, _block_size * _n_gaps));
215 244 : LibmeshPetscCall(createPetscVector(_cmc_sys_Wij_rhs, _block_size * _n_gaps));
216 :
217 : // Energy conservation components
218 244 : LibmeshPetscCall(createPetscMatrix(
219 : _hc_time_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
220 244 : LibmeshPetscCall(createPetscVector(_hc_time_derivative_rhs, _block_size * _n_channels));
221 244 : LibmeshPetscCall(createPetscMatrix(
222 : _hc_advective_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
223 244 : LibmeshPetscCall(createPetscVector(_hc_advective_derivative_rhs, _block_size * _n_channels));
224 244 : LibmeshPetscCall(createPetscMatrix(
225 : _hc_cross_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
226 244 : LibmeshPetscCall(createPetscVector(_hc_cross_derivative_rhs, _block_size * _n_channels));
227 244 : LibmeshPetscCall(createPetscVector(_hc_added_heat_rhs, _block_size * _n_channels));
228 244 : LibmeshPetscCall(
229 : createPetscMatrix(_hc_sys_h_mat, _block_size * _n_channels, _block_size * _n_channels));
230 244 : LibmeshPetscCall(createPetscVector(_hc_sys_h_rhs, _block_size * _n_channels));
231 :
232 244 : if ((_n_blocks == _n_cells) && _implicit_bool)
233 : {
234 0 : mooseError(name(),
235 : ": When implicit number of blocks can't be equal to number of cells. This will "
236 : "cause problems with the subchannel interpolation scheme.");
237 : }
238 244 : }
239 :
240 : void
241 226 : SubChannel1PhaseProblem::initialSetup()
242 : {
243 226 : ExternalProblem::initialSetup();
244 452 : _fp = &getUserObject<SinglePhaseFluidProperties>(getParam<UserObjectName>("fp"));
245 452 : _mdot_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::MASS_FLOW_RATE));
246 452 : _SumWij_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::SUM_CROSSFLOW));
247 452 : _P_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::PRESSURE));
248 452 : _DP_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::PRESSURE_DROP));
249 452 : _h_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::ENTHALPY));
250 452 : _T_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::TEMPERATURE));
251 226 : if (_pin_mesh_exist)
252 : {
253 188 : _Tpin_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::PIN_TEMPERATURE));
254 188 : _Dpin_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::PIN_DIAMETER));
255 : }
256 452 : _rho_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::DENSITY));
257 452 : _mu_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::VISCOSITY));
258 452 : _S_flow_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::SURFACE_AREA));
259 452 : _w_perim_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::WETTED_PERIMETER));
260 452 : _q_prime_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::LINEAR_HEAT_RATE));
261 : _displacement_soln =
262 452 : std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::DISPLACEMENT));
263 226 : if (_duct_mesh_exist)
264 : {
265 : _q_prime_duct_soln =
266 60 : std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::DUCT_LINEAR_HEAT_RATE));
267 60 : _Tduct_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::DUCT_TEMPERATURE));
268 : }
269 226 : }
270 :
271 488 : SubChannel1PhaseProblem::~SubChannel1PhaseProblem()
272 : {
273 244 : PetscErrorCode ierr = cleanUp();
274 244 : if (ierr)
275 0 : mooseError(name(), ": Error in memory cleanup");
276 244 : }
277 :
278 : PetscErrorCode
279 244 : SubChannel1PhaseProblem::cleanUp()
280 : {
281 : PetscFunctionBegin;
282 : // We need to clean up the petsc matrices/vectors
283 : // Mass conservation components
284 244 : LibmeshPetscCall(MatDestroy(&_mc_sumWij_mat));
285 244 : LibmeshPetscCall(VecDestroy(&_Wij_vec));
286 244 : LibmeshPetscCall(VecDestroy(&_prod));
287 244 : LibmeshPetscCall(VecDestroy(&_prodp));
288 244 : LibmeshPetscCall(MatDestroy(&_mc_axial_convection_mat));
289 244 : LibmeshPetscCall(VecDestroy(&_mc_axial_convection_rhs));
290 :
291 : // Axial momentum conservation components
292 244 : LibmeshPetscCall(MatDestroy(&_amc_turbulent_cross_flows_mat));
293 244 : LibmeshPetscCall(VecDestroy(&_amc_turbulent_cross_flows_rhs));
294 244 : LibmeshPetscCall(MatDestroy(&_amc_time_derivative_mat));
295 244 : LibmeshPetscCall(VecDestroy(&_amc_time_derivative_rhs));
296 244 : LibmeshPetscCall(MatDestroy(&_amc_advective_derivative_mat));
297 244 : LibmeshPetscCall(VecDestroy(&_amc_advective_derivative_rhs));
298 244 : LibmeshPetscCall(MatDestroy(&_amc_cross_derivative_mat));
299 244 : LibmeshPetscCall(VecDestroy(&_amc_cross_derivative_rhs));
300 244 : LibmeshPetscCall(MatDestroy(&_amc_friction_force_mat));
301 244 : LibmeshPetscCall(VecDestroy(&_amc_friction_force_rhs));
302 244 : LibmeshPetscCall(VecDestroy(&_amc_gravity_rhs));
303 244 : LibmeshPetscCall(MatDestroy(&_amc_pressure_force_mat));
304 244 : LibmeshPetscCall(VecDestroy(&_amc_pressure_force_rhs));
305 244 : LibmeshPetscCall(MatDestroy(&_amc_sys_mdot_mat));
306 244 : LibmeshPetscCall(VecDestroy(&_amc_sys_mdot_rhs));
307 :
308 : // Lateral momentum conservation components
309 244 : LibmeshPetscCall(MatDestroy(&_cmc_time_derivative_mat));
310 244 : LibmeshPetscCall(VecDestroy(&_cmc_time_derivative_rhs));
311 244 : LibmeshPetscCall(MatDestroy(&_cmc_advective_derivative_mat));
312 244 : LibmeshPetscCall(VecDestroy(&_cmc_advective_derivative_rhs));
313 244 : LibmeshPetscCall(MatDestroy(&_cmc_friction_force_mat));
314 244 : LibmeshPetscCall(VecDestroy(&_cmc_friction_force_rhs));
315 244 : LibmeshPetscCall(MatDestroy(&_cmc_pressure_force_mat));
316 244 : LibmeshPetscCall(VecDestroy(&_cmc_pressure_force_rhs));
317 244 : LibmeshPetscCall(MatDestroy(&_cmc_sys_Wij_mat));
318 244 : LibmeshPetscCall(VecDestroy(&_cmc_sys_Wij_rhs));
319 :
320 : // Energy conservation components
321 244 : LibmeshPetscCall(MatDestroy(&_hc_time_derivative_mat));
322 244 : LibmeshPetscCall(VecDestroy(&_hc_time_derivative_rhs));
323 244 : LibmeshPetscCall(MatDestroy(&_hc_advective_derivative_mat));
324 244 : LibmeshPetscCall(VecDestroy(&_hc_advective_derivative_rhs));
325 244 : LibmeshPetscCall(MatDestroy(&_hc_cross_derivative_mat));
326 244 : LibmeshPetscCall(VecDestroy(&_hc_cross_derivative_rhs));
327 244 : LibmeshPetscCall(VecDestroy(&_hc_added_heat_rhs));
328 244 : LibmeshPetscCall(MatDestroy(&_hc_sys_h_mat));
329 244 : LibmeshPetscCall(VecDestroy(&_hc_sys_h_rhs));
330 :
331 244 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
332 : }
333 :
334 : bool
335 292 : SubChannel1PhaseProblem::solverSystemConverged(const unsigned int)
336 : {
337 292 : return _converged;
338 : }
339 :
340 : PetscScalar
341 270568548 : SubChannel1PhaseProblem::computeInterpolationCoefficients(PetscScalar Peclet)
342 : {
343 270568548 : switch (_interpolation_scheme)
344 : {
345 : case 0: // upwind interpolation
346 : return 1.0;
347 0 : case 1: // downwind interpolation
348 0 : return 0.0;
349 146501892 : case 2: // central_difference interpolation
350 146501892 : return 0.5;
351 39959568 : case 3: // exponential interpolation (Peclet limited)
352 39959568 : return ((Peclet - 1.0) * std::exp(Peclet) + 1) / (Peclet * (std::exp(Peclet) - 1.) + 1e-10);
353 0 : default:
354 0 : mooseError(name(),
355 : ": Interpolation scheme should be a string: upwind, downwind, central_difference, "
356 : "exponential");
357 : }
358 : }
359 :
360 : PetscScalar
361 229901076 : SubChannel1PhaseProblem::computeInterpolatedValue(PetscScalar topValue,
362 : PetscScalar botValue,
363 : PetscScalar Peclet)
364 : {
365 229901076 : PetscScalar alpha = computeInterpolationCoefficients(Peclet);
366 229901076 : return alpha * botValue + (1.0 - alpha) * topValue;
367 : }
368 :
369 : void
370 376 : SubChannel1PhaseProblem::computeWijFromSolve(int iblock)
371 : {
372 376 : unsigned int last_node = (iblock + 1) * _block_size;
373 376 : unsigned int first_node = iblock * _block_size + 1;
374 : // Initial guess, port crossflow of block (iblock) into a vector that will act as my initial guess
375 376 : libMesh::DenseVector<Real> solution_seed(_n_gaps * _block_size, 0.0);
376 15396 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
377 : {
378 330740 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
379 : {
380 315720 : int i = _n_gaps * (iz - first_node) + i_gap; // column wise transfer
381 315720 : solution_seed(i) = _Wij(i_gap, iz);
382 : }
383 : }
384 :
385 : // Solving the combined lateral momentum equation for Wij using a PETSc solver and update vector
386 : // root
387 376 : libMesh::DenseVector<Real> root(_n_gaps * _block_size, 0.0);
388 376 : LibmeshPetscCall(petscSnesSolver(iblock, solution_seed, root));
389 :
390 : // Assign the solution to the cross-flow matrix
391 : int i = 0;
392 15396 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
393 : {
394 330740 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
395 : {
396 315720 : _Wij(i_gap, iz) = root(i);
397 315720 : i++;
398 : }
399 : }
400 376 : }
401 :
402 : void
403 43328 : SubChannel1PhaseProblem::computeSumWij(int iblock)
404 : {
405 43328 : unsigned int last_node = (iblock + 1) * _block_size;
406 43328 : unsigned int first_node = iblock * _block_size + 1;
407 : // Add to solution vector if explicit
408 43328 : if (!_implicit_bool)
409 : {
410 510658 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
411 : {
412 10569050 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
413 : {
414 10084680 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
415 : Real sumWij = 0.0;
416 : // Calculate sum of crossflow into channel i from channels j around i
417 : unsigned int counter = 0;
418 41899320 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
419 : {
420 31814640 : sumWij += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz);
421 31814640 : counter++;
422 : }
423 : // The net crossflow coming out of cell i [kg/sec]
424 10084680 : _SumWij_soln->set(node_out, sumWij);
425 : }
426 : }
427 : }
428 : // Add to matrix if explicit
429 : else
430 : {
431 225216 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
432 : {
433 208176 : unsigned int iz_ind = iz - first_node;
434 8752920 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
435 : {
436 : // Calculate sum of crossflow into channel i from channels j around i
437 : unsigned int counter = 0;
438 34684488 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
439 : {
440 26139744 : PetscInt row = i_ch + _n_channels * iz_ind;
441 26139744 : PetscInt col = i_gap + _n_gaps * iz_ind;
442 26139744 : PetscScalar value = _subchannel_mesh.getCrossflowSign(i_ch, counter);
443 26139744 : LibmeshPetscCall(MatSetValues(_mc_sumWij_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
444 26139744 : counter++;
445 : }
446 : }
447 : }
448 17040 : LibmeshPetscCall(MatAssemblyBegin(_mc_sumWij_mat, MAT_FINAL_ASSEMBLY));
449 17040 : LibmeshPetscCall(MatAssemblyEnd(_mc_sumWij_mat, MAT_FINAL_ASSEMBLY));
450 17040 : if (_segregated_bool)
451 : {
452 : Vec loc_prod;
453 : Vec loc_Wij;
454 14346 : LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &loc_prod));
455 14346 : LibmeshPetscCall(VecDuplicate(_Wij_vec, &loc_Wij));
456 14346 : LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
457 : loc_Wij, _Wij, first_node, last_node, _n_gaps));
458 14346 : LibmeshPetscCall(MatMult(_mc_sumWij_mat, loc_Wij, loc_prod));
459 14346 : LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
460 : loc_prod, *_SumWij_soln, first_node, last_node, _n_channels));
461 14346 : LibmeshPetscCall(VecDestroy(&loc_prod));
462 14346 : LibmeshPetscCall(VecDestroy(&loc_Wij));
463 : }
464 : }
465 43328 : }
466 :
467 : void
468 43328 : SubChannel1PhaseProblem::computeMdot(int iblock)
469 : {
470 43328 : unsigned int last_node = (iblock + 1) * _block_size;
471 43328 : unsigned int first_node = iblock * _block_size + 1;
472 43328 : if (!_implicit_bool)
473 : {
474 510658 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
475 : {
476 484370 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
477 10569050 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
478 : {
479 10084680 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
480 10084680 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
481 10084680 : auto volume = dz * (*_S_flow_soln)(node_in);
482 10084680 : auto time_term = _TR * ((*_rho_soln)(node_out)-_rho_soln->old(node_out)) * volume / _dt;
483 : // Wij positive out of i into j;
484 10084680 : auto mdot_out = (*_mdot_soln)(node_in) - (*_SumWij_soln)(node_out)-time_term;
485 10084680 : if (mdot_out < 0)
486 : {
487 0 : _console << "Wij = : " << _Wij << "\n";
488 0 : mooseError(name(),
489 : " : Calculation of negative mass flow mdot_out = : ",
490 : mdot_out,
491 : " Axial Level= : ",
492 : iz,
493 : " - Implicit solves are required for recirculating flow.");
494 : }
495 10084680 : _mdot_soln->set(node_out, mdot_out); // kg/sec
496 : }
497 : }
498 : }
499 : else
500 : {
501 225216 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
502 : {
503 208176 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
504 208176 : auto iz_ind = iz - first_node;
505 8752920 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
506 : {
507 8544744 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
508 8544744 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
509 8544744 : auto volume = dz * (*_S_flow_soln)(node_in);
510 :
511 : // Adding time derivative to the RHS
512 8544744 : auto time_term = _TR * ((*_rho_soln)(node_out)-_rho_soln->old(node_out)) * volume / _dt;
513 8544744 : PetscInt row_vec = i_ch + _n_channels * iz_ind;
514 8544744 : PetscScalar value_vec = -1.0 * time_term;
515 8544744 : LibmeshPetscCall(
516 : VecSetValues(_mc_axial_convection_rhs, 1, &row_vec, &value_vec, INSERT_VALUES));
517 :
518 : // Imposing bottom boundary condition or adding of diagonal elements
519 8544744 : if (iz == first_node)
520 : {
521 672012 : PetscScalar value_vec = (*_mdot_soln)(node_in);
522 672012 : PetscInt row_vec = i_ch + _n_channels * iz_ind;
523 672012 : LibmeshPetscCall(
524 : VecSetValues(_mc_axial_convection_rhs, 1, &row_vec, &value_vec, ADD_VALUES));
525 : }
526 : else
527 : {
528 7872732 : PetscInt row = i_ch + _n_channels * iz_ind;
529 7872732 : PetscInt col = i_ch + _n_channels * (iz_ind - 1);
530 7872732 : PetscScalar value = -1.0;
531 7872732 : LibmeshPetscCall(
532 : MatSetValues(_mc_axial_convection_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
533 : }
534 :
535 : // Adding diagonal elements
536 8544744 : PetscInt row = i_ch + _n_channels * iz_ind;
537 8544744 : PetscInt col = i_ch + _n_channels * iz_ind;
538 8544744 : PetscScalar value = 1.0;
539 8544744 : LibmeshPetscCall(
540 : MatSetValues(_mc_axial_convection_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
541 :
542 : // Adding cross flows RHS
543 8544744 : if (_segregated_bool)
544 : {
545 4384260 : PetscScalar value_vec_2 = -1.0 * (*_SumWij_soln)(node_out);
546 4384260 : PetscInt row_vec_2 = i_ch + _n_channels * iz_ind;
547 4384260 : LibmeshPetscCall(
548 : VecSetValues(_mc_axial_convection_rhs, 1, &row_vec_2, &value_vec_2, ADD_VALUES));
549 : }
550 : }
551 : }
552 17040 : LibmeshPetscCall(MatAssemblyBegin(_mc_axial_convection_mat, MAT_FINAL_ASSEMBLY));
553 17040 : LibmeshPetscCall(MatAssemblyEnd(_mc_axial_convection_mat, MAT_FINAL_ASSEMBLY));
554 :
555 17040 : if (_segregated_bool)
556 : {
557 : KSP ksploc;
558 : PC pc;
559 : Vec sol;
560 14346 : LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &sol));
561 14346 : LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
562 14346 : LibmeshPetscCall(KSPSetOperators(ksploc, _mc_axial_convection_mat, _mc_axial_convection_mat));
563 14346 : LibmeshPetscCall(KSPGetPC(ksploc, &pc));
564 14346 : LibmeshPetscCall(PCSetType(pc, PCJACOBI));
565 14346 : LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
566 14346 : LibmeshPetscCall(KSPSetFromOptions(ksploc));
567 14346 : LibmeshPetscCall(KSPSolve(ksploc, _mc_axial_convection_rhs, sol));
568 14346 : LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
569 : sol, *_mdot_soln, first_node, last_node, _n_channels));
570 14346 : LibmeshPetscCall(VecZeroEntries(_mc_axial_convection_rhs));
571 14346 : LibmeshPetscCall(KSPDestroy(&ksploc));
572 14346 : LibmeshPetscCall(VecDestroy(&sol));
573 : }
574 : }
575 43328 : }
576 :
577 : void
578 43328 : SubChannel1PhaseProblem::computeDP(int iblock)
579 : {
580 43328 : unsigned int last_node = (iblock + 1) * _block_size;
581 43328 : unsigned int first_node = iblock * _block_size + 1;
582 43328 : if (!_implicit_bool)
583 : {
584 510658 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
585 : {
586 484370 : auto k_grid = _subchannel_mesh.getKGrid();
587 484370 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
588 10569050 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
589 : {
590 10084680 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
591 10084680 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
592 10084680 : auto rho_in = (*_rho_soln)(node_in);
593 10084680 : auto rho_out = (*_rho_soln)(node_out);
594 10084680 : auto mu_in = (*_mu_soln)(node_in);
595 10084680 : auto S = (*_S_flow_soln)(node_in);
596 10084680 : auto w_perim = (*_w_perim_soln)(node_in);
597 : // hydraulic diameter in the i direction
598 10084680 : auto Dh_i = 4.0 * S / w_perim;
599 10084680 : auto time_term = _TR * ((*_mdot_soln)(node_out)-_mdot_soln->old(node_out)) * dz / _dt -
600 10084680 : dz * 2.0 * (*_mdot_soln)(node_out) * (rho_out - _rho_soln->old(node_out)) /
601 10084680 : rho_in / _dt;
602 : auto mass_term1 =
603 10084680 : std::pow((*_mdot_soln)(node_out), 2.0) * (1.0 / S / rho_out - 1.0 / S / rho_in);
604 10084680 : auto mass_term2 = -2.0 * (*_mdot_soln)(node_out) * (*_SumWij_soln)(node_out) / S / rho_in;
605 : auto crossflow_term = 0.0;
606 : auto turbulent_term = 0.0;
607 : unsigned int counter = 0;
608 41899320 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
609 : {
610 31814640 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
611 : unsigned int ii_ch = chans.first;
612 : unsigned int jj_ch = chans.second;
613 31814640 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
614 31814640 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
615 31814640 : auto * node_out_i = _subchannel_mesh.getChannelNode(ii_ch, iz);
616 31814640 : auto * node_out_j = _subchannel_mesh.getChannelNode(jj_ch, iz);
617 31814640 : auto rho_i = (*_rho_soln)(node_in_i);
618 31814640 : auto rho_j = (*_rho_soln)(node_in_j);
619 31814640 : auto Si = (*_S_flow_soln)(node_in_i);
620 31814640 : auto Sj = (*_S_flow_soln)(node_in_j);
621 : Real u_star = 0.0;
622 : // figure out donor axial velocity
623 31814640 : if (_Wij(i_gap, iz) > 0.0)
624 15051936 : u_star = (*_mdot_soln)(node_out_i) / Si / rho_i;
625 : else
626 16762704 : u_star = (*_mdot_soln)(node_out_j) / Sj / rho_j;
627 :
628 31814640 : crossflow_term +=
629 31814640 : _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * u_star;
630 :
631 31814640 : turbulent_term += _WijPrime(i_gap, iz) * (2 * (*_mdot_soln)(node_out) / rho_in / S -
632 31814640 : (*_mdot_soln)(node_out_j) / Sj / rho_j -
633 31814640 : (*_mdot_soln)(node_out_i) / Si / rho_i);
634 31814640 : counter++;
635 : }
636 10084680 : turbulent_term *= _CT;
637 10084680 : auto Re = (((*_mdot_soln)(node_in) / S) * Dh_i / mu_in);
638 10084680 : _friction_args.Re = Re;
639 10084680 : _friction_args.i_ch = i_ch;
640 10084680 : _friction_args.S = S;
641 10084680 : _friction_args.w_perim = w_perim;
642 10084680 : auto fi = computeFrictionFactor(_friction_args);
643 : /// Upwind local form loss
644 : auto ki = 0.0;
645 10084680 : if ((*_mdot_soln)(node_out) >= 0)
646 10084680 : ki = k_grid[i_ch][iz - 1];
647 : else
648 0 : ki = k_grid[i_ch][iz];
649 10084680 : auto friction_term = (fi * dz / Dh_i + ki) * 0.5 *
650 10084680 : (*_mdot_soln)(node_out)*std::abs((*_mdot_soln)(node_out)) /
651 10084680 : (S * (*_rho_soln)(node_out));
652 10084680 : auto gravity_term = _g_grav * (*_rho_soln)(node_out)*dz * S;
653 10084680 : auto DP = std::pow(S, -1.0) * (time_term + mass_term1 + mass_term2 + crossflow_term +
654 10084680 : turbulent_term + friction_term + gravity_term); // Pa
655 10084680 : _DP_soln->set(node_out, DP);
656 : }
657 484370 : }
658 : }
659 : else
660 : {
661 17040 : LibmeshPetscCall(MatZeroEntries(_amc_time_derivative_mat));
662 17040 : LibmeshPetscCall(MatZeroEntries(_amc_advective_derivative_mat));
663 17040 : LibmeshPetscCall(MatZeroEntries(_amc_cross_derivative_mat));
664 17040 : LibmeshPetscCall(MatZeroEntries(_amc_friction_force_mat));
665 17040 : LibmeshPetscCall(VecZeroEntries(_amc_time_derivative_rhs));
666 17040 : LibmeshPetscCall(VecZeroEntries(_amc_advective_derivative_rhs));
667 17040 : LibmeshPetscCall(VecZeroEntries(_amc_cross_derivative_rhs));
668 17040 : LibmeshPetscCall(VecZeroEntries(_amc_friction_force_rhs));
669 17040 : LibmeshPetscCall(VecZeroEntries(_amc_gravity_rhs));
670 17040 : LibmeshPetscCall(MatZeroEntries(_amc_sys_mdot_mat));
671 17040 : LibmeshPetscCall(VecZeroEntries(_amc_sys_mdot_rhs));
672 225216 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
673 : {
674 208176 : auto k_grid = _subchannel_mesh.getKGrid();
675 208176 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
676 208176 : auto iz_ind = iz - first_node;
677 8752920 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
678 : {
679 : // inlet and outlet nodes
680 8544744 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
681 8544744 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
682 :
683 : // interpolation weight coefficient
684 : PetscScalar Pe = 0.5;
685 8544744 : if (_interpolation_scheme == 3)
686 : {
687 : // Compute the Peclet number
688 977544 : auto S_in = (*_S_flow_soln)(node_in);
689 977544 : auto S_out = (*_S_flow_soln)(node_out);
690 977544 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
691 977544 : auto w_perim_in = (*_w_perim_soln)(node_in);
692 977544 : auto w_perim_out = (*_w_perim_soln)(node_out);
693 977544 : auto w_perim_interp = this->computeInterpolatedValue(w_perim_out, w_perim_in, 0.5);
694 : auto mdot_loc =
695 977544 : this->computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), 0.5);
696 977544 : auto mu_in = (*_mu_soln)(node_in);
697 977544 : auto mu_out = (*_mu_soln)(node_out);
698 977544 : auto mu_interp = this->computeInterpolatedValue(mu_out, mu_in, 0.5);
699 977544 : auto Dh_i = 4.0 * S_interp / w_perim_interp;
700 : // Compute friction
701 977544 : auto Re = ((mdot_loc / S_interp) * Dh_i / mu_interp);
702 977544 : _friction_args.Re = Re;
703 977544 : _friction_args.i_ch = i_ch;
704 977544 : _friction_args.S = S_interp;
705 977544 : _friction_args.w_perim = w_perim_interp;
706 977544 : auto fi = computeFrictionFactor(_friction_args);
707 : /// Upwind local form loss
708 : auto ki = 0.0;
709 977544 : if ((*_mdot_soln)(node_out) >= 0)
710 977544 : ki = k_grid[i_ch][iz - 1];
711 : else
712 0 : ki = k_grid[i_ch][iz];
713 977544 : Pe = 1.0 / ((fi * dz / Dh_i + ki) * 0.5) * mdot_loc / std::abs(mdot_loc);
714 : }
715 8544744 : auto alpha = computeInterpolationCoefficients(Pe);
716 :
717 : // inlet, outlet, and interpolated density
718 8544744 : auto rho_in = (*_rho_soln)(node_in);
719 8544744 : auto rho_out = (*_rho_soln)(node_out);
720 8544744 : auto rho_interp = computeInterpolatedValue(rho_out, rho_in, Pe);
721 :
722 : // inlet, outlet, and interpolated viscosity
723 8544744 : auto mu_in = (*_mu_soln)(node_in);
724 8544744 : auto mu_out = (*_mu_soln)(node_out);
725 8544744 : auto mu_interp = computeInterpolatedValue(mu_out, mu_in, Pe);
726 :
727 : // inlet, outlet, and interpolated axial surface area
728 8544744 : auto S_in = (*_S_flow_soln)(node_in);
729 8544744 : auto S_out = (*_S_flow_soln)(node_out);
730 8544744 : auto S_interp = computeInterpolatedValue(S_out, S_in, Pe);
731 :
732 : // inlet, outlet, and interpolated wetted perimeter
733 8544744 : auto w_perim_in = (*_w_perim_soln)(node_in);
734 8544744 : auto w_perim_out = (*_w_perim_soln)(node_out);
735 8544744 : auto w_perim_interp = computeInterpolatedValue(w_perim_out, w_perim_in, Pe);
736 :
737 : // hydraulic diameter in the i direction
738 8544744 : auto Dh_i = 4.0 * S_interp / w_perim_interp;
739 :
740 : /// Time derivative term
741 8544744 : if (iz == first_node)
742 : {
743 672012 : PetscScalar value_vec_tt = -1.0 * _TR * alpha * (*_mdot_soln)(node_in)*dz / _dt;
744 672012 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
745 672012 : LibmeshPetscCall(
746 : VecSetValues(_amc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
747 : }
748 : else
749 : {
750 7872732 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
751 7872732 : PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
752 7872732 : PetscScalar value_tt = _TR * alpha * dz / _dt;
753 7872732 : LibmeshPetscCall(MatSetValues(
754 : _amc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
755 : }
756 :
757 : // Adding diagonal elements
758 8544744 : PetscInt row_tt = i_ch + _n_channels * iz_ind;
759 8544744 : PetscInt col_tt = i_ch + _n_channels * iz_ind;
760 8544744 : PetscScalar value_tt = _TR * (1.0 - alpha) * dz / _dt;
761 8544744 : LibmeshPetscCall(MatSetValues(
762 : _amc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
763 :
764 : // Adding RHS elements
765 : PetscScalar mdot_old_interp =
766 8544744 : computeInterpolatedValue(_mdot_soln->old(node_out), _mdot_soln->old(node_in), Pe);
767 8544744 : PetscScalar value_vec_tt = _TR * mdot_old_interp * dz / _dt;
768 8544744 : PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
769 8544744 : LibmeshPetscCall(
770 : VecSetValues(_amc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
771 :
772 : /// Advective derivative term
773 8544744 : if (iz == first_node)
774 : {
775 672012 : PetscScalar value_vec_at = std::pow((*_mdot_soln)(node_in), 2.0) / (S_in * rho_in);
776 672012 : PetscInt row_vec_at = i_ch + _n_channels * iz_ind;
777 672012 : LibmeshPetscCall(VecSetValues(
778 : _amc_advective_derivative_rhs, 1, &row_vec_at, &value_vec_at, ADD_VALUES));
779 : }
780 : else
781 : {
782 7872732 : PetscInt row_at = i_ch + _n_channels * iz_ind;
783 7872732 : PetscInt col_at = i_ch + _n_channels * (iz_ind - 1);
784 7872732 : PetscScalar value_at = -1.0 * std::abs((*_mdot_soln)(node_in)) / (S_in * rho_in);
785 7872732 : LibmeshPetscCall(MatSetValues(
786 : _amc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
787 : }
788 :
789 : // Adding diagonal elements
790 8544744 : PetscInt row_at = i_ch + _n_channels * iz_ind;
791 8544744 : PetscInt col_at = i_ch + _n_channels * iz_ind;
792 8544744 : PetscScalar value_at = std::abs((*_mdot_soln)(node_out)) / (S_out * rho_out);
793 8544744 : LibmeshPetscCall(MatSetValues(
794 : _amc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
795 :
796 : /// Cross derivative term
797 : unsigned int counter = 0;
798 : unsigned int cross_index = iz; // iz-1;
799 34684488 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
800 : {
801 26139744 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
802 : unsigned int ii_ch = chans.first;
803 : unsigned int jj_ch = chans.second;
804 26139744 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
805 26139744 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
806 26139744 : auto * node_out_i = _subchannel_mesh.getChannelNode(ii_ch, iz);
807 26139744 : auto * node_out_j = _subchannel_mesh.getChannelNode(jj_ch, iz);
808 : auto rho_i =
809 26139744 : computeInterpolatedValue((*_rho_soln)(node_out_i), (*_rho_soln)(node_in_i), Pe);
810 : auto rho_j =
811 26139744 : computeInterpolatedValue((*_rho_soln)(node_out_j), (*_rho_soln)(node_in_j), Pe);
812 : auto S_i =
813 26139744 : computeInterpolatedValue((*_S_flow_soln)(node_out_i), (*_S_flow_soln)(node_in_i), Pe);
814 : auto S_j =
815 26139744 : computeInterpolatedValue((*_S_flow_soln)(node_out_j), (*_S_flow_soln)(node_in_j), Pe);
816 : auto u_star = 0.0;
817 : // figure out donor axial velocity
818 26139744 : if (_Wij(i_gap, cross_index) > 0.0)
819 : {
820 14823540 : if (iz == first_node)
821 : {
822 1233240 : u_star = (*_mdot_soln)(node_in_i) / S_i / rho_i;
823 2466480 : PetscScalar value_vec_ct = -1.0 * alpha *
824 1233240 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
825 1233240 : _Wij(i_gap, cross_index) * u_star;
826 1233240 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
827 1233240 : LibmeshPetscCall(VecSetValues(
828 : _amc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
829 : }
830 : else
831 : {
832 13590300 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
833 13590300 : _Wij(i_gap, cross_index) / S_i / rho_i;
834 13590300 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
835 13590300 : PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
836 13590300 : LibmeshPetscCall(MatSetValues(
837 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
838 : }
839 14823540 : PetscScalar value_ct = (1.0 - alpha) *
840 14823540 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
841 14823540 : _Wij(i_gap, cross_index) / S_i / rho_i;
842 14823540 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
843 14823540 : PetscInt col_ct = ii_ch + _n_channels * iz_ind;
844 14823540 : LibmeshPetscCall(MatSetValues(
845 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
846 : }
847 11316204 : else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
848 : {
849 10378680 : if (iz == first_node)
850 : {
851 825144 : u_star = (*_mdot_soln)(node_in_j) / S_j / rho_j;
852 1650288 : PetscScalar value_vec_ct = -1.0 * alpha *
853 825144 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
854 825144 : _Wij(i_gap, cross_index) * u_star;
855 825144 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
856 825144 : LibmeshPetscCall(VecSetValues(
857 : _amc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
858 : }
859 : else
860 : {
861 9553536 : PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
862 9553536 : _Wij(i_gap, cross_index) / S_j / rho_j;
863 9553536 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
864 9553536 : PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
865 9553536 : LibmeshPetscCall(MatSetValues(
866 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
867 : }
868 10378680 : PetscScalar value_ct = (1.0 - alpha) *
869 10378680 : _subchannel_mesh.getCrossflowSign(i_ch, counter) *
870 10378680 : _Wij(i_gap, cross_index) / S_j / rho_j;
871 10378680 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
872 10378680 : PetscInt col_ct = jj_ch + _n_channels * iz_ind;
873 10378680 : LibmeshPetscCall(MatSetValues(
874 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
875 : }
876 :
877 26139744 : if (iz == first_node)
878 : {
879 2083536 : PetscScalar value_vec_ct = -2.0 * alpha * (*_mdot_soln)(node_in)*_CT *
880 2083536 : _WijPrime(i_gap, cross_index) / (rho_interp * S_interp);
881 2083536 : value_vec_ct += alpha * (*_mdot_soln)(node_in_j)*_CT * _WijPrime(i_gap, cross_index) /
882 2083536 : (rho_j * S_j);
883 2083536 : value_vec_ct += alpha * (*_mdot_soln)(node_in_i)*_CT * _WijPrime(i_gap, cross_index) /
884 2083536 : (rho_i * S_i);
885 2083536 : PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
886 2083536 : LibmeshPetscCall(
887 : VecSetValues(_amc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
888 : }
889 : else
890 : {
891 : PetscScalar value_center_ct =
892 24056208 : 2.0 * alpha * _CT * _WijPrime(i_gap, cross_index) / (rho_interp * S_interp);
893 24056208 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
894 24056208 : PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
895 24056208 : LibmeshPetscCall(MatSetValues(
896 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
897 :
898 : PetscScalar value_left_ct =
899 24056208 : -1.0 * alpha * _CT * _WijPrime(i_gap, cross_index) / (rho_j * S_j);
900 24056208 : row_ct = i_ch + _n_channels * iz_ind;
901 24056208 : col_ct = jj_ch + _n_channels * (iz_ind - 1);
902 24056208 : LibmeshPetscCall(MatSetValues(
903 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
904 :
905 : PetscScalar value_right_ct =
906 24056208 : -1.0 * alpha * _CT * _WijPrime(i_gap, cross_index) / (rho_i * S_i);
907 24056208 : row_ct = i_ch + _n_channels * iz_ind;
908 24056208 : col_ct = ii_ch + _n_channels * (iz_ind - 1);
909 24056208 : LibmeshPetscCall(MatSetValues(
910 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
911 : }
912 :
913 : PetscScalar value_center_ct =
914 26139744 : 2.0 * (1.0 - alpha) * _CT * _WijPrime(i_gap, cross_index) / (rho_interp * S_interp);
915 26139744 : PetscInt row_ct = i_ch + _n_channels * iz_ind;
916 26139744 : PetscInt col_ct = i_ch + _n_channels * iz_ind;
917 26139744 : LibmeshPetscCall(MatSetValues(
918 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
919 :
920 : PetscScalar value_left_ct =
921 26139744 : -1.0 * (1.0 - alpha) * _CT * _WijPrime(i_gap, cross_index) / (rho_j * S_j);
922 26139744 : row_ct = i_ch + _n_channels * iz_ind;
923 26139744 : col_ct = jj_ch + _n_channels * iz_ind;
924 26139744 : LibmeshPetscCall(MatSetValues(
925 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
926 :
927 : PetscScalar value_right_ct =
928 26139744 : -1.0 * (1.0 - alpha) * _CT * _WijPrime(i_gap, cross_index) / (rho_i * S_i);
929 26139744 : row_ct = i_ch + _n_channels * iz_ind;
930 26139744 : col_ct = ii_ch + _n_channels * iz_ind;
931 26139744 : LibmeshPetscCall(MatSetValues(
932 : _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
933 26139744 : counter++;
934 : }
935 :
936 : /// Friction term
937 : PetscScalar mdot_interp =
938 8544744 : computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), Pe);
939 8544744 : auto Re = ((mdot_interp / S_interp) * Dh_i / mu_interp);
940 8544744 : _friction_args.Re = Re;
941 8544744 : _friction_args.i_ch = i_ch;
942 8544744 : _friction_args.S = S_interp;
943 8544744 : _friction_args.w_perim = w_perim_interp;
944 8544744 : auto fi = computeFrictionFactor(_friction_args);
945 : /// Upwind local form loss
946 : auto ki = 0.0;
947 8544744 : if ((*_mdot_soln)(node_out) >= 0)
948 8544744 : ki = k_grid[i_ch][iz - 1];
949 : else
950 0 : ki = k_grid[i_ch][iz];
951 8544744 : auto coef = (fi * dz / Dh_i + ki) * 0.5 * std::abs((*_mdot_soln)(node_out)) /
952 8544744 : (S_interp * rho_interp);
953 8544744 : if (iz == first_node)
954 : {
955 672012 : PetscScalar value_vec = -1.0 * alpha * coef * (*_mdot_soln)(node_in);
956 672012 : PetscInt row_vec = i_ch + _n_channels * iz_ind;
957 672012 : LibmeshPetscCall(
958 : VecSetValues(_amc_friction_force_rhs, 1, &row_vec, &value_vec, ADD_VALUES));
959 : }
960 : else
961 : {
962 7872732 : PetscInt row = i_ch + _n_channels * iz_ind;
963 7872732 : PetscInt col = i_ch + _n_channels * (iz_ind - 1);
964 7872732 : PetscScalar value = alpha * coef;
965 7872732 : LibmeshPetscCall(
966 : MatSetValues(_amc_friction_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
967 : }
968 :
969 : // Adding diagonal elements
970 8544744 : PetscInt row = i_ch + _n_channels * iz_ind;
971 8544744 : PetscInt col = i_ch + _n_channels * iz_ind;
972 8544744 : PetscScalar value = (1.0 - alpha) * coef;
973 8544744 : LibmeshPetscCall(
974 : MatSetValues(_amc_friction_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
975 :
976 : /// Gravity force
977 8544744 : PetscScalar value_vec = -1.0 * _g_grav * rho_interp * dz * S_interp;
978 8544744 : PetscInt row_vec = i_ch + _n_channels * iz_ind;
979 8544744 : LibmeshPetscCall(VecSetValues(_amc_gravity_rhs, 1, &row_vec, &value_vec, ADD_VALUES));
980 : }
981 208176 : }
982 : /// Assembling system
983 17040 : LibmeshPetscCall(MatZeroEntries(_amc_sys_mdot_mat));
984 17040 : LibmeshPetscCall(VecZeroEntries(_amc_sys_mdot_rhs));
985 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
986 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
987 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
988 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
989 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
990 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
991 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_friction_force_mat, MAT_FINAL_ASSEMBLY));
992 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_friction_force_mat, MAT_FINAL_ASSEMBLY));
993 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
994 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
995 : // Matrix
996 : #if !PETSC_VERSION_LESS_THAN(3, 15, 0)
997 17040 : LibmeshPetscCall(
998 : MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
999 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1000 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1001 17040 : LibmeshPetscCall(
1002 : MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
1003 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1004 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1005 17040 : LibmeshPetscCall(
1006 : MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_cross_derivative_mat, UNKNOWN_NONZERO_PATTERN));
1007 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1008 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1009 17040 : LibmeshPetscCall(
1010 : MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_friction_force_mat, UNKNOWN_NONZERO_PATTERN));
1011 : #else
1012 : LibmeshPetscCall(
1013 : MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1014 : LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1015 : LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1016 : LibmeshPetscCall(
1017 : MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1018 : LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1019 : LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1020 : LibmeshPetscCall(
1021 : MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1022 : LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1023 : LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1024 : LibmeshPetscCall(
1025 : MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_friction_force_mat, DIFFERENT_NONZERO_PATTERN));
1026 : #endif
1027 17040 : LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1028 17040 : LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
1029 : // RHS
1030 17040 : LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_time_derivative_rhs));
1031 17040 : LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_advective_derivative_rhs));
1032 17040 : LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_cross_derivative_rhs));
1033 17040 : LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_friction_force_rhs));
1034 17040 : LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_gravity_rhs));
1035 17040 : if (_segregated_bool)
1036 : {
1037 : // Assembly the matrix system
1038 14346 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
1039 : _prod, *_mdot_soln, first_node, last_node, _n_channels));
1040 : Vec ls;
1041 14346 : LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &ls));
1042 14346 : LibmeshPetscCall(MatMult(_amc_sys_mdot_mat, _prod, ls));
1043 14346 : LibmeshPetscCall(VecAXPY(ls, -1.0, _amc_sys_mdot_rhs));
1044 : PetscScalar * xx;
1045 14346 : LibmeshPetscCall(VecGetArray(ls, &xx));
1046 131796 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1047 : {
1048 117450 : auto iz_ind = iz - first_node;
1049 4501710 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1050 : {
1051 : // Setting nodes
1052 4384260 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1053 4384260 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
1054 :
1055 : // inlet, outlet, and interpolated axial surface area
1056 4384260 : auto S_in = (*_S_flow_soln)(node_in);
1057 4384260 : auto S_out = (*_S_flow_soln)(node_out);
1058 4384260 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
1059 :
1060 : // Setting solutions
1061 4384260 : if (S_interp != 0)
1062 : {
1063 4384260 : auto DP = std::pow(S_interp, -1.0) * xx[iz_ind * _n_channels + i_ch];
1064 4384260 : _DP_soln->set(node_out, DP);
1065 : }
1066 : else
1067 : {
1068 : auto DP = 0.0;
1069 0 : _DP_soln->set(node_out, DP);
1070 : }
1071 : }
1072 : }
1073 14346 : LibmeshPetscCall(VecDestroy(&ls));
1074 : }
1075 : }
1076 43328 : }
1077 :
1078 : void
1079 43328 : SubChannel1PhaseProblem::computeP(int iblock)
1080 : {
1081 43328 : unsigned int last_node = (iblock + 1) * _block_size;
1082 43328 : unsigned int first_node = iblock * _block_size + 1;
1083 43328 : if (!_implicit_bool)
1084 : {
1085 26288 : if (!_staggered_pressure_bool)
1086 : {
1087 454030 : for (unsigned int iz = last_node; iz > first_node - 1; iz--)
1088 : {
1089 : // Calculate pressure in the inlet of the cell assuming known outlet
1090 8664290 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1091 : {
1092 8231400 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
1093 8231400 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1094 : // update Pressure solution
1095 8231400 : _P_soln->set(node_in, (*_P_soln)(node_out) + (*_DP_soln)(node_out));
1096 : }
1097 : }
1098 : }
1099 : else
1100 : {
1101 56628 : for (unsigned int iz = last_node; iz > first_node - 1; iz--)
1102 : {
1103 : // Calculate pressure in the inlet of the cell assuming known outlet
1104 1904760 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1105 : {
1106 1853280 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
1107 1853280 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1108 : // update Pressure solution
1109 : // Note: assuming uniform axial discretization in the curren code
1110 : // We will need to update this later if we allow non-uniform refinements in the axial
1111 : // direction
1112 : PetscScalar Pe = 0.5;
1113 1853280 : auto alpha = computeInterpolationCoefficients(Pe);
1114 1853280 : if (iz == last_node)
1115 : {
1116 185328 : _P_soln->set(node_in, (*_P_soln)(node_out) + (*_DP_soln)(node_out) / 2.0);
1117 : }
1118 : else
1119 : {
1120 1667952 : _P_soln->set(node_in,
1121 1667952 : (*_P_soln)(node_out) + (1.0 - alpha) * (*_DP_soln)(node_out) +
1122 1667952 : alpha * (*_DP_soln)(node_in));
1123 : }
1124 : }
1125 : }
1126 : }
1127 : }
1128 : else
1129 : {
1130 17040 : if (!_staggered_pressure_bool)
1131 : {
1132 16980 : LibmeshPetscCall(VecZeroEntries(_amc_pressure_force_rhs));
1133 224556 : for (unsigned int iz = last_node; iz > first_node - 1; iz--)
1134 : {
1135 207576 : auto iz_ind = iz - first_node;
1136 : // Calculate pressure in the inlet of the cell assuming known outlet
1137 8730720 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1138 : {
1139 8523144 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
1140 8523144 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1141 :
1142 : // inlet, outlet, and interpolated axial surface area
1143 8523144 : auto S_in = (*_S_flow_soln)(node_in);
1144 8523144 : auto S_out = (*_S_flow_soln)(node_out);
1145 8523144 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
1146 :
1147 : // Creating matrix of coefficients
1148 8523144 : PetscInt row = i_ch + _n_channels * iz_ind;
1149 8523144 : PetscInt col = i_ch + _n_channels * iz_ind;
1150 8523144 : PetscScalar value = -1.0 * S_interp;
1151 8523144 : LibmeshPetscCall(
1152 : MatSetValues(_amc_pressure_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
1153 :
1154 8523144 : if (iz == last_node)
1155 : {
1156 669852 : PetscScalar value = -1.0 * (*_P_soln)(node_out)*S_interp;
1157 669852 : PetscInt row = i_ch + _n_channels * iz_ind;
1158 669852 : LibmeshPetscCall(VecSetValues(_amc_pressure_force_rhs, 1, &row, &value, ADD_VALUES));
1159 : }
1160 : else
1161 : {
1162 7853292 : PetscInt row = i_ch + _n_channels * iz_ind;
1163 7853292 : PetscInt col = i_ch + _n_channels * (iz_ind + 1);
1164 7853292 : PetscScalar value = 1.0 * S_interp;
1165 7853292 : LibmeshPetscCall(
1166 : MatSetValues(_amc_pressure_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
1167 : }
1168 :
1169 8523144 : if (_segregated_bool)
1170 : {
1171 4384260 : auto dp_out = (*_DP_soln)(node_out);
1172 4384260 : PetscScalar value_v = -1.0 * dp_out * S_interp;
1173 4384260 : PetscInt row_v = i_ch + _n_channels * iz_ind;
1174 4384260 : LibmeshPetscCall(
1175 : VecSetValues(_amc_pressure_force_rhs, 1, &row_v, &value_v, ADD_VALUES));
1176 : }
1177 : }
1178 : }
1179 : // Solving pressure problem
1180 16980 : LibmeshPetscCall(MatAssemblyBegin(_amc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
1181 16980 : LibmeshPetscCall(MatAssemblyEnd(_amc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
1182 16980 : if (_segregated_bool)
1183 : {
1184 : KSP ksploc;
1185 : PC pc;
1186 : Vec sol;
1187 14346 : LibmeshPetscCall(VecDuplicate(_amc_pressure_force_rhs, &sol));
1188 14346 : LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
1189 14346 : LibmeshPetscCall(KSPSetOperators(ksploc, _amc_pressure_force_mat, _amc_pressure_force_mat));
1190 14346 : LibmeshPetscCall(KSPGetPC(ksploc, &pc));
1191 14346 : LibmeshPetscCall(PCSetType(pc, PCJACOBI));
1192 14346 : LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
1193 14346 : LibmeshPetscCall(KSPSetFromOptions(ksploc));
1194 14346 : LibmeshPetscCall(KSPSolve(ksploc, _amc_pressure_force_rhs, sol));
1195 : PetscScalar * xx;
1196 14346 : LibmeshPetscCall(VecGetArray(sol, &xx));
1197 : // update Pressure solution
1198 131796 : for (unsigned int iz = last_node; iz > first_node - 1; iz--)
1199 : {
1200 117450 : auto iz_ind = iz - first_node;
1201 4501710 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1202 : {
1203 4384260 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1204 4384260 : PetscScalar value = xx[iz_ind * _n_channels + i_ch];
1205 4384260 : _P_soln->set(node_in, value);
1206 : }
1207 : }
1208 14346 : LibmeshPetscCall(VecZeroEntries(_amc_pressure_force_rhs));
1209 14346 : LibmeshPetscCall(KSPDestroy(&ksploc));
1210 14346 : LibmeshPetscCall(VecDestroy(&sol));
1211 : }
1212 : }
1213 : else
1214 : {
1215 60 : LibmeshPetscCall(VecZeroEntries(_amc_pressure_force_rhs));
1216 660 : for (unsigned int iz = last_node; iz > first_node - 1; iz--)
1217 : {
1218 600 : auto iz_ind = iz - first_node;
1219 : // Calculate pressure in the inlet of the cell assuming known outlet
1220 22200 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1221 : {
1222 21600 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
1223 21600 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1224 :
1225 : // inlet, outlet, and interpolated axial surface area
1226 21600 : auto S_in = (*_S_flow_soln)(node_in);
1227 21600 : auto S_out = (*_S_flow_soln)(node_out);
1228 21600 : auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
1229 :
1230 : // Creating matrix of coefficients
1231 21600 : PetscInt row = i_ch + _n_channels * iz_ind;
1232 21600 : PetscInt col = i_ch + _n_channels * iz_ind;
1233 21600 : PetscScalar value = -1.0 * S_interp;
1234 21600 : LibmeshPetscCall(
1235 : MatSetValues(_amc_pressure_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
1236 :
1237 21600 : if (iz == last_node)
1238 : {
1239 2160 : PetscScalar value = -1.0 * (*_P_soln)(node_out)*S_interp;
1240 2160 : PetscInt row = i_ch + _n_channels * iz_ind;
1241 2160 : LibmeshPetscCall(VecSetValues(_amc_pressure_force_rhs, 1, &row, &value, ADD_VALUES));
1242 :
1243 2160 : auto dp_out = (*_DP_soln)(node_out);
1244 2160 : PetscScalar value_v = -1.0 * dp_out / 2.0 * S_interp;
1245 2160 : PetscInt row_v = i_ch + _n_channels * iz_ind;
1246 2160 : LibmeshPetscCall(
1247 : VecSetValues(_amc_pressure_force_rhs, 1, &row_v, &value_v, ADD_VALUES));
1248 : }
1249 : else
1250 : {
1251 19440 : PetscInt row = i_ch + _n_channels * iz_ind;
1252 19440 : PetscInt col = i_ch + _n_channels * (iz_ind + 1);
1253 19440 : PetscScalar value = 1.0 * S_interp;
1254 19440 : LibmeshPetscCall(
1255 : MatSetValues(_amc_pressure_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
1256 :
1257 19440 : if (_segregated_bool)
1258 : {
1259 0 : auto dp_in = (*_DP_soln)(node_in);
1260 0 : auto dp_out = (*_DP_soln)(node_out);
1261 0 : auto dp_interp = computeInterpolatedValue(dp_out, dp_in, 0.5);
1262 0 : PetscScalar value_v = -1.0 * dp_interp * S_interp;
1263 0 : PetscInt row_v = i_ch + _n_channels * iz_ind;
1264 0 : LibmeshPetscCall(
1265 : VecSetValues(_amc_pressure_force_rhs, 1, &row_v, &value_v, ADD_VALUES));
1266 : }
1267 : }
1268 : }
1269 : }
1270 : // Solving pressure problem
1271 60 : LibmeshPetscCall(MatAssemblyBegin(_amc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
1272 60 : LibmeshPetscCall(MatAssemblyEnd(_amc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
1273 60 : if (_verbose_subchannel)
1274 60 : _console << "Block: " << iblock << " - Axial momentum pressure force matrix assembled"
1275 60 : << std::endl;
1276 :
1277 60 : if (_segregated_bool)
1278 : {
1279 : KSP ksploc;
1280 : PC pc;
1281 : Vec sol;
1282 0 : LibmeshPetscCall(VecDuplicate(_amc_pressure_force_rhs, &sol));
1283 0 : LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
1284 0 : LibmeshPetscCall(KSPSetOperators(ksploc, _amc_pressure_force_mat, _amc_pressure_force_mat));
1285 0 : LibmeshPetscCall(KSPGetPC(ksploc, &pc));
1286 0 : LibmeshPetscCall(PCSetType(pc, PCJACOBI));
1287 0 : LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
1288 0 : LibmeshPetscCall(KSPSetFromOptions(ksploc));
1289 0 : LibmeshPetscCall(KSPSolve(ksploc, _amc_pressure_force_rhs, sol));
1290 : PetscScalar * xx;
1291 0 : LibmeshPetscCall(VecGetArray(sol, &xx));
1292 : // update Pressure solution
1293 0 : for (unsigned int iz = last_node; iz > first_node - 1; iz--)
1294 : {
1295 0 : auto iz_ind = iz - first_node;
1296 0 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1297 : {
1298 0 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1299 0 : PetscScalar value = xx[iz_ind * _n_channels + i_ch];
1300 0 : _P_soln->set(node_in, value);
1301 : }
1302 : }
1303 0 : LibmeshPetscCall(VecZeroEntries(_amc_pressure_force_rhs));
1304 0 : LibmeshPetscCall(KSPDestroy(&ksploc));
1305 0 : LibmeshPetscCall(VecDestroy(&sol));
1306 : }
1307 : }
1308 : }
1309 43328 : }
1310 :
1311 : Real
1312 3846420 : SubChannel1PhaseProblem::computeAddedHeatDuct(unsigned int i_ch, unsigned int iz)
1313 : {
1314 3846420 : if (_duct_mesh_exist)
1315 : {
1316 69300 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
1317 69300 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
1318 : {
1319 29700 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
1320 29700 : auto * node_in_chan = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1321 29700 : auto * node_out_chan = _subchannel_mesh.getChannelNode(i_ch, iz);
1322 29700 : auto * node_in_duct = _subchannel_mesh.getDuctNodeFromChannel(node_in_chan);
1323 29700 : auto * node_out_duct = _subchannel_mesh.getDuctNodeFromChannel(node_out_chan);
1324 29700 : auto heat_rate_in = (*_q_prime_duct_soln)(node_in_duct);
1325 29700 : auto heat_rate_out = (*_q_prime_duct_soln)(node_out_duct);
1326 29700 : return 0.5 * (heat_rate_in + heat_rate_out) * dz / _outer_channels;
1327 : }
1328 : else
1329 : {
1330 : return 0.0;
1331 : }
1332 : }
1333 : else
1334 : {
1335 : return 0;
1336 : }
1337 : }
1338 :
1339 : void
1340 2930 : SubChannel1PhaseProblem::computeT(int iblock)
1341 : {
1342 2930 : unsigned int last_node = (iblock + 1) * _block_size;
1343 2930 : unsigned int first_node = iblock * _block_size + 1;
1344 100876 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1345 : {
1346 4269830 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1347 : {
1348 4171884 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
1349 4171884 : _T_soln->set(node, _fp->T_from_p_h((*_P_soln)(node) + _P_out, (*_h_soln)(node)));
1350 : }
1351 : }
1352 2930 : }
1353 :
1354 : void
1355 3058 : SubChannel1PhaseProblem::computeRho(int iblock)
1356 : {
1357 3058 : unsigned int last_node = (iblock + 1) * _block_size;
1358 3058 : unsigned int first_node = iblock * _block_size + 1;
1359 3058 : if (iblock == 0)
1360 : {
1361 130852 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1362 : {
1363 127794 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
1364 127794 : _rho_soln->set(node, _fp->rho_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node)));
1365 : }
1366 : }
1367 108684 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1368 : {
1369 4438790 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1370 : {
1371 4333164 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
1372 4333164 : _rho_soln->set(node, _fp->rho_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node)));
1373 : }
1374 : }
1375 3058 : }
1376 :
1377 : void
1378 3058 : SubChannel1PhaseProblem::computeMu(int iblock)
1379 : {
1380 3058 : unsigned int last_node = (iblock + 1) * _block_size;
1381 3058 : unsigned int first_node = iblock * _block_size + 1;
1382 3058 : if (iblock == 0)
1383 : {
1384 130852 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1385 : {
1386 127794 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
1387 127794 : _mu_soln->set(node, _fp->mu_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node)));
1388 : }
1389 : }
1390 108684 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1391 : {
1392 4438790 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1393 : {
1394 4333164 : auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
1395 4333164 : _mu_soln->set(node, _fp->mu_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node)));
1396 : }
1397 : }
1398 3058 : }
1399 :
1400 : void
1401 43328 : SubChannel1PhaseProblem::computeWijResidual(int iblock)
1402 : {
1403 : // Cross flow residual
1404 43328 : if (!_implicit_bool)
1405 : {
1406 26288 : unsigned int last_node = (iblock + 1) * _block_size;
1407 26288 : unsigned int first_node = iblock * _block_size + 1;
1408 26288 : const Real & pitch = _subchannel_mesh.getPitch();
1409 510658 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1410 : {
1411 484370 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
1412 16391690 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
1413 : {
1414 15907320 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
1415 : unsigned int i_ch = chans.first;
1416 : unsigned int j_ch = chans.second;
1417 15907320 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1418 15907320 : auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
1419 15907320 : auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
1420 15907320 : auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
1421 15907320 : auto rho_i = (*_rho_soln)(node_in_i);
1422 15907320 : auto rho_j = (*_rho_soln)(node_in_j);
1423 15907320 : auto Si = (*_S_flow_soln)(node_in_i);
1424 15907320 : auto Sj = (*_S_flow_soln)(node_in_j);
1425 15907320 : auto Sij = dz * _subchannel_mesh.getGapWidth(iz, i_gap);
1426 15907320 : auto Lij = pitch;
1427 : // total local form loss in the ij direction
1428 15907320 : auto friction_term = _kij * _Wij(i_gap, iz) * std::abs(_Wij(i_gap, iz));
1429 15907320 : auto DPij = (*_P_soln)(node_in_i) - (*_P_soln)(node_in_j);
1430 : // Figure out donor cell density
1431 : auto rho_star = 0.0;
1432 15907320 : if (_Wij(i_gap, iz) > 0.0)
1433 : rho_star = rho_i;
1434 8381352 : else if (_Wij(i_gap, iz) < 0.0)
1435 : rho_star = rho_j;
1436 : else
1437 1352862 : rho_star = (rho_i + rho_j) / 2.0;
1438 : auto mass_term_out =
1439 15907320 : (*_mdot_soln)(node_out_i) / Si / rho_i + (*_mdot_soln)(node_out_j) / Sj / rho_j;
1440 : auto mass_term_in =
1441 15907320 : (*_mdot_soln)(node_in_i) / Si / rho_i + (*_mdot_soln)(node_in_j) / Sj / rho_j;
1442 15907320 : auto term_out = Sij * rho_star * (Lij / dz) * mass_term_out * _Wij(i_gap, iz);
1443 15907320 : auto term_in = Sij * rho_star * (Lij / dz) * mass_term_in * _Wij(i_gap, iz - 1);
1444 15907320 : auto inertia_term = term_out - term_in;
1445 15907320 : auto pressure_term = 2 * std::pow(Sij, 2.0) * DPij * rho_star;
1446 : auto time_term =
1447 15907320 : _TR * 2.0 * (_Wij(i_gap, iz) - _Wij_old(i_gap, iz)) * Lij * Sij * rho_star / _dt;
1448 :
1449 15907320 : _Wij_residual_matrix(i_gap, iz - 1 - iblock * _block_size) =
1450 15907320 : time_term + friction_term + inertia_term - pressure_term;
1451 : }
1452 : }
1453 : }
1454 : else
1455 : {
1456 : // Initializing to zero the elements of the lateral momentum assembly
1457 17040 : LibmeshPetscCall(MatZeroEntries(_cmc_time_derivative_mat));
1458 17040 : LibmeshPetscCall(MatZeroEntries(_cmc_advective_derivative_mat));
1459 17040 : LibmeshPetscCall(MatZeroEntries(_cmc_friction_force_mat));
1460 17040 : LibmeshPetscCall(MatZeroEntries(_cmc_pressure_force_mat));
1461 17040 : LibmeshPetscCall(VecZeroEntries(_cmc_time_derivative_rhs));
1462 17040 : LibmeshPetscCall(VecZeroEntries(_cmc_advective_derivative_rhs));
1463 17040 : LibmeshPetscCall(VecZeroEntries(_cmc_friction_force_rhs));
1464 17040 : LibmeshPetscCall(VecZeroEntries(_cmc_pressure_force_rhs));
1465 17040 : LibmeshPetscCall(MatZeroEntries(_cmc_sys_Wij_mat));
1466 17040 : LibmeshPetscCall(VecZeroEntries(_cmc_sys_Wij_rhs));
1467 17040 : unsigned int last_node = (iblock + 1) * _block_size;
1468 17040 : unsigned int first_node = iblock * _block_size + 1;
1469 17040 : const Real & pitch = _subchannel_mesh.getPitch();
1470 225216 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1471 : {
1472 208176 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
1473 208176 : auto iz_ind = iz - first_node;
1474 13278048 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
1475 : {
1476 13069872 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
1477 : unsigned int i_ch = chans.first;
1478 : unsigned int j_ch = chans.second;
1479 13069872 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1480 13069872 : auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
1481 13069872 : auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
1482 13069872 : auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
1483 :
1484 : // inlet, outlet, and interpolated densities
1485 13069872 : auto rho_i_in = (*_rho_soln)(node_in_i);
1486 13069872 : auto rho_i_out = (*_rho_soln)(node_out_i);
1487 13069872 : auto rho_i_interp = computeInterpolatedValue(rho_i_out, rho_i_in, 0.5);
1488 13069872 : auto rho_j_in = (*_rho_soln)(node_in_j);
1489 13069872 : auto rho_j_out = (*_rho_soln)(node_out_j);
1490 13069872 : auto rho_j_interp = computeInterpolatedValue(rho_j_out, rho_j_in, 0.5);
1491 :
1492 : // inlet, outlet, and interpolated areas
1493 13069872 : auto S_i_in = (*_S_flow_soln)(node_in_i);
1494 13069872 : auto S_i_out = (*_S_flow_soln)(node_out_i);
1495 13069872 : auto S_j_in = (*_S_flow_soln)(node_in_j);
1496 13069872 : auto S_j_out = (*_S_flow_soln)(node_out_j);
1497 :
1498 : // Cross-sectional gap area
1499 13069872 : auto Sij = dz * _subchannel_mesh.getGapWidth(iz, i_gap);
1500 13069872 : auto Lij = pitch;
1501 :
1502 : // Figure out donor cell density
1503 : auto rho_star = 0.0;
1504 13069872 : if (_Wij(i_gap, iz) > 0.0)
1505 : rho_star = rho_i_interp;
1506 5658102 : else if (_Wij(i_gap, iz) < 0.0)
1507 : rho_star = rho_j_interp;
1508 : else
1509 468762 : rho_star = (rho_i_interp + rho_j_interp) / 2.0;
1510 :
1511 : // Assembling time derivative
1512 13069872 : PetscScalar time_factor = _TR * Lij * Sij * rho_star / _dt;
1513 13069872 : PetscInt row_td = i_gap + _n_gaps * iz_ind;
1514 13069872 : PetscInt col_td = i_gap + _n_gaps * iz_ind;
1515 13069872 : PetscScalar value_td = time_factor;
1516 13069872 : LibmeshPetscCall(MatSetValues(
1517 : _cmc_time_derivative_mat, 1, &row_td, 1, &col_td, &value_td, INSERT_VALUES));
1518 13069872 : PetscScalar value_td_rhs = time_factor * _Wij_old(i_gap, iz);
1519 13069872 : LibmeshPetscCall(
1520 : VecSetValues(_cmc_time_derivative_rhs, 1, &row_td, &value_td_rhs, INSERT_VALUES));
1521 :
1522 : // Assembling inertial term
1523 : PetscScalar Pe = 0.5;
1524 13069872 : auto alpha = computeInterpolationCoefficients(Pe);
1525 13069872 : auto mass_term_out = (*_mdot_soln)(node_out_i) / S_i_out / rho_i_out +
1526 13069872 : (*_mdot_soln)(node_out_j) / S_j_out / rho_j_out;
1527 13069872 : auto mass_term_in = (*_mdot_soln)(node_in_i) / S_i_in / rho_i_in +
1528 13069872 : (*_mdot_soln)(node_in_j) / S_j_in / rho_j_in;
1529 13069872 : auto term_out = Sij * rho_star * (Lij / dz) * mass_term_out / 2.0;
1530 13069872 : auto term_in = Sij * rho_star * (Lij / dz) * mass_term_in / 2.0;
1531 13069872 : if (iz == first_node)
1532 : {
1533 1041768 : PetscInt row_ad = i_gap + _n_gaps * iz_ind;
1534 1041768 : PetscScalar value_ad = term_in * alpha * _Wij(i_gap, iz - 1);
1535 1041768 : LibmeshPetscCall(
1536 : VecSetValues(_cmc_advective_derivative_rhs, 1, &row_ad, &value_ad, ADD_VALUES));
1537 :
1538 1041768 : PetscInt col_ad = i_gap + _n_gaps * iz_ind;
1539 1041768 : value_ad = -1.0 * term_in * (1.0 - alpha) + term_out * alpha;
1540 1041768 : LibmeshPetscCall(MatSetValues(
1541 : _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
1542 :
1543 1041768 : col_ad = i_gap + _n_gaps * (iz_ind + 1);
1544 1041768 : value_ad = term_out * (1.0 - alpha);
1545 1041768 : LibmeshPetscCall(MatSetValues(
1546 : _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
1547 : }
1548 12028104 : else if (iz == last_node)
1549 : {
1550 1041768 : PetscInt row_ad = i_gap + _n_gaps * iz_ind;
1551 1041768 : PetscInt col_ad = i_gap + _n_gaps * (iz_ind - 1);
1552 1041768 : PetscScalar value_ad = -1.0 * term_in * alpha;
1553 1041768 : LibmeshPetscCall(MatSetValues(
1554 : _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
1555 :
1556 1041768 : col_ad = i_gap + _n_gaps * iz_ind;
1557 1041768 : value_ad = -1.0 * term_in * (1.0 - alpha) + term_out * alpha;
1558 1041768 : LibmeshPetscCall(MatSetValues(
1559 : _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
1560 :
1561 1041768 : value_ad = -1.0 * term_out * (1.0 - alpha) * _Wij(i_gap, iz);
1562 1041768 : LibmeshPetscCall(
1563 : VecSetValues(_cmc_advective_derivative_rhs, 1, &row_ad, &value_ad, ADD_VALUES));
1564 : }
1565 : else
1566 : {
1567 10986336 : PetscInt row_ad = i_gap + _n_gaps * iz_ind;
1568 10986336 : PetscInt col_ad = i_gap + _n_gaps * (iz_ind - 1);
1569 10986336 : PetscScalar value_ad = -1.0 * term_in * alpha;
1570 10986336 : LibmeshPetscCall(MatSetValues(
1571 : _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
1572 :
1573 10986336 : col_ad = i_gap + _n_gaps * iz_ind;
1574 10986336 : value_ad = -1.0 * term_in * (1.0 - alpha) + term_out * alpha;
1575 10986336 : LibmeshPetscCall(MatSetValues(
1576 : _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
1577 :
1578 10986336 : col_ad = i_gap + _n_gaps * (iz_ind + 1);
1579 10986336 : value_ad = term_out * (1.0 - alpha);
1580 10986336 : LibmeshPetscCall(MatSetValues(
1581 : _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
1582 : }
1583 : // Assembling friction force
1584 13069872 : PetscInt row_ff = i_gap + _n_gaps * iz_ind;
1585 13069872 : PetscInt col_ff = i_gap + _n_gaps * iz_ind;
1586 13069872 : PetscScalar value_ff = _kij * std::abs(_Wij(i_gap, iz)) / 2.0;
1587 13069872 : LibmeshPetscCall(MatSetValues(
1588 : _cmc_friction_force_mat, 1, &row_ff, 1, &col_ff, &value_ff, INSERT_VALUES));
1589 :
1590 : // Assembling pressure force
1591 13069872 : alpha = computeInterpolationCoefficients(Pe);
1592 :
1593 13069872 : if (!_staggered_pressure_bool)
1594 : {
1595 13033872 : PetscScalar pressure_factor = std::pow(Sij, 2.0) * rho_star;
1596 13033872 : PetscInt row_pf = i_gap + _n_gaps * iz_ind;
1597 13033872 : PetscInt col_pf = i_ch + _n_channels * iz_ind;
1598 13033872 : PetscScalar value_pf = -1.0 * alpha * pressure_factor;
1599 13033872 : LibmeshPetscCall(
1600 : MatSetValues(_cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
1601 13033872 : col_pf = j_ch + _n_channels * iz_ind;
1602 13033872 : value_pf = alpha * pressure_factor;
1603 13033872 : LibmeshPetscCall(
1604 : MatSetValues(_cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
1605 :
1606 13033872 : if (iz == last_node)
1607 : {
1608 1038168 : PetscInt row_pf = i_gap + _n_gaps * iz_ind;
1609 1038168 : PetscScalar value_pf = (1.0 - alpha) * pressure_factor * (*_P_soln)(node_out_i);
1610 1038168 : LibmeshPetscCall(
1611 : VecSetValues(_cmc_pressure_force_rhs, 1, &row_pf, &value_pf, ADD_VALUES));
1612 1038168 : value_pf = -1.0 * (1.0 - alpha) * pressure_factor * (*_P_soln)(node_out_j);
1613 1038168 : LibmeshPetscCall(
1614 : VecSetValues(_cmc_pressure_force_rhs, 1, &row_pf, &value_pf, ADD_VALUES));
1615 : }
1616 : else
1617 : {
1618 11995704 : row_pf = i_gap + _n_gaps * iz_ind;
1619 11995704 : col_pf = i_ch + _n_channels * (iz_ind + 1);
1620 11995704 : value_pf = -1.0 * (1.0 - alpha) * pressure_factor;
1621 11995704 : LibmeshPetscCall(MatSetValues(
1622 : _cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
1623 11995704 : col_pf = j_ch + _n_channels * (iz_ind + 1);
1624 11995704 : value_pf = (1.0 - alpha) * pressure_factor;
1625 11995704 : LibmeshPetscCall(MatSetValues(
1626 : _cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
1627 : }
1628 : }
1629 : else
1630 : {
1631 36000 : PetscScalar pressure_factor = std::pow(Sij, 2.0) * rho_star;
1632 36000 : PetscInt row_pf = i_gap + _n_gaps * iz_ind;
1633 36000 : PetscInt col_pf = i_ch + _n_channels * iz_ind;
1634 36000 : PetscScalar value_pf = -1.0 * pressure_factor;
1635 36000 : LibmeshPetscCall(
1636 : MatSetValues(_cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
1637 36000 : col_pf = j_ch + _n_channels * iz_ind;
1638 36000 : value_pf = pressure_factor;
1639 36000 : LibmeshPetscCall(
1640 : MatSetValues(_cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
1641 : }
1642 : }
1643 : }
1644 : /// Assembling system
1645 17040 : LibmeshPetscCall(MatZeroEntries(_cmc_sys_Wij_mat));
1646 17040 : LibmeshPetscCall(VecZeroEntries(_cmc_sys_Wij_rhs));
1647 17040 : LibmeshPetscCall(MatAssemblyBegin(_cmc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
1648 17040 : LibmeshPetscCall(MatAssemblyEnd(_cmc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
1649 17040 : LibmeshPetscCall(MatAssemblyBegin(_cmc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
1650 17040 : LibmeshPetscCall(MatAssemblyEnd(_cmc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
1651 17040 : LibmeshPetscCall(MatAssemblyBegin(_cmc_friction_force_mat, MAT_FINAL_ASSEMBLY));
1652 17040 : LibmeshPetscCall(MatAssemblyEnd(_cmc_friction_force_mat, MAT_FINAL_ASSEMBLY));
1653 17040 : LibmeshPetscCall(MatAssemblyBegin(_cmc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
1654 17040 : LibmeshPetscCall(MatAssemblyEnd(_cmc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
1655 17040 : LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1656 17040 : LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1657 : // Matrix
1658 : #if !PETSC_VERSION_LESS_THAN(3, 15, 0)
1659 17040 : LibmeshPetscCall(
1660 : MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
1661 17040 : LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1662 17040 : LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1663 17040 : LibmeshPetscCall(
1664 : MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
1665 17040 : LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1666 17040 : LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1667 17040 : LibmeshPetscCall(
1668 : MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_friction_force_mat, UNKNOWN_NONZERO_PATTERN));
1669 : #else
1670 : LibmeshPetscCall(
1671 : MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1672 : LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1673 : LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1674 : LibmeshPetscCall(
1675 : MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
1676 : LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1677 : LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1678 : LibmeshPetscCall(
1679 : MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_friction_force_mat, DIFFERENT_NONZERO_PATTERN));
1680 : #endif
1681 17040 : LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1682 17040 : LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
1683 : // RHS
1684 17040 : LibmeshPetscCall(VecAXPY(_cmc_sys_Wij_rhs, 1.0, _cmc_time_derivative_rhs));
1685 17040 : LibmeshPetscCall(VecAXPY(_cmc_sys_Wij_rhs, 1.0, _cmc_advective_derivative_rhs));
1686 17040 : LibmeshPetscCall(VecAXPY(_cmc_sys_Wij_rhs, 1.0, _cmc_friction_force_rhs));
1687 :
1688 17040 : if (_segregated_bool)
1689 : {
1690 : // Assembly the matrix system
1691 : Vec sol_holder_P;
1692 14346 : LibmeshPetscCall(createPetscVector(sol_holder_P, _block_size * _n_gaps));
1693 : Vec sol_holder_W;
1694 14346 : LibmeshPetscCall(createPetscVector(sol_holder_W, _block_size * _n_gaps));
1695 14346 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
1696 : _prodp, *_P_soln, iblock * _block_size, (iblock + 1) * _block_size - 1, _n_channels));
1697 14346 : LibmeshPetscCall(MatMult(_cmc_sys_Wij_mat, _Wij_vec, sol_holder_W));
1698 14346 : LibmeshPetscCall(VecAXPY(sol_holder_W, -1.0, _cmc_sys_Wij_rhs));
1699 14346 : LibmeshPetscCall(MatMult(_cmc_pressure_force_mat, _prodp, sol_holder_P));
1700 14346 : LibmeshPetscCall(VecAXPY(sol_holder_P, -1.0, _cmc_pressure_force_rhs));
1701 14346 : LibmeshPetscCall(VecAXPY(sol_holder_W, 1.0, sol_holder_P));
1702 : PetscScalar * xx;
1703 14346 : LibmeshPetscCall(VecGetArray(sol_holder_W, &xx));
1704 131796 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1705 : {
1706 117450 : auto iz_ind = iz - first_node;
1707 7164450 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
1708 : {
1709 7047000 : _Wij_residual_matrix(i_gap, iz - 1 - iblock * _block_size) = xx[iz_ind * _n_gaps + i_gap];
1710 : }
1711 : }
1712 14346 : LibmeshPetscCall(VecDestroy(&sol_holder_P));
1713 14346 : LibmeshPetscCall(VecDestroy(&sol_holder_W));
1714 : }
1715 : }
1716 43328 : }
1717 :
1718 : void
1719 46022 : SubChannel1PhaseProblem::computeWijPrime(int iblock)
1720 : {
1721 46022 : unsigned int last_node = (iblock + 1) * _block_size;
1722 46022 : unsigned int first_node = iblock * _block_size + 1;
1723 829294 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1724 : {
1725 783272 : auto dz = _z_grid[iz] - _z_grid[iz - 1];
1726 35783336 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
1727 : {
1728 35000064 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
1729 : unsigned int i_ch = chans.first;
1730 : unsigned int j_ch = chans.second;
1731 35000064 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
1732 35000064 : auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
1733 35000064 : auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
1734 35000064 : auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
1735 35000064 : auto Si_in = (*_S_flow_soln)(node_in_i);
1736 35000064 : auto Sj_in = (*_S_flow_soln)(node_in_j);
1737 35000064 : auto Si_out = (*_S_flow_soln)(node_out_i);
1738 35000064 : auto Sj_out = (*_S_flow_soln)(node_out_j);
1739 35000064 : auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
1740 35000064 : auto Sij = dz * gap;
1741 : auto avg_massflux =
1742 35000064 : 0.5 * (((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
1743 35000064 : ((*_mdot_soln)(node_out_i) + (*_mdot_soln)(node_out_j)) / (Si_out + Sj_out));
1744 35000064 : auto beta = computeBeta(i_gap, iz, /*enthalpy=*/false);
1745 :
1746 35000064 : if (!_implicit_bool)
1747 : {
1748 15907320 : _WijPrime(i_gap, iz) = beta * avg_massflux * Sij;
1749 : }
1750 : else
1751 : {
1752 19092744 : auto iz_ind = iz - first_node;
1753 19092744 : PetscScalar base_value = beta * 0.5 * Sij;
1754 :
1755 : // Bottom values
1756 19092744 : if (iz == first_node)
1757 : {
1758 1222776 : PetscScalar value_tl = -1.0 * base_value / (Si_in + Sj_in) *
1759 1222776 : ((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j));
1760 1222776 : PetscInt row = i_gap + _n_gaps * iz_ind;
1761 1222776 : LibmeshPetscCall(
1762 : VecSetValues(_amc_turbulent_cross_flows_rhs, 1, &row, &value_tl, INSERT_VALUES));
1763 : }
1764 : else
1765 : {
1766 17869968 : PetscScalar value_tl = base_value / (Si_in + Sj_in);
1767 17869968 : PetscInt row = i_gap + _n_gaps * iz_ind;
1768 :
1769 17869968 : PetscInt col_ich = i_ch + _n_channels * (iz_ind - 1);
1770 17869968 : LibmeshPetscCall(MatSetValues(
1771 : _amc_turbulent_cross_flows_mat, 1, &row, 1, &col_ich, &value_tl, INSERT_VALUES));
1772 :
1773 17869968 : PetscInt col_jch = j_ch + _n_channels * (iz_ind - 1);
1774 17869968 : LibmeshPetscCall(MatSetValues(
1775 : _amc_turbulent_cross_flows_mat, 1, &row, 1, &col_jch, &value_tl, INSERT_VALUES));
1776 : }
1777 :
1778 : // Top values
1779 19092744 : PetscScalar value_bl = base_value / (Si_out + Sj_out);
1780 19092744 : PetscInt row = i_gap + _n_gaps * iz_ind;
1781 :
1782 19092744 : PetscInt col_ich = i_ch + _n_channels * iz_ind;
1783 19092744 : LibmeshPetscCall(MatSetValues(
1784 : _amc_turbulent_cross_flows_mat, 1, &row, 1, &col_ich, &value_bl, INSERT_VALUES));
1785 :
1786 19092744 : PetscInt col_jch = j_ch + _n_channels * iz_ind;
1787 19092744 : LibmeshPetscCall(MatSetValues(
1788 : _amc_turbulent_cross_flows_mat, 1, &row, 1, &col_jch, &value_bl, INSERT_VALUES));
1789 : }
1790 : }
1791 : }
1792 :
1793 46022 : if (_implicit_bool)
1794 : {
1795 19734 : LibmeshPetscCall(MatAssemblyBegin(_amc_turbulent_cross_flows_mat, MAT_FINAL_ASSEMBLY));
1796 19734 : LibmeshPetscCall(MatAssemblyEnd(_amc_turbulent_cross_flows_mat, MAT_FINAL_ASSEMBLY));
1797 :
1798 : /// Update turbulent crossflow
1799 : Vec loc_prod;
1800 : Vec loc_Wij;
1801 19734 : LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &loc_prod));
1802 19734 : LibmeshPetscCall(VecDuplicate(_Wij_vec, &loc_Wij));
1803 19734 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
1804 : loc_prod, *_mdot_soln, first_node, last_node, _n_channels));
1805 19734 : LibmeshPetscCall(MatMult(_amc_turbulent_cross_flows_mat, loc_prod, loc_Wij));
1806 19734 : LibmeshPetscCall(VecAXPY(loc_Wij, -1.0, _amc_turbulent_cross_flows_rhs));
1807 19734 : LibmeshPetscCall(populateDenseFromVector<libMesh::DenseMatrix<Real>>(
1808 : loc_Wij, _WijPrime, first_node, last_node, _n_gaps));
1809 19734 : LibmeshPetscCall(VecDestroy(&loc_prod));
1810 19734 : LibmeshPetscCall(VecDestroy(&loc_Wij));
1811 : }
1812 46022 : }
1813 :
1814 : libMesh::DenseVector<Real>
1815 40634 : SubChannel1PhaseProblem::residualFunction(int iblock, libMesh::DenseVector<Real> solution)
1816 : {
1817 40634 : unsigned int last_node = (iblock + 1) * _block_size;
1818 40634 : unsigned int first_node = iblock * _block_size + 1;
1819 40634 : libMesh::DenseVector<Real> Wij_residual_vector(_n_gaps * _block_size, 0.0);
1820 : // Assign the solution to the cross-flow matrix
1821 : int i = 0;
1822 642454 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1823 : {
1824 23556140 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
1825 : {
1826 22954320 : _Wij(i_gap, iz) = solution(i);
1827 22954320 : i++;
1828 : }
1829 : }
1830 :
1831 : // Calculating sum of crossflows
1832 40634 : computeSumWij(iblock);
1833 : // Solving axial flux
1834 40634 : computeMdot(iblock);
1835 : // Calculation of turbulent Crossflow
1836 40634 : computeWijPrime(iblock);
1837 : // Solving for Pressure Drop
1838 40634 : computeDP(iblock);
1839 : // Solving for pressure
1840 40634 : computeP(iblock);
1841 : // Populating lateral crossflow residual matrix
1842 40634 : computeWijResidual(iblock);
1843 :
1844 : // Turn the residual matrix into a residual vector
1845 642454 : for (unsigned int iz = 0; iz < _block_size; iz++)
1846 : {
1847 23556140 : for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
1848 : {
1849 22954320 : int i = _n_gaps * iz + i_gap; // column wise transfer
1850 22954320 : Wij_residual_vector(i) = _Wij_residual_matrix(i_gap, iz);
1851 : }
1852 : }
1853 40634 : return Wij_residual_vector;
1854 : }
1855 :
1856 : PetscErrorCode
1857 376 : SubChannel1PhaseProblem::petscSnesSolver(int iblock,
1858 : const libMesh::DenseVector<Real> & solution,
1859 : libMesh::DenseVector<Real> & root)
1860 : {
1861 : SNES snes;
1862 : KSP ksp;
1863 : PC pc;
1864 : Vec x, r;
1865 : PetscMPIInt size;
1866 : PetscScalar * xx;
1867 :
1868 : PetscFunctionBegin;
1869 376 : PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1870 376 : if (size > 1)
1871 0 : SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Example is only for sequential runs");
1872 376 : LibmeshPetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
1873 376 : LibmeshPetscCall(VecCreate(PETSC_COMM_WORLD, &x));
1874 376 : LibmeshPetscCall(VecSetSizes(x, PETSC_DECIDE, _block_size * _n_gaps));
1875 376 : LibmeshPetscCall(VecSetFromOptions(x));
1876 376 : LibmeshPetscCall(VecDuplicate(x, &r));
1877 :
1878 : #if PETSC_VERSION_LESS_THAN(3, 13, 0)
1879 : LibmeshPetscCall(PetscOptionsSetValue(PETSC_NULL, "-snes_mf", PETSC_NULL));
1880 : #else
1881 376 : LibmeshPetscCall(SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE));
1882 : #endif
1883 : Ctx ctx;
1884 376 : ctx.iblock = iblock;
1885 376 : ctx.schp = this;
1886 376 : LibmeshPetscCall(SNESSetFunction(snes, r, formFunction, &ctx));
1887 376 : LibmeshPetscCall(SNESGetKSP(snes, &ksp));
1888 376 : LibmeshPetscCall(KSPGetPC(ksp, &pc));
1889 376 : LibmeshPetscCall(PCSetType(pc, PCNONE));
1890 376 : LibmeshPetscCall(KSPSetTolerances(ksp, _rtol, _atol, _dtol, _maxit));
1891 376 : LibmeshPetscCall(SNESSetFromOptions(snes));
1892 376 : LibmeshPetscCall(VecGetArray(x, &xx));
1893 316096 : for (unsigned int i = 0; i < _block_size * _n_gaps; i++)
1894 : {
1895 315720 : xx[i] = solution(i);
1896 : }
1897 376 : LibmeshPetscCall(VecRestoreArray(x, &xx));
1898 :
1899 376 : LibmeshPetscCall(SNESSolve(snes, NULL, x));
1900 376 : LibmeshPetscCall(VecGetArray(x, &xx));
1901 316096 : for (unsigned int i = 0; i < _block_size * _n_gaps; i++)
1902 315720 : root(i) = xx[i];
1903 :
1904 376 : LibmeshPetscCall(VecRestoreArray(x, &xx));
1905 376 : LibmeshPetscCall(VecDestroy(&x));
1906 376 : LibmeshPetscCall(VecDestroy(&r));
1907 376 : LibmeshPetscCall(SNESDestroy(&snes));
1908 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
1909 : }
1910 :
1911 : PetscErrorCode
1912 2694 : SubChannel1PhaseProblem::implicitPetscSolve(int iblock)
1913 : {
1914 : bool lag_block_thermal_solve = true;
1915 : Vec b_nest, x_nest; /* approx solution, RHS, exact solution */
1916 : Mat A_nest; /* linear system matrix */
1917 : KSP ksp; /* linear solver context */
1918 : PC pc; /* preconditioner context */
1919 :
1920 : PetscFunctionBegin;
1921 2694 : PetscInt Q = _monolithic_thermal_bool ? 4 : 3;
1922 2694 : std::vector<Mat> mat_array(Q * Q);
1923 2694 : std::vector<Vec> vec_array(Q);
1924 :
1925 : /// Initializing flags
1926 : bool _axial_mass_flow_tight_coupling = true;
1927 : bool _pressure_axial_momentum_tight_coupling = true;
1928 : bool _pressure_cross_momentum_tight_coupling = true;
1929 2694 : unsigned int first_node = iblock * _block_size + 1;
1930 2694 : unsigned int last_node = (iblock + 1) * _block_size;
1931 :
1932 : /// Assembling matrices
1933 : // Computing sum of crossflows with previous iteration
1934 2694 : computeSumWij(iblock);
1935 : // Assembling axial flux matrix
1936 2694 : computeMdot(iblock);
1937 : // Computing turbulent crossflow with previous step axial mass flows
1938 2694 : computeWijPrime(iblock);
1939 : // Assembling for Pressure Drop matrix
1940 2694 : computeDP(iblock);
1941 : // Assembling pressure matrix
1942 2694 : computeP(iblock);
1943 : // Assembling cross fluxes matrix
1944 2694 : computeWijResidual(iblock);
1945 : // If monolithic solve - Assembling enthalpy matrix
1946 2694 : if (_monolithic_thermal_bool)
1947 978 : computeh(iblock);
1948 :
1949 2694 : if (_verbose_subchannel)
1950 : {
1951 2028 : _console << "Starting nested system." << std::endl;
1952 2028 : _console << "Number of simultaneous variables: " << Q << std::endl;
1953 : }
1954 : // Mass conservation
1955 : PetscInt field_num = 0;
1956 2694 : LibmeshPetscCall(
1957 : MatDuplicate(_mc_axial_convection_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 0]));
1958 2694 : LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY));
1959 2694 : LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY));
1960 2694 : mat_array[Q * field_num + 1] = NULL;
1961 : if (_axial_mass_flow_tight_coupling)
1962 : {
1963 2694 : LibmeshPetscCall(MatDuplicate(_mc_sumWij_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 2]));
1964 2694 : LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY));
1965 2694 : LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY));
1966 : }
1967 : else
1968 : {
1969 : mat_array[Q * field_num + 2] = NULL;
1970 : }
1971 2694 : if (_monolithic_thermal_bool)
1972 : {
1973 978 : mat_array[Q * field_num + 3] = NULL;
1974 : }
1975 2694 : LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &vec_array[field_num]));
1976 2694 : LibmeshPetscCall(VecCopy(_mc_axial_convection_rhs, vec_array[field_num]));
1977 : if (!_axial_mass_flow_tight_coupling)
1978 : {
1979 : Vec sumWij_loc;
1980 : LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &sumWij_loc));
1981 : LibmeshPetscCall(VecSet(sumWij_loc, 0.0));
1982 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
1983 : {
1984 : auto iz_ind = iz - first_node;
1985 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
1986 : {
1987 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
1988 :
1989 : PetscScalar value_vec_2 = -1.0 * (*_SumWij_soln)(node_out);
1990 : PetscInt row_vec_2 = i_ch + _n_channels * iz_ind;
1991 : LibmeshPetscCall(VecSetValues(sumWij_loc, 1, &row_vec_2, &value_vec_2, ADD_VALUES));
1992 : }
1993 : }
1994 : LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, sumWij_loc));
1995 : LibmeshPetscCall(VecDestroy(&sumWij_loc));
1996 : }
1997 :
1998 2694 : if (_verbose_subchannel)
1999 2028 : _console << "Mass ok." << std::endl;
2000 : // Axial momentum conservation
2001 : field_num = 1;
2002 : if (_pressure_axial_momentum_tight_coupling)
2003 : {
2004 2694 : LibmeshPetscCall(
2005 : MatDuplicate(_amc_sys_mdot_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 0]));
2006 2694 : LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY));
2007 2694 : LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 0], MAT_FINAL_ASSEMBLY));
2008 : }
2009 : else
2010 : {
2011 : mat_array[Q * field_num + 0] = NULL;
2012 : }
2013 2694 : LibmeshPetscCall(
2014 : MatDuplicate(_amc_pressure_force_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 1]));
2015 2694 : LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY));
2016 2694 : LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY));
2017 2694 : mat_array[Q * field_num + 2] = NULL;
2018 2694 : if (_monolithic_thermal_bool)
2019 : {
2020 978 : mat_array[Q * field_num + 3] = NULL;
2021 : }
2022 2694 : LibmeshPetscCall(VecDuplicate(_amc_pressure_force_rhs, &vec_array[field_num]));
2023 2694 : LibmeshPetscCall(VecCopy(_amc_pressure_force_rhs, vec_array[field_num]));
2024 : if (_pressure_axial_momentum_tight_coupling)
2025 : {
2026 2694 : LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _amc_sys_mdot_rhs));
2027 : }
2028 : else
2029 : {
2030 : unsigned int last_node = (iblock + 1) * _block_size;
2031 : unsigned int first_node = iblock * _block_size + 1;
2032 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2033 : _prod, *_mdot_soln, first_node, last_node, _n_channels));
2034 : Vec ls;
2035 : LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &ls));
2036 : LibmeshPetscCall(MatMult(_amc_sys_mdot_mat, _prod, ls));
2037 : LibmeshPetscCall(VecAXPY(ls, -1.0, _amc_sys_mdot_rhs));
2038 : LibmeshPetscCall(VecAXPY(vec_array[field_num], -1.0, ls));
2039 : LibmeshPetscCall(VecDestroy(&ls));
2040 : }
2041 :
2042 2694 : if (_verbose_subchannel)
2043 2028 : _console << "Lin mom OK." << std::endl;
2044 :
2045 : // Cross momentum conservation
2046 : field_num = 2;
2047 2694 : mat_array[Q * field_num + 0] = NULL;
2048 : if (_pressure_cross_momentum_tight_coupling)
2049 : {
2050 2694 : LibmeshPetscCall(
2051 : MatDuplicate(_cmc_pressure_force_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 1]));
2052 2694 : LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY));
2053 2694 : LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 1], MAT_FINAL_ASSEMBLY));
2054 : }
2055 : else
2056 : {
2057 : mat_array[Q * field_num + 1] = NULL;
2058 : }
2059 :
2060 2694 : LibmeshPetscCall(MatDuplicate(_cmc_sys_Wij_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 2]));
2061 2694 : LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY));
2062 2694 : LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 2], MAT_FINAL_ASSEMBLY));
2063 2694 : if (_monolithic_thermal_bool)
2064 : {
2065 978 : mat_array[Q * field_num + 3] = NULL;
2066 : }
2067 :
2068 2694 : LibmeshPetscCall(VecDuplicate(_cmc_sys_Wij_rhs, &vec_array[field_num]));
2069 2694 : LibmeshPetscCall(VecCopy(_cmc_sys_Wij_rhs, vec_array[field_num]));
2070 : if (_pressure_cross_momentum_tight_coupling)
2071 : {
2072 2694 : LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _cmc_pressure_force_rhs));
2073 : }
2074 : else
2075 : {
2076 : Vec sol_holder_P;
2077 : LibmeshPetscCall(createPetscVector(sol_holder_P, _block_size * _n_gaps));
2078 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2079 : _prodp, *_P_soln, iblock * _block_size, (iblock + 1) * _block_size - 1, _n_channels));
2080 :
2081 : LibmeshPetscCall(MatMult(_cmc_pressure_force_mat, _prodp, sol_holder_P));
2082 : LibmeshPetscCall(VecAXPY(sol_holder_P, -1.0, _cmc_pressure_force_rhs));
2083 : LibmeshPetscCall(VecScale(sol_holder_P, 1.0));
2084 : LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, sol_holder_P));
2085 : LibmeshPetscCall(VecDestroy(&sol_holder_P));
2086 : }
2087 :
2088 2694 : if (_verbose_subchannel)
2089 2028 : _console << "Cross mom ok." << std::endl;
2090 :
2091 : // Energy conservation
2092 2694 : if (_monolithic_thermal_bool)
2093 : {
2094 : field_num = 3;
2095 978 : mat_array[Q * field_num + 0] = NULL;
2096 978 : mat_array[Q * field_num + 1] = NULL;
2097 978 : mat_array[Q * field_num + 2] = NULL;
2098 978 : LibmeshPetscCall(MatDuplicate(_hc_sys_h_mat, MAT_COPY_VALUES, &mat_array[Q * field_num + 3]));
2099 : if (lag_block_thermal_solve)
2100 : {
2101 978 : LibmeshPetscCall(MatZeroEntries(mat_array[Q * field_num + 3]));
2102 978 : LibmeshPetscCall(MatShift(mat_array[Q * field_num + 3], 1.0));
2103 : }
2104 978 : LibmeshPetscCall(MatAssemblyBegin(mat_array[Q * field_num + 3], MAT_FINAL_ASSEMBLY));
2105 978 : LibmeshPetscCall(MatAssemblyEnd(mat_array[Q * field_num + 3], MAT_FINAL_ASSEMBLY));
2106 978 : LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &vec_array[field_num]));
2107 978 : LibmeshPetscCall(VecCopy(_hc_sys_h_rhs, vec_array[field_num]));
2108 : if (lag_block_thermal_solve)
2109 : {
2110 978 : LibmeshPetscCall(VecZeroEntries(vec_array[field_num]));
2111 978 : LibmeshPetscCall(VecShift(vec_array[field_num], 1.0));
2112 : }
2113 :
2114 978 : if (_verbose_subchannel)
2115 810 : _console << "Energy ok." << std::endl;
2116 : }
2117 :
2118 : // Relaxing linear system
2119 : // Weaker relaxation
2120 : if (true)
2121 : {
2122 : // Estimating cross-flow resistances to achieve realizable solves
2123 2694 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2124 : _prod, *_mdot_soln, first_node, last_node, _n_channels));
2125 : Vec mdot_estimate;
2126 2694 : LibmeshPetscCall(createPetscVector(mdot_estimate, _block_size * _n_channels));
2127 : Vec pmat_diag;
2128 2694 : LibmeshPetscCall(createPetscVector(pmat_diag, _block_size * _n_channels));
2129 : Vec p_estimate;
2130 2694 : LibmeshPetscCall(createPetscVector(p_estimate, _block_size * _n_channels));
2131 : Vec unity_vec;
2132 2694 : LibmeshPetscCall(createPetscVector(unity_vec, _block_size * _n_channels));
2133 2694 : LibmeshPetscCall(VecSet(unity_vec, 1.0));
2134 : Vec sol_holder_P;
2135 2694 : LibmeshPetscCall(createPetscVector(sol_holder_P, _block_size * _n_gaps));
2136 : Vec diag_Wij_loc;
2137 2694 : LibmeshPetscCall(createPetscVector(diag_Wij_loc, _block_size * _n_gaps));
2138 : Vec Wij_estimate;
2139 2694 : LibmeshPetscCall(createPetscVector(Wij_estimate, _block_size * _n_gaps));
2140 : Vec unity_vec_Wij;
2141 2694 : LibmeshPetscCall(createPetscVector(unity_vec_Wij, _block_size * _n_gaps));
2142 2694 : LibmeshPetscCall(VecSet(unity_vec_Wij, 1.0));
2143 : Vec _Wij_loc_vec;
2144 2694 : LibmeshPetscCall(createPetscVector(_Wij_loc_vec, _block_size * _n_gaps));
2145 : Vec _Wij_old_loc_vec;
2146 2694 : LibmeshPetscCall(createPetscVector(_Wij_old_loc_vec, _block_size * _n_gaps));
2147 2694 : LibmeshPetscCall(MatMult(mat_array[Q], _prod, mdot_estimate));
2148 2694 : LibmeshPetscCall(MatGetDiagonal(mat_array[Q + 1], pmat_diag));
2149 2694 : LibmeshPetscCall(VecAXPY(pmat_diag, 1e-10, unity_vec));
2150 2694 : LibmeshPetscCall(VecPointwiseDivide(p_estimate, mdot_estimate, pmat_diag));
2151 2694 : LibmeshPetscCall(MatMult(mat_array[2 * Q + 1], p_estimate, sol_holder_P));
2152 2694 : LibmeshPetscCall(VecAXPY(sol_holder_P, -1.0, _cmc_pressure_force_rhs));
2153 2694 : LibmeshPetscCall(MatGetDiagonal(mat_array[2 * Q + 2], diag_Wij_loc));
2154 2694 : LibmeshPetscCall(VecAXPY(diag_Wij_loc, 1e-10, unity_vec_Wij));
2155 2694 : LibmeshPetscCall(VecPointwiseDivide(Wij_estimate, sol_holder_P, diag_Wij_loc));
2156 : Vec sumWij_loc;
2157 2694 : LibmeshPetscCall(createPetscVector(sumWij_loc, _block_size * _n_channels));
2158 93420 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
2159 : {
2160 90726 : auto iz_ind = iz - first_node;
2161 4251210 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
2162 : {
2163 4160484 : PetscScalar sumWij = 0.0;
2164 : unsigned int counter = 0;
2165 16206228 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
2166 : {
2167 12045744 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
2168 : unsigned int i_ch_loc = chans.first;
2169 12045744 : PetscInt row_vec = i_ch_loc + _n_channels * iz_ind;
2170 : PetscScalar loc_Wij_value;
2171 12045744 : LibmeshPetscCall(VecGetValues(sol_holder_P, 1, &row_vec, &loc_Wij_value));
2172 12045744 : sumWij += _subchannel_mesh.getCrossflowSign(i_ch, counter) * loc_Wij_value;
2173 12045744 : counter++;
2174 : }
2175 4160484 : PetscInt row_vec = i_ch + _n_channels * iz_ind;
2176 4160484 : LibmeshPetscCall(VecSetValues(sumWij_loc, 1, &row_vec, &sumWij, INSERT_VALUES));
2177 : }
2178 : }
2179 :
2180 : PetscScalar min_mdot;
2181 2694 : LibmeshPetscCall(VecAbs(_prod));
2182 2694 : LibmeshPetscCall(VecMin(_prod, NULL, &min_mdot));
2183 2694 : if (_verbose_subchannel)
2184 2028 : _console << "Minimum estimated mdot: " << min_mdot << std::endl;
2185 :
2186 2694 : LibmeshPetscCall(VecAbs(sumWij_loc));
2187 2694 : LibmeshPetscCall(VecMax(sumWij_loc, NULL, &_max_sumWij));
2188 2694 : _max_sumWij = std::max(1e-10, _max_sumWij);
2189 2694 : if (_verbose_subchannel)
2190 2028 : _console << "Maximum estimated Wij: " << _max_sumWij << std::endl;
2191 :
2192 2694 : LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
2193 : _Wij_loc_vec, _Wij, first_node, last_node, _n_gaps));
2194 2694 : LibmeshPetscCall(VecAbs(_Wij_loc_vec));
2195 2694 : LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
2196 : _Wij_old_loc_vec, _Wij_old, first_node, last_node, _n_gaps));
2197 2694 : LibmeshPetscCall(VecAbs(_Wij_old_loc_vec));
2198 2694 : LibmeshPetscCall(VecAXPY(_Wij_loc_vec, -1.0, _Wij_old_loc_vec));
2199 : PetscScalar relax_factor;
2200 2694 : LibmeshPetscCall(VecAbs(_Wij_loc_vec));
2201 : #if !PETSC_VERSION_LESS_THAN(3, 16, 0)
2202 2694 : LibmeshPetscCall(VecMean(_Wij_loc_vec, &relax_factor));
2203 : #else
2204 : VecSum(_Wij_loc_vec, &relax_factor);
2205 : relax_factor /= _block_size * _n_gaps;
2206 : #endif
2207 2694 : relax_factor = relax_factor / _max_sumWij + 0.5;
2208 2694 : if (_verbose_subchannel)
2209 2028 : _console << "Relax base value: " << relax_factor << std::endl;
2210 :
2211 : PetscScalar resistance_relaxation = 0.9;
2212 2694 : _added_K = _max_sumWij / min_mdot;
2213 2694 : if (_verbose_subchannel)
2214 2028 : _console << "New cross resistance: " << _added_K << std::endl;
2215 2694 : _added_K = (_added_K * resistance_relaxation + (1.0 - resistance_relaxation) * _added_K_old) *
2216 : relax_factor;
2217 2694 : if (_verbose_subchannel)
2218 2028 : _console << "Relaxed cross resistance: " << _added_K << std::endl;
2219 2694 : if (_added_K < 10 && _added_K >= 1.0)
2220 540 : _added_K = 1.0; //(1.0 - resistance_relaxation);
2221 2694 : if (_added_K < 1.0 && _added_K >= 0.1)
2222 630 : _added_K = 0.5;
2223 2694 : if (_added_K < 0.1 && _added_K >= 0.01)
2224 1104 : _added_K = 1. / 3.;
2225 2694 : if (_added_K < 1e-2 && _added_K >= 1e-3)
2226 0 : _added_K = 0.1;
2227 : if (_added_K < 1e-3)
2228 : _added_K = 1.0 * _added_K;
2229 2694 : if (_verbose_subchannel)
2230 2028 : _console << "Actual added cross resistance: " << _added_K << std::endl;
2231 2694 : LibmeshPetscCall(VecScale(unity_vec_Wij, _added_K));
2232 2694 : _added_K_old = _added_K;
2233 :
2234 : // Adding cross resistances
2235 2694 : LibmeshPetscCall(MatDiagonalSet(mat_array[2 * Q + 2], unity_vec_Wij, ADD_VALUES));
2236 2694 : LibmeshPetscCall(VecDestroy(&mdot_estimate));
2237 2694 : LibmeshPetscCall(VecDestroy(&pmat_diag));
2238 2694 : LibmeshPetscCall(VecDestroy(&unity_vec));
2239 2694 : LibmeshPetscCall(VecDestroy(&p_estimate));
2240 2694 : LibmeshPetscCall(VecDestroy(&sol_holder_P));
2241 2694 : LibmeshPetscCall(VecDestroy(&diag_Wij_loc));
2242 2694 : LibmeshPetscCall(VecDestroy(&unity_vec_Wij));
2243 2694 : LibmeshPetscCall(VecDestroy(&Wij_estimate));
2244 2694 : LibmeshPetscCall(VecDestroy(&sumWij_loc));
2245 2694 : LibmeshPetscCall(VecDestroy(&_Wij_loc_vec));
2246 2694 : LibmeshPetscCall(VecDestroy(&_Wij_old_loc_vec));
2247 :
2248 : // Auto-computing relaxation factors
2249 : PetscScalar relaxation_factor_mdot, relaxation_factor_P, relaxation_factor_Wij;
2250 2694 : relaxation_factor_mdot = 1.0;
2251 2694 : relaxation_factor_P = 1.0; // std::exp(-5.0);
2252 2694 : relaxation_factor_Wij = 0.1;
2253 :
2254 2694 : if (_verbose_subchannel)
2255 : {
2256 2028 : _console << "Relax mdot: " << relaxation_factor_mdot << std::endl;
2257 2028 : _console << "Relax P: " << relaxation_factor_P << std::endl;
2258 2028 : _console << "Relax Wij: " << relaxation_factor_Wij << std::endl;
2259 : }
2260 :
2261 : PetscInt field_num = 0;
2262 : Vec diag_mdot;
2263 2694 : LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_mdot));
2264 2694 : LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_mdot));
2265 2694 : LibmeshPetscCall(VecScale(diag_mdot, 1.0 / relaxation_factor_mdot));
2266 2694 : LibmeshPetscCall(
2267 : MatDiagonalSet(mat_array[Q * field_num + field_num], diag_mdot, INSERT_VALUES));
2268 2694 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2269 : _prod, *_mdot_soln, first_node, last_node, _n_channels));
2270 2694 : LibmeshPetscCall(VecScale(diag_mdot, (1.0 - relaxation_factor_mdot)));
2271 2694 : LibmeshPetscCall(VecPointwiseMult(_prod, _prod, diag_mdot));
2272 2694 : LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _prod));
2273 2694 : LibmeshPetscCall(VecDestroy(&diag_mdot));
2274 :
2275 2694 : if (_verbose_subchannel)
2276 2028 : _console << "mdot relaxed" << std::endl;
2277 :
2278 : field_num = 1;
2279 : Vec diag_P;
2280 2694 : LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_P));
2281 2694 : LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_P));
2282 2694 : LibmeshPetscCall(VecScale(diag_P, 1.0 / relaxation_factor_P));
2283 2694 : LibmeshPetscCall(MatDiagonalSet(mat_array[Q * field_num + field_num], diag_P, INSERT_VALUES));
2284 2694 : if (_verbose_subchannel)
2285 2028 : _console << "Mat assembled" << std::endl;
2286 2694 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2287 : _prod, *_P_soln, first_node, last_node, _n_channels));
2288 2694 : LibmeshPetscCall(VecScale(diag_P, (1.0 - relaxation_factor_P)));
2289 2694 : LibmeshPetscCall(VecPointwiseMult(_prod, _prod, diag_P));
2290 2694 : LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _prod));
2291 2694 : LibmeshPetscCall(VecDestroy(&diag_P));
2292 :
2293 2694 : if (_verbose_subchannel)
2294 2028 : _console << "P relaxed" << std::endl;
2295 :
2296 : field_num = 2;
2297 : Vec diag_Wij;
2298 2694 : LibmeshPetscCall(VecDuplicate(vec_array[field_num], &diag_Wij));
2299 2694 : LibmeshPetscCall(MatGetDiagonal(mat_array[Q * field_num + field_num], diag_Wij));
2300 2694 : LibmeshPetscCall(VecScale(diag_Wij, 1.0 / relaxation_factor_Wij));
2301 2694 : LibmeshPetscCall(MatDiagonalSet(mat_array[Q * field_num + field_num], diag_Wij, INSERT_VALUES));
2302 2694 : LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
2303 : _Wij_vec, _Wij, first_node, last_node, _n_gaps));
2304 2694 : LibmeshPetscCall(VecScale(diag_Wij, (1.0 - relaxation_factor_Wij)));
2305 2694 : LibmeshPetscCall(VecPointwiseMult(_Wij_vec, _Wij_vec, diag_Wij));
2306 2694 : LibmeshPetscCall(VecAXPY(vec_array[field_num], 1.0, _Wij_vec));
2307 2694 : LibmeshPetscCall(VecDestroy(&diag_Wij));
2308 :
2309 2694 : if (_verbose_subchannel)
2310 2028 : _console << "Wij relaxed" << std::endl;
2311 : }
2312 2694 : if (_verbose_subchannel)
2313 2028 : _console << "Linear solver relaxed." << std::endl;
2314 :
2315 : // Creating nested matrices
2316 2694 : LibmeshPetscCall(MatCreateNest(PETSC_COMM_WORLD, Q, NULL, Q, NULL, mat_array.data(), &A_nest));
2317 2694 : LibmeshPetscCall(VecCreateNest(PETSC_COMM_WORLD, Q, NULL, vec_array.data(), &b_nest));
2318 2694 : if (_verbose_subchannel)
2319 2028 : _console << "Nested system created." << std::endl;
2320 :
2321 : /// Setting up linear solver
2322 : // Creating linear solver
2323 2694 : LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
2324 2694 : LibmeshPetscCall(KSPSetType(ksp, KSPFGMRES));
2325 : // Setting KSP operators
2326 2694 : LibmeshPetscCall(KSPSetOperators(ksp, A_nest, A_nest));
2327 : // Set KSP and PC options
2328 2694 : LibmeshPetscCall(KSPGetPC(ksp, &pc));
2329 2694 : LibmeshPetscCall(PCSetType(pc, PCFIELDSPLIT));
2330 2694 : LibmeshPetscCall(KSPSetTolerances(ksp, _rtol, _atol, _dtol, _maxit));
2331 : // Splitting fields
2332 2694 : std::vector<IS> rows(Q);
2333 : // IS rows[Q];
2334 : PetscInt M = 0;
2335 2694 : LibmeshPetscCall(MatNestGetISs(A_nest, rows.data(), NULL));
2336 11754 : for (PetscInt j = 0; j < Q; ++j)
2337 : {
2338 : IS expand1;
2339 9060 : LibmeshPetscCall(ISDuplicate(rows[M], &expand1));
2340 9060 : M += 1;
2341 9060 : LibmeshPetscCall(PCFieldSplitSetIS(pc, NULL, expand1));
2342 9060 : LibmeshPetscCall(ISDestroy(&expand1));
2343 : }
2344 2694 : if (_verbose_subchannel)
2345 2028 : _console << "Linear solver assembled." << std::endl;
2346 :
2347 : /// Solving
2348 2694 : LibmeshPetscCall(VecDuplicate(b_nest, &x_nest));
2349 2694 : LibmeshPetscCall(VecSet(x_nest, 0.0));
2350 2694 : LibmeshPetscCall(KSPSolve(ksp, b_nest, x_nest));
2351 :
2352 : /// Destroying solver elements
2353 2694 : LibmeshPetscCall(VecDestroy(&b_nest));
2354 2694 : LibmeshPetscCall(MatDestroy(&A_nest));
2355 2694 : LibmeshPetscCall(KSPDestroy(&ksp));
2356 33786 : for (PetscInt i = 0; i < Q * Q; i++)
2357 : {
2358 31092 : LibmeshPetscCall(MatDestroy(&mat_array[i]));
2359 : }
2360 11754 : for (PetscInt i = 0; i < Q; i++)
2361 : {
2362 9060 : LibmeshPetscCall(VecDestroy(&vec_array[i]));
2363 : }
2364 2694 : if (_verbose_subchannel)
2365 2028 : _console << "Solver elements destroyed." << std::endl;
2366 :
2367 : /// Recovering the solutions
2368 : Vec sol_mdot, sol_p, sol_Wij;
2369 2694 : if (_verbose_subchannel)
2370 2028 : _console << "Vectors created." << std::endl;
2371 : PetscInt num_vecs;
2372 : Vec * loc_vecs;
2373 2694 : LibmeshPetscCall(VecNestGetSubVecs(x_nest, &num_vecs, &loc_vecs));
2374 2694 : if (_verbose_subchannel)
2375 2028 : _console << "Starting extraction." << std::endl;
2376 2694 : LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &sol_mdot));
2377 2694 : LibmeshPetscCall(VecCopy(loc_vecs[0], sol_mdot));
2378 2694 : LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &sol_p));
2379 2694 : LibmeshPetscCall(VecCopy(loc_vecs[1], sol_p));
2380 2694 : LibmeshPetscCall(VecDuplicate(_cmc_sys_Wij_rhs, &sol_Wij));
2381 2694 : LibmeshPetscCall(VecCopy(loc_vecs[2], sol_Wij));
2382 2694 : if (_verbose_subchannel)
2383 2028 : _console << "Getting individual field solutions from coupled solver." << std::endl;
2384 :
2385 : /// Assigning the solutions to arrays
2386 : PetscScalar * sol_p_array;
2387 2694 : LibmeshPetscCall(VecGetArray(sol_p, &sol_p_array));
2388 : PetscScalar * sol_Wij_array;
2389 2694 : LibmeshPetscCall(VecGetArray(sol_Wij, &sol_Wij_array));
2390 :
2391 : /// Populating Mass flow
2392 2694 : LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
2393 : sol_mdot, *_mdot_soln, first_node, last_node, _n_channels));
2394 :
2395 : /// Populating Pressure
2396 93420 : for (unsigned int iz = last_node; iz > first_node - 1; iz--)
2397 : {
2398 90726 : auto iz_ind = iz - first_node;
2399 4251210 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
2400 : {
2401 4160484 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
2402 4160484 : PetscScalar value = sol_p_array[iz_ind * _n_channels + i_ch];
2403 4160484 : _P_soln->set(node_in, value);
2404 : }
2405 : }
2406 :
2407 : /// Populating Crossflow
2408 2694 : LibmeshPetscCall(populateDenseFromVector<libMesh::DenseMatrix<Real>>(
2409 : sol_Wij, _Wij, first_node, last_node - 1, _n_gaps));
2410 :
2411 : /// Populating Enthalpy
2412 2694 : if (_monolithic_thermal_bool)
2413 : {
2414 : if (lag_block_thermal_solve)
2415 : {
2416 : KSP ksploc;
2417 : PC pc;
2418 : Vec sol;
2419 978 : LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol));
2420 978 : LibmeshPetscCall(KSPCreate(PETSC_COMM_WORLD, &ksploc));
2421 978 : LibmeshPetscCall(KSPSetOperators(ksploc, _hc_sys_h_mat, _hc_sys_h_mat));
2422 978 : LibmeshPetscCall(KSPGetPC(ksploc, &pc));
2423 978 : LibmeshPetscCall(PCSetType(pc, PCJACOBI));
2424 978 : LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
2425 978 : LibmeshPetscCall(KSPSetFromOptions(ksploc));
2426 978 : LibmeshPetscCall(KSPSolve(ksploc, _hc_sys_h_rhs, sol));
2427 : PetscScalar * xx;
2428 978 : LibmeshPetscCall(VecGetArray(sol, &xx));
2429 30798 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
2430 : {
2431 29820 : auto iz_ind = iz - first_node;
2432 1275420 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
2433 : {
2434 1245600 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
2435 1245600 : auto h_out = xx[iz_ind * _n_channels + i_ch];
2436 1245600 : if (h_out < 0)
2437 : {
2438 0 : mooseError(name(),
2439 : " : Calculation of negative Enthalpy h_out = : ",
2440 : h_out,
2441 : " Axial Level= : ",
2442 : iz);
2443 : }
2444 1245600 : _h_soln->set(node_out, h_out);
2445 : }
2446 : }
2447 978 : LibmeshPetscCall(KSPDestroy(&ksploc));
2448 978 : LibmeshPetscCall(VecDestroy(&sol));
2449 : }
2450 : else
2451 : {
2452 : Vec sol_h;
2453 : LibmeshPetscCall(VecDuplicate(_hc_sys_h_rhs, &sol_h));
2454 : LibmeshPetscCall(VecCopy(loc_vecs[3], sol_h));
2455 : PetscScalar * sol_h_array;
2456 : LibmeshPetscCall(VecGetArray(sol_h, &sol_h_array));
2457 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
2458 : {
2459 : auto iz_ind = iz - first_node;
2460 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
2461 : {
2462 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
2463 : auto h_out = sol_h_array[iz_ind * _n_channels + i_ch];
2464 : if (h_out < 0)
2465 : {
2466 : mooseError(name(),
2467 : " : Calculation of negative Enthalpy h_out = : ",
2468 : h_out,
2469 : " Axial Level= : ",
2470 : iz);
2471 : }
2472 : _h_soln->set(node_out, h_out);
2473 : }
2474 : }
2475 : LibmeshPetscCall(VecDestroy(&sol_h));
2476 : }
2477 : }
2478 :
2479 : /// Populating sum_Wij
2480 2694 : LibmeshPetscCall(MatMult(_mc_sumWij_mat, sol_Wij, _prod));
2481 2694 : LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
2482 : _prod, *_SumWij_soln, first_node, last_node, _n_channels));
2483 :
2484 : Vec sumWij_loc;
2485 2694 : LibmeshPetscCall(createPetscVector(sumWij_loc, _block_size * _n_channels));
2486 2694 : LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
2487 : _prod, *_SumWij_soln, first_node, last_node, _n_channels));
2488 :
2489 2694 : LibmeshPetscCall(VecAbs(_prod));
2490 2694 : LibmeshPetscCall(VecMax(_prod, NULL, &_max_sumWij_new));
2491 2694 : if (_verbose_subchannel)
2492 2028 : _console << "Maximum estimated Wij new: " << _max_sumWij_new << std::endl;
2493 2694 : _correction_factor = _max_sumWij_new / _max_sumWij;
2494 2694 : if (_verbose_subchannel)
2495 2028 : _console << "Correction factor: " << _correction_factor << std::endl;
2496 2694 : if (_verbose_subchannel)
2497 2028 : _console << "Solutions assigned to MOOSE variables." << std::endl;
2498 :
2499 : /// Destroying arrays
2500 2694 : LibmeshPetscCall(VecDestroy(&x_nest));
2501 2694 : LibmeshPetscCall(VecDestroy(&sol_mdot));
2502 2694 : LibmeshPetscCall(VecDestroy(&sol_p));
2503 2694 : LibmeshPetscCall(VecDestroy(&sol_Wij));
2504 2694 : LibmeshPetscCall(VecDestroy(&sumWij_loc));
2505 2694 : if (_verbose_subchannel)
2506 2028 : _console << "Solutions destroyed." << std::endl;
2507 :
2508 2694 : PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
2509 2694 : }
2510 :
2511 : void
2512 292 : SubChannel1PhaseProblem::externalSolve()
2513 : {
2514 292 : _console << "Executing subchannel solver\n";
2515 292 : _dt = (isTransient() ? dt() : _one);
2516 292 : _TR = isTransient();
2517 292 : initializeSolution();
2518 292 : if (_verbose_subchannel)
2519 150 : _console << "Solution initialized" << std::endl;
2520 292 : Real P_error = 1.0;
2521 292 : unsigned int P_it = 0;
2522 : unsigned int P_it_max;
2523 :
2524 292 : if (_segregated_bool)
2525 148 : P_it_max = 20 * _n_blocks;
2526 : else
2527 144 : P_it_max = 100;
2528 :
2529 292 : if ((_n_blocks == 1) && (_segregated_bool))
2530 148 : P_it_max = 5;
2531 :
2532 1674 : while ((P_error > _P_tol && P_it < P_it_max))
2533 : {
2534 1382 : P_it += 1;
2535 1382 : if (P_it == P_it_max && _n_blocks != 1)
2536 : {
2537 0 : _console << "Reached maximum number of axial pressure iterations" << std::endl;
2538 0 : _converged = false;
2539 : }
2540 1382 : _console << "Solving Outer Iteration : " << P_it << std::endl;
2541 1382 : auto P_L2norm_old_axial = _P_soln->L2norm();
2542 2764 : for (unsigned int iblock = 0; iblock < _n_blocks; iblock++)
2543 : {
2544 1382 : int last_level = (iblock + 1) * _block_size;
2545 1382 : int first_level = iblock * _block_size + 1;
2546 1382 : Real T_block_error = 1.0;
2547 : auto T_it = 0;
2548 1382 : _console << "Solving Block: " << iblock << " From first level: " << first_level
2549 1382 : << " to last level: " << last_level << std::endl;
2550 4452 : while (T_block_error > _T_tol && T_it < _T_maxit)
2551 : {
2552 3070 : T_it += 1;
2553 3070 : if (T_it == _T_maxit)
2554 : {
2555 0 : _console << "Reached maximum number of temperature iterations for block: " << iblock
2556 0 : << std::endl;
2557 0 : _converged = false;
2558 : }
2559 3070 : auto T_L2norm_old_block = _T_soln->L2norm();
2560 :
2561 3070 : if (_segregated_bool)
2562 : {
2563 376 : computeWijFromSolve(iblock);
2564 :
2565 376 : if (_compute_power)
2566 : {
2567 248 : computeh(iblock);
2568 248 : computeT(iblock);
2569 : }
2570 : }
2571 : else
2572 : {
2573 2694 : LibmeshPetscCall(implicitPetscSolve(iblock));
2574 2694 : computeWijPrime(iblock);
2575 2694 : if (_verbose_subchannel)
2576 2028 : _console << "Done with main solve." << std::endl;
2577 2694 : if (_monolithic_thermal_bool)
2578 : {
2579 : // Enthalpy is already solved from the monolithic solve
2580 978 : computeT(iblock);
2581 : }
2582 : else
2583 : {
2584 1716 : if (_verbose_subchannel)
2585 1218 : _console << "Starting thermal solve." << std::endl;
2586 1716 : if (_compute_power)
2587 : {
2588 1704 : computeh(iblock);
2589 1704 : computeT(iblock);
2590 : }
2591 1716 : if (_verbose_subchannel)
2592 1218 : _console << "Done with thermal solve." << std::endl;
2593 : }
2594 : }
2595 :
2596 3070 : if (_verbose_subchannel)
2597 2118 : _console << "Start updating thermophysical properties." << std::endl;
2598 :
2599 3070 : if (_compute_density)
2600 3058 : computeRho(iblock);
2601 :
2602 3070 : if (_compute_viscosity)
2603 3058 : computeMu(iblock);
2604 :
2605 3070 : if (_verbose_subchannel)
2606 2118 : _console << "Done updating thermophysical properties." << std::endl;
2607 :
2608 : // We must do a global assembly to make sure data is parallel consistent before we do things
2609 : // like compute L2 norms
2610 3070 : _aux->solution().close();
2611 :
2612 3070 : auto T_L2norm_new = _T_soln->L2norm();
2613 3070 : T_block_error =
2614 3070 : std::abs((T_L2norm_new - T_L2norm_old_block) / (T_L2norm_old_block + 1E-14));
2615 3070 : _console << "T_block_error: " << T_block_error << std::endl;
2616 : }
2617 : }
2618 1382 : auto P_L2norm_new_axial = _P_soln->L2norm();
2619 1382 : P_error =
2620 1382 : std::abs((P_L2norm_new_axial - P_L2norm_old_axial) / (P_L2norm_old_axial + _P_out + 1E-14));
2621 1382 : _console << "P_error :" << P_error << std::endl;
2622 1382 : if (_verbose_subchannel)
2623 : {
2624 888 : _console << "Iteration: " << P_it << std::endl;
2625 888 : _console << "Maximum iterations: " << P_it_max << std::endl;
2626 : }
2627 : }
2628 : // update old crossflow matrix
2629 292 : _Wij_old = _Wij;
2630 292 : _console << "Finished executing subchannel solver\n";
2631 :
2632 : /// Assigning temperature to the fuel pins
2633 292 : if (_pin_mesh_exist)
2634 : {
2635 98 : _console << "Commencing calculation of Pin surface temperature \n";
2636 3034 : for (unsigned int i_pin = 0; i_pin < _n_pins; i_pin++)
2637 : {
2638 60066 : for (unsigned int iz = 0; iz < _n_cells + 1; ++iz)
2639 : {
2640 57130 : const auto * pin_node = _subchannel_mesh.getPinNode(i_pin, iz);
2641 : Real sumTemp = 0.0;
2642 : Real rod_counter = 0.0;
2643 : // Calculate sum of pin surface temperatures that the channels around the pin see
2644 346562 : for (auto i_ch : _subchannel_mesh.getPinChannels(i_pin))
2645 : {
2646 289432 : const auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
2647 289432 : auto mu = (*_mu_soln)(node);
2648 289432 : auto S = (*_S_flow_soln)(node);
2649 289432 : auto w_perim = (*_w_perim_soln)(node);
2650 289432 : auto Dh_i = 4.0 * S / w_perim;
2651 289432 : auto Re = (((*_mdot_soln)(node) / S) * Dh_i / mu);
2652 289432 : auto k = _fp->k_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
2653 289432 : auto cp = _fp->cp_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
2654 289432 : auto Pr = (*_mu_soln)(node)*cp / k;
2655 289432 : auto Nu = 0.023 * std::pow(Re, 0.8) * std::pow(Pr, 0.4);
2656 289432 : auto hw = Nu * k / Dh_i;
2657 289432 : if ((*_Dpin_soln)(pin_node) <= 0)
2658 0 : mooseError("Dpin should not be null or negative when computing pin powers: ",
2659 0 : (*_Dpin_soln)(pin_node));
2660 289432 : sumTemp +=
2661 289432 : (*_q_prime_soln)(pin_node) / ((*_Dpin_soln)(pin_node)*M_PI * hw) + (*_T_soln)(node);
2662 289432 : rod_counter += 1.0;
2663 : }
2664 57130 : if (rod_counter > 0)
2665 57130 : _Tpin_soln->set(pin_node, sumTemp / rod_counter);
2666 : else
2667 0 : mooseError("Pin was not found for pin index: " + std::to_string(i_pin));
2668 : }
2669 : }
2670 : }
2671 :
2672 : /// Assigning temperatures to duct
2673 292 : if (_duct_mesh_exist)
2674 : {
2675 30 : _console << "Commencing calculation of duct surface temperature " << std::endl;
2676 30 : auto duct_nodes = _subchannel_mesh.getDuctNodes();
2677 6186 : for (Node * dn : duct_nodes)
2678 : {
2679 6156 : auto * node_chan = _subchannel_mesh.getChannelNodeFromDuct(dn);
2680 6156 : auto mu = (*_mu_soln)(node_chan);
2681 6156 : auto S = (*_S_flow_soln)(node_chan);
2682 6156 : auto w_perim = (*_w_perim_soln)(node_chan);
2683 6156 : auto Dh_i = 4.0 * S / w_perim;
2684 6156 : auto Re = (((*_mdot_soln)(node_chan) / S) * Dh_i / mu);
2685 6156 : auto k = _fp->k_from_p_T((*_P_soln)(node_chan) + _P_out, (*_T_soln)(node_chan));
2686 6156 : auto cp = _fp->cp_from_p_T((*_P_soln)(node_chan) + _P_out, (*_T_soln)(node_chan));
2687 6156 : auto Pr = (*_mu_soln)(node_chan)*cp / k;
2688 6156 : auto Nu = 0.023 * std::pow(Re, 0.8) * std::pow(Pr, 0.4);
2689 6156 : auto hw = Nu * k / Dh_i;
2690 6156 : auto T_chan = (*_q_prime_duct_soln)(dn) / (_subchannel_mesh.getPinDiameter() * M_PI * hw) +
2691 6156 : (*_T_soln)(node_chan);
2692 6156 : _Tduct_soln->set(dn, T_chan);
2693 : }
2694 30 : }
2695 292 : _aux->solution().close();
2696 292 : _aux->update();
2697 :
2698 : Real power_in = 0.0;
2699 : Real power_out = 0.0;
2700 292 : Real Total_surface_area = 0.0;
2701 : Real Total_wetted_perimeter = 0.0;
2702 292 : Real mass_flow_in = 0.0;
2703 : Real mass_flow_out = 0.0;
2704 10740 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
2705 : {
2706 10448 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, 0);
2707 10448 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, _n_cells);
2708 10448 : Total_surface_area += (*_S_flow_soln)(node_in);
2709 10448 : Total_wetted_perimeter += (*_w_perim_soln)(node_in);
2710 10448 : power_in += (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
2711 10448 : power_out += (*_mdot_soln)(node_out) * (*_h_soln)(node_out);
2712 10448 : mass_flow_in += (*_mdot_soln)(node_in);
2713 10448 : mass_flow_out += (*_mdot_soln)(node_out);
2714 : }
2715 292 : auto h_bulk_out = power_out / mass_flow_out;
2716 292 : auto T_bulk_out = _fp->T_from_p_h(_P_out, h_bulk_out);
2717 :
2718 292 : Real bulk_Dh = 4.0 * Total_surface_area / Total_wetted_perimeter;
2719 292 : Real inlet_mu = (*_mu_soln)(_subchannel_mesh.getChannelNode(0, 0));
2720 292 : Real bulk_Re = mass_flow_in * bulk_Dh / (inlet_mu * Total_surface_area);
2721 292 : if (_verbose_subchannel)
2722 : {
2723 150 : _console << " ======================================= " << std::endl;
2724 150 : _console << " ======== Subchannel Print Outs ======== " << std::endl;
2725 150 : _console << " ======================================= " << std::endl;
2726 150 : _console << "Total flow area :" << Total_surface_area << " m^2" << std::endl;
2727 150 : _console << "Assembly hydraulic diameter :" << bulk_Dh << " m" << std::endl;
2728 150 : _console << "Assembly Re number :" << bulk_Re << " [-]" << std::endl;
2729 150 : _console << "Bulk coolant temperature at outlet :" << T_bulk_out << " K" << std::endl;
2730 150 : _console << "Power added to coolant is : " << power_out - power_in << " Watt" << std::endl;
2731 150 : _console << "Mass flow rate in is : " << mass_flow_in << " kg/sec" << std::endl;
2732 150 : _console << "Mass balance is : " << mass_flow_out - mass_flow_in << " kg/sec" << std::endl;
2733 150 : _console << "User defined outlet pressure is : " << _P_out << " Pa" << std::endl;
2734 150 : _console << " ======================================= " << std::endl;
2735 : }
2736 :
2737 292 : if (MooseUtils::absoluteFuzzyLessEqual((power_out - power_in), -1.0))
2738 0 : mooseWarning(
2739 0 : "Energy conservation equation might not be solved correctly, Power added to coolant: " +
2740 0 : std::to_string(power_out - power_in) + " Watt ");
2741 292 : }
2742 :
2743 : void
2744 584 : SubChannel1PhaseProblem::syncSolutions(Direction /*direction*/)
2745 : {
2746 584 : }
|
