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