1 //============================================================================
2 // Name : ProblemCoreFlows
5 // Copyright : CEA Saclay 2014
6 // Description : Generic class for PDEs problems
7 //============================================================================
8 /* A ProblemCoreFlows class */
10 /*! \class ProblemCoreFlows ProblemCoreFlows.hxx "ProblemCoreFlows.hxx"
11 * \brief Generic class for thermal hydraulics problems
12 * \details Common functions to CoreFlows models
15 #ifndef PROBLEMCOREFLOWS_HXX_
16 #define PROBLEMCOREFLOWS_HXX_
33 #include "CdmathException.hxx"
37 //! enumeration linearSolver
38 /*! the linearSolver can be GMRES or BiCGStab (see Petsc documentation) */
41 GMRES,/**< linearSolver is GMRES */
42 BCGS,/**< linearSolver is BiCGSstab */
43 CG/**< linearSolver is CG */
46 //! enumeration preconditioner
47 /*! the preconditioner can be ILU or LU (see Petsc documentation) */
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)*/
57 //! enumeration saveFormat
58 /*! the numerical results are saved using MED, VTK or CSV format */
61 MED,/**< MED format is used */
62 VTK,/**< VTK format is used */
63 CSV/**< CSV format is used */
66 //! enumeration TimeScheme
67 /*! The numerical method can be Explicit or Implicit */
70 Explicit,/**< Explicit numerical scheme */
71 Implicit/**< Implicit numerical scheme */
74 class ProblemCoreFlows
77 //! Constructeur par défaut
79 virtual ~ProblemCoreFlows();
81 // -*-*-*- Gestion du calcul (interface ICoCo -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
84 * \brief Alloue la mémoire et vérifie que le maillage et les conditions limites/initiales sont bien définis
85 * \Details c'est une fonction virtuelle pure, on la surcharge dans le problem Fluide .
88 virtual void initialize()=0;
91 * \brief vide la mémoire et enregistre le résultat final
92 * \Details on la surcharge dans le problem fluid
95 virtual void terminate()=0;
97 /** \fn computeTimeStep
98 * \brief Propose un pas de temps pour le calcul
99 * \details Pour proposer un pas de temps la fonction nécessite de discrétiser les opérateurs
100 * convection, diffusion et sources. Pour chacun on emploie la condition cfl.
101 * En cas de problème durant ce calcul (exemple t=tmax), la fonction renvoie stop=true sinon stop = false
102 * Cest une fonction virtuelle pure, on la surcharge dans Problème Fulide
103 * @param Stop booléen correspond a True si ya un problème dans le code (exemple t=tmax) False sinon
104 * \return dt le pas de temps
106 virtual double computeTimeStep(bool & stop)=0;
109 * \brief Enregistre le nouveau pas de temps dt et l'intègre aux opérateurs
110 * discrets (convection, diffusion, sources)
111 * \details c'est une fonction virtuelle pure, on la surcharge dans le problemfluid
112 * @param dt est le nouvel pas de temps (double)
113 * \return false si dt <0 et True sinon
115 virtual bool initTimeStep(double dt)=0;
117 /** \fn solveTimeStep
118 * \brief calcule les valeurs inconnues au pas de temps +1 .
119 * \details c'est une fonction virtuelle ,
121 * \return Renvoie false en cas de problème durant le calcul (valeurs non physiques..)
123 virtual bool solveTimeStep();//
125 /** \fn validateTimeStep
126 * \brief Valide le calcule au temps courant
127 * \details met à jour le temps présent t=t+dt, sauvegarde les champs inconnus
128 * et reinitialise dt à 0, teste la stationnarité .
129 * c'est une fonction virtuel , on la surchage dans chacun des modèles
131 * \return Renvoie false en cas de problème durant le calcul
133 virtual void validateTimeStep()=0;//
135 /** \fn abortTimeStep
136 * \brief efface les inconnues calculées par solveTimeStep() et reinitialise dt à 0
137 * \details c'est une fonction virtuelle pure, elle est surchargée dans ProblemFluid, TransportEquation et DiffusionEquation
139 virtual void abortTimeStep()=0;
142 * \brief Vérifie tous les paramètres et lance le code en supposant que la cfl ne changera pas
143 * \details Renvoie "false" en cas de problème durant le calcul
144 * C'est une fonction virtuelle pure, on la surcharge dans problème Fluide
146 * \return Renvoie "false" en cas de problème durant le calcul
148 virtual bool run();//
150 /** \fn iterateTimeStep
151 * \brief Calcul d'une sous-itération du pas de temps en cours, typiquement une itération de Newton pour un schéma implicite
152 * \details Deux paramètres booléen (converged et ok) sont retournés
153 * 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)
154 * ok vaut true si le calcul n'a pas rencontré d'erreur ou de problème particulier dans la mise à jour des variables physiques
155 * \param [in] bool, passage par reférence.
158 virtual bool iterateTimeStep(bool &converged) = 0; //??
161 * \brief vérifie la stationnairité du problème .
162 * \details Renvoie "true" si le problème atteint l'état stationnaire et "false" sinon
166 virtual bool isStationary() const;
169 * \brief Calcule la valeur du temps courant .
172 * \param [out] double
174 virtual double presentTime() const;
178 virtual vector<string> getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales)
179 virtual Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc)
180 virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale
181 virtual vector<string> getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement
182 virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement
185 Field getUnknownField() const;
187 //paramètres du calcul -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
189 /** \fn setPresentTime
190 * \brief met à jour _time (le temps courant du calcul)
195 void setPresentTime (double time);
197 /** \fn setMaxNbOfTimeStep
198 * \brief met à jour _maxNbOfTimeStep ( le nombre maximum d'itération du calcul )
203 void setMaxNbOfTimeStep(int maxNbOfTimeStep);
206 * \brief met à jour _timeMax (Le temps maximum du calcul)
211 void setTimeMax(double timeMax);
214 * \brief met à jour la _CFL
219 void setCFL(double cfl);
222 * \brief met à jour _precision (la précision du calcule)
227 void setPrecision(double precision);
229 /** \fn setInitialField
230 * \brief sets the initial field
235 void setInitialField(const Field &VV);
237 /** \fn setInitialField
238 * \brief sets the initial field from a field in a med file
240 * \param [in] string : the file name
241 * \param [in] string : the field name
242 * \param [in] int : the time step number
245 void setInitialField(string fileName, string fieldName, int timeStepNumber);
247 /** \fn setInitialFieldConstant
248 * \brief sets a constant initial field on a mesh stored in a med file
250 * \param [in] string : the file name
251 * \param [in] vector<double> : the value in each cell
254 void setInitialFieldConstant(string fileName, const vector<double> Vconstant);
256 /** \fn setInitialFieldConstant
257 * \brief sets a constant initial field
263 void setInitialFieldConstant(const Mesh& M, const Vector Vconstant);
265 /** \fn setInitialFieldConstant
266 * \brief sets a constant initial field
269 * \param [in] vector<double>
272 void setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant);
274 /** \fn setInitialFieldConstant
275 * \brief sets a constant initial field
277 * \param [in] int the space dimension
278 * \param [in] vector<double> the value in each cell
279 * \param [in] double the lowest value in the x direction
280 * \param [in] double the highest value in the x direction
281 * \param [in] string name of the left boundary
282 * \param [in] string name of the right boundary
283 * \param [in] double the lowest value in the y direction
284 * \param [in] double the highest value in the y direction
285 * \param [in] string name of the back boundary
286 * \param [in] string name of the front boundary
287 * \param [in] double the lowest value in the z direction
288 * \param [in] double the highest value in the z direction
289 * \param [in] string name of the bottom boundary
290 * \param [in] string name of the top boundary
293 void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
294 double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
295 double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
297 /** \fn setInitialFieldStepFunction
298 * \brief sets a step function initial field (Riemann problem)
303 * \param [in] double position of the discontinuity on one of the three axis
304 * \param [in] int direction (axis carrying the discontinuity) : 0 for x, 1 for y, 2 for z
307 void setInitialFieldStepFunction(const Mesh M, const Vector Vleft, const Vector Vright, double disc_pos, int direction=0);
309 /** \fn setInitialFieldStepFunction
310 * \brief sets a constant initial field
312 * \param [in] int the space dimension
313 * \param [in] vector<double> the value left of the discontinuity
314 * \param [in] vector<double> the value right of the discontinuity
315 * \param [in] double the position of the discontinuity in the x direction
316 * \param [in] double the lowest value in the x direction
317 * \param [in] double the highest value in the x direction
318 * \param [in] string name of the left boundary
319 * \param [in] string name of the right boundary
320 * \param [in] double the lowest value in the y direction
321 * \param [in] double the highest value in the y direction
322 * \param [in] string name of the back boundary
323 * \param [in] string name of the front boundary
324 * \param [in] double the lowest value in the z direction
325 * \param [in] double the highest value in the z direction
326 * \param [in] string name of the bottom boundary
327 * \param [in] string name of the top boundary
330 void setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
331 double xmin, double xmax,int nx, string leftSide, string rightSide,
332 double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
333 double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
335 /** \fn setInitialFieldSphericalStepFunction
336 * \brief sets a step function initial field with value Vin inside the ball with radius Radius and Vout outside
339 * \param [in] Vector Vin, value inside the ball
340 * \param [in] Vector Vout, value outside the ball
341 * \param [in] double radius of the ball
342 * \param [in] Vector Center, coordinates of the ball center
345 void setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double Radius, Vector Center);
348 * \brief renvoie _time (le temps courant du calcul)
351 * \param [out] double
355 /** \fn getNbTimeStep
356 * \brief renvoie _nbTimeStep le Numéro d'itération courant
359 * \param [out] unsigned
361 unsigned getNbTimeStep();
364 * \brief renvoie la _CFL
367 * \param [out] double
372 * \brief renvoie _precision (la précision du calcul)
375 * \param [out] double
377 double getPrecision();
380 * \brief renvoie _Mesh (le maillage du problème)
388 * \brief met à jour _fileName le nom du fichier
393 void setFileName(string fileName);
396 * \brief met à jour _FreqSave (la fréquence du sauvgarde de la solution)
401 void setFreqSave(int freqSave);
404 * \brief sauvgarde les données dans des fichiers MED ou VTK
405 * \details c'est une fonction virtuelle pure , on la surcharge
406 * dans chacun des modèles
409 virtual void save() = 0;
411 /** \fn getLinearSolver
412 * \brief renvoie _ksptype (le type du solveur linéaire utilisé)
415 * \param [out] string
417 string getLinearSolver() {
421 /** \fn getNumberOfVariables
422 * \brief le nombre d'inconnues du problème
425 * \return renvoie _nVar (le nombre d'inconnues du problème)
427 int getNumberOfVariables(){
431 /** \fn setWellBalancedCorrection
432 * \brief include a well balanced correction to treat stiff source terms
433 * @param boolean that is true if a well balanced correction should be applied
435 void setWellBalancedCorrection(bool wellBalancedCorr){
436 _wellBalancedCorrection=wellBalancedCorr;
439 /** \fn setLinearSolver
440 * \brief sets the linear solver and preconditioner
441 * \details virtual function overloaded by intanciable classes
442 * @param kspType linear solver type (GMRES or BICGSTAB)
443 * @param pcType preconditioner (ILU,LU or NONE)
445 void setLinearSolver(linearSolver solverName, preconditioner pcType);
447 /** \fn setNewtonSolver
448 * \brief set the Newton algorithm parameters
449 * \param [in] int maximum number of newton iterations
450 * \param [in] double precision required for the convergence of the newton scheme
453 void setNewtonSolver(double precision,int iterations=20)
455 _maxNewtonIts=iterations;
456 _precision_Newton=precision;
459 /** \fn displayConditionNumber
460 * \brief display the condition number of the preconditioned linear systems
462 void displayConditionNumber(bool display=true){
463 _conditionNumber=display;
466 /** \fn setSaveFileFormat
467 * \brief sets the numerical results file format (MED, VTK or CSV)
469 * \param [in] saveFormat
472 void setSaveFileFormat(saveFormat saveFileFormat){
473 _saveFormat=saveFileFormat;
476 //Couplages Thermohydraulique-thermique-neutronique *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
478 /** \fn setHeatPowerField
479 * \brief set the heat power field (variable in space)
484 void setHeatPowerField(Field heatPower){
485 _heatPowerField=heatPower;
486 _heatPowerFieldSet=true;
489 /** \fn setHeatPowerField
490 * \brief set the heat power field (variable in space)
492 * \param [in] string fileName (including file path)
493 * \param [in] string fieldName
496 void setHeatPowerField(string fileName, string fieldName){
497 _heatPowerField=Field(fileName, CELLS,fieldName);
498 _heatPowerFieldSet=true;
501 /** \fn setHeatSource
502 * \brief met à jour la puissance thermique ( _phi )
507 void setHeatSource(double phi){
511 /** \fn getHeatPowerField
512 * \brief renvoie le champs ?? ( _heatPowerField )
517 Field getHeatPowerField(){
518 return _heatPowerField;
521 /** \fn setRodTemperatureField ??
527 void setRodTemperatureField(Field rodTemperature){
528 _rodTemperatureField=rodTemperature;
529 _rodTemperatureFieldSet=true;
532 /** \fn setRodTemperature ??
538 void setRodTemperature(double rodTemp){
539 _rodTemperature=rodTemp;
542 /** \fn getRodTemperatureField
548 virtual Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
549 return _rodTemperatureField;
552 /** \fn setHeatTransfertCoeff
553 * \brief set the heat transfert coefficient for heat exchange between fluid and solid
558 void setHeatTransfertCoeff(double heatTransfertCoeff){
559 _heatTransfertCoeff=heatTransfertCoeff;
563 * \brief met à jour les paramètres de l'affichage
569 void setVerbose(bool verbose, bool system=false)
576 double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
577 std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
578 std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
579 Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
581 // some supplementary functions
583 /** \fn displayMatrix
584 * \brief displays a matrix of size "size x size" for profiling
585 * @param matrix is a pointer of size "size"
586 * @param size, size of the matrix
587 * @param name, string, name or description of the matrix
588 * @return displays the matrix on the terminal
590 void displayMatrix(double *matrix, int size, string name);
592 /** \fn displayMatrix
593 * \brief displays a vector of size "size" for profiling
594 * @param vector is a pointer of size "size"
595 * @param size, size of the vector
596 * @param name, string, name or description of the vector
597 * @return displays the vector on the terminal
599 void displayVector(double *vector, int size, string name);
601 /** \fn getTimeScheme
602 * \brief returns the time scheme name
604 * \param [out] enum TimeScheme (explicit or implicit)
606 TimeScheme getTimeScheme();
608 /** \fn setNumericalScheme
609 * \brief sets the numerical method ( explicit vs implicit )
611 * \param [in] TimeScheme
614 void setTimeScheme( TimeScheme method);
619 int _Ndim;//space dimension
620 int _nVar;//Number of equations to sole
621 int _Nmailles;//number of cells
622 int _Nnodes;//number of nodes
623 int _Nfaces;//number of faces
624 int _neibMaxNb;//maximum number of neighbours around a cell
625 int _neibMaxNbNodes;/* maximum number of nodes around a node */
635 double _maxvp;//valeur propre max pour calcul cfl
636 double _minl;//minimum cell diameter
638 /** boolean used to specify that a well balanced correction should be used */
639 bool _wellBalancedCorrection;
640 TimeScheme _timeScheme;
642 //Linear solver and petsc
648 int _maxPetscIts;//nombre maximum d'iteration gmres autorise au cours d'une resolution de systeme lineaire
649 int _PetscIts;//the number of iterations of the linear solver
650 int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
652 Mat _A;//Linear system matrix
653 Vec _b;//Linear system right hand side
654 double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
655 bool _conditionNumber;//computes an estimate of the condition number
657 //simulation monitoring variables
659 bool _initialDataSet;
660 bool _initializedMemory;
661 bool _restartWithNewTimeScheme;
662 bool _restartWithNewFileName;
663 double _timeMax,_time;
664 int _maxNbOfTimeStep,_nbTimeStep;
666 double _precision_Newton;
667 double _erreur_rel;//norme(Un+1-Un)
668 string _fileName;//name of the calculation
670 ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
672 //Heat transfert variables
673 Field _heatPowerField, _rodTemperatureField;
674 bool _heatPowerFieldSet, _rodTemperatureFieldSet;
675 double _heatTransfertCoeff;
676 double _heatSource, _rodTemperature;
677 double _hsatv, _hsatl;//all models appart from DiffusionEquation will need this
680 bool _verbose, _system;
682 string _path;//path to execution directory used for saving results
683 saveFormat _saveFormat;//file saving format : MED, VTK or CSV
686 #endif /* PROBLEMCOREFLOWS_HXX_ */