SubChannel1PhaseProblem.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "SubChannel1PhaseProblem.h"
#include "SystemBase.h"
#include "libmesh/petsc_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include <iostream>
#include <cmath>
#include "AuxiliarySystem.h"
#include "SCM.h"
#include "SinglePhaseFluidProperties.h"
#include "SCMFrictionClosureBase.h"
#include "SCMHTCClosureBase.h"
#include "SCMMixingClosureBase.h"
#include "TransientBase.h"
#include "ImplicitEuler.h"
 
struct Ctx
{
  int iblock;
  SubChannel1PhaseProblem * schp;
};
 
PetscErrorCode
formFunction(SNES, Vec x, Vec f, void * ctx)
{
  const PetscScalar * xx;
  PetscScalar * ff;
  PetscInt size;
 
  PetscFunctionBegin;
  Ctx * cc = static_cast<Ctx *>(ctx);
  LibmeshPetscCallQ(VecGetSize(x, &size));
 
  libMesh::DenseVector<Real> solution_seed(size, 0.0);
  LibmeshPetscCallQ(VecGetArrayRead(x, &xx));
  for (PetscInt i = 0; i < size; i++)
    solution_seed(i) = xx[i];
 
  LibmeshPetscCallQ(VecRestoreArrayRead(x, &xx));
 
  libMesh::DenseVector<Real> Wij_residual_vector =
      cc->schp->residualFunction(cc->iblock, solution_seed);
 
  LibmeshPetscCallQ(VecGetArray(f, &ff));
  for (int i = 0; i < size; i++)
    ff[i] = Wij_residual_vector(i);
 
  LibmeshPetscCallQ(VecRestoreArray(f, &ff));
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}
 
InputParameters
SubChannel1PhaseProblem::validParams()
{
  // Enumerations
  MooseEnum schemes("upwind downwind central_difference exponential", "central_difference");
  MooseEnum gravity_direction("counter_flow co_flow none", "counter_flow");
 
  // Inputs
  InputParameters params = ExternalProblem::validParams();
  params += PostprocessorInterface::validParams();
  params.addClassDescription("Base class of the subchannel solvers");
  params.addParam<unsigned int>("n_blocks", 1, "The number of blocks in the axial direction");
  params.addParam<Real>("P_tol", 1e-6, "Pressure tolerance");
  params.addRangeCheckedParam<int>(
      "P_maxit",
      0,
      "P_maxit >= 0",
      "Maximum number of pressure iterations; zero selects the solver's automatic limit");
  params.addParam<Real>("T_tol", 1e-6, "Temperature tolerance");
  params.addParam<int>("T_maxit", 100, "Maximum number of iterations for inner temperature loop");
  params.addRangeCheckedParam<Real>(
      "T_relaxation",
      1.0,
      "T_relaxation > 0 & T_relaxation <= 1",
      "Relaxation factor for temperature updates in the inner thermal-hydraulic iteration");
  params.addRangeCheckedParam<unsigned int>(
      "enthalpy_subcycles",
      1,
      "enthalpy_subcycles > 0",
      "Number of enthalpy, temperature, and property updates performed per flow solve");
  params.addRangeCheckedParam<Real>(
      "mass_flow_equation_relaxation",
      1.0,
      "mass_flow_equation_relaxation > 0 & mass_flow_equation_relaxation <= 1",
      "Equation relaxation factor for mass flow rate in the implicit non-segregated solve");
  params.addRangeCheckedParam<Real>(
      "pressure_equation_relaxation",
      1.0,
      "pressure_equation_relaxation > 0 & pressure_equation_relaxation <= 1",
      "Equation relaxation factor for pressure in the implicit non-segregated solve");
  params.addRangeCheckedParam<Real>(
      "crossflow_equation_relaxation",
      0.1,
      "crossflow_equation_relaxation > 0 & crossflow_equation_relaxation <= 1",
      "Equation relaxation factor for crossflow in the implicit non-segregated solve");
  params.addRangeCheckedParam<Real>(
      "mass_flow_relaxation",
      1.0,
      "mass_flow_relaxation > 0 & mass_flow_relaxation <= 1",
      "Post-solve relaxation factor for mass flow rate updates in the implicit non-segregated "
      "solve");
  params.addRangeCheckedParam<Real>(
      "pressure_relaxation",
      1.0,
      "pressure_relaxation > 0 & pressure_relaxation <= 1",
      "Post-solve relaxation factor for pressure updates in the implicit non-segregated solve");
  params.addRangeCheckedParam<Real>(
      "crossflow_relaxation",
      1.0,
      "crossflow_relaxation > 0 & crossflow_relaxation <= 1",
      "Post-solve relaxation factor for crossflow updates in the implicit non-segregated solve");
  params.addParam<PetscReal>("rtol", 1e-6, "Relative tolerance for ksp solver");
  params.addParam<PetscReal>("atol", 1e-6, "Absolute tolerance for ksp solver");
  params.addParam<PetscReal>("dtol", 1e5, "Divergence tolerance or ksp solver");
  params.addParam<PetscInt>("maxit", 1e4, "Maximum number of iterations for ksp solver");
  params.addParam<MooseEnum>(
      "interpolation_scheme",
      schemes,
      "Interpolation scheme used to evaluate values at the cell center. Default is "
      "central_difference");
  params.addParam<MooseEnum>(
      "gravity", gravity_direction, "Direction of gravity. Default is counter_flow");
  params.addParam<bool>(
      "implicit", false, "Boolean to define the use of explicit or implicit solution.");
  params.addDeprecatedParam<bool>(
      "staggered_pressure",
      false,
      "Boolean to define the use of staggered or collocated pressure.",
      "Pressure is a face value and the pressure drop is a cell value, so the collocated "
      "formulation is always used. This parameter is ignored.");
  params.addParam<bool>(
      "segregated", true, "Boolean to define whether to use a segregated solution.");
  params.addParam<bool>(
      "verbose_subchannel", false, "Boolean to print out information related to subchannel solve.");
  params.addRequiredParam<bool>("compute_density", "Flag that enables the calculation of density");
  params.addRequiredParam<bool>("compute_viscosity",
                                "Flag that enables the calculation of viscosity");
  params.addRequiredParam<bool>(
      "compute_power",
      "Flag that informs whether we solve the Enthalpy/Temperature equations or not");
  params.addRequiredParam<PostprocessorName>(
      "P_out",
      "The postprocessor (or scalar) that provides the absolute outlet pressure [Pa]. The solved "
      "pressure variable P is relative to this value.");
  params.addRequiredParam<UserObjectName>("fp", "Fluid properties user object name");
  params.addRequiredParam<UserObjectName>("friction_closure",
                                          "Closure computing the friction factor");
  params.addRequiredParam<UserObjectName>(
      "mixing_closure",
      "Closure computing the turbulent mixing, wire-induced "
      "mixing and sweep flow mixing parameter where applicable");
  params.addParam<UserObjectName>(
      "pin_HTC_closure", "Closure computing HTC on fuel pin (required if pin mesh exists).");
  params.addParam<UserObjectName>("duct_HTC_closure",
                                  "Closure computing HTC on duct (required if duct mesh exists).");
  params.addParam<bool>(
      "full_output", false, "Flag that enables the output of the maximum number of variables.");
  params.addDeprecatedParam<Real>("beta",
                                  "Thermal diffusion coefficient used in turbulent crossflow.",
                                  "Use closure system instead.");
  params.addDeprecatedParam<bool>(
      "constant_beta",
      true,
      "Boolean to define the use of a constant beta or beta correlation (Kim and Chung, 2001)",
      "Use closure system instead.");
 
  params.addParamNamesToGroup(
      "P_tol P_maxit T_tol T_maxit T_relaxation enthalpy_subcycles "
      "mass_flow_equation_relaxation pressure_equation_relaxation crossflow_equation_relaxation "
      "mass_flow_relaxation pressure_relaxation crossflow_relaxation rtol atol dtol maxit",
      "Solver tolerances and iterations");
  params.addParamNamesToGroup("implicit segregated interpolation_scheme", "Solution method");
  params.addParamNamesToGroup("fp friction_closure mixing_closure pin_HTC_closure duct_HTC_closure",
                              "Closures");
  params.addParamNamesToGroup("compute_density compute_viscosity compute_power gravity",
                              "Physics models");
  params.addParamNamesToGroup("verbose_subchannel full_output", "Output");
 
  return params;
}
 
SubChannel1PhaseProblem::SubChannel1PhaseProblem(const InputParameters & params)
  : ExternalProblem(params),
    PostprocessorInterface(this),
    _friction_args(/*i_ch=*/0, /*Re=*/1.0, /*S=*/0.0, /*w_perim=*/0.0),
    _nusselt_args(
        /*Re=*/1.0, /*Pr=*/1.0, std::numeric_limits<unsigned int>::max(), /*iz=*/0, /*i_ch=*/0),
    _P_out(getPostprocessorValue("P_out")),
    _fp(nullptr),
    _subchannel_mesh(SCM::getMesh<SubChannelMesh>(_mesh)),
    _n_blocks(getParam<unsigned int>("n_blocks")),
    _Wij(declareRestartableData<libMesh::DenseMatrix<Real>>("Wij")),
    _Wij_old(declareRestartableData<libMesh::DenseMatrix<Real>>("Wij_old")),
    _g_grav(9.81),
    _kij(_subchannel_mesh.getKij()),
    _one(1.0),
    _compute_density(getParam<bool>("compute_density")),
    _compute_viscosity(getParam<bool>("compute_viscosity")),
    _compute_power(getParam<bool>("compute_power")),
    _pin_mesh_exist(_subchannel_mesh.pinMeshExist()),
    _duct_mesh_exist(_subchannel_mesh.ductMeshExist()),
    _bulk_Re(1.0),
    _bulk_V(1.0),
    _P_tol(getParam<Real>("P_tol")),
    _P_maxit(getParam<int>("P_maxit")),
    _T_tol(getParam<Real>("T_tol")),
    _T_maxit(getParam<int>("T_maxit")),
    _T_relaxation(getParam<Real>("T_relaxation")),
    _enthalpy_subcycles(getParam<unsigned int>("enthalpy_subcycles")),
    _mass_flow_equation_relaxation(getParam<Real>("mass_flow_equation_relaxation")),
    _pressure_equation_relaxation(getParam<Real>("pressure_equation_relaxation")),
    _crossflow_equation_relaxation(getParam<Real>("crossflow_equation_relaxation")),
    _mass_flow_relaxation(getParam<Real>("mass_flow_relaxation")),
    _pressure_relaxation(getParam<Real>("pressure_relaxation")),
    _crossflow_relaxation(getParam<Real>("crossflow_relaxation")),
    _rtol(getParam<PetscReal>("rtol")),
    _atol(getParam<PetscReal>("atol")),
    _dtol(getParam<PetscReal>("dtol")),
    _maxit(getParam<PetscInt>("maxit")),
    _interpolation_scheme(getParam<MooseEnum>("interpolation_scheme")),
    _gravity_direction(getParam<MooseEnum>("gravity")),
    _dir_grav(computeGravityDir(_gravity_direction)),
    _implicit_bool(getParam<bool>("implicit")),
    _segregated_bool(getParam<bool>("segregated")),
    _verbose_subchannel(getParam<bool>("verbose_subchannel")),
    _friction_closure(nullptr),
    _mixing_closure(nullptr),
    _pin_HTC_closure(nullptr),
    _duct_HTC_closure(nullptr),
    _Tpin_soln(nullptr),
    _duct_heat_flux_soln(nullptr),
    _Tduct_soln(nullptr),
    _HTC_soln(nullptr)
{
  if (params.isParamSetByUser("beta") || params.isParamSetByUser("constant_beta"))
    paramError("beta",
               "You are using a deprecated parameter. Please use the mixing_closure system.");
  if (_pin_mesh_exist && !isParamValid("pin_HTC_closure"))
    paramError("pin_HTC_closure", "required when a pin mesh exists.");
  if (_duct_mesh_exist && !isParamValid("duct_HTC_closure"))
    paramError("duct_HTC_closure", "required when a duct mesh exists.");
  if (!_implicit_bool && !_segregated_bool)
    paramError("segregated", "A non-segregated solve requires 'implicit = true'.");
  if (_segregated_bool)
    for (const auto * relaxation_param : {"mass_flow_equation_relaxation",
                                          "pressure_equation_relaxation",
                                          "crossflow_equation_relaxation",
                                          "mass_flow_relaxation",
                                          "pressure_relaxation",
                                          "crossflow_relaxation"})
      if (params.isParamSetByUser(relaxation_param))
        paramError(relaxation_param,
                   "The '",
                   relaxation_param,
                   "' parameter is only used by the implicit non-segregated solve. Set "
                   "'implicit = true' and 'segregated = false' to use it.");
  // NOTE: The four quantities below are 0 for processor_id != 0
  _n_cells = _subchannel_mesh.getNumOfAxialCells();
  _n_gaps = _subchannel_mesh.getNumOfGapsPerLayer();
  _n_pins = _subchannel_mesh.getNumOfPins();
  _n_channels = _subchannel_mesh.getNumOfChannels();
  // NOTE: The four quantities above are 0 for processor_id != 0
  _z_grid = _subchannel_mesh.getZGrid();
  _block_size = _n_cells / _n_blocks;
  // Pressure drop (lives on subchannel nodes)
  _DP.resize(_n_channels, _n_cells + 1);
  _DP.zero();
  // Turbulent crossflow (stuff that live on the gaps)
  if (!_app.isRestarting() && !_app.isRecovering())
  {
    _Wij.resize(_n_gaps, _n_cells + 1);
    _Wij.zero();
    _Wij_old.resize(_n_gaps, _n_cells + 1);
    _Wij_old.zero();
  }
 
  _WijPrime.resize(_n_gaps, _n_cells + 1);
  _WijPrime.zero();
  _Wij_residual_matrix.resize(_n_gaps, _block_size);
  _Wij_residual_matrix.zero();
  _converged = true;
 
  // Mass conservation components
  LibmeshPetscCall(
      createPetscMatrix(_mc_sumWij_mat, _block_size * _n_channels, _block_size * _n_gaps));
  LibmeshPetscCall(createPetscVector(_Wij_vec, _block_size * _n_gaps));
  LibmeshPetscCall(createPetscVector(_prod, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_prodp, _block_size * _n_channels));
  LibmeshPetscCall(createPetscMatrix(
      _mc_axial_convection_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_mc_axial_convection_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscMatrix(
      _mc_density_pressure_mat, _block_size * _n_channels, _block_size * _n_channels));
 
  // Axial momentum conservation components
  LibmeshPetscCall(createPetscMatrix(
      _amc_turbulent_cross_flows_mat, _block_size * _n_gaps, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_amc_turbulent_cross_flows_rhs, _block_size * _n_gaps));
  LibmeshPetscCall(createPetscMatrix(
      _amc_time_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_amc_time_derivative_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscMatrix(
      _amc_advective_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_amc_advective_derivative_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscMatrix(
      _amc_cross_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_amc_cross_derivative_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscMatrix(
      _amc_friction_force_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_amc_friction_force_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_amc_gravity_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscMatrix(
      _amc_pressure_force_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_amc_pressure_force_rhs, _block_size * _n_channels));
  LibmeshPetscCall(
      createPetscMatrix(_amc_sys_mdot_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_amc_sys_mdot_rhs, _block_size * _n_channels));
 
  // Lateral momentum conservation components
  LibmeshPetscCall(
      createPetscMatrix(_cmc_time_derivative_mat, _block_size * _n_gaps, _block_size * _n_gaps));
  LibmeshPetscCall(createPetscVector(_cmc_time_derivative_rhs, _block_size * _n_gaps));
  LibmeshPetscCall(createPetscMatrix(
      _cmc_advective_derivative_mat, _block_size * _n_gaps, _block_size * _n_gaps));
  LibmeshPetscCall(createPetscVector(_cmc_advective_derivative_rhs, _block_size * _n_gaps));
  LibmeshPetscCall(
      createPetscMatrix(_cmc_friction_force_mat, _block_size * _n_gaps, _block_size * _n_gaps));
  LibmeshPetscCall(createPetscVector(_cmc_friction_force_rhs, _block_size * _n_gaps));
  LibmeshPetscCall(
      createPetscMatrix(_cmc_pressure_force_mat, _block_size * _n_gaps, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_cmc_pressure_force_rhs, _block_size * _n_gaps));
  LibmeshPetscCall(
      createPetscMatrix(_cmc_sys_Wij_mat, _block_size * _n_gaps, _block_size * _n_gaps));
  LibmeshPetscCall(createPetscVector(_cmc_sys_Wij_rhs, _block_size * _n_gaps));
 
  // Energy conservation components
  LibmeshPetscCall(createPetscMatrix(
      _hc_time_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_hc_time_derivative_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscMatrix(
      _hc_advective_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_hc_advective_derivative_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscMatrix(
      _hc_cross_derivative_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_hc_cross_derivative_rhs, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_hc_added_heat_rhs, _block_size * _n_channels));
  LibmeshPetscCall(
      createPetscMatrix(_hc_sys_h_mat, _block_size * _n_channels, _block_size * _n_channels));
  LibmeshPetscCall(createPetscVector(_hc_sys_h_rhs, _block_size * _n_channels));
 
  if ((_n_blocks == _n_cells) && _implicit_bool)
  {
    mooseError(name(),
               ": When implicit number of blocks can't be equal to number of cells. This will "
               "cause problems with the subchannel interpolation scheme.");
  }
}
 
void
SubChannel1PhaseProblem::initialSetup()
{
  // Create variables for output and storage. This comes before the base class setup, which runs
  // the user objects executed on initial, such as SCMPinPowerPostprocessor, that read them.
  _mdot_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::MASS_FLOW_RATE));
  _SumWij_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::SUM_CROSSFLOW));
  _P_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::PRESSURE));
  if (getParam<bool>("full_output"))
  {
    _DP_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::PRESSURE_DROP));
    _ff_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::FRICTION_FACTOR));
  }
  _h_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::ENTHALPY));
  _T_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::TEMPERATURE));
  if (_pin_mesh_exist)
  {
    _Tpin_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::PIN_TEMPERATURE));
    _Dpin_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::PIN_DIAMETER));
    _HTC_soln =
        std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::HEAT_TRANSFER_COEFFICIENT));
  }
  _rho_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::DENSITY));
  _mu_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::VISCOSITY));
  _S_flow_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::SURFACE_AREA));
  _w_perim_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::WETTED_PERIMETER));
  _q_prime_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::LINEAR_HEAT_RATE));
  _displacement_soln =
      std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::DISPLACEMENT));
  if (_duct_mesh_exist)
  {
    _duct_heat_flux_soln =
        std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::DUCT_HEAT_FLUX));
    _Tduct_soln = std::make_unique<SolutionHandle>(getVariable(0, SubChannelApp::DUCT_TEMPERATURE));
  }
 
  ExternalProblem::initialSetup();
 
  _fp = &getUserObject<SinglePhaseFluidProperties>(getParam<UserObjectName>("fp"));
  _friction_closure =
      &getUserObject<SCMFrictionClosureBase>(getParam<UserObjectName>("friction_closure"));
  _mixing_closure =
      &getUserObject<SCMMixingClosureBase>(getParam<UserObjectName>("mixing_closure"));
 
  _CT = _mixing_closure->getCT();
 
  if (_pin_mesh_exist)
    _pin_HTC_closure =
        &getUserObject<SCMHTCClosureBase>(getParam<UserObjectName>("pin_HTC_closure"));
  if (_duct_mesh_exist)
    _duct_HTC_closure =
        &getUserObject<SCMHTCClosureBase>(getParam<UserObjectName>("duct_HTC_closure"));
}
 
