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(MPI_Comm comm = MPI_COMM_WORLD);
  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 This function is added because we have not been able yet to swig properly the enum EntityType. It is replaced by an integer.
 255          * \param [in] string : the file name
 256          * \param [in] string : the field name
 257          * \param [in] int : the time step number
 258          * \param [in] int : int corresponding to the enum CELLS, NODES or FACES
 259          * \param [out] void
 260          *  */
 261         void setInitialField(string fileName, string fieldName, int timeStepNumber, int field_support_type);
 262
 263         /** \fn setInitialField
 264          * \brief sets the initial field from a field in a med file
 265          * \details
 266          * \param [in] string : the file name
 267          * \param [in] string : the field name
 268          * \param [in] int : the time step number
 269          * \param [in] EntityType : CELLS, NODES or FACES
 270          * \param [out] void
 271          *  */
 272         void setInitialField(string fileName, string fieldName, int timeStepNumber, EntityType typeField = CELLS);
 273
 274         /** \fn setInitialFieldConstant
 275          * \brief sets a constant initial field on a mesh stored in a med file
 276          * \details
 277          * \param [in] string : the file name
 278          * \param [in] vector<double> : the value in each cell
 279          * \param [in] EntityType : CELLS, NODES or FACES
 280          * \param [out] void
 281          *  */
 282         void setInitialFieldConstant(string fileName, const vector<double> Vconstant, EntityType typeField = CELLS);
 283
 284         /** \fn setInitialFieldConstant
 285          * \brief sets a constant initial field
 286          * \details
 287          * \param [in] Mesh
 288          * \param [in] Vector
 289          * \param [in] EntityType : CELLS, NODES or FACES
 290          * \param [out] void
 291          *  */
 292         void setInitialFieldConstant(const Mesh& M, const Vector Vconstant, EntityType typeField = CELLS);
 293
 294         /** \fn setInitialFieldConstant
 295          * \brief sets a constant initial field
 296          * \details
 297          * \param [in] Mesh
 298          * \param [in] vector<double>
 299          * \param [in] EntityType : CELLS, NODES or FACES
 300          * \param [out] void
 301          *  */
 302         void setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant, EntityType typeField = CELLS);
 303
 304         /** \fn setInitialFieldConstant
 305          * \brief sets a constant initial field
 306          * \details
 307          * \param [in] int the space dimension
 308          * \param [in] vector<double> the value in each cell
 309          * \param [in] double the lowest value in the x direction
 310          * \param [in] double the highest value in the x direction
 311          * \param [in] string name of the left boundary
 312          * \param [in] string name of the right boundary
 313          * \param [in] double the lowest value in the y direction
 314          * \param [in] double the highest value in the y direction
 315          * \param [in] string name of the back boundary
 316          * \param [in] string name of the front boundary
 317          * \param [in] double the lowest value in the z direction
 318          * \param [in] double the highest value in the z direction
 319          * \param [in] string name of the bottom boundary
 320          * \param [in] string name of the top boundary
 321          * \param [in] EntityType : CELLS, NODES or FACES
 322          * \param [out] void
 323          *  */
 324         void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
 325                         double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
 326                         double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="", EntityType typeField = CELLS);
 327
 328         /** \fn setInitialFieldConstant
 329          * \brief sets a constant initial field
 330          * \details This function is added because we have not been able yet to swig roperly the enum EntityType. It is replaced by an integer.
 331          * \param [in] int the space dimension
 332          * \param [in] vector<double> the value in each cell
 333          * \param [in] double the lowest value in the x direction
 334          * \param [in] double the highest value in the x direction
 335          * \param [in] string name of the left boundary
 336          * \param [in] string name of the right boundary
 337          * \param [in] double the lowest value in the y direction
 338          * \param [in] double the highest value in the y direction
 339          * \param [in] string name of the back boundary
 340          * \param [in] string name of the front boundary
 341          * \param [in] double the lowest value in the z direction
 342          * \param [in] double the highest value in the z direction
 343          * \param [in] string name of the bottom boundary
 344          * \param [in] string name of the top boundary
 345          * \param [in] integer corresponding to the field support enum : CELLS, NODES or FACES
 346          * \param [out] void
 347          *  */
 348         void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
 349                         double ymin, double ymax, int ny, string backSide, string frontSide,
 350                         double zmin, double zmax, int nz, string bottomSide, string topSide, int type_of_field );
 351
 352         /** \fn setInitialFieldStepFunction
 353          * \brief sets a step function initial field (Riemann problem)
 354          * \details
 355          * \param [in] Mesh
 356          * \param [in] Vector
 357          * \param [in] Vector
 358          * \param [in] double position of the discontinuity on one of the three axis
 359          * \param [in] int direction (axis carrying the discontinuity) : 0 for x, 1 for y, 2 for z
 360          * \param [in] EntityType : CELLS, NODES or FACES
 361          * \param [out] void
 362          *  */
 363         void setInitialFieldStepFunction(const Mesh M, const Vector Vleft, const Vector Vright, double disc_pos, int direction=0, EntityType typeField = CELLS);
 364
 365         /** \fn setInitialFieldStepFunction
 366          * \brief sets a constant initial field
 367          * \details
 368          * \param [in] int the space dimension
 369          * \param [in] vector<double> the value left of the discontinuity
 370          * \param [in] vector<double> the value right of the discontinuity
 371          * \param [in] double the position of the discontinuity in the x direction
 372          * \param [in] double the lowest value in the x direction
 373          * \param [in] double the highest value in the x direction
 374          * \param [in] string name of the left boundary
 375          * \param [in] string name of the right boundary
 376          * \param [in] double the lowest value in the y direction
 377          * \param [in] double the highest value in the y direction
 378          * \param [in] string name of the back boundary
 379          * \param [in] string name of the front boundary
 380          * \param [in] double the lowest value in the z direction
 381          * \param [in] double the highest value in the z direction
 382          * \param [in] string name of the bottom boundary
 383          * \param [in] string name of the top boundary
 384          * \param [in] EntityType : CELLS, NODES or FACES
 385          * \param [out] void
 386          *  */
 387         void setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
 388                         double xmin, double xmax,int nx, string leftSide, string rightSide,
 389                         double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
 390                         double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="", EntityType typeField = CELLS);
 391
 392         /** \fn setInitialFieldSphericalStepFunction
 393          * \brief sets a step function initial field with value Vin inside the ball with radius Radius and Vout outside
 394          * \details
 395          * \param [in] Mesh
 396          * \param [in] Vector Vin, value inside the ball
 397          * \param [in] Vector Vout, value outside the ball
 398          * \param [in] double radius of the ball
 399          * \param [in] Vector Center, coordinates of the ball center
 400          * \param [in] EntityType : CELLS, NODES or FACES
 401          * \param [out] void
 402          *  */
 403         void setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double Radius, Vector Center, EntityType typeField = CELLS);
 404
 405         /** \fn getTime
 406          * \brief renvoie _time (le temps courant du calcul)
 407          * \details
 408          * \param [in] void
 409          * \param [out] double
 410          *  */
 411         double getTime();
 412
 413         /** \fn getNbTimeStep
 414          * \brief renvoie _nbTimeStep le Numéro d'itération courant
 415          * \details
 416          * \param [in] void
 417          * \param [out] unsigned
 418          *  */
 419         unsigned getNbTimeStep();
 420
 421         /** \fn getCFL
 422          * \brief renvoie la _CFL
 423          * \details
 424          * \param [in] void
 425          * \param [out] double
 426          *  */
 427         double getCFL();
 428
 429         /** \fn getPrecision
 430          * \brief renvoie _precision (la précision du calcul)
 431          * \details
 432          * \param [in] void
 433          * \param [out] double
 434          *  */
 435         double getPrecision();
 436
 437         /** \fn getMesh
 438          * \brief renvoie _Mesh (le maillage du problème)
 439          * \details
 440          * \param [in] void
 441          * \param [out] Mesh
 442          *  */
 443         Mesh getMesh();
 444
 445         /** \fn setFileName
 446          * \brief met à jour _fileName le nom du fichier
 447          * \details
 448          * \param [in]  string
 449          * \param [out] void
 450          *  */
 451         void setFileName(string fileName);
 452
 453         /** \fn setFreqSave
 454          * \brief met à jour _FreqSave (la fréquence du sauvgarde de la solution)
 455          * \details
 456          * \param [in] double
 457          * \param [out] void
 458          *  */
 459         void setFreqSave(int freqSave);
 460
 461         /** \fn save
 462          * \brief sauvgarde les données dans des fichiers MED ou VTK
 463          * \details c'est une fonction virtuelle pure , on la surcharge
 464          * dans chacun des modèles
 465          * @param  void
 466          */
 467         virtual void save() = 0;
 468
 469         /** \fn getLinearSolver
 470          * \brief renvoie _ksptype (le type du solveur linéaire utilisé)
 471          * \details
 472          * \param [in] void
 473          * \param [out] string
 474          *  */
 475         string getLinearSolver() {
 476                 return _ksptype;
 477         };
 478
 479         /** \fn getNonLinearSolver
 480          * \brief renvoie _nonLinearSolver (le type du solveur de Newton utilisé)
 481          * \details
 482          * \param [in] void
 483          * \param [out] string
 484          *  */
 485         nonLinearSolver getNonLinearSolver() {
 486                 return _nonLinearSolver;
 487         };
 488
 489         /** \fn getNumberOfVariables
 490          * \brief le nombre d'inconnues du problème
 491          * \details
 492          * @param void
 493          * \return renvoie _nVar (le nombre d'inconnues du problème)
 494          *  */
 495         int getNumberOfVariables(){
 496                 return _nVar;
 497         };
 498
 499         /** \fn setWellBalancedCorrection
 500          * \brief include a well balanced correction to treat stiff source terms
 501          * @param boolean that is true if a well balanced correction should be applied
 502          * */
 503         void setWellBalancedCorrection(bool wellBalancedCorr){
 504                 _wellBalancedCorrection=wellBalancedCorr;
 505         }
 506
 507         /** \fn setLinearSolver
 508          * \brief sets the linear solver and preconditioner
 509          * \details virtual function overloaded by intanciable classes
 510          * @param kspType linear solver type (GMRES or BICGSTAB)
 511          * @param pcType preconditioner (ILU,LU or NONE)
 512          */
 513         void setLinearSolver(linearSolver solverName, preconditioner pcType, double maxIts=50);
 514
 515         /** \fn setNewtonSolver
 516          * \brief sets the Newton type algorithm for solving the nonlinear algebraic system arising from the discretisation of the PDE
 517          * \param [in] double : precision required for the convergence of the newton scheme
 518          * \param [in] int : maximum number of newton iterations
 519          * \param [in] nonLinearSolver : the algorithm to be used to solve the nonlinear system
 520          * \param [out] void
 521          *  */
 522         void setNewtonSolver(double precision, int iterations=20, nonLinearSolver solverName=Newton_SOLVERLAB);
 523
 524         /** \fn displayConditionNumber
 525          * \brief display the condition number of the preconditioned linear systems
 526          */
 527         void displayConditionNumber(bool display=true){
 528                 _conditionNumber=display;
 529         }
 530
 531         /** \fn setSaveFileFormat
 532          * \brief sets the numerical results file format (MED, VTK or CSV)
 533          * \details
 534          * \param [in] saveFormat
 535          * \param [out] void
 536          *  */
 537         void setSaveFileFormat(saveFormat saveFileFormat){
 538                 _saveFormat=saveFileFormat;
 539         }
 540
 541         /** \fn setResultDirectory
 542          * \brief sets the directory where the results will be saved
 543          * \details
 544          * \param [in] resultsPath
 545          * \param [out] void
 546          *  */
 547         void setResultDirectory(string resultsPath){
 548                 _path=resultsPath;
 549         }
 550
 551         /** \fn getTimeScheme
 552          * \brief returns the  time scheme name
 553          * \param [in] void
 554          * \param [out] enum TimeScheme (explicit or implicit)
 555          *  */
 556         TimeScheme getTimeScheme();
 557
 558         /** \fn setNumericalScheme
 559          * \brief sets the numerical method ( explicit vs implicit )
 560          * \details
 561          * \param [in] TimeScheme
 562          * \param [out] void
 563          *  */
 564         void setTimeScheme( TimeScheme method);
 565
 566         //Couplages Thermohydraulique-thermique-neutronique *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
 567
 568         /** \fn setHeatPowerField
 569          * \brief set the heat power field (variable in space)
 570          * \details
 571          * \param [in] Field
 572          * \param [out] void
 573          *  */
 574         void setHeatPowerField(Field heatPower){
 575                 _heatPowerField=heatPower;
 576                 _heatPowerFieldSet=true;
 577         }
 578
 579         /** \fn setHeatPowerField
 580          * \brief set the heat power field (variable in space)
 581          * \details
 582          * \param [in] string fileName (including file path)
 583          * \param [in] string fieldName
 584          * \param [out] void
 585          *  */
 586         void setHeatPowerField(string fileName, string fieldName){
 587                 _heatPowerField=Field(fileName, CELLS,fieldName);
 588                 _heatPowerFieldSet=true;
 589         }
 590
 591         /** \fn setHeatSource
 592          * \brief met à jour la puissance thermique ( _phi )
 593          * \details
 594          * \param [in] double
 595          * \param [out] void
 596          *  */
 597         void setHeatSource(double phi){
 598                 _heatSource=phi;
 599         }
 600
 601         /** \fn getHeatPowerField
 602          * \brief renvoie le champs ?? ( _heatPowerField )
 603          * \details
 604          * \param [in] void
 605          * \param [out] Field
 606          *  */
 607         Field getHeatPowerField(){
 608                 return _heatPowerField;
 609         }
 610
 611         /** \fn setRodTemperatureField ??
 612          * \brief
 613          * \details
 614          * \param [in] Field
 615          * \param [out] void
 616          *  */
 617         void setRodTemperatureField(Field rodTemperature){
 618                 _rodTemperatureField=rodTemperature;
 619                 _rodTemperatureFieldSet=true;
 620         }
 621
 622         /** \fn setRodTemperature ??
 623          * \brief
 624          * \details
 625          * \param [in] double
 626          * \param [out] void
 627          *  */
 628         void setRodTemperature(double rodTemp){
 629                 _rodTemperature=rodTemp;
 630         }
 631
 632         /** \fn getRodTemperatureField
 633          * \brief
 634          * \details
 635          * \param [in] void
 636          * \param [out] Field
 637          *  */
 638         virtual Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
 639                 return _rodTemperatureField;
 640         }
 641
 642         /** \fn setHeatTransfertCoeff
 643          * \brief set the heat transfert coefficient for heat exchange between fluid and solid
 644          * \details
 645          * \param [in] double
 646          * \param [out] void
 647          *  */
 648         void setHeatTransfertCoeff(double heatTransfertCoeff){
 649                 _heatTransfertCoeff=heatTransfertCoeff;
 650         }
 651
 652         /** \fn setVerbose
 653          * \brief Updates display options
 654          * \details
 655          * \param [in] bool
 656          * \param [in] bool
 657          * \param [out] void
 658          *  */
 659         void setVerbose(bool verbose,  bool system=false)
 660         {
 661                 _verbose = verbose;
 662                 _system = system;
 663         };
 664
 665     //Spectral analysis
 666     double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
 667     std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 668     std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 669     Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 670     std::vector< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
 671     std::vector< Vector > getSingularVectors(int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
 672
 673         //  some supplementary functions
 674
 675         /** \fn displayMatrix
 676          * \brief displays a matrix of size "size x size" for profiling
 677          * @param  matrix is a pointer of size "size"
 678          * @param size, size of the matrix
 679          * @param name, string, name or description of the matrix
 680          * @return displays the matrix on the terminal
 681          *  */
 682         static void displayMatrix(double *matrix, int size, string name="Vector coefficients :");
 683
 684         /** \fn displayMatrix
 685          * \brief displays a vector of size "size" for profiling
 686          * @param  vector is a pointer of size "size"
 687          * @param size, size of the vector
 688          * @param name, string, name or description of the vector
 689          * @return displays the vector on the terminal
 690          *  */
 691         static void displayVector(double *vector, int size, string name="Matrix coefficients :");
 692
 693 protected :
 694
 695         int _Ndim;//space dimension
 696         int _nVar;//Number of equations to sole
 697         int _Nmailles;//number of cells
 698         int _Nnodes;//number of nodes
 699         int _Nfaces;//number of faces
 700         int _neibMaxNb;//maximum number of neighbours around a cell
 701         int _neibMaxNbNodes;/* maximum number of nodes around a node */
 702         Mesh _mesh;
 703         Field _perimeters;
 704
 705         //Main unknown field
 706         Field _VV;
 707
 708         //Numerical method
 709         double _dt;
 710         double _cfl;
 711         double _maxvp;//valeur propre max pour calcul cfl
 712         double _minl;//minimum cell diameter
 713     bool _FECalculation;
 714         /** boolean used to specify that a well balanced correction should be used */
 715         bool _wellBalancedCorrection;
 716         TimeScheme _timeScheme;
 717
 718         //Linear solver and petsc
 719         KSP _ksp;
 720         KSPType _ksptype;
 721         PC _pc;
 722         PCType _pctype;
 723         string _pc_hypre;
 724         nonLinearSolver _nonLinearSolver;
 725         int _maxPetscIts;//nombre maximum d'iteration gmres autorise au cours d'une resolution de systeme lineaire
 726         int _PetscIts;//the number of iterations of the linear solver
 727         int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
 728         int _NEWTON_its;
 729         Mat _A;//Linear system matrix
 730         Vec _b;//Linear system right hand side
 731         double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
 732         bool _conditionNumber;//computes an estimate of the condition number
 733
 734         //simulation monitoring variables
 735         bool _isStationary;
 736         bool _initialDataSet;
 737         bool _initializedMemory;
 738         bool _restartWithNewTimeScheme;
 739         bool _restartWithNewFileName;
 740         double _timeMax,_time;
 741         int _maxNbOfTimeStep,_nbTimeStep;
 742         double _precision;
 743         double _precision_Newton;
 744         double _erreur_rel;//norme(Un+1-Un)
 745         string _fileName;//name of the calculation
 746         int _freqSave;
 747         ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
 748
 749         //Heat transfert variables
 750         Field _heatPowerField, _rodTemperatureField;
 751         bool _heatPowerFieldSet, _rodTemperatureFieldSet;
 752         double _heatTransfertCoeff;
 753         double _heatSource, _rodTemperature;
 754         double _hsatv, _hsatl;//all models appart from DiffusionEquation will need this
 755
 756         //Display variables
 757         bool _verbose, _system;
 758
 759         string _path;//path to execution directory used for saving results
 760         saveFormat _saveFormat;//file saving format : MED, VTK or CSV
 761
 762         //MPI variables
 763         PetscMPIInt    _size;        /* size of communicator */
 764         PetscMPIInt    _rank;        /* processor rank */
 765
 766
 767 };
 768
 769 #endif /* PROBLEMCOREFLOWS_HXX_ */