Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : #include "ExternalProblem.h"
      12             : #include "SubChannelApp.h"
      13             : #include "InterWrapperMesh.h"
      14             : #include "SolutionHandle.h"
      15             : #include "SinglePhaseFluidProperties.h"
      16             : #include <petscdm.h>
      17             : #include <petscdmda.h>
      18             : #include <petscksp.h>
      19             : #include <petscsys.h>
      20             : #include <petscvec.h>
      21             : #include <petscsnes.h>
      22             : 
      23             : class InterWrapper1PhaseProblem;
      24             : 
      25             : /**
      26             :  * Base class for the 1-phase steady-state/transient interwrapper solver.
      27             :  */
      28             : class InterWrapper1PhaseProblem : public ExternalProblem
      29             : {
      30             : public:
      31             :   InterWrapper1PhaseProblem(const InputParameters & params);
      32             :   virtual ~InterWrapper1PhaseProblem();
      33             : 
      34             :   virtual void externalSolve() override;
      35             :   virtual void syncSolutions(Direction direction) override;
      36             :   virtual bool solverSystemConverged(const unsigned int) override;
      37             :   virtual void initialSetup() override;
      38             : 
      39             : protected:
      40             :   /// Standard return structure for reusing in implicit/explicit formulations
      41             :   struct StructPetscMatVec
      42             :   {
      43             :     Mat A;
      44             :     Vec x;
      45             :   };
      46             : 
      47             :   /// Returns friction factor
      48             :   virtual Real computeFrictionFactor(Real Re) = 0;
      49             :   /// Computes diversion crossflow per gap for block with index iblock
      50             :   /// Block is a partition of the whole domain
      51             :   virtual void computeWijFromSolve(int iblock);
      52             :   /// Computes net diversion crossflow per channel for block iblock
      53             :   virtual void computeSumWij(int iblock);
      54             :   /// Computes mass flow per channel for block iblock
      55             :   virtual void computeMdot(int iblock);
      56             :   /// Computes turbulent crossflow per gap for block iblock
      57             :   virtual void computeWijPrime(int iblock);
      58             :   /// Computes Pressure Drop per channel for block iblock
      59             :   virtual void computeDP(int iblock);
      60             :   /// Computes Pressure per channel for block iblock
      61             :   virtual void computeP(int iblock);
      62             :   /// Computes Enthalpy per channel for block iblock
      63             :   virtual void computeh(int iblock);
      64             :   /// Computes Temperature per channel for block iblock
      65             :   virtual void computeT(int iblock);
      66             :   /// Computes Density per channel for block iblock
      67             :   virtual void computeRho(int iblock);
      68             :   /// Computes Viscosity per channel for block iblock
      69             :   virtual void computeMu(int iblock);
      70             :   /// Computes cross fluxes for block iblock
      71             :   virtual void computeWij(int iblock);
      72             :   /// Computes added heat for channel i_ch and cell iz
      73             :   virtual Real computeAddedHeat(unsigned int i_ch, unsigned int iz);
      74             :   /// Computes Residual per gap for block iblock
      75             :   virtual libMesh::DenseVector<Real> residualFunction(int iblock,
      76             :                                                       libMesh::DenseVector<Real> solution);
      77             :   /// Computes solution of nonlinear equation using snes
      78             :   virtual PetscErrorCode petscSnesSolver(int iblock,
      79             :                                          const libMesh::DenseVector<Real> & solution,
      80             :                                          libMesh::DenseVector<Real> & root);
      81             :   friend PetscErrorCode formFunctionIW(SNES snes, Vec x, Vec f, void * info);
      82             : 
      83             :   /// Computes implicit solve using PetSc
      84             :   virtual PetscErrorCode implicitPetscSolve(int iblock);
      85             :   /// Function to initialize the solution fields
      86             :   virtual void initializeSolution();
      87             : 
      88             :   /// Functions that computes the interpolation scheme given the Peclet number
      89             :   virtual PetscScalar computeInterpolationCoefficients(PetscScalar Peclet = 0.0);
      90             :   virtual PetscScalar
      91             :   computeInterpolatedValue(PetscScalar topValue, PetscScalar botValue, PetscScalar Peclet = 0.0);
      92             :   PetscErrorCode cleanUp();
      93             :   const InterWrapperMesh & _subchannel_mesh;
      94             :   /// number of axial blocks
      95             :   unsigned int _n_blocks;
      96             :   libMesh::DenseMatrix<Real> & _Wij;
      97             :   libMesh::DenseMatrix<Real> _Wij_old;
      98             :   libMesh::DenseMatrix<Real> _WijPrime;
      99             :   libMesh::DenseMatrix<Real> _Wij_residual_matrix;
     100             :   const Real _g_grav;
     101             :   const Real & _kij;
     102             :   const unsigned int _n_cells;
     103             :   const unsigned int _n_gaps;
     104             :   const unsigned int _n_assemblies;
     105             :   const unsigned int _n_channels;
     106             :   const unsigned int _block_size;
     107             :   /// axial location of nodes
     108             :   std::vector<Real> _z_grid;
     109             :   Real _one;
     110             :   /// Flag that activates or deactivates the transient parts of the equations we solve by multiplication
     111             :   Real _TR;
     112             :   /// Flag that activates or deactivates the calculation of density
     113             :   const bool _compute_density;
     114             :   /// Flag that activates or deactivates the calculation of viscosity
     115             :   const bool _compute_viscosity;
     116             :   /// Flag that informs if we need to solve the Enthalpy/Temperature equations or not
     117             :   const bool _compute_power;
     118             :   /// Flag that informs if there is a pin mesh or not
     119             :   const bool _pin_mesh_exist;
     120             :   /// Variable that informs whether we exited external solve with a converged solution or not
     121             :   bool _converged;
     122             :   /// Time step
     123             :   const Real & _dt;
     124             :   /// Outlet Pressure
     125             :   const Real & _P_out;
     126             :   /// Thermal diffusion coefficient used in turbulent crossflow
     127             :   const Real & _beta;
     128             :   /// Turbulent modeling parameter used in axial momentum equation
     129             :   const Real & _CT;
     130             :   /// Convergence tolerance for the pressure loop in external solve
     131             :   const Real & _P_tol;
     132             :   /// Convergence tolerance for the temperature loop in external solve
     133             :   const Real & _T_tol;
     134             :   /// Maximum iterations for the inner temperature loop
     135             :   const int & _T_maxit;
     136             :   /// The relative convergence tolerance, (relative decrease) for the ksp linear solver
     137             :   const PetscReal & _rtol;
     138             :   /// The absolute convergence tolerance for the ksp linear solver
     139             :   const PetscReal & _atol;
     140             :   /// The divergence tolerance for the ksp linear solver
     141             :   const PetscReal & _dtol;
     142             :   /// The maximum number of iterations to use for the ksp linear solver
     143             :   const PetscInt & _maxit;
     144             :   /// The interpolation method used in constructing the systems
     145             :   const MooseEnum _interpolation_scheme;
     146             :   /// Flag to define the usage of a implicit or explicit solution
     147             :   const bool _implicit_bool;
     148             :   /// Flag to define the usage of staggered or collocated pressure
     149             :   const bool _staggered_pressure_bool;
     150             :   /// Segregated solve
     151             :   const bool _segregated_bool;
     152             :   /// Thermal monolithic bool
     153             :   const bool _monolithic_thermal_bool;
     154             : 
     155             :   /// Solutions handles and link to TH tables properties
     156             :   const SinglePhaseFluidProperties * _fp;
     157             :   std::unique_ptr<SolutionHandle> _mdot_soln;
     158             :   std::unique_ptr<SolutionHandle> _SumWij_soln;
     159             :   std::unique_ptr<SolutionHandle> _P_soln;
     160             :   std::unique_ptr<SolutionHandle> _DP_soln;
     161             :   std::unique_ptr<SolutionHandle> _h_soln;
     162             :   std::unique_ptr<SolutionHandle> _T_soln;
     163             :   std::unique_ptr<SolutionHandle> _Tpin_soln;
     164             :   std::unique_ptr<SolutionHandle> _rho_soln;
     165             :   std::unique_ptr<SolutionHandle> _mu_soln;
     166             :   std::unique_ptr<SolutionHandle> _S_flow_soln;
     167             :   std::unique_ptr<SolutionHandle> _w_perim_soln;
     168             :   std::unique_ptr<SolutionHandle> _q_prime_soln;
     169             : 
     170             :   /// Petsc Functions
     171             :   virtual PetscErrorCode createPetscVector(Vec & v, PetscInt n);
     172             :   virtual PetscErrorCode createPetscMatrix(Mat & M, PetscInt n, PetscInt m);
     173             :   template <class T>
     174             :   PetscErrorCode populateVectorFromDense(Vec & x,
     175             :                                          const T & solution,
     176             :                                          const unsigned int first_axial_level,
     177             :                                          const unsigned int last_axial_level,
     178             :                                          const unsigned int cross_dimension);
     179             :   template <class T>
     180             :   PetscErrorCode populateDenseFromVector(const Vec & x,
     181             :                                          T & solution,
     182             :                                          const unsigned int first_axial_level,
     183             :                                          const unsigned int last_axial_level,
     184             :                                          const unsigned int cross_dimension);
     185             :   template <class T>
     186             :   PetscErrorCode populateVectorFromHandle(Vec & x,
     187             :                                           const T & solution,
     188             :                                           const unsigned int first_axial_level,
     189             :                                           const unsigned int last_axial_level,
     190             :                                           const unsigned int cross_dimension);
     191             :   template <class T>
     192             :   PetscErrorCode populateSolutionChan(const Vec & x,
     193             :                                       T & solution,
     194             :                                       const unsigned int first_axial_level,
     195             :                                       const unsigned int last_axial_level,
     196             :                                       const unsigned int cross_dimension);
     197             : 
     198             :   template <class T>
     199             :   PetscErrorCode populateSolutionGap(const Vec & x,
     200             :                                      T & solution,
     201             :                                      const unsigned int first_axial_level,
     202             :                                      const unsigned int last_axial_level,
     203             :                                      const unsigned int cross_dimension);
     204             : 
     205             :   //// Matrices and vectors to be used in implicit assembly
     206             : 
     207             :   /// Mass conservation
     208             :   /// Mass conservation - sum of cross fluxes
     209             :   Mat _mc_sumWij_mat;
     210             :   Vec _Wij_vec;
     211             :   Vec _prod;
     212             :   Vec _prodp;
     213             :   /// Mass conservation - axial convection
     214             :   Mat _mc_axial_convection_mat;
     215             :   Vec _mc_axial_convection_rhs;
     216             :   /// Mass conservation - density time derivative
     217             :   /// No implicit matrix
     218             : 
     219             :   /// Axial momentum
     220             :   /// Axial momentum conservation - compute turbulent cross fluxes
     221             :   Mat _amc_turbulent_cross_flows_mat;
     222             :   Vec _amc_turbulent_cross_flows_rhs;
     223             :   /// Axial momentum conservation - time derivative
     224             :   Mat _amc_time_derivative_mat;
     225             :   Vec _amc_time_derivative_rhs;
     226             :   /// Axial momentum conservation - advective (Eulerian) derivative
     227             :   Mat _amc_advective_derivative_mat;
     228             :   Vec _amc_advective_derivative_rhs;
     229             :   /// Axial momentum conservation - cross flux derivative
     230             :   Mat _amc_cross_derivative_mat;
     231             :   Vec _amc_cross_derivative_rhs;
     232             :   /// Axial momentum conservation - friction force
     233             :   Mat _amc_friction_force_mat;
     234             :   Vec _amc_friction_force_rhs;
     235             :   /// Axial momentum conservation - buoyancy force
     236             :   /// No implicit matrix
     237             :   Vec _amc_gravity_rhs;
     238             :   /// Axial momentum conservation - pressure force
     239             :   Mat _amc_pressure_force_mat;
     240             :   Vec _amc_pressure_force_rhs;
     241             :   /// Axial momentum system matrix
     242             :   Mat _amc_sys_mdot_mat;
     243             :   Vec _amc_sys_mdot_rhs;
     244             : 
     245             :   /// Cross momentum
     246             :   /// Cross momentum conservation - time derivative
     247             :   Mat _cmc_time_derivative_mat;
     248             :   Vec _cmc_time_derivative_rhs;
     249             :   /// Cross momentum conservation - advective (Eulerian) derivative
     250             :   Mat _cmc_advective_derivative_mat;
     251             :   Vec _cmc_advective_derivative_rhs;
     252             :   /// Cross momentum conservation - friction force
     253             :   Mat _cmc_friction_force_mat;
     254             :   Vec _cmc_friction_force_rhs;
     255             :   /// Cross momentum conservation - pressure force
     256             :   Mat _cmc_pressure_force_mat;
     257             :   Vec _cmc_pressure_force_rhs;
     258             :   /// Lateral momentum system matrix
     259             :   Mat _cmc_sys_Wij_mat;
     260             :   Vec _cmc_sys_Wij_rhs;
     261             :   Vec _cmc_Wij_channel_dummy;
     262             : 
     263             :   /// Enthalpy
     264             :   /// Enthalpy conservation - time derivative
     265             :   Mat _hc_time_derivative_mat;
     266             :   Vec _hc_time_derivative_rhs;
     267             :   /// Enthalpy conservation - advective (Eulerian) derivative;
     268             :   Mat _hc_advective_derivative_mat;
     269             :   Vec _hc_advective_derivative_rhs;
     270             :   /// Enthalpy conservation - cross flux derivative
     271             :   Mat _hc_cross_derivative_mat;
     272             :   Vec _hc_cross_derivative_rhs;
     273             :   /// Enthalpy conservation - source and sink
     274             :   Vec _hc_added_heat_rhs;
     275             :   /// System matrices
     276             :   Mat _hc_sys_h_mat;
     277             :   Vec _hc_sys_h_rhs;
     278             :   /// No implicit matrix
     279             :   PetscInt _global_counter = 0;
     280             : 
     281             :   /// Added resistances for monolithic convergence
     282             :   PetscScalar _added_K = 0.0;
     283             :   PetscScalar _added_K_old = 1000.0;
     284             :   PetscScalar _max_sumWij;
     285             :   PetscScalar _max_sumWij_new;
     286             :   PetscScalar _correction_factor = 1.0;
     287             : 
     288             : public:
     289             :   static InputParameters validParams();
     290             : };
     291             : 
     292             : template <class T>
     293             : PetscErrorCode
     294         708 : InterWrapper1PhaseProblem::populateDenseFromVector(const Vec & x,
     295             :                                                    T & loc_solution,
     296             :                                                    const unsigned int first_axial_level,
     297             :                                                    const unsigned int last_axial_level,
     298             :                                                    const unsigned int cross_dimension)
     299             : {
     300             :   PetscErrorCode ierr;
     301             :   PetscScalar * xx;
     302             : 
     303             :   PetscFunctionBegin;
     304         708 :   ierr = VecGetArray(x, &xx);
     305         708 :   CHKERRQ(ierr);
     306        4812 :   for (unsigned int iz = first_axial_level; iz < last_axial_level + 1; iz++)
     307             :   {
     308        4104 :     unsigned int iz_ind = iz - first_axial_level;
     309      250344 :     for (unsigned int i_l = 0; i_l < cross_dimension; i_l++)
     310             :     {
     311      246240 :       loc_solution(i_l, iz) = xx[iz_ind * cross_dimension + i_l];
     312             :     }
     313             :   }
     314         708 :   ierr = VecRestoreArray(x, &xx);
     315         708 :   CHKERRQ(ierr);
     316             :   PetscFunctionReturn(LIBMESH_PETSC_SUCCESS);
     317             : }