void
SubChannel1PhaseProblem::computeBulkReynoldsNumber()
{
  if (processor_id() == 0)
  {
    Real viscosity_in = 0.0;
    Real volumetric_flow_in = 0.0;
    Real mass_flow_in = 0.0;
    for (const auto i_ch : make_range(_n_channels))
    {
      auto * node_in = _subchannel_mesh.getChannelNode(i_ch, 0);
      const Real mdot_in = (*_mdot_soln)(node_in);
      const Real rho_in = (*_rho_soln)(node_in);
      viscosity_in += mdot_in * (*_mu_soln)(node_in);
      volumetric_flow_in += mdot_in / rho_in;
      mass_flow_in += mdot_in;
    }
    const Real bulk_Dh = _subchannel_mesh.getAssemblyHydraulicDiameter();
    const Real flow_area = _subchannel_mesh.getAssemblyFlowArea();
    const Real inlet_mu = viscosity_in / mass_flow_in;
    _bulk_Re = mass_flow_in * bulk_Dh / (inlet_mu * flow_area);
    _bulk_V = volumetric_flow_in / flow_area;
  }
  // Only rank 0 owns the subchannel DOFs (SingleRankPartitioner), so share the result
  _communicator.broadcast(_bulk_Re);
  _communicator.broadcast(_bulk_V);
 
  if (MooseUtils::absoluteFuzzyEqual(_bulk_Re, 0.0))
    mooseError("The computed bulk Reynolds number is zero.");
 
  if (MooseUtils::absoluteFuzzyEqual(_bulk_V, 0.0))
    mooseError("The computed bulk velocity is zero.");
 
  _mixing_closure->computeBulkMixingParameters();
}
 
void
SubChannel1PhaseProblem::detectDeformation()
{
  const Real tol = libMesh::TOLERANCE;
  const auto pin_diameter = _subchannel_mesh.getPinDiameter();
 
  if (_pin_mesh_exist)
  {
    for (unsigned int iz = 0; iz < _n_cells + 1; iz++)
      for (unsigned int i_pin = 0; i_pin < _n_pins; i_pin++)
      {
        auto * node = _subchannel_mesh.getPinNode(i_pin, iz);
        const Real Dpin = (*_Dpin_soln)(node);
        if (std::abs(Dpin) <= tol)
          mooseError("Dpin is zero at node ",
                     node->id(),
                     ". You must initialize Dpin to a non-zero value.");
        if (std::abs(Dpin - pin_diameter) > tol)
          _deformation = true;
      }
  }
 
  for (unsigned int iz = 0; iz < _n_cells + 1 && !_deformation; iz++)
    for (unsigned int i_ch = 0; i_ch < _n_channels && !_deformation; i_ch++)
    {
      auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
      auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
 
      if ((subch_type == EChannelType::CORNER || subch_type == EChannelType::EDGE) &&
          std::abs((*_displacement_soln)(node)) > tol)
        _deformation = true;
    }
}
 
SubChannel1PhaseProblem::~SubChannel1PhaseProblem()
{
  PetscErrorCode ierr = cleanUp();
  if (ierr)
    mooseError(name(), ": Error in memory cleanup");
}
 
PetscErrorCode
SubChannel1PhaseProblem::cleanUp()
{
  PetscFunctionBegin;
  // We need to clean up the petsc matrices/vectors
  // Mass conservation components
  LibmeshPetscCall(MatDestroy(&_mc_sumWij_mat));
  LibmeshPetscCall(VecDestroy(&_Wij_vec));
  LibmeshPetscCall(VecDestroy(&_prod));
  LibmeshPetscCall(VecDestroy(&_prodp));
  LibmeshPetscCall(MatDestroy(&_mc_axial_convection_mat));
  LibmeshPetscCall(VecDestroy(&_mc_axial_convection_rhs));
  LibmeshPetscCall(MatDestroy(&_mc_density_pressure_mat));
 
  // Axial momentum conservation components
  LibmeshPetscCall(MatDestroy(&_amc_turbulent_cross_flows_mat));
  LibmeshPetscCall(VecDestroy(&_amc_turbulent_cross_flows_rhs));
  LibmeshPetscCall(MatDestroy(&_amc_time_derivative_mat));
  LibmeshPetscCall(VecDestroy(&_amc_time_derivative_rhs));
  LibmeshPetscCall(MatDestroy(&_amc_advective_derivative_mat));
  LibmeshPetscCall(VecDestroy(&_amc_advective_derivative_rhs));
  LibmeshPetscCall(MatDestroy(&_amc_cross_derivative_mat));
  LibmeshPetscCall(VecDestroy(&_amc_cross_derivative_rhs));
  LibmeshPetscCall(MatDestroy(&_amc_friction_force_mat));
  LibmeshPetscCall(VecDestroy(&_amc_friction_force_rhs));
  LibmeshPetscCall(VecDestroy(&_amc_gravity_rhs));
  LibmeshPetscCall(MatDestroy(&_amc_pressure_force_mat));
  LibmeshPetscCall(VecDestroy(&_amc_pressure_force_rhs));
  LibmeshPetscCall(MatDestroy(&_amc_sys_mdot_mat));
  LibmeshPetscCall(VecDestroy(&_amc_sys_mdot_rhs));
 
  // Lateral momentum conservation components
  LibmeshPetscCall(MatDestroy(&_cmc_time_derivative_mat));
  LibmeshPetscCall(VecDestroy(&_cmc_time_derivative_rhs));
  LibmeshPetscCall(MatDestroy(&_cmc_advective_derivative_mat));
  LibmeshPetscCall(VecDestroy(&_cmc_advective_derivative_rhs));
  LibmeshPetscCall(MatDestroy(&_cmc_friction_force_mat));
  LibmeshPetscCall(VecDestroy(&_cmc_friction_force_rhs));
  LibmeshPetscCall(MatDestroy(&_cmc_pressure_force_mat));
  LibmeshPetscCall(VecDestroy(&_cmc_pressure_force_rhs));
  LibmeshPetscCall(MatDestroy(&_cmc_sys_Wij_mat));
  LibmeshPetscCall(VecDestroy(&_cmc_sys_Wij_rhs));
 
  // Energy conservation components
  LibmeshPetscCall(MatDestroy(&_hc_time_derivative_mat));
  LibmeshPetscCall(VecDestroy(&_hc_time_derivative_rhs));
  LibmeshPetscCall(MatDestroy(&_hc_advective_derivative_mat));
  LibmeshPetscCall(VecDestroy(&_hc_advective_derivative_rhs));
  LibmeshPetscCall(MatDestroy(&_hc_cross_derivative_mat));
  LibmeshPetscCall(VecDestroy(&_hc_cross_derivative_rhs));
  LibmeshPetscCall(VecDestroy(&_hc_added_heat_rhs));
  LibmeshPetscCall(MatDestroy(&_hc_sys_h_mat));
  LibmeshPetscCall(VecDestroy(&_hc_sys_h_rhs));
 
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}
 
bool
SubChannel1PhaseProblem::solverSystemConverged(const unsigned int)
{
  return _converged;
}
 
PetscScalar
SubChannel1PhaseProblem::computeInterpolationCoefficients(PetscScalar Peclet)
{
  switch (_interpolation_scheme)
  {
    case 0: // upwind interpolation
      return 1.0;
    case 1: // downwind interpolation
      return 0.0;
    case 2: // central_difference interpolation
      return 0.5;
    case 3: // exponential interpolation (Peclet limited)
    {
      // Weight of the bottom value, alpha = 1 - 1 / Pe + 1 / (exp(Pe) - 1), from the exponential
      // profile of the steady 1D advection-diffusion equation. It is evaluated with exp(-|Pe|)
      // so that it stays finite for large |Pe|; the original form overflows to NaN once exp(Pe)
      // exceeds the double range. For |Pe| below 1e-6 the series limit 0.5 + Pe / 12 is used
      const auto Pe = Peclet;
      if (std::abs(Pe) < 1e-6)
        return 0.5 + Pe / 12.0;
      const auto e = std::exp(-std::abs(Pe));
      const auto alpha_pos = 1.0 - 1.0 / std::abs(Pe) + e / (1.0 - e);
      return Pe > 0 ? alpha_pos : 1.0 - alpha_pos;
    }
    default:
      mooseError(name(),
                 ": Interpolation scheme should be a string: upwind, downwind, central_difference, "
                 "exponential");
  }
}
 
PetscScalar
SubChannel1PhaseProblem::computeInterpolatedValue(PetscScalar topValue,
                                                  PetscScalar botValue,
                                                  PetscScalar Peclet)
{
  PetscScalar alpha = computeInterpolationCoefficients(Peclet);
  return alpha * botValue + (1.0 - alpha) * topValue;
}
 
void
SubChannel1PhaseProblem::computeWijFromSolve(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  // Initial guess, port crossflow of block (iblock) into a vector that will act as my initial guess
  libMesh::DenseVector<Real> solution_seed(_n_gaps * _block_size, 0.0);
  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
  {
    for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
    {
      int i = _n_gaps * (iz - first_node) + i_gap; // column wise transfer
      solution_seed(i) = _Wij(i_gap, iz);
    }
  }
 
  // Solving the combined lateral momentum equation for Wij using a PETSc solver and update vector
  // root
  libMesh::DenseVector<Real> root(_n_gaps * _block_size, 0.0);
  LibmeshPetscCall(petscSnesSolver(iblock, solution_seed, root));
 
  // Assign the solution to the cross-flow matrix
  int i = 0;
  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
  {
    for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
    {
      _Wij(i_gap, iz) = root(i);
      i++;
    }
  }
}
 
void
SubChannel1PhaseProblem::computeSumWij(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  // Add to solution vector if explicit
  if (!_implicit_bool)
  {
    for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    {
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
        Real sumWij = 0.0;
        // Calculate sum of crossflow into channel i from channels j around i
        unsigned int counter = 0;
        for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
        {
          sumWij += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz);
          counter++;
        }
        // The net crossflow coming out of cell i [kg/sec]
        _SumWij_soln->set(node_out, sumWij);
      }
    }
  }
  // Add to matrix if implicit
  else
  {
    LibmeshPetscCall(MatZeroEntries(_mc_density_pressure_mat));
    for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    {
      unsigned int iz_ind = iz - first_node;
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        // Calculate sum of crossflow into channel i from channels j around i
        unsigned int counter = 0;
        for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
        {
          PetscInt row = i_ch + _n_channels * iz_ind;
          PetscInt col = i_gap + _n_gaps * iz_ind;
          PetscScalar value = _subchannel_mesh.getCrossflowSign(i_ch, counter);
          LibmeshPetscCall(MatSetValues(_mc_sumWij_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
          counter++;
        }
      }
    }
    LibmeshPetscCall(MatAssemblyBegin(_mc_sumWij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_mc_sumWij_mat, MAT_FINAL_ASSEMBLY));
    if (_segregated_bool)
    {
      Vec loc_prod;
      Vec loc_Wij;
      LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &loc_prod));
      LibmeshPetscCall(VecDuplicate(_Wij_vec, &loc_Wij));
      LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
          loc_Wij, _Wij, first_node, last_node, _n_gaps));
      LibmeshPetscCall(MatMult(_mc_sumWij_mat, loc_Wij, loc_prod));
      LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
          loc_prod, *_SumWij_soln, first_node, last_node, _n_channels));
      LibmeshPetscCall(VecDestroy(&loc_prod));
      LibmeshPetscCall(VecDestroy(&loc_Wij));
    }
  }
}
 
