1 //============================================================================
2 // Name : Problème CoreFlows
5 // Copyright : CEA Saclay 2014
6 // Description : Classe générique des problèmes de Thermohydraulique
7 // Il s'agit d'une classe virtuelle pure (non instanciable qui regroupe
8 // les fonctionalités communes à tous les codes thermohydraulique.
9 // Le but est d'avoir une interface commune rendant
10 // plus lisibles et plus evolutifs les codes dévelopés avec CDMATH
11 //============================================================================
13 #include "ProblemCoreFlows.hxx"
14 #include "SparseMatrixPetsc.hxx"
21 ProblemCoreFlows::ProblemCoreFlows(MPI_Comm comm)
23 /* Initialisation of PETSC */
24 //check if PETSC is already initialised
25 PetscBool petscInitialized;
26 PetscInitialized(&petscInitialized);
28 {//check if MPI is already initialised
30 MPI_Initialized(&mpiInitialized);
32 PETSC_COMM_WORLD = comm;
33 PetscInitialize(NULL,NULL,0,0);//Note this is ok if MPI has been been initialised independently from PETSC
35 MPI_Comm_rank(PETSC_COMM_WORLD,&_mpi_rank);
36 MPI_Comm_size(PETSC_COMM_WORLD,&_mpi_size);
37 PetscPrintf(PETSC_COMM_WORLD,"\n Simulation on %d processors\n",_mpi_size);//Prints to standard out, only from the first processor in the communicator. Calls from other processes are ignored.
38 PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Processor [%d] ready for action\n",_mpi_rank);//Prints synchronized output from several processors. Output of the first processor is followed by that of the second, etc.
39 PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
43 PetscPrintf(PETSC_COMM_WORLD,"---- More than one processor detected : running a parallel simulation ----\n");
44 PetscPrintf(PETSC_COMM_WORLD,"---- Limited parallelism : input and output fields remain sequential ----\n");
45 PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrix operations are done in parallel thanks to PETSc----\n");
46 PetscPrintf(PETSC_COMM_WORLD,"---- Processor %d is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n",_mpi_rank);
49 /* Numerical parameter */
57 _precision_Newton=_precision;
60 _stationaryMode=false;
62 _wellBalancedCorrection=false;
64 _nonLinearSolver = Newton_SOLVERLAB;
65 _conditionNumber=false;
67 _runLogFile=new ofstream;
69 /* Monitoring of simulation */
70 _restartWithNewTimeScheme=false;
71 _restartWithNewFileName=false;
72 _fileName = "myCoreFlowsProblem";
83 /* Memory and restart */
84 _initialDataSet=false;
85 _initializedMemory=false;
87 /* PETSc and linear solver parameters */
88 _MaxIterLinearSolver=0;//During several newton iterations, stores the max petssc interations
92 _PetscIts=0;//the number of iterations of the linear solver
93 _ksptype = (char*)&KSPGMRES;
94 _pctype = (char*)&PCILU;
96 /* Physical parameters */
97 _heatPowerFieldSet=false;
98 _heatTransfertCoeff=0;
101 //extracting current directory
102 char result[ PATH_MAX ];
103 getcwd(result, PATH_MAX );
104 _path=string( result );
108 TimeScheme ProblemCoreFlows::getTimeScheme()
113 void ProblemCoreFlows::setTimeScheme(TimeScheme timeScheme)
115 if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
116 _restartWithNewTimeScheme=true;
117 _timeScheme = timeScheme;
120 bool ProblemCoreFlows::isStationary() const
122 return _isStationary;
125 double ProblemCoreFlows::presentTime() const
129 void ProblemCoreFlows::setTimeMax(double timeMax){
132 void ProblemCoreFlows::resetTime (double time)
136 void ProblemCoreFlows::setMaxNbOfTimeStep(int maxNbOfTimeStep){
137 _maxNbOfTimeStep = maxNbOfTimeStep;
139 void ProblemCoreFlows::setCFL(double cfl)
143 void ProblemCoreFlows::setPrecision(double precision)
145 _precision=precision;
147 void ProblemCoreFlows::setInitialField(const Field &VV)
149 if(_Ndim != VV.getSpaceDimension()){
150 *_runLogFile<<"ProblemCoreFlows::setInitialField: mesh has incorrect space dimension"<<endl;
151 _runLogFile->close();
152 throw CdmathException("ProblemCoreFlows::setInitialField: mesh has incorrect space dimension");
154 if(_nVar!=VV.getNumberOfComponents())
156 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect number of components"<<endl;
157 _runLogFile->close();
158 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect number of components");
160 if(_FECalculation && VV.getTypeOfField()!=NODES)
162 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Elements method"<<endl;
163 _runLogFile->close();
164 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Elements method");
166 else if(!_FECalculation && VV.getTypeOfField()==NODES)
168 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on cells or faces for the Finite Volumes method"<<endl;
169 _runLogFile->close();
170 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on cells or faces for the Finite Volumes method");
174 _VV.setName("SOLVERLAB results");
178 _initialDataSet=true;
181 _Nmailles = _mesh.getNumberOfCells();
182 _Nnodes = _mesh.getNumberOfNodes();
183 _Nfaces = _mesh.getNumberOfFaces();
184 _perimeters=Field("Perimeters", CELLS, _mesh,1);
186 // find _minl (delta x) and maximum nb of neibourghs
188 int nbNeib,indexFace;
193 PetscPrintf(PETSC_COMM_SELF,"Processor %d Computing cell perimeters and mesh minimal diameter\n", _mpi_rank);
195 //Compute the maximum number of neighbours for nodes or cells
196 if(VV.getTypeOfField()==NODES)
197 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
200 _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
202 //Compute Delta x and the cell perimeters
203 for (int i=0; i<_mesh.getNumberOfCells(); i++){
204 Ci = _mesh.getCell(i);
207 for (int k=0 ; k<Ci.getNumberOfFaces() ; k++){
208 indexFace=Ci.getFacesId()[k];
209 Fk = _mesh.getFace(indexFace);
210 _minl = min(_minl,Ci.getMeasure()/Fk.getMeasure());
211 _perimeters(i)+=Fk.getMeasure();
214 _minl = min(_minl,Ci.getMeasure());
215 _perimeters(i)=Ci.getNumberOfFaces();
219 cout<<_perimeters<<endl;
222 /*** MPI distribution of parameters ***/
223 MPI_Allreduce(&_initialDataSet, &_initialDataSet, 1, MPIU_BOOL, MPI_LOR, PETSC_COMM_WORLD);
225 int nbVoisinsMax;//Mettre en attribut ?
227 MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD);
228 MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
229 nbVoisinsMax = _neibMaxNbCells;
232 MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD);
233 MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
234 nbVoisinsMax = _neibMaxNbNodes;
236 _d_nnz = (nbVoisinsMax+1)*_nVar;
237 _o_nnz = nbVoisinsMax *_nVar;
240 //Function needed because swig of enum EntityType fails
241 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type)
243 if(_FECalculation && field_support_type!= NODES)
244 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
248 switch(field_support_type)
251 VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
254 VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
257 VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
260 std::ostringstream message;
261 message << "Error ProblemCoreFlows::setInitialField \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
262 throw CdmathException(message.str().c_str());
267 //Function needed because swig of enum EntityType fails
268 void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
269 double ymin, double ymax, int ny, string backSide, string frontSide,
270 double zmin, double zmax, int nz, string bottomSide, string topSide, int field_support_type)
272 if(_FECalculation && field_support_type!= NODES)
273 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
275 EntityType typeField;
277 switch(field_support_type)
289 std::ostringstream message;
290 message << "Error ProblemCoreFlows::setInitialField \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
291 throw CdmathException(message.str().c_str());
294 setInitialFieldConstant( nDim, Vconstant, xmin, xmax, nx, leftSide, rightSide, ymin, ymax, ny, backSide, frontSide, zmin, zmax, nz, bottomSide, topSide, typeField);
296 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, EntityType typeField)
298 if(_FECalculation && typeField!= NODES)
299 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
301 Field VV(fileName, typeField, fieldName, timeStepNumber, order, meshLevel);
305 void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector<double> Vconstant, EntityType typeField)
307 if(_FECalculation && typeField!= NODES)
308 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
310 Field VV("SOLVERLAB results", typeField, M, Vconstant.size());
312 for (int j = 0; j < VV.getNumberOfElements(); j++) {
313 for (int i=0; i< VV.getNumberOfComponents(); i++)
314 VV(j,i) = Vconstant[i];
319 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const Vector Vconstant, EntityType typeField)
321 if(_FECalculation && typeField!= NODES)
322 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
324 Field VV("SOLVERLAB results", typeField, M, Vconstant.getNumberOfRows());
326 for (int j = 0; j < VV.getNumberOfElements(); j++) {
327 for (int i=0; i< VV.getNumberOfComponents(); i++)
328 VV(j,i) = Vconstant(i);
332 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant, EntityType typeField)
334 if(_FECalculation && typeField!= NODES)
335 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
337 Field VV("SOLVERLAB results", typeField, M, Vconstant.size());
339 for (int j = 0; j < VV.getNumberOfElements(); j++) {
340 for (int i=0; i< VV.getNumberOfComponents(); i++)
341 VV(j,i) = Vconstant[i];
345 void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
346 double ymin, double ymax, int ny, string backSide, string frontSide,
347 double zmin, double zmax, int nz, string bottomSide, string topSide, EntityType typeField)
351 M=Mesh(xmin,xmax,nx);
353 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
355 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
357 PetscPrintf(PETSC_COMM_WORLD,"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3\n");
358 *_runLogFile<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
359 _runLogFile->close();
360 throw CdmathException("Space dimension nDim should be between 1 and 3");
363 M.setGroupAtPlan(xmax,0,_precision,rightSide);
364 M.setGroupAtPlan(xmin,0,_precision,leftSide);
366 M.setGroupAtPlan(ymax,1,_precision,frontSide);
367 M.setGroupAtPlan(ymin,1,_precision,backSide);
370 M.setGroupAtPlan(zmax,2,_precision,topSide);
371 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
374 setInitialFieldConstant(M, Vconstant, typeField);
376 void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, int direction, EntityType typeField)
378 if(_FECalculation && typeField!= NODES)
379 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
381 if (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
383 *_runLogFile<<"ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes"<<endl;
384 _runLogFile->close();
385 throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes");
387 Field VV("SOLVERLAB results", typeField, M, VV_Left.getNumberOfRows());
389 double component_value;
391 for (int j = 0; j < VV.getNumberOfElements(); j++)
394 component_value=VV.getElementComponent(j, direction);
397 PetscPrintf(PETSC_COMM_WORLD,"Error : space dimension is %d, direction asked is \%d \n",M.getSpaceDimension(),direction);
398 _runLogFile->close();
399 throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: direction should be an integer between 0 and 2");
402 for (int i=0; i< VV.getNumberOfComponents(); i++)
403 if (component_value< disc_pos )
404 VV(j,i) = VV_Left[i];
406 VV(j,i) = VV_Right[i];
410 void ProblemCoreFlows::setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
411 double xmin, double xmax, int nx, string leftSide, string rightSide,
412 double ymin, double ymax, int ny, string backSide, string frontSide,
413 double zmin, double zmax, int nz, string bottomSide, string topSide, EntityType typeField)
417 M=Mesh(xmin,xmax,nx);
419 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
421 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
424 _runLogFile->close();
425 throw CdmathException("Space dimension nDim should be between 1 and 3");
428 M.setGroupAtPlan(xmax,0,_precision,rightSide);
429 M.setGroupAtPlan(xmin,0,_precision,leftSide);
431 M.setGroupAtPlan(ymax,1,_precision,frontSide);
432 M.setGroupAtPlan(ymin,1,_precision,backSide);
435 M.setGroupAtPlan(zmax,2,_precision,topSide);
436 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
438 Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
439 for(int i=0;i<VV_Left.size(); i++){
440 V_Left(i)=VV_Left[i];
441 V_Right(i)=VV_Right[i];
443 setInitialFieldStepFunction(M, V_Left, V_Right, xstep, typeField);
446 void ProblemCoreFlows::setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double radius, const Vector Center, EntityType typeField)
448 if(_FECalculation && typeField!= NODES)
449 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
451 if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
453 PetscPrintf(PETSC_COMM_WORLD,"Vout.size() = %d, Vin.size()= %d, Center.size() = %d, M.getSpaceDim = %d \n",Vout.size(),Vin.size(),Center.size(), M.getSpaceDimension());
454 throw CdmathException("ProblemCoreFlows::setInitialFieldSphericalStepFunction : Vector size error");
456 int nVar=Vout.size();
457 int spaceDim=M.getSpaceDimension();
458 Field VV("Primitive variables for spherical step function", typeField, M, nVar);
460 Vector currentPoint(spaceDim);
461 for(int i=0;i<VV.getNumberOfElements();i++)
463 currentPoint(0)=VV.getElementComponent(i,0);
466 currentPoint(1)=VV.getElementComponent(i,1);
468 currentPoint(2)=VV.getElementComponent(i,2);
470 if((currentPoint-Center).norm()<radius)
471 for(int j=0;j<nVar;j++)
474 for(int j=0;j<nVar;j++)
480 double ProblemCoreFlows::getTime()
484 unsigned ProblemCoreFlows::getNbTimeStep()
488 double ProblemCoreFlows::getCFL()
492 double ProblemCoreFlows::getPrecision()
496 Mesh ProblemCoreFlows::getMesh()
501 void ProblemCoreFlows::setLinearSolver(linearSolver kspType, preconditioner pcType, double maxIts)
504 // set linear solver algorithm
506 _ksptype = (char*)&KSPGMRES;
507 else if (kspType==CG)
508 _ksptype = (char*)&KSPCG;
509 else if (kspType==BCGS)
510 _ksptype = (char*)&KSPBCGS;
512 PetscPrintf(PETSC_COMM_WORLD,"!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!\n");
513 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
514 _runLogFile->close();
515 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
517 // set preconditioner
519 _pctype = (char*)&PCNONE;
520 else if (pcType ==LU)
521 _pctype = (char*)&PCLU;
522 else if (pcType == ILU)
523 _pctype = (char*)&PCILU;
524 else if (pcType ==CHOLESKY)
525 _pctype = (char*)&PCCHOLESKY;
526 else if (pcType == ICC)
527 _pctype = (char*)&PCICC;
529 PetscPrintf(PETSC_COMM_WORLD,"!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!\n");
530 *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
531 _runLogFile->close();
532 throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
536 void ProblemCoreFlows::setNewtonSolver(double precision, int iterations, nonLinearSolver solverName)
538 _maxNewtonIts=iterations;
539 _precision_Newton=precision;
540 _nonLinearSolver=solverName;
545 // Cette methode lance une execution du ProblemCoreFlows
546 // Elle peut etre utilisee si le probleme n'est couple a aucun autre.
547 // (s'il n'a besoin d'aucun champ d'entree).
548 // Precondition: initialize
549 // Seule la methode terminate peut etre appelée apres
550 bool ProblemCoreFlows::run()
552 if(!_initializedMemory)
554 _runLogFile->close();
555 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
557 bool stop=false; // Does the Problem want to stop (error) ?
558 bool ok; // Is the time interval successfully solved ?
559 _isStationary=false;//in case of a second run with a different physics or cfl
561 PetscPrintf(PETSC_COMM_WORLD,"Running test case %s\n",_fileName.c_str());
563 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
564 *_runLogFile<< "Running test case "<< _fileName<<endl;
567 while(!stop && !_isStationary &&_time<_timeMax && _nbTimeStep<_maxNbOfTimeStep)
569 ok=false; // Is the time interval successfully solved ?
571 // Guess the next time step length
572 _dt=computeTimeStep(stop);
574 PetscPrintf(PETSC_COMM_WORLD,"Failed computing time step %d, time = %.2f, dt= %.2e, stopping calculation",_nbTimeStep,_time,_dt);
575 *_runLogFile << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
578 // Loop on the time interval tries
579 while (!ok && !stop )
581 stop=!initTimeStep(_dt);
582 // Prepare the next time step
584 PetscPrintf(PETSC_COMM_WORLD,"Failed initializing time step %d, time = %.2f, dt= %.2e, stopping calculation",_nbTimeStep,_time,_dt);
585 *_runLogFile << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
588 // Solve the next time step
591 if (!ok) // The resolution failed, try with a new time interval.
593 /* if(_dt>_precision){
594 cout<<"ComputeTimeStep returned _dt="<<_dt<<endl;
595 cout << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with dt/2"<< endl;
596 *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with dt/2"<< endl;
597 double dt=_dt/2;//We chose to divide the time step by 2
598 abortTimeStep();//Cancel the initTimeStep
599 _dt=dt;//new value of time step is previous time step divided by 2 (we do not call computeTimeStep
600 //_cfl*=0.5;//If we change the cfl, we must compute the new time step with computeTimeStep
601 //_dt=computeTimeStep(stop);
604 PetscPrintf(PETSC_COMM_WORLD,"Failed solving time step %d, time = %.2f, dt= %.2e, cfl = %.2f, stopping calculation \n",_nbTimeStep,_time,_dt,_cfl);
605 *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
606 stop=true; // Impossible to solve the next time step, the Problem has given up
610 else // The resolution was successful, validate and go to the next time step.
613 if ((_nbTimeStep-1)%_freqSave ==0){
614 PetscPrintf(PETSC_COMM_WORLD,"Time step = %d, dt = %.2e, time = %.2f, ||Un+1-Un||= %.2e\n\n",_nbTimeStep,_dt,_time,_erreur_rel);
615 *_runLogFile << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl<<endl;
621 PetscPrintf(PETSC_COMM_WORLD,"Stationary state reached\n");
622 *_runLogFile << "Stationary state reached" <<endl;
624 else if(_time>=_timeMax){
625 PetscPrintf(PETSC_COMM_WORLD,"Maximum time %.2f reached\n",_timeMax);
626 *_runLogFile<<"Maximum time "<<_timeMax<<" reached"<<endl;
628 else if(_nbTimeStep>=_maxNbOfTimeStep){
629 PetscPrintf(PETSC_COMM_WORLD,"Maximum number of time steps %d reached\n",_maxNbOfTimeStep);
630 *_runLogFile<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
633 PetscPrintf(PETSC_COMM_WORLD,"Error problem wants to stop!\n");
634 *_runLogFile<<"Error problem wants to stop!"<<endl;
636 PetscPrintf(PETSC_COMM_WORLD,"End of calculation at time t = %.2f and time step number %d\n",_time,_nbTimeStep);
637 *_runLogFile << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
639 _runLogFile->close();
643 void ProblemCoreFlows::displayMatrix(double *matrix, int size, string name)
646 for(int p=0; p<size; p++)
648 for(int q=0; q<size; q++)
650 cout << matrix[p*size+q] << "\t";
657 void ProblemCoreFlows::displayVector(double *vector, int size, string name)
660 for(int q=0; q<size; q++)
662 cout << vector[q] << "\t";
666 void ProblemCoreFlows::setFileName(string fileName){
667 if( _nbTimeStep>0 && fileName!=_fileName)//This is a change of file name during a simulation
668 _restartWithNewFileName=true;
669 _fileName = fileName;
672 void ProblemCoreFlows::setFreqSave(int freqSave){
673 _freqSave = freqSave;
676 bool ProblemCoreFlows::solveTimeStep(){
678 bool converged=false, ok=true;
679 while(!converged && ok && _NEWTON_its < _maxNewtonIts){
680 ok=iterateTimeStep(converged);//resolution du systeme lineaire si schema implicite
682 if(_timeScheme == Implicit && (_nbTimeStep-1)%_freqSave ==0)//To monitor the convergence of the newton scheme
684 PetscPrintf(PETSC_COMM_WORLD," Newton iteration %d, %s iterations : %d maximum variation ||Uk+1-Uk||: %.2e\n",_NEWTON_its,_ksptype,_PetscIts,_erreur_rel);
685 *_runLogFile<< " Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
689 PetscReal sv_max, sv_min;
690 KSPComputeExtremeSingularValues(_ksp, &sv_max, &sv_min);
691 PetscPrintf(PETSC_COMM_WORLD," Singular value max = %.2e, singular value min = %.2e, condition number = %.2e\n",sv_max,sv_min,sv_max/sv_min);
692 *_runLogFile<<" Singular value max = " << sv_max <<", singular value min = " << sv_min <<", condition number = " << sv_max/sv_min <<endl;
698 if(_NEWTON_its >= _maxNewtonIts){
699 PetscPrintf(PETSC_COMM_WORLD,"Maximum number of Newton iterations %d reached\n",_maxNewtonIts);
700 *_runLogFile << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
703 PetscPrintf(PETSC_COMM_WORLD,"iterateTimeStep: solving Newton iteration %d Failed\n",_NEWTON_its);
704 *_runLogFile<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
707 else if(_timeScheme == Implicit && (_nbTimeStep-1)%_freqSave ==0)
709 PetscPrintf(PETSC_COMM_WORLD,"Nombre d'iterations de Newton %d, Nombre max d'iterations %s : %d\n\n",_NEWTON_its, _ksptype, _MaxIterLinearSolver);
711 *_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its << "Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl << endl;
712 _MaxIterLinearSolver = 0;
718 ProblemCoreFlows::~ProblemCoreFlows()
720 PetscBool petscInitialized;
721 PetscInitialized(&petscInitialized);
729 ProblemCoreFlows::getConditionNumber(bool isSingular, double tol) const
731 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
732 return A.getConditionNumber( isSingular, tol);
734 std::vector< double >
735 ProblemCoreFlows::getEigenvalues(int nev, EPSWhich which, double tol) const
737 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
738 return A.getEigenvalues( nev, which, tol);
740 std::vector< Vector >
741 ProblemCoreFlows::getEigenvectors(int nev, EPSWhich which, double tol) const
743 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
744 return A.getEigenvectors( nev, which, tol);
747 ProblemCoreFlows::getEigenvectorsField(int nev, EPSWhich which, double tol) const
749 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
750 MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
753 my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
755 my_eigenfield.setFieldByDataArrayDouble(d);
757 return my_eigenfield;
760 std::vector< double >
761 ProblemCoreFlows::getSingularValues( int nsv, SVDWhich which, double tol) const
763 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
764 return A.getSingularValues( nsv, which, tol);
766 std::vector< Vector >
767 ProblemCoreFlows::getSingularVectors(int nsv, SVDWhich which, double tol) const
769 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
770 return A.getSingularVectors( nsv, which, tol);
774 ProblemCoreFlows::getUnknownField() const
776 if(!_initializedMemory)
778 _runLogFile->close();
779 throw CdmathException("ProblemCoreFlows::getUnknownField() call initialize() first");
785 ProblemCoreFlows::isMEDCoupling64Bits() const
787 #ifdef MEDCOUPLING_USE_64BIT_IDS
795 ProblemCoreFlows::setHeatPowerField(Field heatPower){
797 throw CdmathException("!!!!!!!! ProblemCoreFlows::setHeatPowerField set initial field first");
799 heatPower.getMesh().checkFastEquivalWith(_mesh);
800 _heatPowerField=heatPower;
801 _heatPowerFieldSet=true;
805 ProblemCoreFlows::setHeatPowerField(string fileName, string fieldName, int iteration, int order, int meshLevel){
807 throw CdmathException("!!!!!!!! ProblemCoreFlows::setHeatPowerField set initial field first");
809 _heatPowerField=Field(fileName, CELLS,fieldName, iteration, order, meshLevel);
810 _heatPowerField.getMesh().checkFastEquivalWith(_mesh);
811 _heatPowerFieldSet=true;