SubChannel1PhaseProblem.h

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
 
#pragma once
 
#include "ExternalProblem.h"
#include "PostprocessorInterface.h"
#include "SubChannelApp.h"
#include "QuadSubChannelMesh.h"
#include "SolutionHandle.h"
#include <petscdm.h>
#include <petscdmda.h>
#include <petscksp.h>
#include <petscsys.h>
#include <petscvec.h>
#include <petscsnes.h>
#include <limits>
 
class SinglePhaseFluidProperties;
class SCMFrictionClosureBase;
class SCMHTCClosureBase;
class SCMMixingClosureBase;
 
class SubChannel1PhaseProblem : public ExternalProblem, public PostprocessorInterface
{
public:
  SubChannel1PhaseProblem(const InputParameters & params);
  virtual ~SubChannel1PhaseProblem();
 
  virtual void externalSolve() override;
  virtual void syncSolutions(Direction direction) override;
  virtual bool solverSystemConverged(const unsigned int) override;
  virtual void initialSetup() override;
 
  const SCMHTCClosureBase * getDuctHTCClosure() const { return _duct_HTC_closure; }
  const SCMHTCClosureBase * getPinHTCClosure() const { return _pin_HTC_closure; } // optional
  const SCMFrictionClosureBase * getFrictionClosure() const { return _friction_closure; }
  Real getBulkReynoldsNumber() const { return _bulk_Re; }
 
  struct FrictionStruct
  {
    unsigned int i_ch = 0;
    Real Re = 1.0;
    Real S = 0.0;
    Real w_perim = 0.0;
 
    FrictionStruct() = delete;
    FrictionStruct(unsigned int i_ch_, Real Re_, Real S_, Real w_perim_)
      : i_ch(i_ch_), Re(Re_), S(S_), w_perim(w_perim_)
    {
    }
  } _friction_args;
 
  struct NusseltStruct
  {
    Real Re = 1.0;
    Real Pr = 1.0;
    unsigned int i_pin = std::numeric_limits<unsigned int>::max(); // sentinel (duct) default
    unsigned int iz = 0;
    unsigned int i_ch = 0;
 
    NusseltStruct() = delete;
    NusseltStruct(Real Re_, Real Pr_, unsigned int i_pin_, unsigned int iz_, unsigned int i_ch_)
      : Re(Re_), Pr(Pr_), i_pin(i_pin_), iz(iz_), i_ch(i_ch_)
    {
    }
  } _nusselt_args;
 
  Real getAddedHeatPin(unsigned int i_ch, unsigned int iz) const
  {
    return computeAddedHeatPin(i_ch, iz);
  }
 
  Real getAddedHeatDuct(unsigned int i_ch, unsigned int iz) const
  {
    return computeAddedHeatDuct(i_ch, iz);
  }
 
  const PostprocessorValue & _P_out;
 
  const PostprocessorValue & getOutletPressure() const { return _P_out; }
 
  const SinglePhaseFluidProperties * _fp;
 
  const SinglePhaseFluidProperties * getSinglePhaseFluidProperties() const { return _fp; }
 
protected:
  virtual Real computeAddedHeatPin(unsigned int i_ch, unsigned int iz) const = 0;
 
  virtual Real computeAddedHeatDuct(unsigned int i_ch, unsigned int iz) const;
 
  void computeWijFromSolve(int iblock);
  void computeSumWij(int iblock);
  void computeMdot(int iblock);
  void computeWijPrime(int iblock);
  Real computeMixingParameter(unsigned int i_gap, unsigned int iz) const;
  Real computeSweepFlowMixingParameter(unsigned int i_gap, unsigned int iz) const;
  void computeDP(int iblock);
  void computeP(int iblock);
  virtual void computeh(int iblock) = 0;
  Real computeT(int iblock);
  void computeRho(int iblock);
  void computeMu(int iblock);
  void computeWijResidual(int iblock);
  virtual Real getSubChannelPeripheralDuctWidth(unsigned int i_ch) const = 0;
  libMesh::DenseVector<Real> residualFunction(int iblock, libMesh::DenseVector<Real> solution);
  PetscErrorCode petscSnesSolver(int iblock,
                                 const libMesh::DenseVector<Real> & solution,
                                 libMesh::DenseVector<Real> & root);
  friend PetscErrorCode formFunction(SNES snes, Vec x, Vec f, void * ctx);
 