void
SubChannel1PhaseProblem::computeMdot(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  if (!_implicit_bool)
  {
    for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    {
      auto dz = _z_grid[iz] - _z_grid[iz - 1];
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
        auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
        auto volume = dz * (*_S_flow_soln)(node_in);
        auto time_term = _TR * ((*_rho_soln)(node_out)-_rho_soln->old(node_out)) * volume / _dt;
        // Wij positive out of i into j;
        auto mdot_out = (*_mdot_soln)(node_in) - (*_SumWij_soln)(node_out)-time_term;
        if (mdot_out < 0)
        {
          _console << "Wij = : " << _Wij << "\n";
          mooseError(name(),
                     " : Calculation of negative mass flow mdot_out = : ",
                     mdot_out,
                     " Axial Level= : ",
                     iz,
                     " - Implicit solves are required for recirculating flow.");
        }
        _mdot_soln->set(node_out, mdot_out); // kg/sec
      }
    }
  }
  else
  {
    for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    {
      auto dz = _z_grid[iz] - _z_grid[iz - 1];
      auto iz_ind = iz - first_node;
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
        auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
        auto volume = dz * (*_S_flow_soln)(node_in);
 
        // Adding time derivative to the RHS
        auto time_term = _TR * ((*_rho_soln)(node_out)-_rho_soln->old(node_out)) * volume / _dt;
        PetscInt row_vec = i_ch + _n_channels * iz_ind;
        PetscScalar value_vec = -1.0 * time_term;
        LibmeshPetscCall(
            VecSetValues(_mc_axial_convection_rhs, 1, &row_vec, &value_vec, INSERT_VALUES));
 
        // Linearize transient density with respect to pressure at constant enthalpy. The
        // temperature response at constant h is required even for a rho(T)-only liquid model:
        //
        //   (d rho / d p)_h = (d rho / d p)_T
        //                      - (d rho / d T)_p (d h / d p)_T / (d h / d T)_p.
        //
        // Without this Jacobian contribution, pressure-density coupling remains in the outer
        // fixed-point iteration and scales as O(1 / dt). As the time step decreases, the iteration
        // can become increasingly slow or noncontractive and may require stronger relaxation.
        if (_TR && _compute_density && _compute_power && iz < last_node)
        {
          Real rho, drho_dp_T, drho_dT;
          _fp->rho_from_p_T(
              (*_P_soln)(node_out) + _P_out, (*_T_soln)(node_out), rho, drho_dp_T, drho_dT);
          Real h, dh_dp_T, dh_dT;
          _fp->h_from_p_T((*_P_soln)(node_out) + _P_out, (*_T_soln)(node_out), h, dh_dp_T, dh_dT);
          const Real drho_dp_h = drho_dp_T - drho_dT * dh_dp_T / dh_dT;
          const PetscScalar pressure_coefficient = volume / _dt * drho_dp_h;
          const PetscInt pressure_col = i_ch + _n_channels * (iz_ind + 1);
          LibmeshPetscCall(MatSetValues(_mc_density_pressure_mat,
                                        1,
                                        &row_vec,
                                        1,
                                        &pressure_col,
                                        &pressure_coefficient,
                                        INSERT_VALUES));
          const PetscScalar linearization_rhs = pressure_coefficient * (*_P_soln)(node_out);
          LibmeshPetscCall(
              VecSetValues(_mc_axial_convection_rhs, 1, &row_vec, &linearization_rhs, ADD_VALUES));
        }
 
        // Imposing bottom boundary condition or adding of diagonal elements
        if (iz == first_node)
        {
          PetscScalar value_vec = (*_mdot_soln)(node_in);
          PetscInt row_vec = i_ch + _n_channels * iz_ind;
          LibmeshPetscCall(
              VecSetValues(_mc_axial_convection_rhs, 1, &row_vec, &value_vec, ADD_VALUES));
        }
        else
        {
          PetscInt row = i_ch + _n_channels * iz_ind;
          PetscInt col = i_ch + _n_channels * (iz_ind - 1);
          PetscScalar value = -1.0;
          LibmeshPetscCall(
              MatSetValues(_mc_axial_convection_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
        }
 
        // Adding diagonal elements
        PetscInt row = i_ch + _n_channels * iz_ind;
        PetscInt col = i_ch + _n_channels * iz_ind;
        PetscScalar value = 1.0;
        LibmeshPetscCall(
            MatSetValues(_mc_axial_convection_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
 
        // Adding cross flows RHS
        if (_segregated_bool)
        {
          PetscScalar value_vec_2 = -1.0 * (*_SumWij_soln)(node_out);
          PetscInt row_vec_2 = i_ch + _n_channels * iz_ind;
          LibmeshPetscCall(
              VecSetValues(_mc_axial_convection_rhs, 1, &row_vec_2, &value_vec_2, ADD_VALUES));
        }
      }
    }
    LibmeshPetscCall(MatAssemblyBegin(_mc_axial_convection_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_mc_axial_convection_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_mc_density_pressure_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_mc_density_pressure_mat, MAT_FINAL_ASSEMBLY));
 
    if (_segregated_bool)
    {
      KSP ksploc;
      PC pc;
      Vec sol;
      LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &sol));
      LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksploc));
      LibmeshPetscCall(KSPSetOperators(ksploc, _mc_axial_convection_mat, _mc_axial_convection_mat));
      LibmeshPetscCall(KSPGetPC(ksploc, &pc));
      LibmeshPetscCall(PCSetType(pc, PCJACOBI));
      LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
      LibmeshPetscCall(KSPSetFromOptions(ksploc));
      LibmeshPetscCall(KSPSolve(ksploc, _mc_axial_convection_rhs, sol));
      LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
          sol, *_mdot_soln, first_node, last_node, _n_channels));
      LibmeshPetscCall(VecZeroEntries(_mc_axial_convection_rhs));
      LibmeshPetscCall(KSPDestroy(&ksploc));
      LibmeshPetscCall(VecDestroy(&sol));
    }
  }
}
 
