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