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 "SinglePhaseFluidProperties.h"
#include <petscdm.h>
#include <petscdmda.h>
#include <petscksp.h>
#include <petscsys.h>
#include <petscvec.h>
#include <petscsnes.h>

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;

  virtual Real computeAddedHeatPin(unsigned int i_ch, unsigned int iz) = 0;
  Real computeAddedHeatDuct(unsigned int i_ch, unsigned int iz);

protected:
  struct FrictionStruct
  {
    int i_ch;
    Real Re, S, w_perim;
  } _friction_args;

  virtual Real computeFrictionFactor(FrictionStruct friction_args) = 0;
  void computeWijFromSolve(int iblock);
  void computeSumWij(int iblock);
  void computeMdot(int iblock);
  void computeWijPrime(int iblock);
  virtual Real computeBeta(unsigned int i_gap, unsigned int iz, bool enthalpy) = 0;
  void computeDP(int iblock);
  void computeP(int iblock);
  virtual void computeh(int iblock) = 0;
  void computeT(int iblock);
  void computeRho(int iblock);
  void computeMu(int iblock);
  void computeWijResidual(int iblock);
  virtual Real getSubChannelPeripheralDuctWidth(unsigned int i_ch) = 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;

  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 cleanUp();
  SubChannelMesh & _subchannel_mesh;
  unsigned int _n_blocks;
  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 _dt;
  const PostprocessorValue & _P_out;
  const Real & _CT;
  const Real & _P_tol;
  const Real & _T_tol;
  const int & _T_maxit;
  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 _monolithic_thermal_bool;
  const bool _verbose_subchannel;
  const bool _deformation;

  const SinglePhaseFluidProperties * _fp;
  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;

  inline PetscErrorCode createPetscVector(Vec & v, PetscInt n)
  {
    PetscFunctionBegin;
    LibmeshPetscCall(VecCreate(PETSC_COMM_WORLD, &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_WORLD, &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 _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;

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