  PetscErrorCode implicitPetscSolve(int iblock);
 
  virtual void initializeSolution() = 0;
  void detectDeformation();
  void computeBulkReynoldsNumber();
 
  PetscScalar computeInterpolationCoefficients(PetscScalar Peclet = 0.0);
  PetscScalar
  computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet = 0.0);
 
  Real computeGravityDir(const MooseEnum & dir) const
  {
    switch (dir)
    {
      case 0: // counter_flow
        return 1.0;
      case 1: // co_flow
        return -1.0;
      case 2: // none
        return 0.0;
      default:
        mooseError(name(), ": Invalid gravity direction: expected counter_flow, co_flow, or none");
    }
  }
 
  PetscErrorCode solveAndPopulateEnthalpy(
      Mat A, Vec rhs, unsigned int first_node, unsigned int last_node, const char * ksp_prefix);
 
  PetscErrorCode cleanUp();
  SubChannelMesh & _subchannel_mesh;
  unsigned int _n_blocks;
  libMesh::DenseMatrix<Real> _DP;
  libMesh::DenseMatrix<Real> & _Wij;
  libMesh::DenseMatrix<Real> & _Wij_old;
  libMesh::DenseMatrix<Real> _WijPrime;
  libMesh::DenseMatrix<Real> _Wij_residual_matrix;
  const Real _g_grav;
  const Real & _kij;
  unsigned int _n_cells;
  unsigned int _n_gaps;
  unsigned int _n_pins;
  unsigned int _n_channels;
  unsigned int _block_size;
  std::vector<Real> _z_grid;
  Real _one;
  Real _TR;
  const bool _compute_density;
  const bool _compute_viscosity;
  const bool _compute_power;
  const bool _pin_mesh_exist;
  const bool _duct_mesh_exist;
  bool _converged;
  Real _bulk_Re;
  bool _time_integrator_checked = false;
  Real _dt;
  Real _CT;
  const Real & _P_tol;
  const int & _P_maxit;
  const Real & _T_tol;
  const int & _T_maxit;
  const Real & _T_relaxation;
  const unsigned int & _enthalpy_subcycles;
  const Real & _mass_flow_equation_relaxation;
  const Real & _pressure_equation_relaxation;
  const Real & _crossflow_equation_relaxation;
  const Real & _mass_flow_relaxation;
  const Real & _pressure_relaxation;
  const Real & _crossflow_relaxation;
  const PetscReal & _rtol;
  const PetscReal & _atol;
  const PetscReal & _dtol;
  const PetscInt & _maxit;
  const MooseEnum _interpolation_scheme;
  const MooseEnum _gravity_direction;
  const Real _dir_grav;
  const bool _implicit_bool;
  const bool _staggered_pressure_bool;
  const bool _segregated_bool;
  const bool _verbose_subchannel;
  bool _deformation = false;
  const SCMFrictionClosureBase * _friction_closure;
  const SCMMixingClosureBase * _mixing_closure;
  const SCMHTCClosureBase * _pin_HTC_closure;
  const SCMHTCClosureBase * _duct_HTC_closure;
 
  std::unique_ptr<SolutionHandle> _mdot_soln;
  std::unique_ptr<SolutionHandle> _SumWij_soln;
  std::unique_ptr<SolutionHandle> _P_soln;
  std::unique_ptr<SolutionHandle> _DP_soln;
  std::unique_ptr<SolutionHandle> _h_soln;
  std::unique_ptr<SolutionHandle> _T_soln;
  std::unique_ptr<SolutionHandle> _Tpin_soln;
  std::unique_ptr<SolutionHandle> _Dpin_soln;
  std::unique_ptr<SolutionHandle> _rho_soln;
  std::unique_ptr<SolutionHandle> _mu_soln;
  std::unique_ptr<SolutionHandle> _S_flow_soln;
  std::unique_ptr<SolutionHandle> _w_perim_soln;
  std::unique_ptr<SolutionHandle> _q_prime_soln;
  std::unique_ptr<SolutionHandle> _duct_heat_flux_soln; // Only used for ducted assemblies
  std::unique_ptr<SolutionHandle> _Tduct_soln;          // Only used for ducted assemblies
  std::unique_ptr<SolutionHandle> _displacement_soln;
  std::unique_ptr<SolutionHandle> _ff_soln;
  std::unique_ptr<SolutionHandle> _HTC_soln;
 
  inline PetscErrorCode createPetscVector(Vec & v, PetscInt n)
  {
    PetscFunctionBegin;
    LibmeshPetscCall(VecCreate(PETSC_COMM_SELF, &v));
    LibmeshPetscCall(PetscObjectSetName((PetscObject)v, "Solution"));
    LibmeshPetscCall(VecSetSizes(v, PETSC_DECIDE, n));
    LibmeshPetscCall(VecSetFromOptions(v));
    LibmeshPetscCall(VecZeroEntries(v));
    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }
 
  inline PetscErrorCode createPetscMatrix(Mat & M, PetscInt n, PetscInt m)
  {
    PetscFunctionBegin;
    LibmeshPetscCall(MatCreate(PETSC_COMM_SELF, &M));
    LibmeshPetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, n, m));
    LibmeshPetscCall(MatSetFromOptions(M));
    LibmeshPetscCall(MatSetUp(M));
    PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
  }
 
  template <class T>
  PetscErrorCode populateVectorFromDense(Vec & x,
                                         const T & solution,
                                         const unsigned int first_axial_level,
                                         const unsigned int last_axial_level,
                                         const unsigned int cross_dimension);
  template <class T>
  PetscErrorCode populateDenseFromVector(const Vec & x,
                                         T & solution,
                                         const unsigned int first_axial_level,
                                         const unsigned int last_axial_level,
                                         const unsigned int cross_dimension);
  template <class T>
  PetscErrorCode populateVectorFromHandle(Vec & x,
                                          const T & solution,
                                          const unsigned int first_axial_level,
                                          const unsigned int last_axial_level,
                                          const unsigned int cross_dimension);
  template <class T>
  PetscErrorCode populateSolutionChan(const Vec & x,
                                      T & solution,
                                      const unsigned int first_axial_level,
                                      const unsigned int last_axial_level,
                                      const unsigned int cross_dimension);
 
  Mat _mc_sumWij_mat;
  Vec _Wij_vec;
  Vec _prod;
  Vec _prodp;
  Mat _mc_axial_convection_mat;
  Vec _mc_axial_convection_rhs;
  Mat _mc_density_pressure_mat;
 
  Mat _amc_turbulent_cross_flows_mat;
  Vec _amc_turbulent_cross_flows_rhs;
  Mat _amc_time_derivative_mat;
  Vec _amc_time_derivative_rhs;
  Mat _amc_advective_derivative_mat;
  Vec _amc_advective_derivative_rhs;
  Mat _amc_cross_derivative_mat;
  Vec _amc_cross_derivative_rhs;
  Mat _amc_friction_force_mat;
  Vec _amc_friction_force_rhs;
  Vec _amc_gravity_rhs;
  Mat _amc_pressure_force_mat;
  Vec _amc_pressure_force_rhs;
  Mat _amc_sys_mdot_mat;
  Vec _amc_sys_mdot_rhs;
 
  Mat _cmc_time_derivative_mat;
  Vec _cmc_time_derivative_rhs;
  Mat _cmc_advective_derivative_mat;
  Vec _cmc_advective_derivative_rhs;
  Mat _cmc_friction_force_mat;
  Vec _cmc_friction_force_rhs;
  Mat _cmc_pressure_force_mat;
  Vec _cmc_pressure_force_rhs;
  Mat _cmc_sys_Wij_mat;
  Vec _cmc_sys_Wij_rhs;
 
  Mat _hc_time_derivative_mat;
  Vec _hc_time_derivative_rhs;
  Mat _hc_advective_derivative_mat;
  Vec _hc_advective_derivative_rhs;
  Mat _hc_cross_derivative_mat;
  Vec _hc_cross_derivative_rhs;
  Vec _hc_added_heat_rhs;
  Mat _hc_sys_h_mat;
  Vec _hc_sys_h_rhs;
 
  PetscScalar _added_K = 0.0;
  PetscScalar _added_K_old = 1000.0;
  PetscScalar _max_sumWij;
  PetscScalar _max_sumWij_new;
  PetscScalar _correction_factor = 1.0;
  Real _pressure_fixed_point_error = 1.0;
 