void
SubChannel1PhaseProblem::computeDP(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  if (!_implicit_bool)
  {
    for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    {
      auto k_grid = _subchannel_mesh.getKGrid();
      auto dz = _z_grid[iz] - _z_grid[iz - 1];
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
        auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
        auto rho_in = (*_rho_soln)(node_in);
        auto rho_out = (*_rho_soln)(node_out);
        auto mu_in = (*_mu_soln)(node_in);
        auto S = (*_S_flow_soln)(node_in);
        auto w_perim = (*_w_perim_soln)(node_in);
        auto time_term = _TR * ((*_mdot_soln)(node_out)-_mdot_soln->old(node_out)) * dz / _dt -
                         dz * 2.0 * (*_mdot_soln)(node_out) * (rho_out - _rho_soln->old(node_out)) /
                             rho_in / _dt;
        auto mass_term1 =
            Utility::pow<2>((*_mdot_soln)(node_out)) * (1.0 / S / rho_out - 1.0 / S / rho_in);
        auto mass_term2 = -2.0 * (*_mdot_soln)(node_out) * (*_SumWij_soln)(node_out) / S / rho_in;
        auto crossflow_term = 0.0;
        auto turbulent_term = 0.0;
        unsigned int counter = 0;
        for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
        {
          auto chans = _subchannel_mesh.getGapChannels(i_gap);
          unsigned int ii_ch = chans.first;
          unsigned int jj_ch = chans.second;
          auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
          auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
          auto * node_out_i = _subchannel_mesh.getChannelNode(ii_ch, iz);
          auto * node_out_j = _subchannel_mesh.getChannelNode(jj_ch, iz);
          auto rho_i = (*_rho_soln)(node_in_i);
          auto rho_j = (*_rho_soln)(node_in_j);
          auto Si = (*_S_flow_soln)(node_in_i);
          auto Sj = (*_S_flow_soln)(node_in_j);
          Real u_star = 0.0;
          // figure out donor axial velocity
          if (_Wij(i_gap, iz) > 0.0)
            u_star = (*_mdot_soln)(node_out_i) / Si / rho_i;
          else
            u_star = (*_mdot_soln)(node_out_j) / Sj / rho_j;
 
          crossflow_term +=
              _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * u_star;
 
          turbulent_term += _WijPrime(i_gap, iz) * (2 * (*_mdot_soln)(node_out) / rho_in / S -
                                                    (*_mdot_soln)(node_out_j) / Sj / rho_j -
                                                    (*_mdot_soln)(node_out_i) / Si / rho_i);
          counter++;
        }
        turbulent_term *= _CT;
        auto ki = 0.0;
        if ((*_mdot_soln)(node_out) >= 0)
          ki = k_grid[i_ch][iz - 1];
        else
          ki = k_grid[i_ch][iz];
        // Cell Peclet number of the exponential scheme, from the friction of the cell evaluated at
        // the arithmetic midpoint, as in the implicit branch
        PetscScalar Pe = 0.5;
        if (_interpolation_scheme == 3)
        {
          const auto S_mid = computeInterpolatedValue((*_S_flow_soln)(node_out), S, 0.5);
          const auto w_perim_mid =
              computeInterpolatedValue((*_w_perim_soln)(node_out), w_perim, 0.5);
          const auto mdot_mid =
              computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), 0.5);
          const auto mu_mid = computeInterpolatedValue((*_mu_soln)(node_out), mu_in, 0.5);
          const auto Dh_mid = 4.0 * S_mid / w_perim_mid;
          _friction_args = FrictionStruct(
              i_ch, std::abs(mdot_mid) / S_mid * Dh_mid / mu_mid, S_mid, w_perim_mid);
          const Real ff_mid = _friction_closure->computeFrictionFactor(_friction_args);
          Pe = 1.0 / ((ff_mid * dz / Dh_mid + ki) * 0.5) * mdot_mid / std::abs(mdot_mid);
        }
        // The friction factor and the friction force are cell-center quantities
        const auto S_c = computeInterpolatedValue((*_S_flow_soln)(node_out), S, Pe);
        const auto w_perim_c = computeInterpolatedValue((*_w_perim_soln)(node_out), w_perim, Pe);
        const auto mu_c = computeInterpolatedValue((*_mu_soln)(node_out), mu_in, Pe);
        const auto rho_c = computeInterpolatedValue(rho_out, rho_in, Pe);
        const auto mdot_c =
            computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), Pe);
        const auto Dh_c = 4.0 * S_c / w_perim_c;
        auto Re = std::abs(mdot_c) / S_c * Dh_c / mu_c;
        _friction_args = FrictionStruct(i_ch, Re, S_c, w_perim_c);
        Real ff = _friction_closure->computeFrictionFactor(_friction_args);
        if (_ff_soln)
          _ff_soln->set(node_out, ff);
        auto friction_term =
            (ff * dz / Dh_c + ki) * 0.5 * mdot_c * std::abs(mdot_c) / (S_c * rho_c);
        auto gravity_term = _dir_grav * _g_grav * (*_rho_soln)(node_out)*dz * S;
        auto DP = (1 / S) * (time_term + mass_term1 + mass_term2 + crossflow_term + turbulent_term +
                             friction_term + gravity_term); // Pa
        _DP(i_ch, iz) = DP;
        if (_DP_soln)
          _DP_soln->set(node_out, DP);
      }
    }
  }
  else
  {
    LibmeshPetscCall(MatZeroEntries(_amc_time_derivative_mat));
    LibmeshPetscCall(MatZeroEntries(_amc_advective_derivative_mat));
    LibmeshPetscCall(MatZeroEntries(_amc_cross_derivative_mat));
    LibmeshPetscCall(MatZeroEntries(_amc_friction_force_mat));
    LibmeshPetscCall(VecZeroEntries(_amc_time_derivative_rhs));
    LibmeshPetscCall(VecZeroEntries(_amc_advective_derivative_rhs));
    LibmeshPetscCall(VecZeroEntries(_amc_cross_derivative_rhs));
    LibmeshPetscCall(VecZeroEntries(_amc_friction_force_rhs));
    LibmeshPetscCall(VecZeroEntries(_amc_gravity_rhs));
    LibmeshPetscCall(MatZeroEntries(_amc_sys_mdot_mat));
    LibmeshPetscCall(VecZeroEntries(_amc_sys_mdot_rhs));
    for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    {
      auto k_grid = _subchannel_mesh.getKGrid();
      auto dz = _z_grid[iz] - _z_grid[iz - 1];
      auto iz_ind = iz - first_node;
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        // inlet and outlet nodes
        auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
        auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
 
        // interpolation weight coefficient
        PetscScalar Pe = 0.5;
        if (_interpolation_scheme == 3)
        {
          // Compute the Peclet number
          auto S_in = (*_S_flow_soln)(node_in);
          auto S_out = (*_S_flow_soln)(node_out);
          auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
          auto w_perim_in = (*_w_perim_soln)(node_in);
          auto w_perim_out = (*_w_perim_soln)(node_out);
          auto w_perim_interp = this->computeInterpolatedValue(w_perim_out, w_perim_in, 0.5);
          auto mdot_loc =
              this->computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), 0.5);
          auto mu_in = (*_mu_soln)(node_in);
          auto mu_out = (*_mu_soln)(node_out);
          auto mu_interp = this->computeInterpolatedValue(mu_out, mu_in, 0.5);
          auto Dh_c = 4.0 * S_interp / w_perim_interp;
          // Compute friction factor
          auto Re = ((std::abs(mdot_loc) / S_interp) * Dh_c / mu_interp);
          _friction_args = FrictionStruct(i_ch, Re, S_interp, w_perim_interp);
          Real ff = _friction_closure->computeFrictionFactor(_friction_args);
          if (_ff_soln)
            _ff_soln->set(node_out, ff);
          auto ki = 0.0;
          if ((*_mdot_soln)(node_out) >= 0)
            ki = k_grid[i_ch][iz - 1];
          else
            ki = k_grid[i_ch][iz];
          Pe = 1.0 / ((ff * dz / Dh_c + ki) * 0.5) * mdot_loc / std::abs(mdot_loc);
        }
        auto alpha = computeInterpolationCoefficients(Pe);
 
        // inlet, outlet, and interpolated density
        auto rho_in = (*_rho_soln)(node_in);
        auto rho_out = (*_rho_soln)(node_out);
        auto rho_interp = computeInterpolatedValue(rho_out, rho_in, Pe);
 
        // inlet, outlet, and interpolated viscosity
        auto mu_in = (*_mu_soln)(node_in);
        auto mu_out = (*_mu_soln)(node_out);
        auto mu_interp = computeInterpolatedValue(mu_out, mu_in, Pe);
 
        // inlet, outlet, and interpolated axial surface area
        auto S_in = (*_S_flow_soln)(node_in);
        auto S_out = (*_S_flow_soln)(node_out);
        auto S_interp = computeInterpolatedValue(S_out, S_in, Pe);
 
        // inlet, outlet, and interpolated wetted perimeter
        auto w_perim_in = (*_w_perim_soln)(node_in);
        auto w_perim_out = (*_w_perim_soln)(node_out);
        auto w_perim_interp = computeInterpolatedValue(w_perim_out, w_perim_in, Pe);
 
        // cell-center hydraulic diameter
        auto Dh_c = 4.0 * S_interp / w_perim_interp;
 
        // Keep temporal storage independent of the spatial interpolation scheme. Applying the
        // upwind coefficient to this term puts the storage on the upstream node and leaves the
        // first axial row without transient inertia in the upwind limit.
        PetscInt row_tt = i_ch + _n_channels * iz_ind;
        PetscInt col_tt = i_ch + _n_channels * iz_ind;
        PetscScalar value_tt = _TR * dz / _dt;
        LibmeshPetscCall(MatSetValues(
            _amc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
 
        PetscScalar value_vec_tt = _TR * _mdot_soln->old(node_out) * dz / _dt;
        PetscInt row_vec_tt = i_ch + _n_channels * iz_ind;
        LibmeshPetscCall(
            VecSetValues(_amc_time_derivative_rhs, 1, &row_vec_tt, &value_vec_tt, ADD_VALUES));
 
        if (iz == first_node)
        {
          PetscScalar value_vec_at = Utility::pow<2>((*_mdot_soln)(node_in)) / (S_in * rho_in);
          PetscInt row_vec_at = i_ch + _n_channels * iz_ind;
          LibmeshPetscCall(VecSetValues(
              _amc_advective_derivative_rhs, 1, &row_vec_at, &value_vec_at, ADD_VALUES));
        }
        else
        {
          PetscInt row_at = i_ch + _n_channels * iz_ind;
          PetscInt col_at = i_ch + _n_channels * (iz_ind - 1);
          PetscScalar value_at = -1.0 * (*_mdot_soln)(node_in) / (S_in * rho_in);
          LibmeshPetscCall(MatSetValues(
              _amc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
        }
 
        // Adding diagonal elements
        PetscInt row_at = i_ch + _n_channels * iz_ind;
        PetscInt col_at = i_ch + _n_channels * iz_ind;
        PetscScalar value_at = (*_mdot_soln)(node_out) / (S_out * rho_out);
        LibmeshPetscCall(MatSetValues(
            _amc_advective_derivative_mat, 1, &row_at, 1, &col_at, &value_at, INSERT_VALUES));
 
        unsigned int counter = 0;
        unsigned int cross_index = iz; // iz-1;
        for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
        {
          auto chans = _subchannel_mesh.getGapChannels(i_gap);
          unsigned int ii_ch = chans.first;
          unsigned int jj_ch = chans.second;
          auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
          auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
          auto * node_out_i = _subchannel_mesh.getChannelNode(ii_ch, iz);
          auto * node_out_j = _subchannel_mesh.getChannelNode(jj_ch, iz);
          auto rho_i =
              computeInterpolatedValue((*_rho_soln)(node_out_i), (*_rho_soln)(node_in_i), Pe);
          auto rho_j =
              computeInterpolatedValue((*_rho_soln)(node_out_j), (*_rho_soln)(node_in_j), Pe);
          auto S_i =
              computeInterpolatedValue((*_S_flow_soln)(node_out_i), (*_S_flow_soln)(node_in_i), Pe);
          auto S_j =
              computeInterpolatedValue((*_S_flow_soln)(node_out_j), (*_S_flow_soln)(node_in_j), Pe);
          auto u_star = 0.0;
          // figure out donor axial velocity
          if (_Wij(i_gap, cross_index) > 0.0)
          {
            if (iz == first_node)
            {
              u_star = (*_mdot_soln)(node_in_i) / S_i / rho_i;
              PetscScalar value_vec_ct = -1.0 * alpha *
                                         _subchannel_mesh.getCrossflowSign(i_ch, counter) *
                                         _Wij(i_gap, cross_index) * u_star;
              PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
              LibmeshPetscCall(VecSetValues(
                  _amc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
            }
            else
            {
              PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
                                     _Wij(i_gap, cross_index) / S_i / rho_i;
              PetscInt row_ct = i_ch + _n_channels * iz_ind;
              PetscInt col_ct = ii_ch + _n_channels * (iz_ind - 1);
              LibmeshPetscCall(MatSetValues(
                  _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
            }
            PetscScalar value_ct = (1.0 - alpha) *
                                   _subchannel_mesh.getCrossflowSign(i_ch, counter) *
                                   _Wij(i_gap, cross_index) / S_i / rho_i;
            PetscInt row_ct = i_ch + _n_channels * iz_ind;
            PetscInt col_ct = ii_ch + _n_channels * iz_ind;
            LibmeshPetscCall(MatSetValues(
                _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
          }
          else if (_Wij(i_gap, cross_index) < 0.0) // _Wij=0 operations not necessary
          {
            if (iz == first_node)
            {
              u_star = (*_mdot_soln)(node_in_j) / S_j / rho_j;
              PetscScalar value_vec_ct = -1.0 * alpha *
                                         _subchannel_mesh.getCrossflowSign(i_ch, counter) *
                                         _Wij(i_gap, cross_index) * u_star;
              PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
              LibmeshPetscCall(VecSetValues(
                  _amc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
            }
            else
            {
              PetscScalar value_ct = alpha * _subchannel_mesh.getCrossflowSign(i_ch, counter) *
                                     _Wij(i_gap, cross_index) / S_j / rho_j;
              PetscInt row_ct = i_ch + _n_channels * iz_ind;
              PetscInt col_ct = jj_ch + _n_channels * (iz_ind - 1);
              LibmeshPetscCall(MatSetValues(
                  _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
            }
            PetscScalar value_ct = (1.0 - alpha) *
                                   _subchannel_mesh.getCrossflowSign(i_ch, counter) *
                                   _Wij(i_gap, cross_index) / S_j / rho_j;
            PetscInt row_ct = i_ch + _n_channels * iz_ind;
            PetscInt col_ct = jj_ch + _n_channels * iz_ind;
            LibmeshPetscCall(MatSetValues(
                _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_ct, ADD_VALUES));
          }
 
          if (iz == first_node)
          {
            PetscScalar value_vec_ct = -2.0 * alpha * (*_mdot_soln)(node_in)*_CT *
                                       _WijPrime(i_gap, cross_index) / (rho_interp * S_interp);
            value_vec_ct += alpha * (*_mdot_soln)(node_in_j)*_CT * _WijPrime(i_gap, cross_index) /
                            (rho_j * S_j);
            value_vec_ct += alpha * (*_mdot_soln)(node_in_i)*_CT * _WijPrime(i_gap, cross_index) /
                            (rho_i * S_i);
            PetscInt row_vec_ct = i_ch + _n_channels * iz_ind;
            LibmeshPetscCall(
                VecSetValues(_amc_cross_derivative_rhs, 1, &row_vec_ct, &value_vec_ct, ADD_VALUES));
          }
          else
          {
            PetscScalar value_center_ct =
                2.0 * alpha * _CT * _WijPrime(i_gap, cross_index) / (rho_interp * S_interp);
            PetscInt row_ct = i_ch + _n_channels * iz_ind;
            PetscInt col_ct = i_ch + _n_channels * (iz_ind - 1);
            LibmeshPetscCall(MatSetValues(
                _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
 
            PetscScalar value_left_ct =
                -1.0 * alpha * _CT * _WijPrime(i_gap, cross_index) / (rho_j * S_j);
            row_ct = i_ch + _n_channels * iz_ind;
            col_ct = jj_ch + _n_channels * (iz_ind - 1);
            LibmeshPetscCall(MatSetValues(
                _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
 
            PetscScalar value_right_ct =
                -1.0 * alpha * _CT * _WijPrime(i_gap, cross_index) / (rho_i * S_i);
            row_ct = i_ch + _n_channels * iz_ind;
            col_ct = ii_ch + _n_channels * (iz_ind - 1);
            LibmeshPetscCall(MatSetValues(
                _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
          }
 
          PetscScalar value_center_ct =
              2.0 * (1.0 - alpha) * _CT * _WijPrime(i_gap, cross_index) / (rho_interp * S_interp);
          PetscInt row_ct = i_ch + _n_channels * iz_ind;
          PetscInt col_ct = i_ch + _n_channels * iz_ind;
          LibmeshPetscCall(MatSetValues(
              _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_center_ct, ADD_VALUES));
 
          PetscScalar value_left_ct =
              -1.0 * (1.0 - alpha) * _CT * _WijPrime(i_gap, cross_index) / (rho_j * S_j);
          row_ct = i_ch + _n_channels * iz_ind;
          col_ct = jj_ch + _n_channels * iz_ind;
          LibmeshPetscCall(MatSetValues(
              _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_left_ct, ADD_VALUES));
 
          PetscScalar value_right_ct =
              -1.0 * (1.0 - alpha) * _CT * _WijPrime(i_gap, cross_index) / (rho_i * S_i);
          row_ct = i_ch + _n_channels * iz_ind;
          col_ct = ii_ch + _n_channels * iz_ind;
          LibmeshPetscCall(MatSetValues(
              _amc_cross_derivative_mat, 1, &row_ct, 1, &col_ct, &value_right_ct, ADD_VALUES));
          counter++;
        }
 
        PetscScalar mdot_interp =
            computeInterpolatedValue((*_mdot_soln)(node_out), (*_mdot_soln)(node_in), Pe);
        auto Re = ((std::abs(mdot_interp) / S_interp) * Dh_c / mu_interp);
        _friction_args = FrictionStruct(i_ch, Re, S_interp, w_perim_interp);
        Real ff = _friction_closure->computeFrictionFactor(_friction_args);
        if (_ff_soln)
          _ff_soln->set(node_out, ff);
        auto ki = 0.0;
        if ((*_mdot_soln)(node_out) >= 0)
          ki = k_grid[i_ch][iz - 1];
        else
          ki = k_grid[i_ch][iz];
        auto coef = (ff * dz / Dh_c + ki) * 0.5 * std::abs(mdot_interp) / (S_interp * rho_interp);
        if (iz == first_node)
        {
          PetscScalar value_vec = -1.0 * alpha * coef * (*_mdot_soln)(node_in);
          PetscInt row_vec = i_ch + _n_channels * iz_ind;
          LibmeshPetscCall(
              VecSetValues(_amc_friction_force_rhs, 1, &row_vec, &value_vec, ADD_VALUES));
        }
        else
        {
          PetscInt row = i_ch + _n_channels * iz_ind;
          PetscInt col = i_ch + _n_channels * (iz_ind - 1);
          PetscScalar value = alpha * coef;
          LibmeshPetscCall(
              MatSetValues(_amc_friction_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
        }
 
        // Adding diagonal elements
        PetscInt row = i_ch + _n_channels * iz_ind;
        PetscInt col = i_ch + _n_channels * iz_ind;
        PetscScalar value = (1.0 - alpha) * coef;
        LibmeshPetscCall(
            MatSetValues(_amc_friction_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
 
        PetscScalar value_vec = _dir_grav * -1.0 * _g_grav * rho_interp * dz * S_interp;
        PetscInt row_vec = i_ch + _n_channels * iz_ind;
        LibmeshPetscCall(VecSetValues(_amc_gravity_rhs, 1, &row_vec, &value_vec, ADD_VALUES));
      }
    }
    LibmeshPetscCall(MatZeroEntries(_amc_sys_mdot_mat));
    LibmeshPetscCall(VecZeroEntries(_amc_sys_mdot_rhs));
    LibmeshPetscCall(MatAssemblyBegin(_amc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_amc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_amc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_cross_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_amc_friction_force_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_friction_force_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    // Matrix
#if !PETSC_VERSION_LESS_THAN(3, 15, 0)
    LibmeshPetscCall(
        MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_cross_derivative_mat, UNKNOWN_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_friction_force_mat, UNKNOWN_NONZERO_PATTERN));
#else
    LibmeshPetscCall(
        MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_cross_derivative_mat, DIFFERENT_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_amc_sys_mdot_mat, 1.0, _amc_friction_force_mat, DIFFERENT_NONZERO_PATTERN));
#endif
    LibmeshPetscCall(MatAssemblyBegin(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_sys_mdot_mat, MAT_FINAL_ASSEMBLY));
    // RHS
    LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_time_derivative_rhs));
    LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_advective_derivative_rhs));
    LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_cross_derivative_rhs));
    LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_friction_force_rhs));
    LibmeshPetscCall(VecAXPY(_amc_sys_mdot_rhs, 1.0, _amc_gravity_rhs));
    if (_segregated_bool)
    {
      // Assembly the matrix system
      LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
          _prod, *_mdot_soln, first_node, last_node, _n_channels));
      Vec ls;
      LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &ls));
      LibmeshPetscCall(MatMult(_amc_sys_mdot_mat, _prod, ls));
      LibmeshPetscCall(VecAXPY(ls, -1.0, _amc_sys_mdot_rhs));
      PetscScalar * xx;
      LibmeshPetscCall(VecGetArray(ls, &xx));
      for (unsigned int iz = first_node; iz < last_node + 1; iz++)
      {
        auto iz_ind = iz - first_node;
        for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
        {
          // Setting nodes
          auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
          auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
 
          // inlet, outlet, and interpolated axial surface area
          auto S_in = (*_S_flow_soln)(node_in);
          auto S_out = (*_S_flow_soln)(node_out);
          auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
 
          // Setting solutions
          if (S_interp != 0)
          {
            auto DP = (1 / S_interp) * xx[iz_ind * _n_channels + i_ch];
            _DP(i_ch, iz) = DP;
            if (_DP_soln)
              _DP_soln->set(node_out, DP);
          }
          else
          {
            auto DP = 0.0;
            _DP(i_ch, iz) = DP;
            if (_DP_soln)
              _DP_soln->set(node_out, DP);
          }
        }
      }
      LibmeshPetscCall(VecDestroy(&ls));
    }
  }
}
 
void
SubChannel1PhaseProblem::computeP(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  if (!_implicit_bool)
  {
    for (unsigned int iz = last_node; iz > first_node - 1; iz--)
    {
      // Calculate pressure in the inlet of the cell assuming known outlet
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
        auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
        // update Pressure solution
        _P_soln->set(node_in, (*_P_soln)(node_out) + _DP(i_ch, iz));
      }
    }
  }
  else
  {
    LibmeshPetscCall(VecZeroEntries(_amc_pressure_force_rhs));
    for (unsigned int iz = last_node; iz > first_node - 1; iz--)
    {
      auto iz_ind = iz - first_node;
      // Calculate pressure in the inlet of the cell assuming known outlet
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
        auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
 
        // inlet, outlet, and interpolated axial surface area
        auto S_in = (*_S_flow_soln)(node_in);
        auto S_out = (*_S_flow_soln)(node_out);
        auto S_interp = computeInterpolatedValue(S_out, S_in, 0.5);
 
        // Creating matrix of coefficients
        PetscInt row = i_ch + _n_channels * iz_ind;
        PetscInt col = i_ch + _n_channels * iz_ind;
        PetscScalar value = -1.0 * S_interp;
        LibmeshPetscCall(
            MatSetValues(_amc_pressure_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
 
        // Row iz_ind is the pressure force on cell iz, S_c * (P(iz) - P(iz - 1)), and its unknown
        // is the inlet pressure P(iz - 1). The outlet pressure of the top cell of the block is
        // known (the outlet boundary condition or the block above), so it goes to the right-hand
        // side. For the other cells, P(iz) is the unknown of the next row.
        if (iz == last_node)
        {
          PetscScalar value = -1.0 * (*_P_soln)(node_out)*S_interp;
          PetscInt row = i_ch + _n_channels * iz_ind;
          LibmeshPetscCall(VecSetValues(_amc_pressure_force_rhs, 1, &row, &value, ADD_VALUES));
        }
        else
        {
          PetscInt row = i_ch + _n_channels * iz_ind;
          PetscInt col = i_ch + _n_channels * (iz_ind + 1);
          PetscScalar value = 1.0 * S_interp;
          LibmeshPetscCall(
              MatSetValues(_amc_pressure_force_mat, 1, &row, 1, &col, &value, INSERT_VALUES));
        }
 
        // The segregated solver takes the cell pressure drop DP from the axial momentum solve
        // (computeDP), so each row reduces to P(iz - 1) = P(iz) + DP(iz). The monolithic solver
        // omits this term because the same matrix is the pressure force in its axial momentum
        // equation, which determines P directly.
        if (_segregated_bool)
        {
          auto dp_out = _DP(i_ch, iz);
          PetscScalar value_v = -1.0 * dp_out * S_interp;
          PetscInt row_v = i_ch + _n_channels * iz_ind;
          LibmeshPetscCall(VecSetValues(_amc_pressure_force_rhs, 1, &row_v, &value_v, ADD_VALUES));
        }
      }
    }
    // Solving pressure problem
    LibmeshPetscCall(MatAssemblyBegin(_amc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
    if (_segregated_bool)
    {
      KSP ksploc;
      PC pc;
      Vec sol;
      LibmeshPetscCall(VecDuplicate(_amc_pressure_force_rhs, &sol));
      LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksploc));
      LibmeshPetscCall(KSPSetOperators(ksploc, _amc_pressure_force_mat, _amc_pressure_force_mat));
      LibmeshPetscCall(KSPGetPC(ksploc, &pc));
      LibmeshPetscCall(PCSetType(pc, PCJACOBI));
      LibmeshPetscCall(KSPSetTolerances(ksploc, _rtol, _atol, _dtol, _maxit));
      LibmeshPetscCall(KSPSetFromOptions(ksploc));
      LibmeshPetscCall(KSPSolve(ksploc, _amc_pressure_force_rhs, sol));
      PetscScalar * xx;
      LibmeshPetscCall(VecGetArray(sol, &xx));
      // update Pressure solution
      for (unsigned int iz = last_node; iz > first_node - 1; iz--)
      {
        auto iz_ind = iz - first_node;
        for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
        {
          auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
          PetscScalar value = xx[iz_ind * _n_channels + i_ch];
          _P_soln->set(node_in, value);
        }
      }
      LibmeshPetscCall(VecZeroEntries(_amc_pressure_force_rhs));
      LibmeshPetscCall(KSPDestroy(&ksploc));
      LibmeshPetscCall(VecDestroy(&sol));
    }
  }
}
 
Real
SubChannel1PhaseProblem::computeT(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  std::vector<Real> residual;
  residual.reserve(_block_size * _n_channels);
  Real residual_norm_sq = 0.0;
  Real temperature_norm_sq = 0.0;
  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
  {
    for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
    {
      auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
      const Real T = (*_T_soln)(node);
      const Real T_from_ph = _fp->T_from_p_h((*_P_soln)(node) + _P_out, (*_h_soln)(node));
      residual.push_back(T_from_ph - T);
      residual_norm_sq += Utility::pow<2>(residual.back());
      temperature_norm_sq += Utility::pow<2>(T);
    }
  }
 
  // Set the temperature solution with relaxation if needed
  std::size_t i = 0;
  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
    {
      auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
      _T_soln->set(node, (*_T_soln)(node) + _T_relaxation * residual[i]);
      ++i;
    }
 
  return std::sqrt(residual_norm_sq) / (std::sqrt(temperature_norm_sq) + 1e-14);
}
 
void
SubChannel1PhaseProblem::computeRho(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  if (iblock == 0)
  {
    for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
    {
      auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
      _rho_soln->set(node, _fp->rho_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node)));
    }
  }
  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
  {
    for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
    {
      auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
      _rho_soln->set(node, _fp->rho_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node)));
    }
  }
}
 
void
SubChannel1PhaseProblem::computeMu(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  if (iblock == 0)
  {
    for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
    {
      auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
      _mu_soln->set(node, _fp->mu_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node)));
    }
  }
  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
  {
    for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
    {
      auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
      _mu_soln->set(node, _fp->mu_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node)));
    }
  }
}
 
