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