public:
  static InputParameters validParams();
};
 
template <class T>
PetscErrorCode
SubChannel1PhaseProblem::populateDenseFromVector(const Vec & x,
                                                 T & loc_solution,
                                                 const unsigned int first_axial_level,
                                                 const unsigned int last_axial_level,
                                                 const unsigned int cross_dimension)
{
  PetscScalar * xx;
 
  PetscFunctionBegin;
  LibmeshPetscCall(VecGetArray(x, &xx));
  for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
  {
    unsigned int iz_ind = iz - first_axial_level;
    for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
    {
      loc_solution(i_l, iz) = xx[iz_ind * cross_dimension + i_l];
    }
  }
  LibmeshPetscCall(VecRestoreArray(x, &xx));
 
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}
 
template <class T>
PetscErrorCode
SubChannel1PhaseProblem::populateVectorFromHandle(Vec & x,
                                                  const T & loc_solution,
                                                  const unsigned int first_axial_level,
                                                  const unsigned int last_axial_level,
                                                  const unsigned int cross_dimension)
{
  PetscScalar * xx;
 
  PetscFunctionBegin;
  LibmeshPetscCall(VecGetArray(x, &xx));
  for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
  {
    unsigned int iz_ind = iz - first_axial_level;
    for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
    {
      auto * loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
      xx[iz_ind * cross_dimension + i_l] = loc_solution(loc_node);
    }
  }
  LibmeshPetscCall(VecRestoreArray(x, &xx));
 
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}
 