void
SubChannel1PhaseProblem::computeWijResidual(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  // Cross flow residual
  if (!_implicit_bool)
  {
    const Real & pitch = _subchannel_mesh.getPitch();
    for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    {
      auto dz = _z_grid[iz] - _z_grid[iz - 1];
      for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
      {
        auto chans = _subchannel_mesh.getGapChannels(i_gap);
        unsigned int i_ch = chans.first;
        unsigned int j_ch = chans.second;
        auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
        auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
        auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
        auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
        auto rho_i = (*_rho_soln)(node_in_i);
        auto rho_j = (*_rho_soln)(node_in_j);
        auto Si = (*_S_flow_soln)(node_in_i);
        auto Sj = (*_S_flow_soln)(node_in_j);
        auto Sij = dz * _subchannel_mesh.getGapWidth(iz, i_gap);
        auto Lij = pitch;
        // total local form loss in the ij direction
        auto friction_term = _kij * _Wij(i_gap, iz) * std::abs(_Wij(i_gap, iz));
        auto DPij = (*_P_soln)(node_in_i) - (*_P_soln)(node_in_j);
        // Figure out donor cell density
        auto rho_star = 0.0;
        if (_Wij(i_gap, iz) > 0.0)
          rho_star = rho_i;
        else if (_Wij(i_gap, iz) < 0.0)
          rho_star = rho_j;
        else
          rho_star = (rho_i + rho_j) / 2.0;
        auto mass_term_out =
            (*_mdot_soln)(node_out_i) / (*_S_flow_soln)(node_out_i) / (*_rho_soln)(node_out_i) +
            (*_mdot_soln)(node_out_j) / (*_S_flow_soln)(node_out_j) / (*_rho_soln)(node_out_j);
        auto mass_term_in =
            (*_mdot_soln)(node_in_i) / Si / rho_i + (*_mdot_soln)(node_in_j) / Sj / rho_j;
        auto term_out = Sij * rho_star * (Lij / dz) * mass_term_out * _Wij(i_gap, iz);
        auto term_in = Sij * rho_star * (Lij / dz) * mass_term_in * _Wij(i_gap, iz - 1);
        auto inertia_term = term_out - term_in;
        auto pressure_term = 2 * Utility::pow<2>(Sij) * DPij * rho_star;
        auto time_term =
            _TR * 2.0 * (_Wij(i_gap, iz) - _Wij_old(i_gap, iz)) * Lij * Sij * rho_star / _dt;
 
        _Wij_residual_matrix(i_gap, iz - 1 - iblock * _block_size) =
            time_term + friction_term + inertia_term - pressure_term;
      }
    }
  }
  else
  {
    // Initializing to zero the elements of the lateral momentum assembly
    LibmeshPetscCall(MatZeroEntries(_cmc_time_derivative_mat));
    LibmeshPetscCall(MatZeroEntries(_cmc_advective_derivative_mat));
    LibmeshPetscCall(MatZeroEntries(_cmc_friction_force_mat));
    LibmeshPetscCall(MatZeroEntries(_cmc_pressure_force_mat));
    LibmeshPetscCall(VecZeroEntries(_cmc_time_derivative_rhs));
    LibmeshPetscCall(VecZeroEntries(_cmc_advective_derivative_rhs));
    LibmeshPetscCall(VecZeroEntries(_cmc_friction_force_rhs));
    LibmeshPetscCall(VecZeroEntries(_cmc_pressure_force_rhs));
    LibmeshPetscCall(MatZeroEntries(_cmc_sys_Wij_mat));
    LibmeshPetscCall(VecZeroEntries(_cmc_sys_Wij_rhs));
    const Real & pitch = _subchannel_mesh.getPitch();
    for (unsigned int iz = first_node; iz < last_node + 1; iz++)
    {
      auto dz = _z_grid[iz] - _z_grid[iz - 1];
      auto iz_ind = iz - first_node;
      for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
      {
        auto chans = _subchannel_mesh.getGapChannels(i_gap);
        unsigned int i_ch = chans.first;
        unsigned int j_ch = chans.second;
        auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
        auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
        auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
        auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
 
        // inlet, outlet, and interpolated densities
        auto rho_i_in = (*_rho_soln)(node_in_i);
        auto rho_i_out = (*_rho_soln)(node_out_i);
        auto rho_i_interp = computeInterpolatedValue(rho_i_out, rho_i_in, 0.5);
        auto rho_j_in = (*_rho_soln)(node_in_j);
        auto rho_j_out = (*_rho_soln)(node_out_j);
        auto rho_j_interp = computeInterpolatedValue(rho_j_out, rho_j_in, 0.5);
 
        // inlet, outlet, and interpolated areas
        auto S_i_in = (*_S_flow_soln)(node_in_i);
        auto S_i_out = (*_S_flow_soln)(node_out_i);
        auto S_j_in = (*_S_flow_soln)(node_in_j);
        auto S_j_out = (*_S_flow_soln)(node_out_j);
 
        // Cross-sectional gap area
        auto Sij = dz * _subchannel_mesh.getGapWidth(iz, i_gap);
        auto Lij = pitch;
 
        // Figure out donor cell density
        auto rho_star = 0.0;
        if (_Wij(i_gap, iz) > 0.0)
          rho_star = rho_i_interp;
        else if (_Wij(i_gap, iz) < 0.0)
          rho_star = rho_j_interp;
        else
          rho_star = (rho_i_interp + rho_j_interp) / 2.0;
 
        // Assembling time derivative
        PetscScalar time_factor = _TR * Lij * Sij * rho_star / _dt;
        PetscInt row_td = i_gap + _n_gaps * iz_ind;
        PetscInt col_td = i_gap + _n_gaps * iz_ind;
        PetscScalar value_td = time_factor;
        LibmeshPetscCall(MatSetValues(
            _cmc_time_derivative_mat, 1, &row_td, 1, &col_td, &value_td, INSERT_VALUES));
        PetscScalar value_td_rhs = time_factor * _Wij_old(i_gap, iz);
        LibmeshPetscCall(
            VecSetValues(_cmc_time_derivative_rhs, 1, &row_td, &value_td_rhs, INSERT_VALUES));
 
        // Assembling inertial term
        PetscScalar Pe = 0.5;
        auto alpha = computeInterpolationCoefficients(Pe);
        auto mass_term_out = (*_mdot_soln)(node_out_i) / S_i_out / rho_i_out +
                             (*_mdot_soln)(node_out_j) / S_j_out / rho_j_out;
        auto mass_term_in = (*_mdot_soln)(node_in_i) / S_i_in / rho_i_in +
                            (*_mdot_soln)(node_in_j) / S_j_in / rho_j_in;
        auto term_out = Sij * rho_star * (Lij / dz) * mass_term_out / 2.0;
        auto term_in = Sij * rho_star * (Lij / dz) * mass_term_in / 2.0;
        if (iz == first_node)
        {
          PetscInt row_ad = i_gap + _n_gaps * iz_ind;
          PetscScalar value_ad = term_in * alpha * _Wij(i_gap, iz - 1);
          LibmeshPetscCall(
              VecSetValues(_cmc_advective_derivative_rhs, 1, &row_ad, &value_ad, ADD_VALUES));
 
          PetscInt col_ad = i_gap + _n_gaps * iz_ind;
          value_ad = -1.0 * term_in * (1.0 - alpha) + term_out * alpha;
          LibmeshPetscCall(MatSetValues(
              _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
 
          col_ad = i_gap + _n_gaps * (iz_ind + 1);
          value_ad = term_out * (1.0 - alpha);
          LibmeshPetscCall(MatSetValues(
              _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
        }
        else if (iz == last_node)
        {
          PetscInt row_ad = i_gap + _n_gaps * iz_ind;
          PetscInt col_ad = i_gap + _n_gaps * (iz_ind - 1);
          PetscScalar value_ad = -1.0 * term_in * alpha;
          LibmeshPetscCall(MatSetValues(
              _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
 
          col_ad = i_gap + _n_gaps * iz_ind;
          value_ad = -1.0 * term_in * (1.0 - alpha) + term_out * alpha;
          LibmeshPetscCall(MatSetValues(
              _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
 
          value_ad = -1.0 * term_out * (1.0 - alpha) * _Wij(i_gap, iz);
          LibmeshPetscCall(
              VecSetValues(_cmc_advective_derivative_rhs, 1, &row_ad, &value_ad, ADD_VALUES));
        }
        else
        {
          PetscInt row_ad = i_gap + _n_gaps * iz_ind;
          PetscInt col_ad = i_gap + _n_gaps * (iz_ind - 1);
          PetscScalar value_ad = -1.0 * term_in * alpha;
          LibmeshPetscCall(MatSetValues(
              _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
 
          col_ad = i_gap + _n_gaps * iz_ind;
          value_ad = -1.0 * term_in * (1.0 - alpha) + term_out * alpha;
          LibmeshPetscCall(MatSetValues(
              _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
 
          col_ad = i_gap + _n_gaps * (iz_ind + 1);
          value_ad = term_out * (1.0 - alpha);
          LibmeshPetscCall(MatSetValues(
              _cmc_advective_derivative_mat, 1, &row_ad, 1, &col_ad, &value_ad, INSERT_VALUES));
        }
        // Assembling friction force
        PetscInt row_ff = i_gap + _n_gaps * iz_ind;
        PetscInt col_ff = i_gap + _n_gaps * iz_ind;
        PetscScalar value_ff = _kij * std::abs(_Wij(i_gap, iz)) / 2.0;
        LibmeshPetscCall(MatSetValues(
            _cmc_friction_force_mat, 1, &row_ff, 1, &col_ff, &value_ff, INSERT_VALUES));
 
        // Assembling pressure force
        alpha = computeInterpolationCoefficients(Pe);
 
        PetscScalar pressure_factor = Utility::pow<2>(Sij) * rho_star;
        PetscInt row_pf = i_gap + _n_gaps * iz_ind;
        PetscInt col_pf = i_ch + _n_channels * iz_ind;
        PetscScalar value_pf = -1.0 * alpha * pressure_factor;
        LibmeshPetscCall(
            MatSetValues(_cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
        col_pf = j_ch + _n_channels * iz_ind;
        value_pf = alpha * pressure_factor;
        LibmeshPetscCall(
            MatSetValues(_cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
 
        if (iz == last_node)
        {
          PetscInt row_pf = i_gap + _n_gaps * iz_ind;
          PetscScalar value_pf = (1.0 - alpha) * pressure_factor * (*_P_soln)(node_out_i);
          LibmeshPetscCall(
              VecSetValues(_cmc_pressure_force_rhs, 1, &row_pf, &value_pf, ADD_VALUES));
          value_pf = -1.0 * (1.0 - alpha) * pressure_factor * (*_P_soln)(node_out_j);
          LibmeshPetscCall(
              VecSetValues(_cmc_pressure_force_rhs, 1, &row_pf, &value_pf, ADD_VALUES));
        }
        else
        {
          row_pf = i_gap + _n_gaps * iz_ind;
          col_pf = i_ch + _n_channels * (iz_ind + 1);
          value_pf = -1.0 * (1.0 - alpha) * pressure_factor;
          LibmeshPetscCall(
              MatSetValues(_cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
          col_pf = j_ch + _n_channels * (iz_ind + 1);
          value_pf = (1.0 - alpha) * pressure_factor;
          LibmeshPetscCall(
              MatSetValues(_cmc_pressure_force_mat, 1, &row_pf, 1, &col_pf, &value_pf, ADD_VALUES));
        }
      }
    }
    LibmeshPetscCall(MatZeroEntries(_cmc_sys_Wij_mat));
    LibmeshPetscCall(VecZeroEntries(_cmc_sys_Wij_rhs));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_time_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_advective_derivative_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_friction_force_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_friction_force_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_pressure_force_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    // Matrix
#if !PETSC_VERSION_LESS_THAN(3, 15, 0)
    LibmeshPetscCall(
        MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_time_derivative_mat, UNKNOWN_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_advective_derivative_mat, UNKNOWN_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_friction_force_mat, UNKNOWN_NONZERO_PATTERN));
#else
    LibmeshPetscCall(
        MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_time_derivative_mat, DIFFERENT_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_advective_derivative_mat, DIFFERENT_NONZERO_PATTERN));
    LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(
        MatAXPY(_cmc_sys_Wij_mat, 1.0, _cmc_friction_force_mat, DIFFERENT_NONZERO_PATTERN));
#endif
    LibmeshPetscCall(MatAssemblyBegin(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_cmc_sys_Wij_mat, MAT_FINAL_ASSEMBLY));
    // RHS
    LibmeshPetscCall(VecAXPY(_cmc_sys_Wij_rhs, 1.0, _cmc_time_derivative_rhs));
    LibmeshPetscCall(VecAXPY(_cmc_sys_Wij_rhs, 1.0, _cmc_advective_derivative_rhs));
    LibmeshPetscCall(VecAXPY(_cmc_sys_Wij_rhs, 1.0, _cmc_friction_force_rhs));
 
    if (_segregated_bool)
    {
      // Assembly the matrix system
      Vec sol_holder_P;
      LibmeshPetscCall(createPetscVector(sol_holder_P, _block_size * _n_gaps));
      Vec sol_holder_W;
      LibmeshPetscCall(createPetscVector(sol_holder_W, _block_size * _n_gaps));
      LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
          _prodp, *_P_soln, first_node - 1, last_node - 1, _n_channels));
      LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
          _Wij_vec, _Wij, first_node, last_node, _n_gaps));
      LibmeshPetscCall(MatMult(_cmc_sys_Wij_mat, _Wij_vec, sol_holder_W));
      LibmeshPetscCall(VecAXPY(sol_holder_W, -1.0, _cmc_sys_Wij_rhs));
      LibmeshPetscCall(MatMult(_cmc_pressure_force_mat, _prodp, sol_holder_P));
      LibmeshPetscCall(VecAXPY(sol_holder_P, -1.0, _cmc_pressure_force_rhs));
      LibmeshPetscCall(VecAXPY(sol_holder_W, 1.0, sol_holder_P));
      PetscScalar * xx;
      LibmeshPetscCall(VecGetArray(sol_holder_W, &xx));
      for (unsigned int iz = first_node; iz < last_node + 1; iz++)
      {
        auto iz_ind = iz - first_node;
        for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
        {
          _Wij_residual_matrix(i_gap, iz - 1 - iblock * _block_size) = xx[iz_ind * _n_gaps + i_gap];
        }
      }
      LibmeshPetscCall(VecDestroy(&sol_holder_P));
      LibmeshPetscCall(VecDestroy(&sol_holder_W));
    }
  }
}
 
void
SubChannel1PhaseProblem::computeWijPrime(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  _mixing_closure->computeBlockMixingParameters(first_node, last_node);
  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
  {
    auto dz = _z_grid[iz] - _z_grid[iz - 1];
    for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
    {
      auto chans = _subchannel_mesh.getGapChannels(i_gap);
      unsigned int i_ch = chans.first;
      unsigned int j_ch = chans.second;
      auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
      auto * node_out_i = _subchannel_mesh.getChannelNode(i_ch, iz);
      auto * node_in_j = _subchannel_mesh.getChannelNode(j_ch, iz - 1);
      auto * node_out_j = _subchannel_mesh.getChannelNode(j_ch, iz);
      auto Si_in = (*_S_flow_soln)(node_in_i);
      auto Sj_in = (*_S_flow_soln)(node_in_j);
      auto Si_out = (*_S_flow_soln)(node_out_i);
      auto Sj_out = (*_S_flow_soln)(node_out_j);
      auto gap = _subchannel_mesh.getGapWidth(iz, i_gap);
      auto Sij = dz * gap;
      // Mass flux of the gap: total mass flow rate of the two subchannels over their total flow
      // area, averaged over the inlet and outlet of the cell. This is the area-weighted mean of rho
      // V of the two subchannels; Eq. (30) of Pacio et al. (2022) uses the arithmetic mean rho (V_i
      // + V_j) / 2, which differs when the two subchannels have different flow areas
      auto avg_massflux =
          0.5 * (((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j)) / (Si_in + Sj_in) +
                 ((*_mdot_soln)(node_out_i) + (*_mdot_soln)(node_out_j)) / (Si_out + Sj_out));
      auto beta = computeMixingParameter(i_gap, iz);
 
      if (!_implicit_bool)
      {
        _WijPrime(i_gap, iz) = beta * avg_massflux * Sij;
      }
      else
      {
        auto iz_ind = iz - first_node;
        PetscScalar base_value = beta * 0.5 * Sij;
 
        // Bottom values
        if (iz == first_node)
        {
          PetscScalar value_tl = -1.0 * base_value / (Si_in + Sj_in) *
                                 ((*_mdot_soln)(node_in_i) + (*_mdot_soln)(node_in_j));
          PetscInt row = i_gap + _n_gaps * iz_ind;
          LibmeshPetscCall(
              VecSetValues(_amc_turbulent_cross_flows_rhs, 1, &row, &value_tl, INSERT_VALUES));
        }
        else
        {
          PetscScalar value_tl = base_value / (Si_in + Sj_in);
          PetscInt row = i_gap + _n_gaps * iz_ind;
 
          PetscInt col_ich = i_ch + _n_channels * (iz_ind - 1);
          LibmeshPetscCall(MatSetValues(
              _amc_turbulent_cross_flows_mat, 1, &row, 1, &col_ich, &value_tl, INSERT_VALUES));
 
          PetscInt col_jch = j_ch + _n_channels * (iz_ind - 1);
          LibmeshPetscCall(MatSetValues(
              _amc_turbulent_cross_flows_mat, 1, &row, 1, &col_jch, &value_tl, INSERT_VALUES));
        }
 
        // Top values
        PetscScalar value_bl = base_value / (Si_out + Sj_out);
        PetscInt row = i_gap + _n_gaps * iz_ind;
 
        PetscInt col_ich = i_ch + _n_channels * iz_ind;
        LibmeshPetscCall(MatSetValues(
            _amc_turbulent_cross_flows_mat, 1, &row, 1, &col_ich, &value_bl, INSERT_VALUES));
 
        PetscInt col_jch = j_ch + _n_channels * iz_ind;
        LibmeshPetscCall(MatSetValues(
            _amc_turbulent_cross_flows_mat, 1, &row, 1, &col_jch, &value_bl, INSERT_VALUES));
      }
    }
  }
 
  if (_implicit_bool)
  {
    LibmeshPetscCall(MatAssemblyBegin(_amc_turbulent_cross_flows_mat, MAT_FINAL_ASSEMBLY));
    LibmeshPetscCall(MatAssemblyEnd(_amc_turbulent_cross_flows_mat, MAT_FINAL_ASSEMBLY));
 
    Vec loc_prod;
    Vec loc_Wij;
    LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &loc_prod));
    LibmeshPetscCall(VecDuplicate(_Wij_vec, &loc_Wij));
    LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
        loc_prod, *_mdot_soln, first_node, last_node, _n_channels));
    LibmeshPetscCall(MatMult(_amc_turbulent_cross_flows_mat, loc_prod, loc_Wij));
    LibmeshPetscCall(VecAXPY(loc_Wij, -1.0, _amc_turbulent_cross_flows_rhs));
    LibmeshPetscCall(populateDenseFromVector<libMesh::DenseMatrix<Real>>(
        loc_Wij, _WijPrime, first_node, last_node, _n_gaps));
    LibmeshPetscCall(VecDestroy(&loc_prod));
    LibmeshPetscCall(VecDestroy(&loc_Wij));
  }
}
 
Real
SubChannel1PhaseProblem::computeMixingParameter(unsigned int i_gap, unsigned int iz) const
{
  auto beta = _mixing_closure->computeMixingParameter(i_gap, iz);
  if (!std::isfinite(beta) || beta < 0.0)
    mooseError(name(),
               ": Mixing closure returned invalid beta = ",
               beta,
               " for gap ",
               i_gap,
               " at axial index ",
               iz,
               ". Beta must be finite and non-negative.");
 
  return beta;
}
 
Real
SubChannel1PhaseProblem::computeSweepFlowMixingParameter(unsigned int i_gap, unsigned int iz) const
{
  auto beta = _mixing_closure->computeSweepFlowMixingParameter(i_gap, iz);
  if (!std::isfinite(beta) || beta < 0.0)
    mooseError(name(),
               ": Mixing closure returned invalid sweep-flow coefficient = ",
               beta,
               " for gap ",
               i_gap,
               " at axial index ",
               iz,
               ". sweep-flow coefficient must be finite and non-negative.");
 
  return beta;
}
 
libMesh::DenseVector<Real>
SubChannel1PhaseProblem::residualFunction(int iblock, libMesh::DenseVector<Real> solution)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  libMesh::DenseVector<Real> Wij_residual_vector(_n_gaps * _block_size, 0.0);
  // Assign the solution to the cross-flow matrix
  int i = 0;
  for (unsigned int iz = first_node; iz < last_node + 1; iz++)
  {
    for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
    {
      _Wij(i_gap, iz) = solution(i);
      i++;
    }
  }
 
  // Calculating sum of crossflows
  computeSumWij(iblock);
  // Solving axial flux
  computeMdot(iblock);
  // Calculation of turbulent Crossflow
  computeWijPrime(iblock);
  // Solving for Pressure Drop
  computeDP(iblock);
  // Solving for pressure
  computeP(iblock);
  // Populating lateral crossflow residual matrix
  computeWijResidual(iblock);
 
  // Turn the residual matrix into a residual vector
  for (unsigned int iz = 0; iz < _block_size; iz++)
  {
    for (unsigned int i_gap = 0; i_gap < _n_gaps; i_gap++)
    {
      int i = _n_gaps * iz + i_gap; // column wise transfer
      Wij_residual_vector(i) = _Wij_residual_matrix(i_gap, iz);
    }
  }
  return Wij_residual_vector;
}
 
PetscErrorCode
SubChannel1PhaseProblem::petscSnesSolver(int iblock,
                                         const libMesh::DenseVector<Real> & solution,
                                         libMesh::DenseVector<Real> & root)
{
  SNES snes;
  KSP ksp;
  PC pc;
  Vec x, r;
  PetscScalar * xx;
 
  PetscFunctionBegin;
  LibmeshPetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
  LibmeshPetscCall(VecCreate(PETSC_COMM_SELF, &x));
  LibmeshPetscCall(VecSetSizes(x, PETSC_DECIDE, _block_size * _n_gaps));
  LibmeshPetscCall(VecSetFromOptions(x));
  LibmeshPetscCall(VecDuplicate(x, &r));
 
#if PETSC_VERSION_LESS_THAN(3, 13, 0)
  LibmeshPetscCall(PetscOptionsSetValue(PETSC_NULL, "-snes_mf", PETSC_NULL));
#else
  LibmeshPetscCall(SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE));
#endif
  Ctx ctx;
  ctx.iblock = iblock;
  ctx.schp = this;
  LibmeshPetscCall(SNESSetFunction(snes, r, formFunction, &ctx));
  LibmeshPetscCall(SNESGetKSP(snes, &ksp));
  LibmeshPetscCall(KSPGetPC(ksp, &pc));
  LibmeshPetscCall(PCSetType(pc, PCNONE));
  LibmeshPetscCall(KSPSetTolerances(ksp, _rtol, _atol, _dtol, _maxit));
  LibmeshPetscCall(SNESSetFromOptions(snes));
  LibmeshPetscCall(VecGetArray(x, &xx));
  for (unsigned int i = 0; i < _block_size * _n_gaps; i++)
  {
    xx[i] = solution(i);
  }
  LibmeshPetscCall(VecRestoreArray(x, &xx));
 
  LibmeshPetscCall(SNESSolve(snes, NULL, x));
  // A diverged solve leaves x at its initial guess, which is then used as the crossflow of this
  // iteration. Report it rather than stopping, since the unpreconditioned linear solve can break
  // down near convergence
  SNESConvergedReason reason;
  LibmeshPetscCall(SNESGetConvergedReason(snes, &reason));
  if (reason < 0)
  {
    PetscInt iterations;
    PetscReal residual_norm;
    LibmeshPetscCall(SNESGetIterationNumber(snes, &iterations));
    LibmeshPetscCall(SNESGetFunctionNorm(snes, &residual_norm));
    mooseWarning("Crossflow nonlinear solve did not converge: ",
                 SNESConvergedReasons[reason],
                 " (",
                 static_cast<int>(reason),
                 ") after ",
                 iterations,
                 " iterations; residual norm = ",
                 residual_norm,
                 ". The crossflow of this iteration is not updated.");
  }
  LibmeshPetscCall(VecGetArray(x, &xx));
  for (unsigned int i = 0; i < _block_size * _n_gaps; i++)
    root(i) = xx[i];
 
  LibmeshPetscCall(VecRestoreArray(x, &xx));
  LibmeshPetscCall(VecDestroy(&x));
  LibmeshPetscCall(VecDestroy(&r));
  LibmeshPetscCall(SNESDestroy(&snes));
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}
 
PetscErrorCode
SubChannel1PhaseProblem::solveAndPopulateEnthalpy(
    Mat A, Vec rhs, unsigned int first_node, unsigned int last_node, const char * ksp_prefix)
{
  PetscFunctionBegin;
 
  // Create solution vector with rhs layout
  Vec x = nullptr;
  LibmeshPetscCall(VecDuplicate(rhs, &x));
 
  // KSP setup
  KSP ksp = nullptr;
  PC pc = nullptr;
  LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
  LibmeshPetscCall(KSPSetOperators(ksp, A, A));
  LibmeshPetscCall(KSPGetPC(ksp, &pc));
  LibmeshPetscCall(PCSetType(pc, PCJACOBI));
  LibmeshPetscCall(KSPSetTolerances(ksp, _rtol, _atol, _dtol, _maxit));
  if (ksp_prefix && *ksp_prefix)
    LibmeshPetscCall(KSPSetOptionsPrefix(ksp, ksp_prefix));
  LibmeshPetscCall(KSPSetFromOptions(ksp));
 
  // Solve
  LibmeshPetscCall(KSPSolve(ksp, rhs, x));
  KSPConvergedReason reason;
  LibmeshPetscCall(KSPGetConvergedReason(ksp, &reason));
  if (reason < 0)
  {
    PetscInt iterations;
    PetscReal residual_norm;
    LibmeshPetscCall(KSPGetIterationNumber(ksp, &iterations));
    LibmeshPetscCall(KSPGetResidualNorm(ksp, &residual_norm));
    mooseError(name(),
               ": enthalpy linear solve failed: ",
               KSPConvergedReasons[reason],
               " (",
               static_cast<int>(reason),
               ") after ",
               iterations,
               " iterations; residual norm = ",
               residual_norm,
               ".");
  }
 
  // Scatter to _h_soln with sanity checks
  PetscScalar * xx = nullptr;
  LibmeshPetscCall(VecGetArray(x, &xx));
  for (unsigned int iz = first_node; iz <= last_node; ++iz)
  {
    const unsigned int iz_ind = iz - first_node;
    for (unsigned int i_ch = 0; i_ch < _n_channels; ++i_ch)
    {
      auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
      const PetscScalar h_out = xx[iz_ind * _n_channels + i_ch];
      if (h_out < 0.0)
        mooseError(
            name(), " : Calculation of negative Enthalpy h_out = ", h_out, " Axial Level = ", iz);
      _h_soln->set(node_out, h_out);
    }
  }
  LibmeshPetscCall(VecRestoreArray(x, &xx));
 
  // Cleanup
  LibmeshPetscCall(KSPDestroy(&ksp));
  LibmeshPetscCall(VecDestroy(&x));
 
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}
 
Real
SubChannel1PhaseProblem::computeAddedHeatDuct(unsigned int i_ch, unsigned int iz) const
{
  mooseAssert(iz > 0, "Trapezoidal rule requires starting at index 1 at least");
  if (_duct_mesh_exist)
  {
    auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
    if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
    {
      auto dz = _z_grid[iz] - _z_grid[iz - 1];
      auto * node_in_chan = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
      auto * node_out_chan = _subchannel_mesh.getChannelNode(i_ch, iz);
      auto * node_in_duct = _subchannel_mesh.getDuctNodeFromChannel(node_in_chan);
      auto * node_out_duct = _subchannel_mesh.getDuctNodeFromChannel(node_out_chan);
      auto heat_rate_in = (*_duct_heat_flux_soln)(node_in_duct);
      auto heat_rate_out = (*_duct_heat_flux_soln)(node_out_duct);
      auto width = getSubChannelPeripheralDuctWidth(i_ch);
      return 0.5 * (heat_rate_in + heat_rate_out) * dz * width;
    }
    else
    {
      return 0.0;
    }
  }
  else
  {
    return 0.0;
  }
}
 
PetscErrorCode
SubChannel1PhaseProblem::implicitPetscSolve(int iblock)
{
  PetscFunctionBegin;
  // ---------- helper functions -----------------------------
  auto V = [&](const std::string & s)
  {
    if (_verbose_subchannel)
      _console << s << std::endl;
  };
 
  auto DupMatAssembled = [&](Mat src, Mat * dst)
  {
    if (src)
    {
      LibmeshPetscCall(MatDuplicate(src, MAT_COPY_VALUES, dst));
      LibmeshPetscCall(MatAssemblyBegin(*dst, MAT_FINAL_ASSEMBLY));
      LibmeshPetscCall(MatAssemblyEnd(*dst, MAT_FINAL_ASSEMBLY));
    }
    else
      *dst = NULL;
  };
 
  auto DupVecCopy = [&](Vec src, Vec * dst)
  {
    LibmeshPetscCall(VecDuplicate(src, dst));
    LibmeshPetscCall(VecCopy(src, *dst));
  };
 
  const PetscInt Q = 3; // [mass conservation, axial momentum, cross momentum]
 
  // small indexer
  auto Idx = [&](PetscInt r, PetscInt c) { return Q * r + c; };
 
  // arrays that MUST be declared before lambdas use them
  std::vector<Mat> mat_array(Q * Q, NULL);
  std::vector<Vec> vec_array(Q, NULL);
 
  // generic assembler for one governing equation row in the nested matrix
  auto AssembleEquation = [&](PetscInt f,
                              Mat A0,
                              Mat A1,
                              Mat A2,             // three blocks in row f (can be nullptr)
                              Vec rhs,            // base RHS for equation f
                              Vec rhs_add,        // optional extra RHS to add (can be nullptr)
                              const char * label) // e.g. "Mass", "Lin mom", "Cross mom"
  {
    DupMatAssembled(A0, &mat_array[Idx(f, 0)]);
    DupMatAssembled(A1, &mat_array[Idx(f, 1)]);
    DupMatAssembled(A2, &mat_array[Idx(f, 2)]);
    DupVecCopy(rhs, &vec_array[f]);
    if (rhs_add)
      LibmeshPetscCall(VecAXPY(vec_array[f], 1.0, rhs_add));
    V(std::string(label) + " system assembled");
  };
 
  // -----------------------------------------------------------------------------
  // Helper lambda that applies per-equation under-relaxation by modifying BOTH the
  // matrix block and the RHS for that equation.
  //
  // Specifically, with A_ff for the equation and D = diag(A_ff):
  //   1) Matrix diagonal scaling: D <- D / alpha, then A_ff's diagonal is replaced
  //      with D. For alpha < 1 this increases diagonal dominance (more damping).
  //   2) RHS blending with the previous solution x_old:
  //        rhs_f <- rhs_f + (1 - alpha) * (D / alpha) * x_old
  //      where x_old is provided by the caller via `populate(work)`.
  //
  // Net effect: the solved x satisfies
  //        A_ff x = rhs_f_original + ((1 - alpha)/alpha) * D * (x_old - x),
  // which damps updates toward x_old without changing the converged solution.
  // This equation relaxation is independent of the optional post-solve solution relaxation.
  // -----------------------------------------------------------------------------
  auto relaxEquation =
      [&](Mat diagonal_block, Vec rhs, Vec work, const Real relaxation, auto && populate)
  {
    if (relaxation == 1.0)
      return;
 
    Vec diagonal = nullptr;
    LibmeshPetscCall(VecDuplicate(rhs, &diagonal));
 
    // 1) diagonal_block: diag <- diag / relaxation
    LibmeshPetscCall(MatGetDiagonal(diagonal_block, diagonal));
    LibmeshPetscCall(VecScale(diagonal, 1.0 / relaxation));
    LibmeshPetscCall(MatDiagonalSet(diagonal_block, diagonal, INSERT_VALUES));
 
    // 2) work <- x_old (caller-provided populator)
    LibmeshPetscCall(populate(work));
 
    // 3) rhs += (1 - relaxation) * (diag .* work)
    LibmeshPetscCall(VecScale(diagonal, 1.0 - relaxation));
    LibmeshPetscCall(VecPointwiseMult(work, work, diagonal));
    LibmeshPetscCall(VecAXPY(rhs, 1.0, work));
 
    LibmeshPetscCall(VecDestroy(&diagonal));
  };
 
  // indices
  const unsigned int first_node = iblock * _block_size + 1;
  const unsigned int last_node = (iblock + 1) * _block_size;
 
  // ---------- assemble per-block operators -----------------
  computeSumWij(iblock);
  computeMdot(iblock);
  computeWijPrime(iblock);
  computeDP(iblock);
  computeP(iblock);
  computeWijResidual(iblock);
 
  V("Starting nested system.");
 
  // Populate nested matrix with the individual physics
  // equation 0: Mass conservation
  AssembleEquation(/*f=*/0,
                   /*A0=*/_mc_axial_convection_mat,
                   /*A1=*/_mc_density_pressure_mat,
                   /*A2=*/_mc_sumWij_mat,
                   /*rhs=*/_mc_axial_convection_rhs,
                   /*rhs_add=*/nullptr,
                   /*label=*/"Mass");
 
  // equation 1: Axial momentum conservation
  AssembleEquation(/*f=*/1,
                   /*A0=*/_amc_sys_mdot_mat,
                   /*A1=*/_amc_pressure_force_mat,
                   /*A2=*/nullptr,
                   /*rhs=*/_amc_pressure_force_rhs,
                   /*rhs_add=*/_amc_sys_mdot_rhs,
                   /*label=*/"Lin mom");
 
  // equation 2: Cross momentum conservation
  AssembleEquation(/*f=*/2,
                   /*A0=*/nullptr,
                   /*A1=*/_cmc_pressure_force_mat,
                   /*A2=*/_cmc_sys_Wij_mat,
                   /*rhs=*/_cmc_sys_Wij_rhs,
                   /*rhs_add=*/_cmc_pressure_force_rhs,
                   /*label=*/"Cross mom");
 
  // ========================== Relaxation  ====================
  if (true)
  {
    LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
        _prod, *_mdot_soln, first_node, last_node, _n_channels));
 
    Vec mdot_estimate;
    LibmeshPetscCall(createPetscVector(mdot_estimate, _block_size * _n_channels));
    Vec pmat_diag;
    LibmeshPetscCall(createPetscVector(pmat_diag, _block_size * _n_channels));
    Vec p_estimate;
    LibmeshPetscCall(createPetscVector(p_estimate, _block_size * _n_channels));
    Vec unity_vec;
    LibmeshPetscCall(createPetscVector(unity_vec, _block_size * _n_channels));
    LibmeshPetscCall(VecSet(unity_vec, 1.0));
    Vec sol_holder_P;
    LibmeshPetscCall(createPetscVector(sol_holder_P, _block_size * _n_gaps));
    Vec unity_vec_Wij;
    LibmeshPetscCall(createPetscVector(unity_vec_Wij, _block_size * _n_gaps));
    LibmeshPetscCall(VecSet(unity_vec_Wij, 1.0));
    Vec _Wij_loc_vec;
    LibmeshPetscCall(createPetscVector(_Wij_loc_vec, _block_size * _n_gaps));
    Vec _Wij_old_loc_vec;
    LibmeshPetscCall(createPetscVector(_Wij_old_loc_vec, _block_size * _n_gaps));
 
    // ---- scale estimates ----
    // mdot_estimate = A(1,0) * mdot
    LibmeshPetscCall(MatMult(mat_array[Q /* (1,0) */], _prod, mdot_estimate));
 
    // p_estimate = mdot_est / (diag(A(1,1)) + eps)
    LibmeshPetscCall(MatGetDiagonal(mat_array[Q + 1], pmat_diag));
    LibmeshPetscCall(VecAXPY(pmat_diag, 1e-10, unity_vec));
    LibmeshPetscCall(VecPointwiseDivide(p_estimate, mdot_estimate, pmat_diag));
 
    // sol_holder_P = A(2,1) * p_estimate - rhs_cmc_pressure
    LibmeshPetscCall(MatMult(mat_array[2 * Q + 1], p_estimate, sol_holder_P));
    LibmeshPetscCall(VecAXPY(sol_holder_P, -1.0, _cmc_pressure_force_rhs));
 
    // sumWij_loc from sol_holder_P (accumulate)
    Vec sumWij_loc;
    LibmeshPetscCall(createPetscVector(sumWij_loc, _block_size * _n_channels));
    for (unsigned int iz = first_node; iz <= last_node; ++iz)
    {
      const auto iz_ind = iz - first_node;
      for (unsigned int i_ch = 0; i_ch < _n_channels; ++i_ch)
      {
        PetscScalar sumWij = 0.0;
        unsigned int counter = 0;
        for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
        {
          PetscInt row_vec = i_gap + _n_gaps * iz_ind;
          PetscScalar loc_Wij_value;
          LibmeshPetscCall(VecGetValues(sol_holder_P, 1, &row_vec, &loc_Wij_value));
          sumWij += _subchannel_mesh.getCrossflowSign(i_ch, counter) * loc_Wij_value;
          counter++;
        }
        PetscInt row_vec = i_ch + _n_channels * iz_ind;
        LibmeshPetscCall(VecSetValues(sumWij_loc, 1, &row_vec, &sumWij, INSERT_VALUES));
      }
    }
    LibmeshPetscCall(VecAssemblyBegin(sumWij_loc));
    LibmeshPetscCall(VecAssemblyEnd(sumWij_loc));
 
    // ---- robust scale measurements ----
    PetscScalar min_mdot, sum_mdot, avg_mdot;
    PetscInt n;
    LibmeshPetscCall(VecAbs(_prod));
    LibmeshPetscCall(VecMin(_prod, NULL, &min_mdot));
    LibmeshPetscCall(VecSum(_prod, &sum_mdot)); // Sum all entries
    LibmeshPetscCall(VecGetSize(_prod, &n));    // Number of entries
    avg_mdot = sum_mdot / static_cast<PetscScalar>(n);
 
    V("Average estimated mdot: " + std::to_string(avg_mdot));
    V("Minimum estimated mdot: " + std::to_string(min_mdot));
 
    LibmeshPetscCall(VecAbs(sumWij_loc));
    LibmeshPetscCall(VecMax(sumWij_loc, NULL, &_max_sumWij));
    _max_sumWij = std::max(1e-10, _max_sumWij);
    V("Maximum estimated Wij: " + std::to_string(_max_sumWij));
 
    LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
        _Wij_loc_vec, _Wij, first_node, last_node, _n_gaps));
    LibmeshPetscCall(VecAbs(_Wij_loc_vec));
    LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
        _Wij_old_loc_vec, _Wij_old, first_node, last_node, _n_gaps));
    LibmeshPetscCall(VecAbs(_Wij_old_loc_vec));
    LibmeshPetscCall(VecAXPY(_Wij_loc_vec, -1.0, _Wij_old_loc_vec));
 
    PetscScalar relax_factor;
    LibmeshPetscCall(VecAbs(_Wij_loc_vec));
#if !PETSC_VERSION_LESS_THAN(3, 16, 0)
    LibmeshPetscCall(VecMean(_Wij_loc_vec, &relax_factor));
#else
    VecSum(_Wij_loc_vec, &relax_factor);
    relax_factor /= _block_size * _n_gaps;
#endif
    relax_factor = relax_factor / _max_sumWij + 0.5;
    V("Relax base value: " + std::to_string(relax_factor));
 
    // ---- crossflow resistance inflation ----
    const PetscScalar resistance_relaxation = 0.9;
    _added_K = _max_sumWij / min_mdot;
    V("New cross resistance: " + std::to_string(_added_K));
    _added_K = (_added_K * resistance_relaxation + (1.0 - resistance_relaxation) * _added_K_old) *
               relax_factor;
    V("Relaxed cross resistance: " + std::to_string(_added_K));
 
    // Snap-up lower-bounding
    if (_added_K < 10 && _added_K >= 1.0)
      _added_K = 1.0;
    if (_added_K < 1.0 && _added_K >= 0.1)
      _added_K = 0.5;
    if (_added_K < 0.1 && _added_K >= 0.01)
      _added_K = 1. / 3.;
    if (_added_K < 1e-2 && _added_K >= 1e-3)
      _added_K = 0.1;
    V("Actual added cross resistance: " + std::to_string(_added_K));
    LibmeshPetscCall(VecScale(unity_vec_Wij, _added_K));
    _added_K_old = _added_K;
 
    LibmeshPetscCall(MatDiagonalSet(mat_array[2 * Q + 2], unity_vec_Wij, ADD_VALUES));
 
    // ---- cleanup temp vectors used above ----
    LibmeshPetscCall(VecDestroy(&mdot_estimate));
    LibmeshPetscCall(VecDestroy(&pmat_diag));
    LibmeshPetscCall(VecDestroy(&unity_vec));
    LibmeshPetscCall(VecDestroy(&p_estimate));
    LibmeshPetscCall(VecDestroy(&sol_holder_P));
    LibmeshPetscCall(VecDestroy(&unity_vec_Wij));
    LibmeshPetscCall(VecDestroy(&sumWij_loc));
    LibmeshPetscCall(VecDestroy(&_Wij_loc_vec));
    LibmeshPetscCall(VecDestroy(&_Wij_old_loc_vec));
  }
 
  V("Relax mdot: " + std::to_string(_mass_flow_equation_relaxation));
  V("Relax P: " + std::to_string(_pressure_equation_relaxation));
  V("Relax Wij: " + std::to_string(_crossflow_equation_relaxation));
 
  relaxEquation(mat_array[Idx(0, 0)],
                vec_array[0],
                _prod,
                _mass_flow_equation_relaxation,
                [&](Vec work)
                {
                  return populateVectorFromHandle<SolutionHandle>(
                      work, *_mdot_soln, first_node, last_node, _n_channels);
                });
  V("mdot relaxed");
 
  relaxEquation(mat_array[Idx(1, 1)],
                vec_array[1],
                _prodp,
                _pressure_equation_relaxation,
                [&](Vec work)
                {
                  // The pressure block represents the node below each axial momentum cell because
                  // the outlet pressure is prescribed.
                  return populateVectorFromHandle<SolutionHandle>(
                      work, *_P_soln, first_node - 1, last_node - 1, _n_channels);
                });
  V("P relaxed");
 
  relaxEquation(mat_array[Idx(2, 2)],
                vec_array[2],
                _Wij_vec,
                _crossflow_equation_relaxation,
                [&](Vec work)
                {
                  return populateVectorFromDense<libMesh::DenseMatrix<Real>>(
                      work, _Wij, first_node, last_node, _n_gaps);
                });
  V("Wij relaxed");
  V("Linear solver relaxed");
 
  // ======================== Create and configure KSP =========================
  Mat A_nest;
  Vec b_nest;
  Vec x_nest;
  LibmeshPetscCall(MatCreateNest(PETSC_COMM_SELF, Q, NULL, Q, NULL, mat_array.data(), &A_nest));
  LibmeshPetscCall(VecCreateNest(PETSC_COMM_SELF, Q, NULL, vec_array.data(), &b_nest));
  V("Nested system created");
 
  KSP ksp;
  PC pc;
  LibmeshPetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
  LibmeshPetscCall(KSPSetOptionsPrefix(ksp, "scm_coupled_"));
  LibmeshPetscCall(KSPSetType(ksp, KSPFGMRES));
  LibmeshPetscCall(KSPSetOperators(ksp, A_nest, A_nest));
  LibmeshPetscCall(KSPGetPC(ksp, &pc));
  LibmeshPetscCall(PCSetType(pc, PCFIELDSPLIT));
  LibmeshPetscCall(KSPSetTolerances(ksp, _rtol, _atol, _dtol, _maxit));
 
  // split equations
  std::vector<IS> rows(Q);
  LibmeshPetscCall(MatNestGetISs(A_nest, rows.data(), NULL));
  for (PetscInt j = 0; j < Q; ++j)
  {
    IS part;
    LibmeshPetscCall(ISDuplicate(rows[j], &part));
    LibmeshPetscCall(PCFieldSplitSetIS(pc, NULL, part));
    LibmeshPetscCall(ISDestroy(&part));
  }
  LibmeshPetscCall(KSPSetFromOptions(ksp));
  V("Linear solver assembled");
 
  // ============================== Solve =====================================
  LibmeshPetscCall(VecDuplicate(b_nest, &x_nest));
  LibmeshPetscCall(VecSet(x_nest, 0.0));
  LibmeshPetscCall(KSPSolve(ksp, b_nest, x_nest));
  KSPConvergedReason reason;
  LibmeshPetscCall(KSPGetConvergedReason(ksp, &reason));
  if (reason < 0)
  {
    PetscInt iterations;
    PetscReal residual_norm;
    LibmeshPetscCall(KSPGetIterationNumber(ksp, &iterations));
    LibmeshPetscCall(KSPGetResidualNorm(ksp, &residual_norm));
    mooseError(name(),
               ": coupled mass/momentum linear solve failed: ",
               KSPConvergedReasons[reason],
               " (",
               static_cast<int>(reason),
               ") after ",
               iterations,
               " iterations; residual norm = ",
               residual_norm,
               ".");
  }
 
  // destroy solver containers first
  LibmeshPetscCall(VecDestroy(&b_nest));
  LibmeshPetscCall(MatDestroy(&A_nest));
  LibmeshPetscCall(KSPDestroy(&ksp));
  for (PetscInt i = 0; i < Q * Q; i++)
    LibmeshPetscCall(MatDestroy(&mat_array[i]));
  for (PetscInt i = 0; i < Q; i++)
    LibmeshPetscCall(VecDestroy(&vec_array[i]));
  V("Solver elements destroyed");
 
  // ====================== Extract & scatter the solution =====================
  Vec sol_mdot, sol_p, sol_Wij;
  V("Vectors to hold solution created");
  PetscInt num_vecs;
  Vec * loc_vecs;
  LibmeshPetscCall(VecNestGetSubVecs(x_nest, &num_vecs, &loc_vecs));
  LibmeshPetscCall(VecDuplicate(_mc_axial_convection_rhs, &sol_mdot));
  LibmeshPetscCall(VecCopy(loc_vecs[0], sol_mdot));
  LibmeshPetscCall(VecDuplicate(_amc_sys_mdot_rhs, &sol_p));
  LibmeshPetscCall(VecCopy(loc_vecs[1], sol_p));
  LibmeshPetscCall(VecDuplicate(_cmc_sys_Wij_rhs, &sol_Wij));
  LibmeshPetscCall(VecCopy(loc_vecs[2], sol_Wij));
  V("Solution from coupled solver copied to solution vectors");
 
  // Apply independent post-solve relaxation to the fixed-point solution update.
  auto relaxSolution = [&](Vec solution, Vec old_solution, const Real relaxation)
  {
    LibmeshPetscCall(VecScale(solution, relaxation));
    LibmeshPetscCall(VecAXPY(solution, 1.0 - relaxation, old_solution));
  };
  LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
      _prod, *_mdot_soln, first_node, last_node, _n_channels));
  LibmeshPetscCall(populateVectorFromHandle<SolutionHandle>(
      _prodp, *_P_soln, first_node - 1, last_node - 1, _n_channels));
 
  // Measure the pressure update before post-solve relaxation so changing that relaxation does
  // not change the meaning of the outer pressure tolerance.
  Vec pressure_residual;
  LibmeshPetscCall(VecDuplicate(sol_p, &pressure_residual));
  LibmeshPetscCall(VecCopy(sol_p, pressure_residual));
  LibmeshPetscCall(VecAXPY(pressure_residual, -1.0, _prodp));
 
  const PetscScalar * residual_array;
  const PetscScalar * old_pressure_array;
  PetscInt pressure_size;
  LibmeshPetscCall(VecGetSize(pressure_residual, &pressure_size));
  LibmeshPetscCall(VecGetArrayRead(pressure_residual, &residual_array));
  LibmeshPetscCall(VecGetArrayRead(_prodp, &old_pressure_array));
  Real residual_norm_sq = 0.0;
  Real pressure_norm_sq = 0.0;
  for (PetscInt i = 0; i < pressure_size; ++i)
  {
    residual_norm_sq += Utility::pow<2>(residual_array[i]);
    pressure_norm_sq += Utility::pow<2>(old_pressure_array[i] + _P_out);
  }
  LibmeshPetscCall(VecRestoreArrayRead(pressure_residual, &residual_array));
  LibmeshPetscCall(VecRestoreArrayRead(_prodp, &old_pressure_array));
  LibmeshPetscCall(VecDestroy(&pressure_residual));
  _pressure_fixed_point_error =
      std::max(_pressure_fixed_point_error,
               std::sqrt(residual_norm_sq) / (std::sqrt(pressure_norm_sq) + 1e-14));
 
  relaxSolution(sol_mdot, _prod, _mass_flow_relaxation);
  relaxSolution(sol_p, _prodp, _pressure_relaxation);
  LibmeshPetscCall(populateVectorFromDense<libMesh::DenseMatrix<Real>>(
      _Wij_vec, _Wij, first_node, last_node, _n_gaps));
  relaxSolution(sol_Wij, _Wij_vec, _crossflow_relaxation);
 
  // mass flow
  LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
      sol_mdot, *_mdot_soln, first_node, last_node, _n_channels));
 
  // pressure
  {
    PetscScalar * sol_p_array;
    LibmeshPetscCall(VecGetArray(sol_p, &sol_p_array));
    for (unsigned int iz = last_node; iz > first_node - 1; iz--)
    {
      const auto iz_ind = iz - first_node;
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
        PetscScalar value = sol_p_array[iz_ind * _n_channels + i_ch];
        _P_soln->set(node_in, value);
      }
    }
    LibmeshPetscCall(VecRestoreArray(sol_p, &sol_p_array));
  }
 
  // crossflow dense + sumWij + correction factor
  LibmeshPetscCall(populateDenseFromVector<libMesh::DenseMatrix<Real>>(
      sol_Wij, _Wij, first_node, last_node, _n_gaps));
 
  LibmeshPetscCall(MatMult(_mc_sumWij_mat, sol_Wij, _prod));
  LibmeshPetscCall(populateSolutionChan<SolutionHandle>(
      _prod, *_SumWij_soln, first_node, last_node, _n_channels));
 
  LibmeshPetscCall(VecAbs(_prod));
  LibmeshPetscCall(VecMax(_prod, NULL, &_max_sumWij_new));
  V("Maximum estimated Wij new: " + std::to_string(_max_sumWij_new));
  _correction_factor = _max_sumWij_new / _max_sumWij;
  V("Correction factor: " + std::to_string(_correction_factor));
  V("Solutions assigned to MOOSE variables.");
 
  // cleanup solution objects
  LibmeshPetscCall(VecDestroy(&x_nest));
  LibmeshPetscCall(VecDestroy(&sol_mdot));
  LibmeshPetscCall(VecDestroy(&sol_p));
  LibmeshPetscCall(VecDestroy(&sol_Wij));
  V("Solutions destroyed.");
 
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}
 
