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();
80 // -*-*-*- Gestion du calcul -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
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 //paramètres du calcul -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
187 /** \fn setPresentTime
188 * \brief met à jour _time (le temps courant du calcul)
193 void setPresentTime (double time);
195 /** \fn setMaxNbOfTimeStep
196 * \brief met à jour _maxNbOfTimeStep ( le nombre maximum d'itération du calcul )
201 void setMaxNbOfTimeStep(int maxNbOfTimeStep);
204 * \brief met à jour _timeMax (Le temps maximum du calcul)
209 void setTimeMax(double timeMax);
212 * \brief met à jour la _CFL
217 void setCFL(double cfl);
220 * \brief met à jour _precision (la précision du calcule)
225 void setPrecision(double precision);
227 /** \fn setInitialField
228 * \brief sets the initial field
233 void setInitialField(const Field &VV);
235 /** \fn setInitialField
236 * \brief sets the initial field from a field in a med file
238 * \param [in] string : the file name
239 * \param [in] string : the field name
240 * \param [in] int : the time step number
243 void setInitialField(string fileName, string fieldName, int timeStepNumber);
245 /** \fn setInitialFieldConstant
246 * \brief sets a constant initial field on a mesh stored in a med file
248 * \param [in] string : the file name
249 * \param [in] vector<double> : the value in each cell
252 void setInitialFieldConstant(string fileName, const vector<double> Vconstant);
254 /** \fn setInitialFieldConstant
255 * \brief sets a constant initial field
261 void setInitialFieldConstant(const Mesh& M, const Vector Vconstant);
263 /** \fn setInitialFieldConstant
264 * \brief sets a constant initial field
267 * \param [in] vector<double>
270 void setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant);
272 /** \fn setInitialFieldConstant
273 * \brief sets a constant initial field
275 * \param [in] int the space dimension
276 * \param [in] vector<double> the value in each cell
277 * \param [in] double the lowest value in the x direction
278 * \param [in] double the highest value in the x direction
279 * \param [in] string name of the left boundary
280 * \param [in] string name of the right boundary
281 * \param [in] double the lowest value in the y direction
282 * \param [in] double the highest value in the y direction
283 * \param [in] string name of the back boundary
284 * \param [in] string name of the front boundary
285 * \param [in] double the lowest value in the z direction
286 * \param [in] double the highest value in the z direction
287 * \param [in] string name of the bottom boundary
288 * \param [in] string name of the top boundary
291 void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
292 double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
293 double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
295 /** \fn setInitialFieldStepFunction
296 * \brief sets a step function initial field (Riemann problem)
301 * \param [in] double position of the discontinuity on one of the three axis
302 * \param [in] int direction (axis carrying the discontinuity) : 0 for x, 1 for y, 2 for z
305 void setInitialFieldStepFunction(const Mesh M, const Vector Vleft, const Vector Vright, double disc_pos, int direction=0);
307 /** \fn setInitialFieldStepFunction
308 * \brief sets a constant initial field
310 * \param [in] int the space dimension
311 * \param [in] vector<double> the value left of the discontinuity
312 * \param [in] vector<double> the value right of the discontinuity
313 * \param [in] double the position of the discontinuity in the x direction
314 * \param [in] double the lowest value in the x direction
315 * \param [in] double the highest value in the x direction
316 * \param [in] string name of the left boundary
317 * \param [in] string name of the right boundary
318 * \param [in] double the lowest value in the y direction
319 * \param [in] double the highest value in the y direction
320 * \param [in] string name of the back boundary
321 * \param [in] string name of the front boundary
322 * \param [in] double the lowest value in the z direction
323 * \param [in] double the highest value in the z direction
324 * \param [in] string name of the bottom boundary
325 * \param [in] string name of the top boundary
328 void setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
329 double xmin, double xmax,int nx, string leftSide, string rightSide,
330 double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
331 double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
333 /** \fn setInitialFieldSphericalStepFunction
334 * \brief sets a step function initial field with value Vin inside the ball with radius Radius and Vout outside
337 * \param [in] Vector Vin, value inside the ball
338 * \param [in] Vector Vout, value outside the ball
339 * \param [in] double radius of the ball
340 * \param [in] Vector Center, coordinates of the ball center
343 void setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double Radius, Vector Center);
346 * \brief renvoie _time (le temps courant du calcul)
349 * \param [out] double
353 /** \fn getNbTimeStep
354 * \brief renvoie _nbTimeStep le Numéro d'itération courant
357 * \param [out] unsigned
359 unsigned getNbTimeStep();
362 * \brief renvoie la _CFL
365 * \param [out] double
370 * \brief renvoie _precision (la précision du calcul)
373 * \param [out] double
375 double getPrecision();
378 * \brief renvoie _Mesh (le maillage du problème)
386 * \brief met à jour _fileName le nom du fichier
391 void setFileName(string fileName);
394 * \brief met à jour _FreqSave (la fréquence du sauvgarde de la solution)
399 void setFreqSave(int freqSave);
402 * \brief sauvgarde les données dans des fichiers MED ou VTK
403 * \details c'est une fonction virtuelle pure , on la surcharge
404 * dans chacun des modèles
407 virtual void save() = 0;
409 /** \fn getLinearSolver
410 * \brief renvoie _ksptype (le type du solveur linéaire utilisé)
413 * \param [out] string
415 string getLinearSolver() {
419 /** \fn getNumberOfVariables
420 * \brief le nombre d'inconnues du problème
423 * \return renvoie _nVar (le nombre d'inconnues du problème)
425 int getNumberOfVariables(){
429 /** \fn setWellBalancedCorrection
430 * \brief include a well balanced correction to treat stiff source terms
431 * @param boolean that is true if a well balanced correction should be applied
433 void setWellBalancedCorrection(bool wellBalancedCorr){
434 _wellBalancedCorrection=wellBalancedCorr;
437 /** \fn setLinearSolver
438 * \brief sets the linear solver and preconditioner
439 * \details virtual function overloaded by intanciable classes
440 * @param kspType linear solver type (GMRES or BICGSTAB)
441 * @param pcType preconditioner (ILU,LU or NONE)
443 void setLinearSolver(linearSolver solverName, preconditioner pcType);
445 /** \fn setNewtonSolver
446 * \brief set the Newton algorithm parameters
447 * \param [in] int maximum number of newton iterations
448 * \param [in] double precision required for the convergence of the newton scheme
451 void setNewtonSolver(double precision,int iterations=20)
453 _maxNewtonIts=iterations;
454 _precision_Newton=precision;
457 /** \fn displayConditionNumber
458 * \brief display the condition number of the preconditioned linear systems
460 void displayConditionNumber(bool display=true){
461 _conditionNumber=display;
464 /** \fn setSaveFileFormat
465 * \brief sets the numerical results file format (MED, VTK or CSV)
467 * \param [in] saveFormat
470 void setSaveFileFormat(saveFormat saveFileFormat){
471 _saveFormat=saveFileFormat;
474 //Couplages Thermohydraulique-thermique-neutronique *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
476 /** \fn setHeatPowerField
477 * \brief set the heat power field (variable in space)
482 void setHeatPowerField(Field heatPower){
483 _heatPowerField=heatPower;
484 _heatPowerFieldSet=true;
487 /** \fn setHeatPowerField
488 * \brief set the heat power field (variable in space)
490 * \param [in] string fileName (including file path)
491 * \param [in] string fieldName
494 void setHeatPowerField(string fileName, string fieldName){
495 _heatPowerField=Field(fileName, CELLS,fieldName);
496 _heatPowerFieldSet=true;
499 /** \fn setHeatSource
500 * \brief met à jour la puissance thermique ( _phi )
505 void setHeatSource(double phi){
509 /** \fn getHeatPowerField
510 * \brief renvoie le champs ?? ( _heatPowerField )
515 Field getHeatPowerField(){
516 return _heatPowerField;
519 /** \fn setRodTemperatureField ??
525 void setRodTemperatureField(Field rodTemperature){
526 _rodTemperatureField=rodTemperature;
527 _rodTemperatureFieldSet=true;
530 /** \fn setRodTemperature ??
536 void setRodTemperature(double rodTemp){
537 _rodTemperature=rodTemp;
540 /** \fn getRodTemperatureField
546 virtual Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
547 return _rodTemperatureField;
550 /** \fn setHeatTransfertCoeff
551 * \brief set the heat transfert coefficient for heat exchange between fluid and solid
556 void setHeatTransfertCoeff(double heatTransfertCoeff){
557 _heatTransfertCoeff=heatTransfertCoeff;
561 * \brief met à jour les paramètres de l'affichage
567 void setVerbose(bool verbose, bool system=false)
574 double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
575 std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
576 std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
577 Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
579 // some supplementary functions
581 /** \fn displayMatrix
582 * \brief displays a matrix of size "size x size" for profiling
583 * @param matrix is a pointer of size "size"
584 * @param size, size of the matrix
585 * @param name, string, name or description of the matrix
586 * @return displays the matrix on the terminal
588 void displayMatrix(double *matrix, int size, string name);
590 /** \fn displayMatrix
591 * \brief displays a vector of size "size" for profiling
592 * @param vector is a pointer of size "size"
593 * @param size, size of the vector
594 * @param name, string, name or description of the vector
595 * @return displays the vector on the terminal
597 void displayVector(double *vector, int size, string name);
599 /** \fn getTimeScheme
600 * \brief returns the time scheme name
602 * \param [out] enum TimeScheme (explicit or implicit)
604 TimeScheme getTimeScheme();
606 /** \fn setNumericalScheme
607 * \brief sets the numerical method ( explicit vs implicit )
609 * \param [in] TimeScheme
612 void setTimeScheme( TimeScheme method);
617 int _Ndim;//space dimension
618 int _nVar;//Number of equations to sole
619 int _Nmailles;//number of cells
620 int _Nnodes;//number of nodes
621 int _Nfaces;//number of faces
622 int _neibMaxNb;//maximum number of neighbours around a cell
623 int _neibMaxNbNodes;/* maximum number of nodes around a node */
633 double _maxvp;//valeur propre max pour calcul cfl
634 double _minl;//minimum cell diameter
636 /** boolean used to specify that a well balanced correction should be used */
637 bool _wellBalancedCorrection;
638 TimeScheme _timeScheme;
640 //Linear solver and petsc
646 int _maxPetscIts;//nombre maximum d'iteration gmres autorise au cours d'une resolution de systeme lineaire
647 int _PetscIts;//the number of iterations of the linear solver
648 int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
650 Mat _A;//Linear system matrix
651 Vec _b;//Linear system right hand side
652 double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
653 bool _conditionNumber;//computes an estimate of the condition number
655 //simulation monitoring variables
657 bool _initialDataSet;
658 bool _initializedMemory;
659 bool _restartWithNewTimeScheme;
660 bool _restartWithNewFileName;
661 double _timeMax,_time;
662 int _maxNbOfTimeStep,_nbTimeStep;
664 double _precision_Newton;
665 double _erreur_rel;//norme(Un+1-Un)
666 string _fileName;//name of the calculation
668 ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
670 //Heat transfert variables
671 Field _heatPowerField, _rodTemperatureField;
672 bool _heatPowerFieldSet, _rodTemperatureFieldSet;
673 double _heatTransfertCoeff;
674 double _heatSource, _rodTemperature;
675 double _hsatv, _hsatl;//all models appart from DiffusionEquation will need this
678 bool _verbose, _system;
680 string _path;//path to execution directory used for saving results
681 saveFormat _saveFormat;//file saving format : MED, VTK or CSV
684 #endif /* PROBLEMCOREFLOWS_HXX_ */