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