void
SubChannel1PhaseProblem::externalSolve()
{
  _console << "Executing subchannel solver\n";
  _dt = (isTransient() ? dt() : _one);
  _TR = isTransient();
 
  // The subchannel solver hardcodes a first-order backward (implicit) Euler time discretization, so
  // any other time integrator a user selects is silently ignored. Warn once if one is requested.
  if (!_time_integrator_checked)
  {
    _time_integrator_checked = true;
    if (isTransient())
      if (auto * transient = dynamic_cast<TransientBase *>(_app.getExecutioner()))
        for (const auto * ti : transient->getTimeIntegrators())
          if (!dynamic_cast<const ImplicitEuler *>(ti))
            mooseWarning("The subchannel solver always uses implicit (backward) Euler time "
                         "integration; the requested '",
                         ti->type(),
                         "' time integrator is ignored.");
  }
 
  initializeSolution();
  computeBulkReynoldsNumber();
  // Small helper functions to reduce repetition
  // Verbose print helper (no-op unless _verbose_subchannel is true)
  auto V = [&](const std::string & s)
  {
    if (_verbose_subchannel)
      _console << s << std::endl;
  };
  auto saveValues = [&](const SolutionHandle & solution,
                        const unsigned int first_node,
                        const unsigned int last_node)
  {
    std::vector<Real> values;
    values.reserve((last_node - first_node + 1) * _n_channels);
    for (unsigned int iz = first_node; iz <= last_node; ++iz)
      for (unsigned int i_ch = 0; i_ch < _n_channels; ++i_ch)
        values.push_back(solution(_subchannel_mesh.getChannelNode(i_ch, iz)));
    return values;
  };
  auto relativeChange = [&](const SolutionHandle & solution,
                            const std::vector<Real> & old_values,
                            const unsigned int first_node,
                            const unsigned int last_node,
                            const Real reference_offset)
  {
    Real difference_norm_sq = 0.0;
    Real reference_norm_sq = 0.0;
    std::size_t i = 0;
    for (unsigned int iz = first_node; iz <= last_node; ++iz)
      for (unsigned int i_ch = 0; i_ch < _n_channels; ++i_ch)
      {
        const Real value = solution(_subchannel_mesh.getChannelNode(i_ch, iz));
        difference_norm_sq += Utility::pow<2>(value - old_values[i]);
        reference_norm_sq += Utility::pow<2>(old_values[i] + reference_offset);
        ++i;
      }
    return std::sqrt(difference_norm_sq) / (std::sqrt(reference_norm_sq) + 1e-14);
  };
  V("Solution initialized");
  Real P_error = 1.0;
  unsigned int P_it = 0;
  unsigned int P_it_max;
  bool temperature_converged = true;
 
  if (_segregated_bool)
    P_it_max = 20 * _n_blocks;
  else
    P_it_max = 100;
 
  if ((_n_blocks == 1) && (_segregated_bool))
    P_it_max = 5;
 
  if (_P_maxit > 0)
    P_it_max = _P_maxit;
 
  while ((P_error > _P_tol && P_it < P_it_max))
  {
    P_it += 1;
    temperature_converged = true;
    _pressure_fixed_point_error = 0.0;
    _console << "Solving Outer Iteration : " << P_it << std::endl;
    const auto P_old = saveValues(*_P_soln, 0, _n_cells);
    for (unsigned int iblock = 0; iblock < _n_blocks; iblock++)
    {
      int last_level = (iblock + 1) * _block_size;
      int first_level = iblock * _block_size + 1;
      Real T_block_error = 1.0;
      auto T_it = 0;
      _console << "Solving Block: " << iblock << " From first level: " << first_level
               << " to last level: " << last_level << std::endl;
 
      while (T_block_error > _T_tol && T_it < _T_maxit)
      {
        if (processor_id() == 0)
        {
          if (_segregated_bool)
            computeWijFromSolve(iblock);
          else
          {
            LibmeshPetscCall(implicitPetscSolve(iblock));
            computeWijPrime(iblock);
            V("Done with main solve.");
          }
        }
 
        // A value of one preserves the original ordering by refreshing the flow solution before
        // every thermal update. Larger values opt into lagging flow during thermal subcycles.
        for (const auto enthalpy_subcycle : make_range(_enthalpy_subcycles))
        {
          if (T_block_error <= _T_tol || T_it >= _T_maxit)
            break;
 
          V("Enthalpy subcycle: " + std::to_string(enthalpy_subcycle + 1));
          T_it += 1;
          const auto T_old = saveValues(*_T_soln, first_level, last_level);
          // We are only computing quantities on rank 0
          if (processor_id() > 0)
          {
            // Rank zero supplies the error through the maximum reduction below.
            T_block_error = 0.0;
            goto aux_close;
          }
 
          if (_compute_power)
          {
            computeh(iblock);
            T_block_error = computeT(iblock);
            V("Done with thermal solve.");
          }
 
          V("Start updating thermophysical properties.");
          if (_compute_density)
            computeRho(iblock);
          if (_compute_viscosity)
            computeMu(iblock);
          V("Done updating thermophysical properties.");
 
          // We must do a global assembly to make sure data is parallel consistent before we do
          // things like compute L2 norms
        aux_close:
          _aux->solution().close();
 
          if (!_compute_power)
            T_block_error =
                relativeChange(*_T_soln, T_old, first_level, last_level, /*reference_offset=*/0.0);
          _console << "T_block_error: " << T_block_error << std::endl;
 
          // All processes must have the same iteration count
          comm().max(T_block_error);
        }
      }
      const bool block_converged = T_block_error <= _T_tol;
      temperature_converged &= block_converged;
      if (!block_converged)
      {
        _console << "Reached maximum number of temperature iterations for block: " << iblock
                 << std::endl;
      }
    }
    P_error = _segregated_bool ? relativeChange(*_P_soln, P_old, 0, _n_cells, _P_out)
                               : _pressure_fixed_point_error;
    comm().max(P_error);
    _console << "P_error :" << P_error << std::endl;
    V("Iteration:  " + std::to_string(P_it));
    V("Maximum iterations: " + std::to_string(P_it_max));
  }
  // Cache only the final iteration status. Earlier outer iterations may fail their thermal
  // tolerance and subsequently recover.
  const bool pressure_converged = P_error <= _P_tol;
  if (!pressure_converged)
  {
    _console << "Reached maximum number of axial pressure iterations" << std::endl;
  }
  _converged = pressure_converged && temperature_converged;
  // update old crossflow matrix
  _Wij_old = _Wij;
  _console << "Finished executing subchannel solver\n";
 
  // set SumWij at the inlet equal to the one on the first axial level  (for visualization purposes)
  for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
  {
    auto * node_in = _subchannel_mesh.getChannelNode(i_ch, 0);
    auto * node_out = _subchannel_mesh.getChannelNode(i_ch, 1);
    _SumWij_soln->set(node_in, (*_SumWij_soln)(node_out)); // kg/sec
  }
 
  if (_pin_mesh_exist)
  {
    // Assign average HTC to subchannels. This is exact if all pins have the same diameter
    for (unsigned int iz = 0; iz < _n_cells + 1; ++iz)
    {
      for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
      {
        const auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
        auto mu = (*_mu_soln)(node);
        auto S = (*_S_flow_soln)(node);
        auto w_perim = (*_w_perim_soln)(node);
        auto Dh_i = 4.0 * S / w_perim;
        auto Re = ((std::abs((*_mdot_soln)(node)) / S) * Dh_i / mu);
        auto k = _fp->k_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
        auto cp = _fp->cp_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
        auto Pr = (*_mu_soln)(node)*cp / k;
        // Create Friction structure
        _friction_args = FrictionStruct(i_ch, Re, S, w_perim);
 
        Real sumhw = 0.0;
        for (auto i_pin : _subchannel_mesh.getChannelPins(i_ch))
        {
          // Create nusselt number structure
          _nusselt_args = NusseltStruct(Re, Pr, i_pin, iz, i_ch);
 
          // Compute HTC
          sumhw += _pin_HTC_closure->computeHTC(_friction_args, _nusselt_args, k);
        }
 
        // Set HTC
        _HTC_soln->set(node, sumhw / _subchannel_mesh.getChannelPins(i_ch).size());
      }
    }
    _HTC_soln->close();
 
    _console << "Commencing calculation of Pin surface temperature \n";
    for (unsigned int i_pin = 0; i_pin < _n_pins; i_pin++)
    {
      for (unsigned int iz = 0; iz < _n_cells + 1; ++iz)
      {
        const auto * pin_node = _subchannel_mesh.getPinNode(i_pin, iz);
        Real sumTemp = 0.0;
        // Calculate sum of pin surface temperatures that the channels around the pin see
        for (auto i_ch : _subchannel_mesh.getPinChannels(i_pin))
        {
          const auto * node = _subchannel_mesh.getChannelNode(i_ch, iz);
          auto mu = (*_mu_soln)(node);
          auto S = (*_S_flow_soln)(node);
          auto w_perim = (*_w_perim_soln)(node);
          auto Dh_i = 4.0 * S / w_perim;
          auto Re = ((std::abs((*_mdot_soln)(node)) / S) * Dh_i / mu);
          auto k = _fp->k_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
          auto cp = _fp->cp_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
          auto Pr = (*_mu_soln)(node)*cp / k;
          // Create Friction structure
          _friction_args = FrictionStruct(i_ch, Re, S, w_perim);
          // Create nusselt number structure
          _nusselt_args = NusseltStruct(Re, Pr, i_pin, iz, i_ch);
          // Compute HTC
          auto hw = _pin_HTC_closure->computeHTC(_friction_args, _nusselt_args, k);
          // Compute surface temperature contribution from subchannel side
          sumTemp +=
              (*_q_prime_soln)(pin_node) / ((*_Dpin_soln)(pin_node)*M_PI * hw) + (*_T_soln)(node);
        }
        if (_subchannel_mesh.getPinChannels(i_pin).size() > 0)
          _Tpin_soln->set(pin_node, sumTemp / _subchannel_mesh.getPinChannels(i_pin).size());
        else
          mooseError("Pin was not found for pin index:  " + std::to_string(i_pin));
      }
    }
  }
 
  if (_duct_mesh_exist && processor_id() == 0)
  {
    _console << "Commencing calculation of duct surface temperature " << std::endl;
    auto duct_nodes = _subchannel_mesh.getDuctNodes();
    for (Node * dn : duct_nodes)
    {
      auto * node_chan = _subchannel_mesh.getChannelNodeFromDuct(dn);
      auto mu = (*_mu_soln)(node_chan);
      auto S = (*_S_flow_soln)(node_chan);
      auto w_perim = (*_w_perim_soln)(node_chan);
      auto Dh_i = 4.0 * S / w_perim;
      auto Re = ((std::abs((*_mdot_soln)(node_chan)) / S) * Dh_i / mu);
      auto k = _fp->k_from_p_T((*_P_soln)(node_chan) + _P_out, (*_T_soln)(node_chan));
      auto cp = _fp->cp_from_p_T((*_P_soln)(node_chan) + _P_out, (*_T_soln)(node_chan));
      auto Pr = (*_mu_soln)(node_chan)*cp / k;
 
      // Create nusselt number structure (consistent with pin case)
      const libMesh::Point & node_point = *_subchannel_mesh.getChannelNodeFromDuct(dn);
      const unsigned int iz = _subchannel_mesh.getZIndex(node_point);
      const unsigned int i_ch = _subchannel_mesh.channelIndex(node_point);
 
      // Create nusselt number structure
      _nusselt_args = NusseltStruct(Re, Pr, std::numeric_limits<unsigned int>::max(), iz, i_ch);
 
      // Create Friction structure
      _friction_args = FrictionStruct(i_ch, Re, S, w_perim);
 
      // Compute HTC
      auto hw = _duct_HTC_closure->computeHTC(_friction_args, _nusselt_args, k);
 
      // Compute Channel Temperature
      auto T_chan = (*_duct_heat_flux_soln)(dn) / hw + (*_T_soln)(node_chan);
      _Tduct_soln->set(dn, T_chan);
    }
  }
  _aux->solution().close();
  _aux->update();
 
  // Called before the rank-0 early return since it broadcasts to all ranks
  computeBulkReynoldsNumber();
 
  if (processor_id() != 0)
    return;
  Real power_in = 0.0;
  Real power_out = 0.0;
  Real mass_flow_in = 0.0;
  Real mass_flow_out = 0.0;
  for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
  {
    auto * node_in = _subchannel_mesh.getChannelNode(i_ch, 0);
    auto * node_out = _subchannel_mesh.getChannelNode(i_ch, _n_cells);
    const Real mdot_in = (*_mdot_soln)(node_in);
    power_in += mdot_in * (*_h_soln)(node_in);
    power_out += (*_mdot_soln)(node_out) * (*_h_soln)(node_out);
    mass_flow_in += mdot_in;
    mass_flow_out += (*_mdot_soln)(node_out);
  }
  auto h_bulk_out = power_out / mass_flow_out;
  auto T_bulk_out = _fp->T_from_p_h(_P_out, h_bulk_out);
 
  const Real bulk_Dh = _subchannel_mesh.getAssemblyHydraulicDiameter();
  if (_verbose_subchannel)
  {
    _console << " ======================================= " << std::endl;
    _console << " ======== Subchannel Print Outs ======== " << std::endl;
    _console << " ======================================= " << std::endl;
    _console << "Total flow area :" << _subchannel_mesh.getAssemblyFlowArea() << " m^2"
             << std::endl;
    _console << "Assembly hydraulic diameter :" << bulk_Dh << " m" << std::endl;
    _console << "Assembly Re number :" << _bulk_Re << " [-]" << std::endl;
    _console << "Bulk coolant temperature at outlet :" << T_bulk_out << " K" << std::endl;
    _console << "Power added to coolant is : " << power_out - power_in << " Watt" << std::endl;
    _console << "Mass flow rate in is : " << mass_flow_in << " kg/sec" << std::endl;
    _console << "Mass balance is : " << mass_flow_out - mass_flow_in << " kg/sec" << std::endl;
    _console << "User defined outlet pressure is : " << _P_out << " Pa" << std::endl;
    _console << " ======================================= " << std::endl;
  }
 
  if (MooseUtils::absoluteFuzzyLessEqual((power_out - power_in), -1.0))
    mooseWarning(
        "Energy conservation equation might not be solved correctly, Power added to coolant:  " +
        std::to_string(power_out - power_in) + " Watt ");
}
 
void
SubChannel1PhaseProblem::syncSolutions(Direction /*direction*/)
{
}