template <class T>
PetscErrorCode
SubChannel1PhaseProblem::populateVectorFromDense(Vec & x,
                                                 const T & loc_solution,
                                                 const unsigned int first_axial_level,
                                                 const unsigned int last_axial_level,
                                                 const unsigned int cross_dimension)
{
  PetscScalar * xx;
  PetscFunctionBegin;
  LibmeshPetscCall(VecGetArray(x, &xx));
  for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
  {
    unsigned int iz_ind = iz - first_axial_level;
    for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
    {
      xx[iz_ind * cross_dimension + i_l] = loc_solution(i_l, iz);
    }
  }
  LibmeshPetscCall(VecRestoreArray(x, &xx));
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}
 
template <class T>
PetscErrorCode
SubChannel1PhaseProblem::populateSolutionChan(const Vec & x,
                                              T & loc_solution,
                                              const unsigned int first_axial_level,
                                              const unsigned int last_axial_level,
                                              const unsigned int cross_dimension)
{
  PetscScalar * xx;
  PetscFunctionBegin;
  LibmeshPetscCall(VecGetArray(x, &xx));
  Node * loc_node;
  for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
  {
    unsigned int iz_ind = iz - first_axial_level;
    for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
    {
      loc_node = _subchannel_mesh.getChannelNode(i_l, iz);
      loc_solution.set(loc_node, xx[iz_ind * cross_dimension + i_l]);
    }
  }
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}