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         /** \fn setStationaryMode
 192          * \brief Perform the search of a stationary regime
 193          * \details
 194          * \param [in] bool
 195          * \param [out]
 196          *  */
 197         virtual void setStationaryMode(bool stationaryMode){ _stationaryMode=stationaryMode;};
 198
 199         /** \fn getStationaryMode
 200          * \brief Tells if we are seeking a stationary regime
 201          * \details
 202          * \param [in]
 203          * \param [out] bool
 204          *  */
 205         virtual bool getStationaryMode(){return _stationaryMode;};
 206
 207         /** \fn resetTime
 208          * \brief sets the current time (typically to start a new calculation)
 209          * \details
 210          * \param [in] double
 211          * \param [out] void
 212          *  */
 213         void resetTime (double time);
 214
 215         /*
 216         //Coupling interface
 217         virtual void getInputMEDDoubleFieldTemplate(const std::string& name, MEDDoubleField& afield) const;//Renvoie le format de champs attendu (maillage, composantes etc)
 218         virtual vector<string> getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement
 219         virtual void getOutputMEDDoubleField(const std::string& name, MEDDoubleField& afield)//Renvoie un champs pour le postraitement ou le couplage
 220          */
 221
 222         /** Set input fields to prepare the simulation or coupling **/
 223         virtual vector<string> getInputFieldsNames()=0;
 224         virtual void setInputField(const string& nameField, Field& inputField )=0;//supply of a required input field
 225
 226      /*! @brief (Optional) Provide the code with a scalar double data.
 227      *
 228      * See Problem documentation for more details on the time semantic of a scalar value.
 229      *
 230      * @param[in] name name of the scalar value that is given to the code.
 231      * @param[in] val value passed to the code.
 232      * @throws ICoCo::WrongArgument exception if the scalar name ('name' parameter) is invalid.
 233      */
 234     /* virtual void setInputDoubleValue(const std::string& name, const double& val); */
 235
 236     /*! @brief (Optional) Retrieve a scalar double value from the code.
 237      *
 238      * See Problem documentation for more details on the time semantic of a scalar value.
 239      *
 240      * @param[in] name name of the scalar value to be read from the code.
 241      * @return the double value read from the code.
 242      * @throws ICoCo::WrongArgument exception if the scalar name ('name' parameter) is invalid.
 243      */
 244     /* virtual double getOutputDoubleValue(const std::string& name) const; */
 245
 246     /*! @brief (Optional) Similar to setInputDoubleValue() but for an int value.
 247      * @sa setInputDoubleValue()
 248      */
 249     /* virtual void setInputIntValue(const std::string& name, const int& val); */
 250
 251     /*! @brief (Optional) Similar to getOutputDoubleValue() but for an int value.
 252      * @sa getOutputDoubleValue()
 253      */
 254     /* virtual int getOutputIntValue(const std::string& name) const; */
 255
 256     /*! @brief (Optional) Similar to setInputDoubleValue() but for an string value.
 257      * @sa setInputDoubleValue()
 258      */
 259     /* virtual void setInputStringValue(const std::string& name, const std::string& val); */
 260
 261     /*! @brief (Optional) Similar to getOutputDoubleValue() but for an string value.
 262      * @sa getOutputDoubleValue()
 263      */
 264     /* virtual std::string getOutputStringValue(const std::string& name) const; */
 265
 266    /*! @brief Return ICoCo interface major version number.
 267      * @return ICoCo interface major version number (2 at present)
 268      */
 269     static int GetICoCoMajorVersion() { return 2; }
 270
 271     /*! @brief (Optional) Get MEDCoupling major version, if the code was built with MEDCoupling support.
 272      *
 273      * This can be used to assess compatibility between codes when coupling them.
 274      *
 275      * @return the MEDCoupling major version number (typically 7, 8, 9, ...)
 276      */
 277     virtual int getMEDCouplingMajorVersion() const{ return MEDCOUPLING_VERSION_MAJOR; };
 278
 279     /*! @brief (Optional) Indicate whether the code was built with a 64-bits version of MEDCoupling.
 280      *
 281      * Implemented if the code was built with MEDCoupling support.
 282      * This can be used to assess compatibility between codes when coupling them.
 283      *
 284      * @return the MEDCoupling major version number
 285      */
 286     virtual bool isMEDCoupling64Bits() const;
 287
 288     /*! @brief (Optional) Get the list of input scalars accepted by the code.
 289      *
 290      * @return the list of scalar names that represent inputs of the code
 291      * @throws ICoCo::WrongContext exception if called before initialize() or after terminate().
 292      */
 293     /* virtual std::vector<std::string> getInputValuesNames() const; */
 294
 295     /*! @brief (Optional) Get the list of output scalars that can be provided by the code.
 296      *
 297      * @return the list of scalar names that can be returned by the code
 298      * @throws ICoCo::WrongContext exception if called before initialize() or after terminate().
 299      */
 300     /* virtual std::vector<std::string> getOutputValuesNames() const; */
 301
 302         Field getUnknownField() const;
 303
 304         //paramètres du calcul -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
 305
 306         /** \fn setMaxNbOfTimeStep
 307          * \brief met à jour _maxNbOfTimeStep ( le nombre maximum d'itération du calcul )
 308          * \details
 309          * \param [in] int
 310          * \param [out] void
 311          *  */
 312         void setMaxNbOfTimeStep(int maxNbOfTimeStep);
 313
 314         /** \fn setTimeMax
 315          * \brief met à jour _timeMax (Le temps maximum du calcul)
 316          * \details
 317          * \param [in] double
 318          * \param [out] void
 319          *  */
 320         void setTimeMax(double timeMax);
 321
 322         /** \fn setCFL
 323          * \brief met à jour la _CFL
 324          * \details
 325          * \param [in] double
 326          * \param [out] void
 327          *  */
 328         void setCFL(double cfl);
 329
 330         /** \fn setPrecision
 331          * \brief met à jour _precision (la précision du calcule)
 332          * \details
 333          * \param [in] double
 334          * \param [out] void
 335          *  */
 336         void setPrecision(double precision);
 337
 338         /** \fn setInitialField
 339          * \brief sets the initial field
 340          * \details
 341          * \param [in] Field
 342          * \param [out] void
 343          *  */
 344         void setInitialField(const Field &VV);
 345
 346         /** \fn setInitialField
 347          * \brief sets the initial field from a field in a med file.
 348          * \details This function is added because we have not been able yet to swig properly the enum EntityType. It is replaced by an integer.
 349          * \param [in] string : the file name
 350          * \param [in] string : the field name
 351          * \param [in] int : the time step number
 352          * \param [in] int : int corresponding to the enum CELLS, NODES or FACES
 353          * \param [out] void
 354          *  */
 355         void setInitialField(string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type);
 356
 357         /** \fn setInitialField
 358          * \brief sets the initial field from a field in a med file
 359          * \details
 360          * \param [in] string : the file name
 361          * \param [in] string : the field name
 362          * \param [in] int : the time step number
 363          * \param [in] EntityType : CELLS, NODES or FACES
 364          * \param [out] void
 365          *  */
 366         void setInitialField(string fileName, string fieldName, int timeStepNumber, int order = 0, int meshLevel=0, EntityType typeField = CELLS);
 367
 368         /** \fn setInitialFieldConstant
 369          * \brief sets a constant initial field on a mesh stored in a med file
 370          * \details
 371          * \param [in] string : the file name
 372          * \param [in] vector<double> : the value in each cell
 373          * \param [in] EntityType : CELLS, NODES or FACES
 374          * \param [out] void
 375          *  */
 376         void setInitialFieldConstant(string fileName, const vector<double> Vconstant, EntityType typeField = CELLS);
 377
 378         /** \fn setInitialFieldConstant
 379          * \brief sets a constant initial field
 380          * \details
 381          * \param [in] Mesh
 382          * \param [in] Vector
 383          * \param [in] EntityType : CELLS, NODES or FACES
 384          * \param [out] void
 385          *  */
 386         void setInitialFieldConstant(const Mesh& M, const Vector Vconstant, EntityType typeField = CELLS);
 387
 388         /** \fn setInitialFieldConstant
 389          * \brief sets a constant initial field
 390          * \details
 391          * \param [in] Mesh
 392          * \param [in] vector<double>
 393          * \param [in] EntityType : CELLS, NODES or FACES
 394          * \param [out] void
 395          *  */
 396         void setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant, EntityType typeField = CELLS);
 397
 398         /** \fn setInitialFieldConstant
 399          * \brief sets a constant initial field
 400          * \details
 401          * \param [in] int the space dimension
 402          * \param [in] vector<double> the value in each cell
 403          * \param [in] double the lowest value in the x direction
 404          * \param [in] double the highest value in the x direction
 405          * \param [in] string name of the left boundary
 406          * \param [in] string name of the right boundary
 407          * \param [in] double the lowest value in the y direction
 408          * \param [in] double the highest value in the y direction
 409          * \param [in] string name of the back boundary
 410          * \param [in] string name of the front boundary
 411          * \param [in] double the lowest value in the z direction
 412          * \param [in] double the highest value in the z direction
 413          * \param [in] string name of the bottom boundary
 414          * \param [in] string name of the top boundary
 415          * \param [in] EntityType : CELLS, NODES or FACES
 416          * \param [out] void
 417          *  */
 418         void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
 419                         double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
 420                         double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="", EntityType typeField = CELLS);
 421
 422         /** \fn setInitialFieldConstant
 423          * \brief sets a constant initial field
 424          * \details This function is added because we have not been able yet to swig roperly the enum EntityType. It is replaced by an integer.
 425          * \param [in] int the space dimension
 426          * \param [in] vector<double> the value in each cell
 427          * \param [in] double the lowest value in the x direction
 428          * \param [in] double the highest value in the x direction
 429          * \param [in] string name of the left boundary
 430          * \param [in] string name of the right boundary
 431          * \param [in] double the lowest value in the y direction
 432          * \param [in] double the highest value in the y direction
 433          * \param [in] string name of the back boundary
 434          * \param [in] string name of the front boundary
 435          * \param [in] double the lowest value in the z direction
 436          * \param [in] double the highest value in the z direction
 437          * \param [in] string name of the bottom boundary
 438          * \param [in] string name of the top boundary
 439          * \param [in] integer corresponding to the field support enum : CELLS, NODES or FACES
 440          * \param [out] void
 441          *  */
 442         void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
 443                         double ymin, double ymax, int ny, string backSide, string frontSide,
 444                         double zmin, double zmax, int nz, string bottomSide, string topSide, int type_of_field );
 445
 446         /** \fn setInitialFieldStepFunction
 447          * \brief sets a step function initial field (Riemann problem)
 448          * \details
 449          * \param [in] Mesh
 450          * \param [in] Vector
 451          * \param [in] Vector
 452          * \param [in] double position of the discontinuity on one of the three axis
 453          * \param [in] int direction (axis carrying the discontinuity) : 0 for x, 1 for y, 2 for z
 454          * \param [in] EntityType : CELLS, NODES or FACES
 455          * \param [out] void
 456          *  */
 457         void setInitialFieldStepFunction(const Mesh M, const Vector Vleft, const Vector Vright, double disc_pos, int direction=0, EntityType typeField = CELLS);
 458
 459         /** \fn setInitialFieldStepFunction
 460          * \brief sets a constant initial field
 461          * \details
 462          * \param [in] int the space dimension
 463          * \param [in] vector<double> the value left of the discontinuity
 464          * \param [in] vector<double> the value right of the discontinuity
 465          * \param [in] double the position of the discontinuity in the x direction
 466          * \param [in] double the lowest value in the x direction
 467          * \param [in] double the highest value in the x direction
 468          * \param [in] string name of the left boundary
 469          * \param [in] string name of the right boundary
 470          * \param [in] double the lowest value in the y direction
 471          * \param [in] double the highest value in the y direction
 472          * \param [in] string name of the back boundary
 473          * \param [in] string name of the front boundary
 474          * \param [in] double the lowest value in the z direction
 475          * \param [in] double the highest value in the z direction
 476          * \param [in] string name of the bottom boundary
 477          * \param [in] string name of the top boundary
 478          * \param [in] EntityType : CELLS, NODES or FACES
 479          * \param [out] void
 480          *  */
 481         void setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
 482                         double xmin, double xmax,int nx, string leftSide, string rightSide,
 483                         double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
 484                         double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="", EntityType typeField = CELLS);
 485
 486         /** \fn setInitialFieldSphericalStepFunction
 487          * \brief sets a step function initial field with value Vin inside the ball with radius Radius and Vout outside
 488          * \details
 489          * \param [in] Mesh
 490          * \param [in] Vector Vin, value inside the ball
 491          * \param [in] Vector Vout, value outside the ball
 492          * \param [in] double radius of the ball
 493          * \param [in] Vector Center, coordinates of the ball center
 494          * \param [in] EntityType : CELLS, NODES or FACES
 495          * \param [out] void
 496          *  */
 497         void setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double Radius, Vector Center, EntityType typeField = CELLS);
 498
 499         /** \fn getTime
 500          * \brief renvoie _time (le temps courant du calcul)
 501          * \details
 502          * \param [in] void
 503          * \param [out] double
 504          *  */
 505         double getTime();
 506
 507         /** \fn getNbTimeStep
 508          * \brief renvoie _nbTimeStep le Numéro d'itération courant
 509          * \details
 510          * \param [in] void
 511          * \param [out] unsigned
 512          *  */
 513         unsigned getNbTimeStep();
 514
 515         /** \fn getCFL
 516          * \brief renvoie la _CFL
 517          * \details
 518          * \param [in] void
 519          * \param [out] double
 520          *  */
 521         double getCFL();
 522
 523         /** \fn getPrecision
 524          * \brief renvoie _precision (la précision du calcul)
 525          * \details
 526          * \param [in] void
 527          * \param [out] double
 528          *  */
 529         double getPrecision();
 530
 531         /** \fn getMesh
 532          * \brief renvoie _Mesh (le maillage du problème)
 533          * \details
 534          * \param [in] void
 535          * \param [out] Mesh
 536          *  */
 537         Mesh getMesh();
 538
 539         /** \fn setFileName
 540          * \brief met à jour _fileName le nom du fichier
 541          * \details
 542          * \param [in]  string
 543          * \param [out] void
 544          *  */
 545         void setFileName(string fileName);
 546
 547         /** \fn setFreqSave
 548          * \brief met à jour _FreqSave (la fréquence du sauvgarde de la solution)
 549          * \details
 550          * \param [in] double
 551          * \param [out] void
 552          *  */
 553         void setFreqSave(int freqSave);
 554
 555         /** \fn save
 556          * \brief sauvgarde les données dans des fichiers MED ou VTK
 557          * \details c'est une fonction virtuelle pure , on la surcharge
 558          * dans chacun des modèles
 559          * @param  void
 560          */
 561         virtual void save() = 0;
 562
 563         /** \fn getLinearSolver
 564          * \brief renvoie _ksptype (le type du solveur linéaire utilisé)
 565          * \details
 566          * \param [in] void
 567          * \param [out] string
 568          *  */
 569         string getLinearSolver() {
 570                 return _ksptype;
 571         };
 572
 573         /** \fn getNonLinearSolver
 574          * \brief renvoie _nonLinearSolver (le type du solveur de Newton utilisé)
 575          * \details
 576          * \param [in] void
 577          * \param [out] string
 578          *  */
 579         nonLinearSolver getNonLinearSolver() {
 580                 return _nonLinearSolver;
 581         };
 582
 583         /** \fn getNumberOfVariables
 584          * \brief le nombre d'inconnues du problème
 585          * \details
 586          * @param void
 587          * \return renvoie _nVar (le nombre d'inconnues du problème)
 588          *  */
 589         int getNumberOfVariables(){
 590                 return _nVar;
 591         };
 592
 593         /** \fn setWellBalancedCorrection
 594          * \brief include a well balanced correction to treat stiff source terms
 595          * @param boolean that is true if a well balanced correction should be applied
 596          * */
 597         void setWellBalancedCorrection(bool wellBalancedCorr){
 598                 _wellBalancedCorrection=wellBalancedCorr;
 599         }
 600
 601         /** \fn setLinearSolver
 602          * \brief sets the linear solver and preconditioner
 603          * \details virtual function overloaded by intanciable classes
 604          * @param kspType linear solver type (GMRES or BICGSTAB)
 605          * @param pcType preconditioner (ILU,LU or NONE)
 606          */
 607         void setLinearSolver(linearSolver solverName, preconditioner pcType, double maxIts=50);
 608
 609         /** \fn setNewtonSolver
 610          * \brief sets the Newton type algorithm for solving the nonlinear algebraic system arising from the discretisation of the PDE
 611          * \param [in] double : precision required for the convergence of the newton scheme
 612          * \param [in] int : maximum number of newton iterations
 613          * \param [in] nonLinearSolver : the algorithm to be used to solve the nonlinear system
 614          * \param [out] void
 615          *  */
 616         void setNewtonSolver(double precision, int iterations=20, nonLinearSolver solverName=Newton_SOLVERLAB);
 617
 618         /** \fn displayConditionNumber
 619          * \brief display the condition number of the preconditioned linear systems
 620          */
 621         void displayConditionNumber(bool display=true){
 622                 _conditionNumber=display;
 623         }
 624
 625         /** \fn setSaveFileFormat
 626          * \brief sets the numerical results file format (MED, VTK or CSV)
 627          * \details
 628          * \param [in] saveFormat
 629          * \param [out] void
 630          *  */
 631         void setSaveFileFormat(saveFormat saveFileFormat){
 632                 _saveFormat=saveFileFormat;
 633         }
 634
 635         /** \fn setResultDirectory
 636          * \brief sets the directory where the results will be saved
 637          * \details
 638          * \param [in] resultsPath
 639          * \param [out] void
 640          *  */
 641         void setResultDirectory(string resultsPath){
 642                 _path=resultsPath;
 643         }
 644
 645         /** \fn getTimeScheme
 646          * \brief returns the  time scheme name
 647          * \param [in] void
 648          * \param [out] enum TimeScheme (explicit or implicit)
 649          *  */
 650         TimeScheme getTimeScheme();
 651
 652         /** \fn setNumericalScheme
 653          * \brief sets the numerical method ( explicit vs implicit )
 654          * \details
 655          * \param [in] TimeScheme
 656          * \param [out] void
 657          *  */
 658         void setTimeScheme( TimeScheme method);
 659
 660         //Couplages Thermohydraulique-thermique-neutronique *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
 661
 662         /** \fn setHeatPowerField
 663          * \brief set the heat power field (variable in space)
 664          * \details
 665          * \param [in] Field
 666          * \param [out] void
 667          *  */
 668         void setHeatPowerField(Field heatPower);
 669
 670         /** \fn setHeatPowerField
 671          * \brief set the heat power field (variable in space)
 672          * \details
 673          * \param [in] string fileName (including file path)
 674          * \param [in] string fieldName
 675          * \param [out] void
 676          *  */
 677         void setHeatPowerField(string fileName, string fieldName, int iteration = 0, int order = 0, int meshLevel=0);
 678
 679         /** \fn setHeatSource
 680          * \brief sets a constant heat power field
 681          * \details
 682          * \param [in] double
 683          * \param [out] void
 684          *  */
 685         void setHeatSource(double phi){
 686                 _heatSource=phi;
 687                 _isStationary=false;//Source term may be changed after previously reaching a stationary state
 688         }
 689
 690         /** \fn getHeatPowerField
 691          * \brief returns the heat power field
 692          * \details
 693          * \param [in] void
 694          * \param [out] Field
 695          *  */
 696         Field getHeatPowerField(){
 697                 return _heatPowerField;
 698         }
 699
 700         /** \fn setHeatTransfertCoeff
 701          * \brief set the heat transfert coefficient for heat exchange between fluid and solid
 702          * \details
 703          * \param [in] double
 704          * \param [out] void
 705          *  */
 706         void setHeatTransfertCoeff(double heatTransfertCoeff){
 707                 _heatTransfertCoeff=heatTransfertCoeff;
 708                 _isStationary=false;//Source term may be changed after previously reaching a stationary state
 709         }
 710
 711         /** \fn setVerbose
 712          * \brief Updates display options
 713          * \details
 714          * \param [in] bool
 715          * \param [in] bool
 716          * \param [out] void
 717          *  */
 718         void setVerbose(bool verbose,  bool system=false)
 719         {
 720                 _verbose = verbose;
 721                 _system = system;
 722         };
 723
 724     //Spectral analysis
 725     double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
 726     std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 727     std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 728     Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 729     std::vector< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
 730     std::vector< Vector > getSingularVectors(int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
 731
 732         //  some supplementary functions
 733
 734         /** \fn displayMatrix
 735          * \brief displays a matrix of size "size x size" for profiling
 736          * @param  matrix is a pointer of size "size"
 737          * @param size, size of the matrix
 738          * @param name, string, name or description of the matrix
 739          * @return displays the matrix on the terminal
 740          *  */
 741         static void displayMatrix(double *matrix, int size, string name="Vector coefficients :");
 742
 743         /** \fn displayMatrix
 744          * \brief displays a vector of size "size" for profiling
 745          * @param  vector is a pointer of size "size"
 746          * @param size, size of the vector
 747          * @param name, string, name or description of the vector
 748          * @return displays the vector on the terminal
 749          *  */
 750         static void displayVector(double *vector, int size, string name="Matrix coefficients :");
 751
 752 protected :
 753
 754         int _Ndim;//space dimension
 755         int _nVar;//Number of equations to sole
 756         int _Nmailles;//number of cells
 757         int _Nnodes;//number of nodes
 758         int _Nfaces;//number of faces
 759         int _neibMaxNbCells;//maximum number of neighbours around a cell
 760         int _neibMaxNbNodes;/* maximum number of nodes around a node */
 761         Mesh _mesh;
 762         Field _perimeters;
 763
 764         //Main unknown field
 765         Field _VV;
 766
 767         //Numerical method
 768         double _dt;
 769         double _cfl;
 770         double _maxvp;//valeur propre max pour calcul cfl
 771         double _minl;//minimum cell diameter
 772     bool _FECalculation;
 773         /** boolean used to specify that a well balanced correction should be used */
 774         bool _wellBalancedCorrection;
 775         TimeScheme _timeScheme;
 776
 777         //Linear solver and petsc
 778         KSP _ksp;
 779         KSPType _ksptype;
 780         PC _pc;
 781         PCType _pctype;
 782         string _pc_hypre;
 783         nonLinearSolver _nonLinearSolver;
 784         int _maxPetscIts;//nombre maximum d'iteration gmres autorise au cours d'une resolution de systeme lineaire
 785         int _PetscIts;//the number of iterations of the linear solver
 786         int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
 787         int _NEWTON_its;
 788         Mat _A;//Linear system matrix
 789         Vec _b;//Linear system right hand side
 790         double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
 791         bool _conditionNumber;//computes an estimate of the condition number
 792
 793         //simulation monitoring variables
 794         bool _isStationary;
 795         bool _initialDataSet;
 796         bool _initializedMemory;
 797         bool _stationaryMode;//ICoCo V2
 798         bool _restartWithNewTimeScheme;
 799         bool _restartWithNewFileName;
 800         double _timeMax,_time;
 801         int _maxNbOfTimeStep,_nbTimeStep;
 802         double _precision;
 803         double _precision_Newton;
 804         double _erreur_rel;//norme(Un+1-Un)
 805         string _fileName;//name of the calculation
 806         int _freqSave;
 807         ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
 808
 809         //Heat transfert variables
 810         Field _heatPowerField;
 811         bool _heatPowerFieldSet;
 812         double _heatTransfertCoeff;
 813         double _heatSource;
 814         double _hsatv, _hsatl;//all models appart from DiffusionEquation will need this
 815
 816         //Display variables
 817         bool _verbose, _system;
 818
 819         string _path;//path to execution directory used for saving results
 820         saveFormat _saveFormat;//file saving format : MED, VTK or CSV
 821
 822         //MPI related variables
 823         PetscMPIInt    _mpi_size;        /* size of communicator */
 824         PetscMPIInt    _mpi_rank;        /* processor rank */
 825         VecScatter         _scat;                       /* For the distribution of a local vector */
 826         int _globalNbUnknowns, _localNbUnknowns;
 827         int _d_nnz, _o_nnz;                     /* local and "non local" numbers of non zeros corfficients */
 828 };
 829
 830 #endif /* PROBLEMCOREFLOWS_HXX_ */