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