InterWrapper1PhaseProblem.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 "SubChannelApp.h"
#include "InterWrapperMesh.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 InterWrapper1PhaseProblem;

class InterWrapper1PhaseProblem : public ExternalProblem
{
public:
  InterWrapper1PhaseProblem(const InputParameters & params);
  virtual ~InterWrapper1PhaseProblem();

  virtual void externalSolve() override;
  virtual void syncSolutions(Direction direction) override;
  virtual bool solverSystemConverged(const unsigned int) override;
  virtual void initialSetup() override;

protected:
  struct StructPetscMatVec
  {
    Mat A;
    Vec x;
  };

  virtual Real computeFrictionFactor(Real Re) = 0;
  virtual void computeWijFromSolve(int iblock);
  virtual void computeSumWij(int iblock);
  virtual void computeMdot(int iblock);
  virtual void computeWijPrime(int iblock);
  virtual void computeDP(int iblock);
  virtual void computeP(int iblock);
  virtual void computeh(int iblock);
  virtual void computeT(int iblock);
  virtual void computeRho(int iblock);
  virtual void computeMu(int iblock);
  virtual void computeWij(int iblock);
  virtual Real computeAddedHeat(unsigned int i_ch, unsigned int iz);
  virtual libMesh::DenseVector<Real> residualFunction(int iblock,
                                                      libMesh::DenseVector<Real> solution);
  virtual PetscErrorCode petscSnesSolver(int iblock,
                                         const libMesh::DenseVector<Real> & solution,
                                         libMesh::DenseVector<Real> & root);
  friend PetscErrorCode formFunctionIW(SNES snes, Vec x, Vec f, void * info);

  virtual PetscErrorCode implicitPetscSolve(int iblock);
  virtual void initializeSolution();

  virtual PetscScalar computeInterpolationCoefficients(PetscScalar Peclet = 0.0);
  virtual PetscScalar
  computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet = 0.0);
  PetscErrorCode cleanUp();
  const InterWrapperMesh & _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;
  const unsigned int _n_cells;
  const unsigned int _n_gaps;
  const unsigned int _n_assemblies;
  const unsigned int _n_channels;
  const 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;
  bool _converged;
  const Real & _dt;
  const Real & _P_out;
  const Real & _beta;
  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 bool _implicit_bool;
  const bool _staggered_pressure_bool;
  const bool _segregated_bool;
  const bool _monolithic_thermal_bool;

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

  virtual PetscErrorCode createPetscVector(Vec & v, PetscInt n);
  virtual PetscErrorCode createPetscMatrix(Mat & M, PetscInt n, PetscInt m);
  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);

  template <class T>
  PetscErrorCode populateSolutionGap(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;
  Vec _cmc_Wij_channel_dummy;

  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;
  PetscInt _global_counter = 0;

  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
InterWrapper1PhaseProblem::populateDenseFromVector(const Vec & x,
                                                   T & loc_solution,
                                                   const unsigned int first_axial_level,
                                                   const unsigned int last_axial_level,
                                                   const unsigned int cross_dimension)
{
  PetscErrorCode ierr;
  PetscScalar * xx;

  PetscFunctionBegin;
  ierr = VecGetArray(x, &xx);
  CHKERRQ(ierr);
  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];
    }
  }
  ierr = VecRestoreArray(x, &xx);
  CHKERRQ(ierr);
  PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
}