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.addRequiredParam<unsigned int>("n_blocks", "The number of blocks in the axial direction");
  params.addParam<Real>("P_tol", 1e-6, "Pressure tolerance");
  params.addParam<Real>("T_tol", 1e-6, "Temperature tolerance");
  params.addParam<int>("T_maxit", 100, "Maximum number of iterations for inner temperature loop");
  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 for the method. 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.addParam<bool>("staggered_pressure",
                        false,
                        "Boolean to define the use of staggered or collocated pressure.");
  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 T_tol T_maxit rtol atol dtol maxit",
                              "Solver tolerances and iterations");
  params.addParamNamesToGroup("implicit segregated staggered_pressure 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")),
    _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()),
    _P_tol(getParam<Real>("P_tol")),
    _T_tol(getParam<Real>("T_tol")),
    _T_maxit(getParam<int>("T_maxit")),
    _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")),
    _staggered_pressure_bool(getParam<bool>("staggered_pressure")),
    _segregated_bool(getParam<bool>("segregated")),
    _verbose_subchannel(getParam<bool>("verbose_subchannel")),
    _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.");
  // 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));

  // 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()
{
  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();

  // Create variables for output and storage
  _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));
    _pin_HTC_closure =
        &getUserObject<SCMHTCClosureBase>(getParam<UserObjectName>("pin_HTC_closure"));
  }
  _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));
    _duct_HTC_closure =
        &getUserObject<SCMHTCClosureBase>(getParam<UserObjectName>("duct_HTC_closure"));
  }
}

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));

  // 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)
      return ((Peclet - 1.0) * std::exp(Peclet) + 1) / (Peclet * (std::exp(Peclet) - 1.) + 1e-10);
    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
  {
    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));

        // 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));

    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);
        // hydraulic diameter in the i direction
        auto Dh_i = 4.0 * S / w_perim;
        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 Re = (((*_mdot_soln)(node_in) / S) * Dh_i / mu_in);
        _friction_args = FrictionStruct(i_ch, Re, S, w_perim);
        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 friction_term = (ff * dz / Dh_i + ki) * 0.5 *
                             (*_mdot_soln)(node_out)*std::abs((*_mdot_soln)(node_out)) /
                             (S * (*_rho_soln)(node_out));
        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_i = 4.0 * S_interp / w_perim_interp;
          // Compute friction factor
          auto Re = ((mdot_loc / S_interp) * Dh_i / 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_i + 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);

        // hydraulic diameter in the i direction
        auto Dh_i = 4.0 * S_interp / w_perim_interp;

        if (iz == first_node)
        {
          PetscScalar value_vec_tt = -1.0 * _TR * alpha * (*_mdot_soln)(node_in)*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));
        }
        else
        {
          PetscInt row_tt = i_ch + _n_channels * iz_ind;
          PetscInt col_tt = i_ch + _n_channels * (iz_ind - 1);
          PetscScalar value_tt = _TR * alpha * dz / _dt;
          LibmeshPetscCall(MatSetValues(
              _amc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));
        }

        // Adding diagonal elements
        PetscInt row_tt = i_ch + _n_channels * iz_ind;
        PetscInt col_tt = i_ch + _n_channels * iz_ind;
        PetscScalar value_tt = _TR * (1.0 - alpha) * dz / _dt;
        LibmeshPetscCall(MatSetValues(
            _amc_time_derivative_mat, 1, &row_tt, 1, &col_tt, &value_tt, INSERT_VALUES));

        // Adding RHS elements
        PetscScalar mdot_old_interp =
            computeInterpolatedValue(_mdot_soln->old(node_out), _mdot_soln->old(node_in), Pe);
        PetscScalar value_vec_tt = _TR * mdot_old_interp * 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 * std::abs((*_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 = std::abs((*_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 = ((mdot_interp / S_interp) * Dh_i / 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_i + ki) * 0.5 * std::abs((*_mdot_soln)(node_out)) /
                    (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)
  {
    if (!_staggered_pressure_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
    {
      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
          // Note: assuming uniform axial discretization in the curren code
          // We will need to update this later if we allow non-uniform refinements in the axial
          // direction
          PetscScalar Pe = 0.5;
          auto alpha = computeInterpolationCoefficients(Pe);
          if (iz == last_node)
          {
            _P_soln->set(node_in, (*_P_soln)(node_out) + _DP(i_ch, iz) / 2.0);
          }
          else
          {
            _P_soln->set(node_in,
                         (*_P_soln)(node_out) + (1.0 - alpha) * _DP(i_ch, iz) +
                             alpha * _DP(i_ch, iz - 1));
          }
        }
      }
    }
  }
  else
  {
    if (!_staggered_pressure_bool)
    {
      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));

          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));
          }

          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));
      }
    }
    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));

          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));

            auto dp_out = _DP(i_ch, iz);
            PetscScalar value_v = -1.0 * dp_out / 2.0 * S_interp;
            PetscInt row_v = i_ch + _n_channels * iz_ind;
            LibmeshPetscCall(
                VecSetValues(_amc_pressure_force_rhs, 1, &row_v, &value_v, 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));

            if (_segregated_bool)
            {
              auto dp_in = _DP(i_ch, iz - 1);
              auto dp_out = _DP(i_ch, iz);
              auto dp_interp = computeInterpolatedValue(dp_out, dp_in, 0.5);
              PetscScalar value_v = -1.0 * dp_interp * 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 (_verbose_subchannel)
        _console << "Block: " << iblock << " - Axial momentum pressure force matrix assembled"
                 << std::endl;

      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));
      }
    }
  }
}

