CoreFlows/Models/inc/ProblemCoreFlows.hxx

   1 //============================================================================
   2 // Name        : ProblemCoreFlows
   3 // Author      : M. Ndjinga
   4 // Version     :
   5 // Copyright   : CEA Saclay 2014
   6 // Description : Generic class for PDEs problems
   7 //============================================================================
   8 /* A ProblemCoreFlows class */
   9
  10 /*! \class ProblemCoreFlows ProblemCoreFlows.hxx "ProblemCoreFlows.hxx"
  11  *  \brief Generic class for thermal hydraulics problems
  12  *  \details  Common functions to CoreFlows models
  13  */
  14
  15 #ifndef PROBLEMCOREFLOWS_HXX_
  16 #define PROBLEMCOREFLOWS_HXX_
  17
  18 #include <iostream>
  19 #include <fstream>
  20 #include <unistd.h>
  21 #include <vector>
  22 #include <string>
  23 #include <map>
  24
  25 #include <petscksp.h>
  26 #include <slepceps.h>
  27 #include <slepcsvd.h>
  28
  29 #include "Field.hxx"
  30 #include "Mesh.hxx"
  31 #include "Cell.hxx"
  32 #include "Face.hxx"
  33 #include "CdmathException.hxx"
  34
  35 using namespace std;
  36
  37 //! enumeration linearSolver
  38 /*! the linearSolver can be GMRES or BiCGStab (see Petsc documentation) */
  39 enum linearSolver
  40 {
  41         GMRES,/**< linearSolver is GMRES */
  42         BCGS,/**< linearSolver is BiCGSstab */
  43         CG/**< linearSolver is CG */
  44 };
  45
  46 //! enumeration preconditioner
  47 /*! the preconditioner can be ILU or LU  (see Petsc documentation) */
  48 enum preconditioner
  49 {
  50         ILU,/**< preconditioner is ILU(0) */
  51         LU,/**< preconditioner is actually a direct solver (LU factorisation)*/
  52         NONE,/**< no preconditioner used */
  53         ICC,/**< preconditioner is ICC(0) */
  54         CHOLESKY/**< preconditioner is actually a direct solver for symmetric matrices (CHOLESKY factorisation)*/
  55 };
  56
  57 //! enumeration nonLinearSolver
  58 /*! the nonLinearSolver can be Newton_SOLVERLAB or using PETSc, Newton_PETSC_LINESEARCH, Newton_PETSC_TRUSTREGION (see Petsc documentation) */
  59 enum nonLinearSolver
  60 {
  61         Newton_SOLVERLAB,/**< nonLinearSolver is Newton_SOLVERLAB */
  62         Newton_PETSC_LINESEARCH,/**< nonLinearSolver is Newton_PETSC_LINESEARCH */
  63         Newton_PETSC_LINESEARCH_BASIC,/**< nonLinearSolver is Newton_PETSC_LINESEARCH_BASIC */
  64         Newton_PETSC_LINESEARCH_BT,/**< nonLinearSolver is Newton_PETSC_LINESEARCH_BT */
  65         Newton_PETSC_LINESEARCH_SECANT,/**< nonLinearSolver is Newton_PETSC_LINESEARCH_SECANT */
  66         Newton_PETSC_LINESEARCH_NLEQERR,/**< nonLinearSolver is Newton_PETSC_LINESEARCH_LEQERR */
  67         Newton_PETSC_TRUSTREGION,/**< nonLinearSolver is Newton_PETSC_TRUSTREGION */
  68         Newton_PETSC_NGMRES,/**< nonLinearSolver is Newton_PETSC_NGMRES */
  69         Newton_PETSC_ASPIN/**< nonLinearSolver is Newton_PETSC_ASPIN */
  70 };
  71
  72 //! enumeration saveFormat
  73 /*! the numerical results are saved using MED, VTK or CSV format */
  74 enum saveFormat
  75 {
  76         MED,/**< MED format is used  */
  77         VTK,/**< VTK format is used */
  78         CSV/**< CSV format is used */
  79 };
  80
  81 //! enumeration TimeScheme
  82 /*! The numerical method can be Explicit or Implicit  */
  83 enum TimeScheme
  84 {
  85         Explicit,/**<  Explicit numerical scheme */
  86         Implicit/**< Implicit numerical scheme */
  87 };
  88
  89 class ProblemCoreFlows
  90 {
  91 public :
  92         //! Constructeur par défaut
  93         ProblemCoreFlows();
  94         virtual ~ProblemCoreFlows();
  95
  96         // -*-*-*- Gestion du calcul (interface ICoCo) -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
  97
  98         /** \fn initialize
  99          * \brief Alloue la mémoire et vérifie que  le maillage et les conditions limites/initiales sont bien définis
 100          * \Details c'est une fonction virtuelle pure, on la surcharge dans le problem Fluide .
 101          * @param  void
 102          *  */
 103         virtual void initialize()=0;
 104
 105         /** \fn terminate
 106          * \brief vide la mémoire et enregistre le résultat final
 107          * \Details on la surcharge dans le problem fluid
 108          *  @param void
 109          *  */
 110         virtual void terminate()=0;
 111
 112         /** \fn computeTimeStep
 113          * \brief Propose un pas de temps pour le calcul
 114          * \details Pour proposer un pas de temps la fonction nécessite de discrétiser les opérateurs
 115          *  convection, diffusion et sources. Pour chacun on emploie la condition cfl.
 116          *   En cas de problème durant ce calcul (exemple t=tmax), la fonction renvoie stop=true sinon stop = false
 117          *   Cest une fonction virtuelle pure, on la surcharge dans Problème Fulide
 118          *  @param Stop booléen correspond a True si ya un problème dans le code (exemple t=tmax) False sinon
 119          *  \return  dt le pas de temps
 120          *  */
 121         virtual double computeTimeStep(bool & stop)=0;
 122
 123         /** \fn initTimeStep
 124          * \brief Enregistre le nouveau pas de temps dt et l'intègre aux opérateurs
 125          *  discrets (convection, diffusion, sources)
 126          *  \details c'est une fonction virtuelle pure, on la surcharge dans le problemfluid
 127          *  @param  dt est le nouvel pas de temps (double)
 128          *  \return false si dt <0 et True sinon
 129          *                         */
 130         virtual bool initTimeStep(double dt)=0;
 131
 132         /** \fn solveTimeStep
 133          * \brief calcule les valeurs inconnues au pas de temps +1 .
 134          *  \details c'est une fonction virtuelle
 135          *  @param  void
 136          *  \return Renvoie false en cas de problème durant le calcul (valeurs non physiques..)
 137          *  */
 138         virtual bool solveTimeStep();//
 139
 140         /** \fn validateTimeStep
 141          * \brief Valide le calcule au temps courant
 142          * \details met à jour le temps présent t=t+dt, sauvegarde les champs inconnus
 143          * et reinitialise dt à 0, teste la stationnarité .
 144          * c'est une fonction virtuel , on la surchage dans chacun des modèles
 145          * @param  void
 146          * \return  Renvoie false en cas de problème durant le calcul
 147          *  */
 148         virtual void validateTimeStep()=0;//
 149
 150         /** \fn abortTimeStep
 151          * \brief efface les inconnues calculées par solveTimeStep() et reinitialise dt à 0
 152          * \details c'est une fonction virtuelle pure, elle est surchargée dans ProblemFluid, TransportEquation et DiffusionEquation
 153          *  */
 154         virtual void abortTimeStep()=0;
 155
 156         /** \fn run
 157          * \brief Vérifie tous les paramètres et lance le code en supposant que la cfl ne changera pas
 158          * \details  Renvoie "false" en cas de problème durant le calcul
 159          * C'est une fonction virtuelle pure, on la surcharge dans problème Fluide
 160          * @param void
 161          * \return Renvoie "false" en cas de problème durant le calcul
 162          *  */
 163         virtual bool run();//
 164
 165         /** \fn iterateTimeStep
 166          * \brief Calcul d'une sous-itération du pas de temps en cours, typiquement une itération de Newton pour un schéma implicite
 167          * \details Deux paramètres booléen (converged et ok) sont retournés
 168          * converged, Vaut true si on peut passer au pas de temps suivant (shéma explicite ou convergence du schéma de Newton dans le schéma implicite)
 169          * ok  vaut true si le calcul n'a pas rencontré d'erreur ou de problème particulier dans la mise à jour des variables physiques
 170          * \param [in] bool, passage par reférence.
 171          * \param [out] bool
 172          *  */
 173         virtual bool iterateTimeStep(bool &converged) = 0; //??
 174
 175         /** \fn isStationary
 176          * \brief vérifie la stationnairité du problème .
 177          * \details Renvoie "true" si le problème atteint l'état stationnaire et "false" sinon
 178          * \param [in] void
 179          * \param [out] bool
 180          *  */
 181         virtual bool isStationary() const;
 182
 183         /** \fn presentTime
 184          * \brief Calcule la valeur du temps courant .
 185          * \details
 186          * \param [in] void
 187          * \param [out] double
 188          *  */
 189         virtual double presentTime() const;
 190
 191         /*
 192         //Coupling interface
 193         virtual vector<string> getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales)
 194         virtual  Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc)
 195         virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale
 196         virtual vector<string> getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement
 197         virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement
 198          */
 199
 200         Field getUnknownField() const;
 201
 202         //paramètres du calcul -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
 203
 204         /** \fn setPresentTime
 205          * \brief met à jour _time (le temps courant du calcul)
 206          * \details
 207          * \param [in] double
 208          * \param [out] void
 209          *  */
 210         void setPresentTime (double time);
 211
 212         /** \fn setMaxNbOfTimeStep
 213          * \brief met à jour _maxNbOfTimeStep ( le nombre maximum d'itération du calcul )
 214          * \details
 215          * \param [in] int
 216          * \param [out] void
 217          *  */
 218         void setMaxNbOfTimeStep(int maxNbOfTimeStep);
 219
 220         /** \fn setTimeMax
 221          * \brief met à jour _timeMax (Le temps maximum du calcul)
 222          * \details
 223          * \param [in] double
 224          * \param [out] void
 225          *  */
 226         void setTimeMax(double timeMax);
 227
 228         /** \fn setCFL
 229          * \brief met à jour la _CFL
 230          * \details
 231          * \param [in] double
 232          * \param [out] void
 233          *  */
 234         void setCFL(double cfl);
 235
 236         /** \fn setPrecision
 237          * \brief met à jour _precision (la précision du calcule)
 238          * \details
 239          * \param [in] double
 240          * \param [out] void
 241          *  */
 242         void setPrecision(double precision);
 243
 244         /** \fn setInitialField
 245          * \brief sets the initial field
 246          * \details
 247          * \param [in] Field
 248          * \param [out] void
 249          *  */
 250         void setInitialField(const Field &VV);
 251
 252         /** \fn setInitialField
 253          * \brief sets the initial field from a field in a med file
 254          * \details
 255          * \param [in] string : the file name
 256          * \param [in] string : the field name
 257          * \param [in] int : the time step number
 258          * \param [out] void
 259          *  */
 260         void setInitialField(string fileName, string fieldName, int timeStepNumber);
 261
 262         /** \fn setInitialFieldConstant
 263          * \brief sets a constant initial field on a mesh stored in a med file
 264          * \details
 265          * \param [in] string : the file name
 266          * \param [in] vector<double> : the value in each cell
 267          * \param [out] void
 268          *  */
 269         void setInitialFieldConstant(string fileName, const vector<double> Vconstant);
 270
 271         /** \fn setInitialFieldConstant
 272          * \brief sets a constant initial field
 273          * \details
 274          * \param [in] Mesh
 275          * \param [in] Vector
 276          * \param [out] void
 277          *  */
 278         void setInitialFieldConstant(const Mesh& M, const Vector Vconstant);
 279
 280         /** \fn setInitialFieldConstant
 281          * \brief sets a constant initial field
 282          * \details
 283          * \param [in] Mesh
 284          * \param [in] vector<double>
 285          * \param [out] void
 286          *  */
 287         void setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant);
 288
 289         /** \fn setInitialFieldConstant
 290          * \brief sets a constant initial field
 291          * \details
 292          * \param [in] int the space dimension
 293          * \param [in] vector<double> the value in each cell
 294          * \param [in] double the lowest value in the x direction
 295          * \param [in] double the highest value in the x direction
 296          * \param [in] string name of the left boundary
 297          * \param [in] string name of the right boundary
 298          * \param [in] double the lowest value in the y direction
 299          * \param [in] double the highest value in the y direction
 300          * \param [in] string name of the back boundary
 301          * \param [in] string name of the front boundary
 302          * \param [in] double the lowest value in the z direction
 303          * \param [in] double the highest value in the z direction
 304          * \param [in] string name of the bottom boundary
 305          * \param [in] string name of the top boundary
 306          * \param [out] void
 307          *  */
 308         void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
 309                         double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
 310                         double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
 311
 312         /** \fn setInitialFieldStepFunction
 313          * \brief sets a step function initial field (Riemann problem)
 314          * \details
 315          * \param [in] Mesh
 316          * \param [in] Vector
 317          * \param [in] Vector
 318          * \param [in] double position of the discontinuity on one of the three axis
 319          * \param [in] int direction (axis carrying the discontinuity) : 0 for x, 1 for y, 2 for z
 320          * \param [out] void
 321          *  */
 322         void setInitialFieldStepFunction(const Mesh M, const Vector Vleft, const Vector Vright, double disc_pos, int direction=0);
 323
 324         /** \fn setInitialFieldStepFunction
 325          * \brief sets a constant initial field
 326          * \details
 327          * \param [in] int the space dimension
 328          * \param [in] vector<double> the value left of the discontinuity
 329          * \param [in] vector<double> the value right of the discontinuity
 330          * \param [in] double the position of the discontinuity in the x direction
 331          * \param [in] double the lowest value in the x direction
 332          * \param [in] double the highest value in the x direction
 333          * \param [in] string name of the left boundary
 334          * \param [in] string name of the right boundary
 335          * \param [in] double the lowest value in the y direction
 336          * \param [in] double the highest value in the y direction
 337          * \param [in] string name of the back boundary
 338          * \param [in] string name of the front boundary
 339          * \param [in] double the lowest value in the z direction
 340          * \param [in] double the highest value in the z direction
 341          * \param [in] string name of the bottom boundary
 342          * \param [in] string name of the top boundary
 343          * \param [out] void
 344          *  */
 345         void setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
 346                         double xmin, double xmax,int nx, string leftSide, string rightSide,
 347                         double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
 348                         double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
 349
 350         /** \fn setInitialFieldSphericalStepFunction
 351          * \brief sets a step function initial field with value Vin inside the ball with radius Radius and Vout outside
 352          * \details
 353          * \param [in] Mesh
 354          * \param [in] Vector Vin, value inside the ball
 355          * \param [in] Vector Vout, value outside the ball
 356          * \param [in] double radius of the ball
 357          * \param [in] Vector Center, coordinates of the ball center
 358          * \param [out] void
 359          *  */
 360         void setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double Radius, Vector Center);
 361
 362         /** \fn getTime
 363          * \brief renvoie _time (le temps courant du calcul)
 364          * \details
 365          * \param [in] void
 366          * \param [out] double
 367          *  */
 368         double getTime();
 369
 370         /** \fn getNbTimeStep
 371          * \brief renvoie _nbTimeStep le Numéro d'itération courant
 372          * \details
 373          * \param [in] void
 374          * \param [out] unsigned
 375          *  */
 376         unsigned getNbTimeStep();
 377
 378         /** \fn getCFL
 379          * \brief renvoie la _CFL
 380          * \details
 381          * \param [in] void
 382          * \param [out] double
 383          *  */
 384         double getCFL();
 385
 386         /** \fn getPrecision
 387          * \brief renvoie _precision (la précision du calcul)
 388          * \details
 389          * \param [in] void
 390          * \param [out] double
 391          *  */
 392         double getPrecision();
 393
 394         /** \fn getMesh
 395          * \brief renvoie _Mesh (le maillage du problème)
 396          * \details
 397          * \param [in] void
 398          * \param [out] Mesh
 399          *  */
 400         Mesh getMesh();
 401
 402         /** \fn setFileName
 403          * \brief met à jour _fileName le nom du fichier
 404          * \details
 405          * \param [in]  string
 406          * \param [out] void
 407          *  */
 408         void setFileName(string fileName);
 409
 410         /** \fn setFreqSave
 411          * \brief met à jour _FreqSave (la fréquence du sauvgarde de la solution)
 412          * \details
 413          * \param [in] double
 414          * \param [out] void
 415          *  */
 416         void setFreqSave(int freqSave);
 417
 418         /** \fn save
 419          * \brief sauvgarde les données dans des fichiers MED ou VTK
 420          * \details c'est une fonction virtuelle pure , on la surcharge
 421          * dans chacun des modèles
 422          * @param  void
 423          */
 424         virtual void save() = 0;
 425
 426         /** \fn getLinearSolver
 427          * \brief renvoie _ksptype (le type du solveur linéaire utilisé)
 428          * \details
 429          * \param [in] void
 430          * \param [out] string
 431          *  */
 432         string getLinearSolver() {
 433                 return _ksptype;
 434         };
 435
 436         /** \fn getNonLinearSolver
 437          * \brief renvoie _nonLinearSolver (le type du solveur de Newton utilisé)
 438          * \details
 439          * \param [in] void
 440          * \param [out] string
 441          *  */
 442         nonLinearSolver getNonLinearSolver() {
 443                 return _nonLinearSolver;
 444         };
 445
 446         /** \fn getNumberOfVariables
 447          * \brief le nombre d'inconnues du problème
 448          * \details
 449          * @param void
 450          * \return renvoie _nVar (le nombre d'inconnues du problème)
 451          *  */
 452         int getNumberOfVariables(){
 453                 return _nVar;
 454         };
 455
 456         /** \fn setWellBalancedCorrection
 457          * \brief include a well balanced correction to treat stiff source terms
 458          * @param boolean that is true if a well balanced correction should be applied
 459          * */
 460         void setWellBalancedCorrection(bool wellBalancedCorr){
 461                 _wellBalancedCorrection=wellBalancedCorr;
 462         }
 463
 464         /** \fn setLinearSolver
 465          * \brief sets the linear solver and preconditioner
 466          * \details virtual function overloaded by intanciable classes
 467          * @param kspType linear solver type (GMRES or BICGSTAB)
 468          * @param pcType preconditioner (ILU,LU or NONE)
 469          */
 470         void setLinearSolver(linearSolver solverName, preconditioner pcType);
 471
 472         /** \fn setNewtonSolver
 473          * \brief sets the Newton type algorithm for solving the nonlinear algebraic system arising from the discretisation of the PDE
 474          * \param [in] double : precision required for the convergence of the newton scheme
 475          * \param [in] int : maximum number of newton iterations
 476          * \param [in] nonLinearSolver : the algorithm to be used to solve the nonlinear system
 477          * \param [out] void
 478          *  */
 479         void setNewtonSolver(double precision, int iterations=20, nonLinearSolver solverName=Newton_SOLVERLAB);
 480
 481         /** \fn displayConditionNumber
 482          * \brief display the condition number of the preconditioned linear systems
 483          */
 484         void displayConditionNumber(bool display=true){
 485                 _conditionNumber=display;
 486         }
 487
 488         /** \fn setSaveFileFormat
 489          * \brief sets the numerical results file format (MED, VTK or CSV)
 490          * \details
 491          * \param [in] saveFormat
 492          * \param [out] void
 493          *  */
 494         void setSaveFileFormat(saveFormat saveFileFormat){
 495                 _saveFormat=saveFileFormat;
 496         }
 497
 498         /** \fn setResultDirectory
 499          * \brief sets the directory where the results will be saved
 500          * \details
 501          * \param [in] resultsPath
 502          * \param [out] void
 503          *  */
 504         void setResultDirectory(string resultsPath){
 505                 _path=resultsPath;
 506         }
 507
 508         /** \fn getTimeScheme
 509          * \brief returns the  time scheme name
 510          * \param [in] void
 511          * \param [out] enum TimeScheme (explicit or implicit)
 512          *  */
 513         TimeScheme getTimeScheme();
 514
 515         /** \fn setNumericalScheme
 516          * \brief sets the numerical method ( explicit vs implicit )
 517          * \details
 518          * \param [in] TimeScheme
 519          * \param [out] void
 520          *  */
 521         void setTimeScheme( TimeScheme method);
 522
 523         //Couplages Thermohydraulique-thermique-neutronique *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
 524
 525         /** \fn setHeatPowerField
 526          * \brief set the heat power field (variable in space)
 527          * \details
 528          * \param [in] Field
 529          * \param [out] void
 530          *  */
 531         void setHeatPowerField(Field heatPower){
 532                 _heatPowerField=heatPower;
 533                 _heatPowerFieldSet=true;
 534         }
 535
 536         /** \fn setHeatPowerField
 537          * \brief set the heat power field (variable in space)
 538          * \details
 539          * \param [in] string fileName (including file path)
 540          * \param [in] string fieldName
 541          * \param [out] void
 542          *  */
 543         void setHeatPowerField(string fileName, string fieldName){
 544                 _heatPowerField=Field(fileName, CELLS,fieldName);
 545                 _heatPowerFieldSet=true;
 546         }
 547
 548         /** \fn setHeatSource
 549          * \brief met à jour la puissance thermique ( _phi )
 550          * \details
 551          * \param [in] double
 552          * \param [out] void
 553          *  */
 554         void setHeatSource(double phi){
 555                 _heatSource=phi;
 556         }
 557
 558         /** \fn getHeatPowerField
 559          * \brief renvoie le champs ?? ( _heatPowerField )
 560          * \details
 561          * \param [in] void
 562          * \param [out] Field
 563          *  */
 564         Field getHeatPowerField(){
 565                 return _heatPowerField;
 566         }
 567
 568         /** \fn setRodTemperatureField ??
 569          * \brief
 570          * \details
 571          * \param [in] Field
 572          * \param [out] void
 573          *  */
 574         void setRodTemperatureField(Field rodTemperature){
 575                 _rodTemperatureField=rodTemperature;
 576                 _rodTemperatureFieldSet=true;
 577         }
 578
 579         /** \fn setRodTemperature ??
 580          * \brief
 581          * \details
 582          * \param [in] double
 583          * \param [out] void
 584          *  */
 585         void setRodTemperature(double rodTemp){
 586                 _rodTemperature=rodTemp;
 587         }
 588
 589         /** \fn getRodTemperatureField
 590          * \brief
 591          * \details
 592          * \param [in] void
 593          * \param [out] Field
 594          *  */
 595         virtual Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
 596                 return _rodTemperatureField;
 597         }
 598
 599         /** \fn setHeatTransfertCoeff
 600          * \brief set the heat transfert coefficient for heat exchange between fluid and solid
 601          * \details
 602          * \param [in] double
 603          * \param [out] void
 604          *  */
 605         void setHeatTransfertCoeff(double heatTransfertCoeff){
 606                 _heatTransfertCoeff=heatTransfertCoeff;
 607         }
 608
 609         /** \fn setDISPLAY
 610          * \brief met à jour les paramètres de l'affichage
 611          * \details
 612          * \param [in] bool
 613          * \param [in] bool
 614          * \param [out] void
 615          *  */
 616         void setVerbose(bool verbose,  bool system=false)
 617         {
 618                 _verbose = verbose;
 619                 _system = system;
 620         };
 621
 622     //Spectral analysis
 623     double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
 624     std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 625     std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 626     Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 627     std::vector< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
 628     std::vector< Vector > getSingularVectors(int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
 629
 630         //  some supplementary functions
 631
 632         /** \fn displayMatrix
 633          * \brief displays a matrix of size "size x size" for profiling
 634          * @param  matrix is a pointer of size "size"
 635          * @param size, size of the matrix
 636          * @param name, string, name or description of the matrix
 637          * @return displays the matrix on the terminal
 638          *  */
 639         static void displayMatrix(double *matrix, int size, string name="Vector coefficients :");
 640
 641         /** \fn displayMatrix
 642          * \brief displays a vector of size "size" for profiling
 643          * @param  vector is a pointer of size "size"
 644          * @param size, size of the vector
 645          * @param name, string, name or description of the vector
 646          * @return displays the vector on the terminal
 647          *  */
 648         static void displayVector(double *vector, int size, string name="Matrix coefficients :");
 649
 650 protected :
 651
 652         int _Ndim;//space dimension
 653         int _nVar;//Number of equations to sole
 654         int _Nmailles;//number of cells
 655         int _Nnodes;//number of nodes
 656         int _Nfaces;//number of faces
 657         int _neibMaxNb;//maximum number of neighbours around a cell
 658         int _neibMaxNbNodes;/* maximum number of nodes around a node */
 659         Mesh _mesh;
 660         Field _perimeters;
 661
 662         //Main unknown field
 663         Field _VV;
 664
 665         //Numerical method
 666         double _dt;
 667         double _cfl;
 668         double _maxvp;//valeur propre max pour calcul cfl
 669         double _minl;//minimum cell diameter
 670     bool _FECalculation;
 671         /** boolean used to specify that a well balanced correction should be used */
 672         bool _wellBalancedCorrection;
 673         TimeScheme _timeScheme;
 674
 675         //Linear solver and petsc
 676         KSP _ksp;
 677         KSPType _ksptype;
 678         PC _pc;
 679         PCType _pctype;
 680         string _pc_hypre;
 681         nonLinearSolver _nonLinearSolver;
 682         int _maxPetscIts;//nombre maximum d'iteration gmres autorise au cours d'une resolution de systeme lineaire
 683         int _PetscIts;//the number of iterations of the linear solver
 684         int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
 685         int _NEWTON_its;
 686         Mat _A;//Linear system matrix
 687         Vec _b;//Linear system right hand side
 688         double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
 689         bool _conditionNumber;//computes an estimate of the condition number
 690
 691         //simulation monitoring variables
 692         bool _isStationary;
 693         bool _initialDataSet;
 694         bool _initializedMemory;
 695         bool _restartWithNewTimeScheme;
 696         bool _restartWithNewFileName;
 697         double _timeMax,_time;
 698         int _maxNbOfTimeStep,_nbTimeStep;
 699         double _precision;
 700         double _precision_Newton;
 701         double _erreur_rel;//norme(Un+1-Un)
 702         string _fileName;//name of the calculation
 703         int _freqSave;
 704         ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
 705
 706         //Heat transfert variables
 707         Field _heatPowerField, _rodTemperatureField;
 708         bool _heatPowerFieldSet, _rodTemperatureFieldSet;
 709         double _heatTransfertCoeff;
 710         double _heatSource, _rodTemperature;
 711         double _hsatv, _hsatl;//all models appart from DiffusionEquation will need this
 712
 713         //Display variables
 714         bool _verbose, _system;
 715
 716         string _path;//path to execution directory used for saving results
 717         saveFormat _saveFormat;//file saving format : MED, VTK or CSV
 718
 719 };
 720
 721 #endif /* PROBLEMCOREFLOWS_HXX_ */