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