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;
95 _pctype = (char*)&PCNONE;
97 _pctype = (char*)&PCLU;
99 /* Physical parameters */
100 _heatPowerFieldSet=false;
101 _heatTransfertCoeff=0;
104 //extracting current directory
105 char result[ PATH_MAX ];
106 getcwd(result, PATH_MAX );
107 _path=string( result );
111 TimeScheme ProblemCoreFlows::getTimeScheme()
116 void ProblemCoreFlows::setTimeScheme(TimeScheme timeScheme)
118 if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
119 _restartWithNewTimeScheme=true;
120 _timeScheme = timeScheme;
123 bool ProblemCoreFlows::isStationary() const
125 return _isStationary;
128 double ProblemCoreFlows::presentTime() const
132 void ProblemCoreFlows::setTimeMax(double timeMax){
135 void ProblemCoreFlows::resetTime (double time)
139 void ProblemCoreFlows::setMaxNbOfTimeStep(int maxNbOfTimeStep){
140 _maxNbOfTimeStep = maxNbOfTimeStep;
142 void ProblemCoreFlows::setCFL(double cfl)
146 void ProblemCoreFlows::setPrecision(double precision)
148 _precision=precision;
150 void ProblemCoreFlows::setInitialField(const Field &VV)
152 if(_Ndim != VV.getSpaceDimension()){
153 *_runLogFile<<"ProblemCoreFlows::setInitialField: mesh has incorrect space dimension"<<endl;
154 _runLogFile->close();
155 throw CdmathException("ProblemCoreFlows::setInitialField: mesh has incorrect space dimension");
157 if(_nVar!=VV.getNumberOfComponents())
159 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect number of components"<<endl;
160 _runLogFile->close();
161 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect number of components");
163 if(_FECalculation && VV.getTypeOfField()!=NODES)
165 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Elements method"<<endl;
166 _runLogFile->close();
167 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Elements method");
169 else if(!_FECalculation && VV.getTypeOfField()==NODES)
171 *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on cells or faces for the Finite Volumes method"<<endl;
172 _runLogFile->close();
173 throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on cells or faces for the Finite Volumes method");
177 _VV.setName("SOLVERLAB results");
181 _initialDataSet=true;
184 _Nmailles = _mesh.getNumberOfCells();
185 _Nnodes = _mesh.getNumberOfNodes();
186 _Nfaces = _mesh.getNumberOfFaces();
187 _perimeters=Field("Perimeters", CELLS, _mesh,1);
189 // find _minl (delta x) and maximum nb of neibourghs
191 int nbNeib,indexFace;
196 PetscPrintf(PETSC_COMM_SELF,"Processor %d Computing cell perimeters and mesh minimal diameter\n", _mpi_rank);
198 //Compute the maximum number of neighbours for nodes or cells
199 if(VV.getTypeOfField()==NODES)
200 _neibMaxNbNodes=_mesh.getMaxNbNeighbours(NODES);
203 _neibMaxNbCells=_mesh.getMaxNbNeighbours(CELLS);
205 //Compute Delta x and the cell perimeters
206 for (int i=0; i<_mesh.getNumberOfCells(); i++){
207 Ci = _mesh.getCell(i);
210 for (int k=0 ; k<Ci.getNumberOfFaces() ; k++){
211 indexFace=Ci.getFacesId()[k];
212 Fk = _mesh.getFace(indexFace);
213 _minl = min(_minl,Ci.getMeasure()/Fk.getMeasure());
214 _perimeters(i)+=Fk.getMeasure();
217 _minl = min(_minl,Ci.getMeasure());
218 _perimeters(i)=Ci.getNumberOfFaces();
222 cout<<_perimeters<<endl;
225 /*** MPI distribution of parameters ***/
226 MPI_Allreduce(&_initialDataSet, &_initialDataSet, 1, MPIU_BOOL, MPI_LOR, PETSC_COMM_WORLD);
228 int nbVoisinsMax;//Mettre en attribut ?
230 MPI_Bcast(&_Nmailles , 1, MPI_INT, 0, PETSC_COMM_WORLD);
231 MPI_Bcast(&_neibMaxNbCells, 1, MPI_INT, 0, PETSC_COMM_WORLD);
232 nbVoisinsMax = _neibMaxNbCells;
235 MPI_Bcast(&_Nnodes , 1, MPI_INT, 0, PETSC_COMM_WORLD);
236 MPI_Bcast(&_neibMaxNbNodes, 1, MPI_INT, 0, PETSC_COMM_WORLD);
237 nbVoisinsMax = _neibMaxNbNodes;
239 _d_nnz = (nbVoisinsMax+1)*_nVar;
240 _o_nnz = nbVoisinsMax *_nVar;
243 //Function needed because swig of enum EntityType fails
244 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, int field_support_type)
246 if(_FECalculation && field_support_type!= NODES)
247 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
251 switch(field_support_type)
254 VV = Field(fileName, CELLS, fieldName, timeStepNumber, order, meshLevel);
257 VV = Field(fileName, NODES, fieldName, timeStepNumber, order, meshLevel);
260 VV = Field(fileName, FACES, fieldName, timeStepNumber, order, meshLevel);
263 std::ostringstream message;
264 message << "Error ProblemCoreFlows::setInitialField \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
265 throw CdmathException(message.str().c_str());
270 //Function needed because swig of enum EntityType fails
271 void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
272 double ymin, double ymax, int ny, string backSide, string frontSide,
273 double zmin, double zmax, int nz, string bottomSide, string topSide, int field_support_type)
275 if(_FECalculation && field_support_type!= NODES)
276 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
278 EntityType typeField;
280 switch(field_support_type)
292 std::ostringstream message;
293 message << "Error ProblemCoreFlows::setInitialField \n Accepted field support integers are "<< CELLS <<" (for CELLS), "<<NODES<<" (for NODES), and "<< FACES <<" (for FACES)" ;
294 throw CdmathException(message.str().c_str());
297 setInitialFieldConstant( nDim, Vconstant, xmin, xmax, nx, leftSide, rightSide, ymin, ymax, ny, backSide, frontSide, zmin, zmax, nz, bottomSide, topSide, typeField);
299 void ProblemCoreFlows::setInitialField(string fileName, string fieldName, int timeStepNumber, int order, int meshLevel, EntityType typeField)
301 if(_FECalculation && typeField!= NODES)
302 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
304 Field VV(fileName, typeField, fieldName, timeStepNumber, order, meshLevel);
308 void ProblemCoreFlows::setInitialFieldConstant(string fileName, const vector<double> Vconstant, EntityType typeField)
310 if(_FECalculation && typeField!= NODES)
311 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
313 Field VV("SOLVERLAB results", typeField, M, Vconstant.size());
315 for (int j = 0; j < VV.getNumberOfElements(); j++) {
316 for (int i=0; i< VV.getNumberOfComponents(); i++)
317 VV(j,i) = Vconstant[i];
322 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const Vector Vconstant, EntityType typeField)
324 if(_FECalculation && typeField!= NODES)
325 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
327 Field VV("SOLVERLAB results", typeField, M, Vconstant.getNumberOfRows());
329 for (int j = 0; j < VV.getNumberOfElements(); j++) {
330 for (int i=0; i< VV.getNumberOfComponents(); i++)
331 VV(j,i) = Vconstant(i);
335 void ProblemCoreFlows:: setInitialFieldConstant(const Mesh& M, const vector<double> Vconstant, EntityType typeField)
337 if(_FECalculation && typeField!= NODES)
338 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
340 Field VV("SOLVERLAB results", typeField, M, Vconstant.size());
342 for (int j = 0; j < VV.getNumberOfElements(); j++) {
343 for (int i=0; i< VV.getNumberOfComponents(); i++)
344 VV(j,i) = Vconstant[i];
348 void ProblemCoreFlows::setInitialFieldConstant( int nDim, const vector<double> Vconstant, double xmin, double xmax, int nx, string leftSide, string rightSide,
349 double ymin, double ymax, int ny, string backSide, string frontSide,
350 double zmin, double zmax, int nz, string bottomSide, string topSide, EntityType typeField)
354 M=Mesh(xmin,xmax,nx);
356 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
358 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
360 PetscPrintf(PETSC_COMM_WORLD,"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3\n");
361 *_runLogFile<<"ProblemCoreFlows::setInitialFieldConstant: Space dimension nDim should be between 1 and 3"<<endl;
362 _runLogFile->close();
363 throw CdmathException("Space dimension nDim should be between 1 and 3");
366 M.setGroupAtPlan(xmax,0,_precision,rightSide);
367 M.setGroupAtPlan(xmin,0,_precision,leftSide);
369 M.setGroupAtPlan(ymax,1,_precision,frontSide);
370 M.setGroupAtPlan(ymin,1,_precision,backSide);
373 M.setGroupAtPlan(zmax,2,_precision,topSide);
374 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
377 setInitialFieldConstant(M, Vconstant, typeField);
379 void ProblemCoreFlows::setInitialFieldStepFunction(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, int direction, EntityType typeField)
381 if(_FECalculation && typeField!= NODES)
382 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
384 if (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
386 *_runLogFile<<"ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes"<<endl;
387 _runLogFile->close();
388 throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: Vectors VV_Left and VV_Right have different sizes");
390 Field VV("SOLVERLAB results", typeField, M, VV_Left.getNumberOfRows());
392 double component_value;
394 for (int j = 0; j < VV.getNumberOfElements(); j++)
397 component_value=VV.getElementComponent(j, direction);
400 PetscPrintf(PETSC_COMM_WORLD,"Error : space dimension is %d, direction asked is \%d \n",M.getSpaceDimension(),direction);
401 _runLogFile->close();
402 throw CdmathException( "ProblemCoreFlows::setStepFunctionInitialField: direction should be an integer between 0 and 2");
405 for (int i=0; i< VV.getNumberOfComponents(); i++)
406 if (component_value< disc_pos )
407 VV(j,i) = VV_Left[i];
409 VV(j,i) = VV_Right[i];
413 void ProblemCoreFlows::setInitialFieldStepFunction( int nDim, const vector<double> VV_Left, vector<double> VV_Right, double xstep,
414 double xmin, double xmax, int nx, string leftSide, string rightSide,
415 double ymin, double ymax, int ny, string backSide, string frontSide,
416 double zmin, double zmax, int nz, string bottomSide, string topSide, EntityType typeField)
420 M=Mesh(xmin,xmax,nx);
422 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
424 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
427 _runLogFile->close();
428 throw CdmathException("Space dimension nDim should be between 1 and 3");
431 M.setGroupAtPlan(xmax,0,_precision,rightSide);
432 M.setGroupAtPlan(xmin,0,_precision,leftSide);
434 M.setGroupAtPlan(ymax,1,_precision,frontSide);
435 M.setGroupAtPlan(ymin,1,_precision,backSide);
438 M.setGroupAtPlan(zmax,2,_precision,topSide);
439 M.setGroupAtPlan(zmin,2,_precision,bottomSide);
441 Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
442 for(int i=0;i<VV_Left.size(); i++){
443 V_Left(i)=VV_Left[i];
444 V_Right(i)=VV_Right[i];
446 setInitialFieldStepFunction(M, V_Left, V_Right, xstep, typeField);
449 void ProblemCoreFlows::setInitialFieldSphericalStepFunction(const Mesh M, const Vector Vin, const Vector Vout, double radius, const Vector Center, EntityType typeField)
451 if(_FECalculation && typeField!= NODES)
452 cout<<"Warning : finite element simulation should have initial field on nodes!!!"<<endl;
454 if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
456 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());
457 throw CdmathException("ProblemCoreFlows::setInitialFieldSphericalStepFunction : Vector size error");
459 int nVar=Vout.size();
460 int spaceDim=M.getSpaceDimension();
461 Field VV("Primitive variables for spherical step function", typeField, M, nVar);
463 Vector currentPoint(spaceDim);
464 for(int i=0;i<VV.getNumberOfElements();i++)
466 currentPoint(0)=VV.getElementComponent(i,0);
469 currentPoint(1)=VV.getElementComponent(i,1);
471 currentPoint(2)=VV.getElementComponent(i,2);
473 if((currentPoint-Center).norm()<radius)
474 for(int j=0;j<nVar;j++)
477 for(int j=0;j<nVar;j++)
483 double ProblemCoreFlows::getTime()
487 unsigned ProblemCoreFlows::getNbTimeStep()
491 double ProblemCoreFlows::getCFL()
495 double ProblemCoreFlows::getPrecision()
499 Mesh ProblemCoreFlows::getMesh()
504 void ProblemCoreFlows::setLinearSolver(linearSolver kspType, preconditioner pcType, double maxIts)
507 // set linear solver algorithm
509 _ksptype = (char*)&KSPGMRES;
510 else if (kspType==CG)
511 _ksptype = (char*)&KSPCG;
512 else if (kspType==BCGS)
513 _ksptype = (char*)&KSPBCGS;
515 PetscPrintf(PETSC_COMM_WORLD,"!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!\n");
516 *_runLogFile << "!!! Error : only 'GMRES', 'CG' or 'BCGS' is acceptable as a linear solver !!!" << endl;
517 _runLogFile->close();
518 throw CdmathException("!!! Error : only 'GMRES', 'CG' or 'BCGS' algorithm is acceptable !!!");
520 // set preconditioner
522 _pctype = (char*)&PCNONE;
523 else if (pcType ==LU)
524 _pctype = (char*)&PCLU;
525 else if (pcType == ILU)
526 _pctype = (char*)&PCILU;
527 else if (pcType ==CHOLESKY)
528 _pctype = (char*)&PCCHOLESKY;
529 else if (pcType == ICC)
530 _pctype = (char*)&PCICC;
532 PetscPrintf(PETSC_COMM_WORLD,"!!! Error : only 'NOPC', 'LU', 'ILU', 'CHOLESKY' or 'ICC' preconditioners are acceptable !!!\n");
533 *_runLogFile << "!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" << endl;
534 _runLogFile->close();
535 throw CdmathException("!!! Error : only 'NOPC' or 'LU' or 'ILU' preconditioners are acceptable !!!" );
539 void ProblemCoreFlows::setNewtonSolver(double precision, int iterations, nonLinearSolver solverName)
541 _maxNewtonIts=iterations;
542 _precision_Newton=precision;
543 _nonLinearSolver=solverName;
548 // Cette methode lance une execution du ProblemCoreFlows
549 // Elle peut etre utilisee si le probleme n'est couple a aucun autre.
550 // (s'il n'a besoin d'aucun champ d'entree).
551 // Precondition: initialize
552 // Seule la methode terminate peut etre appelée apres
553 bool ProblemCoreFlows::run()
555 if(!_initializedMemory)
557 _runLogFile->close();
558 throw CdmathException("ProblemCoreFlows::run() call initialize() first");
560 bool stop=false; // Does the Problem want to stop (error) ?
561 bool ok; // Is the time interval successfully solved ?
562 _isStationary=false;//in case of a second run with a different physics or cfl
564 PetscPrintf(PETSC_COMM_WORLD,"Running test case %s\n",_fileName.c_str());
566 _runLogFile->open((_fileName+".log").c_str(), ios::out | ios::trunc);;//for creation of a log file to save the history of the simulation
567 *_runLogFile<< "Running test case "<< _fileName<<endl;
570 while(!stop && !_isStationary &&_time<_timeMax && _nbTimeStep<_maxNbOfTimeStep)
572 ok=false; // Is the time interval successfully solved ?
574 // Guess the next time step length
575 _dt=computeTimeStep(stop);
577 PetscPrintf(PETSC_COMM_WORLD,"Failed computing time step %d, time = %.2f, dt= %.2e, stopping calculation",_nbTimeStep,_time,_dt);
578 *_runLogFile << "Failed computing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
581 // Loop on the time interval tries
582 while (!ok && !stop )
584 stop=!initTimeStep(_dt);
585 // Prepare the next time step
587 PetscPrintf(PETSC_COMM_WORLD,"Failed initializing time step %d, time = %.2f, dt= %.2e, stopping calculation",_nbTimeStep,_time,_dt);
588 *_runLogFile << "Failed initializing time step "<<_nbTimeStep<<", time = " << _time <<", dt= "<<_dt<<", stopping calculation"<< endl;
591 // Solve the next time step
594 if (!ok) // The resolution failed, try with a new time interval.
596 /* if(_dt>_precision){
597 cout<<"ComputeTimeStep returned _dt="<<_dt<<endl;
598 cout << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with dt/2"<< endl;
599 *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", time = " << _time <<" _dt= "<<_dt<<", cfl= "<<_cfl<<", trying again with dt/2"<< endl;
600 double dt=_dt/2;//We chose to divide the time step by 2
601 abortTimeStep();//Cancel the initTimeStep
602 _dt=dt;//new value of time step is previous time step divided by 2 (we do not call computeTimeStep
603 //_cfl*=0.5;//If we change the cfl, we must compute the new time step with computeTimeStep
604 //_dt=computeTimeStep(stop);
607 PetscPrintf(PETSC_COMM_WORLD,"Failed solving time step %d, time = %.2f, dt= %.2e, cfl = %.2f, stopping calculation \n",_nbTimeStep,_time,_dt,_cfl);
608 *_runLogFile << "Failed solving time step "<<_nbTimeStep<<", _time = " << _time<<" _dt= "<<_dt<<", cfl= "<<_cfl <<", stopping calculation"<< endl;
609 stop=true; // Impossible to solve the next time step, the Problem has given up
613 else // The resolution was successful, validate and go to the next time step.
616 if ((_nbTimeStep-1)%_freqSave ==0){
617 PetscPrintf(PETSC_COMM_WORLD,"Time step = %d, dt = %.2e, time = %.2f, ||Un+1-Un||= %.2e\n\n",_nbTimeStep,_dt,_time,_erreur_rel);
618 *_runLogFile << "Time step = "<< _nbTimeStep << ", dt = "<< _dt <<", time = "<<_time << ", ||Un+1-Un||= "<<_erreur_rel<<endl<<endl;
624 PetscPrintf(PETSC_COMM_WORLD,"Stationary state reached\n");
625 *_runLogFile << "Stationary state reached" <<endl;
627 else if(_time>=_timeMax){
628 PetscPrintf(PETSC_COMM_WORLD,"Maximum time %.2f reached\n",_timeMax);
629 *_runLogFile<<"Maximum time "<<_timeMax<<" reached"<<endl;
631 else if(_nbTimeStep>=_maxNbOfTimeStep){
632 PetscPrintf(PETSC_COMM_WORLD,"Maximum number of time steps %d reached\n",_maxNbOfTimeStep);
633 *_runLogFile<<"Maximum number of time steps "<<_maxNbOfTimeStep<<" reached"<<endl;
636 PetscPrintf(PETSC_COMM_WORLD,"Error problem wants to stop!\n");
637 *_runLogFile<<"Error problem wants to stop!"<<endl;
639 PetscPrintf(PETSC_COMM_WORLD,"End of calculation at time t = %.2f and time step number %d\n",_time,_nbTimeStep);
640 *_runLogFile << "End of calculation time t= " << _time << " at time step number "<< _nbTimeStep << endl;
642 _runLogFile->close();
646 void ProblemCoreFlows::displayMatrix(double *matrix, int size, string name)
649 for(int p=0; p<size; p++)
651 for(int q=0; q<size; q++)
653 cout << matrix[p*size+q] << "\t";
660 void ProblemCoreFlows::displayVector(double *vector, int size, string name)
663 for(int q=0; q<size; q++)
665 cout << vector[q] << "\t";
669 void ProblemCoreFlows::setFileName(string fileName){
670 if( _nbTimeStep>0 && fileName!=_fileName)//This is a change of file name during a simulation
671 _restartWithNewFileName=true;
672 _fileName = fileName;
675 void ProblemCoreFlows::setFreqSave(int freqSave){
676 _freqSave = freqSave;
679 bool ProblemCoreFlows::solveTimeStep(){
681 bool converged=false, ok=true;
682 while(!converged && ok && _NEWTON_its < _maxNewtonIts){
683 ok=iterateTimeStep(converged);//resolution du systeme lineaire si schema implicite
685 if(_timeScheme == Implicit && (_nbTimeStep-1)%_freqSave ==0)//To monitor the convergence of the newton scheme
687 PetscPrintf(PETSC_COMM_WORLD," Newton iteration %d, %s iterations : %d maximum variation ||Uk+1-Uk||: %.2e\n",_NEWTON_its,_ksptype,_PetscIts,_erreur_rel);
688 *_runLogFile<< " Newton iteration " << _NEWTON_its<< ", "<< _ksptype << " iterations : " << _PetscIts<< " maximum variation ||Uk+1-Uk||: " << _erreur_rel << endl;
692 PetscReal sv_max, sv_min;
693 KSPComputeExtremeSingularValues(_ksp, &sv_max, &sv_min);
694 PetscPrintf(PETSC_COMM_WORLD," Singular value max = %.2e, singular value min = %.2e, condition number = %.2e\n",sv_max,sv_min,sv_max/sv_min);
695 *_runLogFile<<" Singular value max = " << sv_max <<", singular value min = " << sv_min <<", condition number = " << sv_max/sv_min <<endl;
701 if(_NEWTON_its >= _maxNewtonIts){
702 PetscPrintf(PETSC_COMM_WORLD,"Maximum number of Newton iterations %d reached\n",_maxNewtonIts);
703 *_runLogFile << "Maximum number of Newton iterations "<<_maxNewtonIts<<" reached"<< endl;
706 PetscPrintf(PETSC_COMM_WORLD,"iterateTimeStep: solving Newton iteration %d Failed\n",_NEWTON_its);
707 *_runLogFile<<"iterateTimeStep: solving Newton iteration "<<_NEWTON_its<<" Failed"<<endl;
710 else if(_timeScheme == Implicit && (_nbTimeStep-1)%_freqSave ==0)
712 PetscPrintf(PETSC_COMM_WORLD,"Nombre d'iterations de Newton %d, Nombre max d'iterations %s : %d\n\n",_NEWTON_its, _ksptype, _MaxIterLinearSolver);
714 *_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its << "Nombre max d'iterations "<< _ksptype << " : " << _MaxIterLinearSolver << endl << endl;
715 _MaxIterLinearSolver = 0;
721 ProblemCoreFlows::~ProblemCoreFlows()
723 PetscBool petscInitialized;
724 PetscInitialized(&petscInitialized);
732 ProblemCoreFlows::getConditionNumber(bool isSingular, double tol) const
734 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
735 return A.getConditionNumber( isSingular, tol);
737 std::vector< double >
738 ProblemCoreFlows::getEigenvalues(int nev, EPSWhich which, double tol) const
740 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
741 return A.getEigenvalues( nev, which, tol);
743 std::vector< Vector >
744 ProblemCoreFlows::getEigenvectors(int nev, EPSWhich which, double tol) const
746 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
747 return A.getEigenvectors( nev, which, tol);
750 ProblemCoreFlows::getEigenvectorsField(int nev, EPSWhich which, double tol) const
752 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
753 MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
756 my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
758 my_eigenfield.setFieldByDataArrayDouble(d);
760 return my_eigenfield;
763 std::vector< double >
764 ProblemCoreFlows::getSingularValues( int nsv, SVDWhich which, double tol) const
766 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
767 return A.getSingularValues( nsv, which, tol);
769 std::vector< Vector >
770 ProblemCoreFlows::getSingularVectors(int nsv, SVDWhich which, double tol) const
772 SparseMatrixPetsc A = SparseMatrixPetsc(_A);
773 return A.getSingularVectors( nsv, which, tol);
777 ProblemCoreFlows::getUnknownField() const
779 if(!_initializedMemory)
781 _runLogFile->close();
782 throw CdmathException("ProblemCoreFlows::getUnknownField() call initialize() first");
788 ProblemCoreFlows::isMEDCoupling64Bits() const
790 #ifdef MEDCOUPLING_USE_64BIT_IDS