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