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 nonLinearSolver
58 /*! the nonLinearSolver can be Newton_SOLVERLAB or using PETSc, Newton_PETSC_LINESEARCH, Newton_PETSC_TRUSTREGION (see Petsc documentation) */
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 */
72 //! enumeration saveFormat
73 /*! the numerical results are saved using MED, VTK or CSV format */
76 MED,/**< MED format is used */
77 VTK,/**< VTK format is used */
78 CSV/**< CSV format is used */
81 //! enumeration TimeScheme
82 /*! The numerical method can be Explicit or Implicit */
85 Explicit,/**< Explicit numerical scheme */
86 Implicit/**< Implicit numerical scheme */
89 class ProblemCoreFlows
92 //! Constructeur par défaut
94 virtual ~ProblemCoreFlows();
96 // -*-*-*- Gestion du calcul (interface ICoCo) -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
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 .
103 virtual void initialize()=0;
106 * \brief vide la mémoire et enregistre le résultat final
107 * \Details on la surcharge dans le problem fluid
110 virtual void terminate()=0;
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
121 virtual double computeTimeStep(bool & stop)=0;
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
130 virtual bool initTimeStep(double dt)=0;
132 /** \fn solveTimeStep
133 * \brief calcule les valeurs inconnues au pas de temps +1 .
134 * \details c'est une fonction virtuelle
136 * \return Renvoie false en cas de problème durant le calcul (valeurs non physiques..)
138 virtual bool solveTimeStep();//
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
146 * \return Renvoie false en cas de problème durant le calcul
148 virtual void validateTimeStep()=0;//
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
154 virtual void abortTimeStep()=0;
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
161 * \return Renvoie "false" en cas de problème durant le calcul
163 virtual bool run();//
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.
173 virtual bool iterateTimeStep(bool &converged) = 0; //??
176 * \brief vérifie la stationnairité du problème .
177 * \details Renvoie "true" si le problème atteint l'état stationnaire et "false" sinon
181 virtual bool isStationary() const;
184 * \brief Calcule la valeur du temps courant .
187 * \param [out] double
189 virtual double presentTime() const;
193 virtual vector<string> getInputFieldsNames()=0 ;//Renvoie les noms des champs dont le problème a besoin (données initiales)
194 virtual Field& getInputFieldTemplate(const string& name)=0;//Renvoie le format de champs attendu (maillage, composantes etc)
195 virtual void setInputField(const string& name, const Field& afield)=0;//enregistre les valeurs d'une donnée initiale
196 virtual vector<string> getOutputFieldsNames()=0 ;//liste tous les champs que peut fournir le code pour le postraitement
197 virtual Field& getOutputField(const string& nameField )=0;//Renvoie un champs pour le postraitement
200 Field getUnknownField() const;
202 //paramètres du calcul -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
204 /** \fn setPresentTime
205 * \brief met à jour _time (le temps courant du calcul)
210 void setPresentTime (double time);
212 /** \fn setMaxNbOfTimeStep
213 * \brief met à jour _maxNbOfTimeStep ( le nombre maximum d'itération du calcul )
218 void setMaxNbOfTimeStep(int maxNbOfTimeStep);
221 * \brief met à jour _timeMax (Le temps maximum du calcul)
226 void setTimeMax(double timeMax);
229 * \brief met à jour la _CFL
234 void setCFL(double cfl);
237 * \brief met à jour _precision (la précision du calcule)
242 void setPrecision(double precision);
244 /** \fn setInitialField
245 * \brief sets the initial field
250 void setInitialField(const Field &VV);
252 /** \fn setInitialField
253 * \brief sets the initial field from a field in a med file
255 * \param [in] string : the file name
256 * \param [in] string : the field name
257 * \param [in] int : the time step number
260 void setInitialField(string fileName, string fieldName, int timeStepNumber);
262 /** \fn setInitialFieldConstant
263 * \brief sets a constant initial field on a mesh stored in a med file
265 * \param [in] string : the file name
266 * \param [in] vector<double> : the value in each cell
269 void setInitialFieldConstant(string fileName, const vector<double> Vconstant);
271 /** \fn setInitialFieldConstant
272 * \brief sets a constant initial field
278 void setInitialFieldConstant(const Mesh& M, const Vector Vconstant);
280 /** \fn setInitialFieldConstant
281 * \brief sets a constant initial field
284 * \param [in] vector<double>
287 void setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant);
289 /** \fn setInitialFieldConstant
290 * \brief sets a constant initial field
292 * \param [in] int the space dimension
293 * \param [in] vector<double> the value in each cell
294 * \param [in] double the lowest value in the x direction
295 * \param [in] double the highest value in the x direction
296 * \param [in] string name of the left boundary
297 * \param [in] string name of the right boundary
298 * \param [in] double the lowest value in the y direction
299 * \param [in] double the highest value in the y direction
300 * \param [in] string name of the back boundary
301 * \param [in] string name of the front boundary
302 * \param [in] double the lowest value in the z direction
303 * \param [in] double the highest value in the z direction
304 * \param [in] string name of the bottom boundary
305 * \param [in] string name of the top boundary
308 void setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax,int nx, string leftSide, string rightSide,
309 double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
310 double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
312 /** \fn setInitialFieldStepFunction
313 * \brief sets a step function initial field (Riemann problem)
318 * \param [in] double position of the discontinuity on one of the three axis
319 * \param [in] int direction (axis carrying the discontinuity) : 0 for x, 1 for y, 2 for z
322 void setInitialFieldStepFunction(const Mesh M, const Vector Vleft, const Vector Vright, double disc_pos, int direction=0);
324 /** \fn setInitialFieldStepFunction
325 * \brief sets a constant initial field
327 * \param [in] int the space dimension
328 * \param [in] vector<double> the value left of the discontinuity
329 * \param [in] vector<double> the value right of the discontinuity
330 * \param [in] double the position of the discontinuity in the x direction
331 * \param [in] double the lowest value in the x direction
332 * \param [in] double the highest value in the x direction
333 * \param [in] string name of the left boundary
334 * \param [in] string name of the right boundary
335 * \param [in] double the lowest value in the y direction
336 * \param [in] double the highest value in the y direction
337 * \param [in] string name of the back boundary
338 * \param [in] string name of the front boundary
339 * \param [in] double the lowest value in the z direction
340 * \param [in] double the highest value in the z direction
341 * \param [in] string name of the bottom boundary
342 * \param [in] string name of the top boundary
345 void setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
346 double xmin, double xmax,int nx, string leftSide, string rightSide,
347 double ymin=0, double ymax=0, int ny=0, string backSide="", string frontSide="",
348 double zmin=0, double zmax=0, int nz=0, string bottomSide="", string topSide="");
350 /** \fn setInitialFieldSphericalStepFunction
351 * \brief sets a step function initial field with value Vin inside the ball with radius Radius and Vout outside
354 * \param [in] Vector Vin, value inside the ball
355 * \param [in] Vector Vout, value outside the ball
356 * \param [in] double radius of the ball
357 * \param [in] Vector Center, coordinates of the ball center
360 void setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double Radius, Vector Center);
363 * \brief renvoie _time (le temps courant du calcul)
366 * \param [out] double
370 /** \fn getNbTimeStep
371 * \brief renvoie _nbTimeStep le Numéro d'itération courant
374 * \param [out] unsigned
376 unsigned getNbTimeStep();
379 * \brief renvoie la _CFL
382 * \param [out] double
387 * \brief renvoie _precision (la précision du calcul)
390 * \param [out] double
392 double getPrecision();
395 * \brief renvoie _Mesh (le maillage du problème)
403 * \brief met à jour _fileName le nom du fichier
408 void setFileName(string fileName);
411 * \brief met à jour _FreqSave (la fréquence du sauvgarde de la solution)
416 void setFreqSave(int freqSave);
419 * \brief sauvgarde les données dans des fichiers MED ou VTK
420 * \details c'est une fonction virtuelle pure , on la surcharge
421 * dans chacun des modèles
424 virtual void save() = 0;
426 /** \fn getLinearSolver
427 * \brief renvoie _ksptype (le type du solveur linéaire utilisé)
430 * \param [out] string
432 string getLinearSolver() {
436 /** \fn getNonLinearSolver
437 * \brief renvoie _nonLinearSolver (le type du solveur de Newton utilisé)
440 * \param [out] string
442 nonLinearSolver getNonLinearSolver() {
443 return _nonLinearSolver;
446 /** \fn getNumberOfVariables
447 * \brief le nombre d'inconnues du problème
450 * \return renvoie _nVar (le nombre d'inconnues du problème)
452 int getNumberOfVariables(){
456 /** \fn setWellBalancedCorrection
457 * \brief include a well balanced correction to treat stiff source terms
458 * @param boolean that is true if a well balanced correction should be applied
460 void setWellBalancedCorrection(bool wellBalancedCorr){
461 _wellBalancedCorrection=wellBalancedCorr;
464 /** \fn setLinearSolver
465 * \brief sets the linear solver and preconditioner
466 * \details virtual function overloaded by intanciable classes
467 * @param kspType linear solver type (GMRES or BICGSTAB)
468 * @param pcType preconditioner (ILU,LU or NONE)
470 void setLinearSolver(linearSolver solverName, preconditioner pcType, double maxIts=50);
472 /** \fn setNewtonSolver
473 * \brief sets the Newton type algorithm for solving the nonlinear algebraic system arising from the discretisation of the PDE
474 * \param [in] double : precision required for the convergence of the newton scheme
475 * \param [in] int : maximum number of newton iterations
476 * \param [in] nonLinearSolver : the algorithm to be used to solve the nonlinear system
479 void setNewtonSolver(double precision, int iterations=20, nonLinearSolver solverName=Newton_SOLVERLAB);
481 /** \fn displayConditionNumber
482 * \brief display the condition number of the preconditioned linear systems
484 void displayConditionNumber(bool display=true){
485 _conditionNumber=display;
488 /** \fn setSaveFileFormat
489 * \brief sets the numerical results file format (MED, VTK or CSV)
491 * \param [in] saveFormat
494 void setSaveFileFormat(saveFormat saveFileFormat){
495 _saveFormat=saveFileFormat;
498 /** \fn setResultDirectory
499 * \brief sets the directory where the results will be saved
501 * \param [in] resultsPath
504 void setResultDirectory(string resultsPath){
508 /** \fn getTimeScheme
509 * \brief returns the time scheme name
511 * \param [out] enum TimeScheme (explicit or implicit)
513 TimeScheme getTimeScheme();
515 /** \fn setNumericalScheme
516 * \brief sets the numerical method ( explicit vs implicit )
518 * \param [in] TimeScheme
521 void setTimeScheme( TimeScheme method);
523 //Couplages Thermohydraulique-thermique-neutronique *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
525 /** \fn setHeatPowerField
526 * \brief set the heat power field (variable in space)
531 void setHeatPowerField(Field heatPower){
532 _heatPowerField=heatPower;
533 _heatPowerFieldSet=true;
536 /** \fn setHeatPowerField
537 * \brief set the heat power field (variable in space)
539 * \param [in] string fileName (including file path)
540 * \param [in] string fieldName
543 void setHeatPowerField(string fileName, string fieldName){
544 _heatPowerField=Field(fileName, CELLS,fieldName);
545 _heatPowerFieldSet=true;
548 /** \fn setHeatSource
549 * \brief met à jour la puissance thermique ( _phi )
554 void setHeatSource(double phi){
558 /** \fn getHeatPowerField
559 * \brief renvoie le champs ?? ( _heatPowerField )
564 Field getHeatPowerField(){
565 return _heatPowerField;
568 /** \fn setRodTemperatureField ??
574 void setRodTemperatureField(Field rodTemperature){
575 _rodTemperatureField=rodTemperature;
576 _rodTemperatureFieldSet=true;
579 /** \fn setRodTemperature ??
585 void setRodTemperature(double rodTemp){
586 _rodTemperature=rodTemp;
589 /** \fn getRodTemperatureField
595 virtual Field& getRodTemperatureField(){ // ?? je ne retrouve pas cet attribut dans le file.cxx
596 return _rodTemperatureField;
599 /** \fn setHeatTransfertCoeff
600 * \brief set the heat transfert coefficient for heat exchange between fluid and solid
605 void setHeatTransfertCoeff(double heatTransfertCoeff){
606 _heatTransfertCoeff=heatTransfertCoeff;
610 * \brief met à jour les paramètres de l'affichage
616 void setVerbose(bool verbose, bool system=false)
623 double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
624 std::vector< double > getEigenvalues (int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
625 std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
626 Field getEigenvectorsField(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
627 std::vector< double > getSingularValues( int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
628 std::vector< Vector > getSingularVectors(int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
630 // some supplementary functions
632 /** \fn displayMatrix
633 * \brief displays a matrix of size "size x size" for profiling
634 * @param matrix is a pointer of size "size"
635 * @param size, size of the matrix
636 * @param name, string, name or description of the matrix
637 * @return displays the matrix on the terminal
639 static void displayMatrix(double *matrix, int size, string name="Vector coefficients :");
641 /** \fn displayMatrix
642 * \brief displays a vector of size "size" for profiling
643 * @param vector is a pointer of size "size"
644 * @param size, size of the vector
645 * @param name, string, name or description of the vector
646 * @return displays the vector on the terminal
648 static void displayVector(double *vector, int size, string name="Matrix coefficients :");
652 int _Ndim;//space dimension
653 int _nVar;//Number of equations to sole
654 int _Nmailles;//number of cells
655 int _Nnodes;//number of nodes
656 int _Nfaces;//number of faces
657 int _neibMaxNb;//maximum number of neighbours around a cell
658 int _neibMaxNbNodes;/* maximum number of nodes around a node */
668 double _maxvp;//valeur propre max pour calcul cfl
669 double _minl;//minimum cell diameter
671 /** boolean used to specify that a well balanced correction should be used */
672 bool _wellBalancedCorrection;
673 TimeScheme _timeScheme;
675 //Linear solver and petsc
681 nonLinearSolver _nonLinearSolver;
682 int _maxPetscIts;//nombre maximum d'iteration gmres autorise au cours d'une resolution de systeme lineaire
683 int _PetscIts;//the number of iterations of the linear solver
684 int _maxNewtonIts;//nombre maximum d'iteration de Newton autorise au cours de la resolution d'un pas de temps
686 Mat _A;//Linear system matrix
687 Vec _b;//Linear system right hand side
688 double _MaxIterLinearSolver;//nombre maximum d'iteration gmres obtenu au cours par les resolution de systemes lineaires au cours d'un pas de tmeps
689 bool _conditionNumber;//computes an estimate of the condition number
691 //simulation monitoring variables
693 bool _initialDataSet;
694 bool _initializedMemory;
695 bool _restartWithNewTimeScheme;
696 bool _restartWithNewFileName;
697 double _timeMax,_time;
698 int _maxNbOfTimeStep,_nbTimeStep;
700 double _precision_Newton;
701 double _erreur_rel;//norme(Un+1-Un)
702 string _fileName;//name of the calculation
704 ofstream * _runLogFile;//for creation of a log file to save the history of the simulation
706 //Heat transfert variables
707 Field _heatPowerField, _rodTemperatureField;
708 bool _heatPowerFieldSet, _rodTemperatureFieldSet;
709 double _heatTransfertCoeff;
710 double _heatSource, _rodTemperature;
711 double _hsatv, _hsatl;//all models appart from DiffusionEquation will need this
714 bool _verbose, _system;
716 string _path;//path to execution directory used for saving results
717 saveFormat _saveFormat;//file saving format : MED, VTK or CSV
721 #endif /* PROBLEMCOREFLOWS_HXX_ */