void
SubChannel1PhaseProblem::computeT(int iblock)
{
  const unsigned int last_node = (iblock + 1) * _block_size;
  const unsigned int first_node = iblock * _block_size + 1;
  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, _fp->T_from_p_h((*_P_soln)(node) + _P_out, (*_h_soln)(node)));
    }
  }
}

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);

        if (!_staggered_pressure_bool)
        {
          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));
          }
        }
        else
        {
          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 * 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 = 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;
  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;
      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));
  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));

  // 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.
  // -----------------------------------------------------------------------------
  auto RelaxEquation =
      [&](Mat A_ff, Vec rhs_f, Vec like_vec, Vec work, PetscScalar alpha, auto && populate)
  {
    Vec d = nullptr;
    LibmeshPetscCall(VecDuplicate(like_vec, &d));

    // 1) A_ff: diag <- diag / alpha
    LibmeshPetscCall(MatGetDiagonal(A_ff, d));
    LibmeshPetscCall(VecScale(d, 1.0 / alpha));
    LibmeshPetscCall(MatDiagonalSet(A_ff, d, INSERT_VALUES));

    // 2) work <- x_old (caller-provided populator)
    LibmeshPetscCall(populate(work));

    // 3) rhs_f += (1 - alpha) * (diag .* work)
    LibmeshPetscCall(VecScale(d, (1.0 - alpha)));
    LibmeshPetscCall(VecPointwiseMult(work, work, d));
    LibmeshPetscCall(VecAXPY(rhs_f, 1.0, work));

    LibmeshPetscCall(VecDestroy(&d));
  };

  // 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=*/nullptr,
                   /*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))
        {
          auto chans = _subchannel_mesh.getGapChannels(i_gap);
          unsigned int i_ch_loc = chans.first;
          PetscInt row_vec = i_ch_loc + _n_channels * 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;
    LibmeshPetscCall(VecAbs(_prod));
    LibmeshPetscCall(VecMin(_prod, NULL, &min_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));

    // ---- per-equation under-relaxation ----
    const PetscScalar relaxation_factor_mdot = 1.0;
    const PetscScalar relaxation_factor_P = 1.0;
    const PetscScalar relaxation_factor_Wij = 0.1;

    V("Relax mdot: " + std::to_string(relaxation_factor_mdot));
    V("Relax P: " + std::to_string(relaxation_factor_P));
    V("Relax Wij: " + std::to_string(relaxation_factor_Wij));

    // mdot
    RelaxEquation(mat_array[Idx(0, 0)],
                  vec_array[0],
                  vec_array[0],
                  _prod,
                  relaxation_factor_mdot,
                  [&](Vec dst)
                  {
                    return populateVectorFromHandle<SolutionHandle>(
                        dst, *_mdot_soln, first_node, last_node, _n_channels);
                  });
    V("mdot relaxed");

    // pressure
    RelaxEquation(mat_array[Idx(1, 1)],
                  vec_array[1],
                  vec_array[1],
                  _prod,
                  relaxation_factor_P,
                  [&](Vec dst)
                  {
                    return populateVectorFromHandle<SolutionHandle>(
                        dst, *_P_soln, first_node, last_node, _n_channels);
                  });
    V("P relaxed");

    // crossflow
    RelaxEquation(mat_array[Idx(2, 2)],
                  vec_array[2],
                  vec_array[2],
                  _Wij_vec,
                  relaxation_factor_Wij,
                  [&](Vec dst)
                  {
                    return populateVectorFromDense<libMesh::DenseMatrix<Real>>(
                        dst, _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(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));
  }
  V("Linear solver assembled");

  // ============================== Solve =====================================
  LibmeshPetscCall(VecDuplicate(b_nest, &x_nest));
  LibmeshPetscCall(VecSet(x_nest, 0.0));
  LibmeshPetscCall(KSPSolve(ksp, b_nest, x_nest));

  // 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");

  // 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();
  // 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;
  };
  V("Solution initialized");
  Real P_error = 1.0;
  unsigned int P_it = 0;
  unsigned int P_it_max;

  if (_segregated_bool)
    P_it_max = 20 * _n_blocks;
  else
    P_it_max = 100;

  if ((_n_blocks == 1) && (_segregated_bool))
    P_it_max = 5;

  while ((P_error > _P_tol && P_it < P_it_max))
  {
    P_it += 1;
    if (P_it == P_it_max && _n_blocks != 1)
    {
      _console << "Reached maximum number of axial pressure iterations" << std::endl;
      _converged = false;
    }
    _console << "Solving Outer Iteration : " << P_it << std::endl;
    auto P_L2norm_old_axial = _P_soln->L2norm();
    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)
      {
        T_it += 1;
        if (T_it == _T_maxit)
        {
          _console << "Reached maximum number of temperature iterations for block: " << iblock
                   << std::endl;
          _converged = false;
        }
        auto T_L2norm_old_block = _T_soln->L2norm();
        // We are only computing quantities on rank 0
        if (processor_id() > 0)
          goto aux_close;

        if (_segregated_bool)
        {
          computeWijFromSolve(iblock);
          if (_compute_power)
          {
            computeh(iblock);
            computeT(iblock);
          }
        }
        else
        {
          LibmeshPetscCall(implicitPetscSolve(iblock));
          computeWijPrime(iblock);
          V("Done with main solve.");
          if (_compute_power)
          {
            computeh(iblock);
            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();

        auto T_L2norm_new = _T_soln->L2norm();
        T_block_error =
            std::abs((T_L2norm_new - T_L2norm_old_block) / (T_L2norm_old_block + 1E-14));
        _console << "T_block_error: " << T_block_error << std::endl;

        // All processes must have the same iteration count
        comm().max(T_block_error);
      }
    }
    auto P_L2norm_new_axial = _P_soln->L2norm();
    P_error =
        std::abs((P_L2norm_new_axial - P_L2norm_old_axial) / (P_L2norm_old_axial + _P_out + 1E-14));
    _console << "P_error :" << P_error << std::endl;
    V("Iteration:  " + std::to_string(P_it));
    V("Maximum iterations: " + std::to_string(P_it_max));
  }
  // 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 = (((*_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 = (((*_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 = (((*_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();

  if (processor_id() != 0)
    return;
  Real power_in = 0.0;
  Real power_out = 0.0;
  Real viscosity_in = 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);
    viscosity_in += mdot_in * (*_mu_soln)(node_in);
    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);

  Real bulk_Dh = _subchannel_mesh.getAssemblyHydraulicDiameter();
  Real inlet_mu = viscosity_in / mass_flow_in;
  Real bulk_Re = mass_flow_in * bulk_Dh / (inlet_mu * _subchannel_mesh.getAssemblyFlowArea());
  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*/)
{
}