1 #include "TransportEquation.hxx"
8 TransportEquation::TransportEquation(phase fluid, pressureMagnitude pEstimate,vector<double> vitesseTransport){
9 if(pEstimate==around1bar300KTransport){
11 if(fluid==GasPhase){//Nitrogen pressure 1 bar and temperature 27°C
12 _href=3.11e5; //nitrogen enthalpy at 1 bar and 300K
13 _cpref=1041;//nitrogen specific heat at constant pressure 1 bar and 300K
14 //saturation data for nitrogen at 1 bar and 77K
15 _hsatv=0.77e5;//nitrogen vapour enthalpy at saturation at 1 bar
16 _hsatl=-1.22e5;//nitrogen liquid enthalpy at saturation at 1 bar
17 _rhosatv=4.556;//nitrogen vapour density at saturation at 1 bar
18 _rhosatl=806.6;//nitrogen liquid density at saturation at 1 bar
21 //Water at pressure 1 bar and temperature 27°C
22 _href=1.127e5; //water enthalpy at 1 bar and 300K
23 _cpref=4181;//water specific heat at 1 bar and 300K
24 //saturation data for water at 1 bar and 373K
25 _hsatv=2.675e6;//Gas enthalpy at saturation at 1 bar
26 _hsatl=4.175e5;//water enthalpy at saturation at 1 bar
27 _rhosatv=0.6;//Gas density at saturation at 1 bar
28 _rhosatl=958;//water density at saturation at 1 bar
31 else{//around155bars600K
34 _href=2.675e6; //Gas enthalpy at 155 bars and 618K
35 _cpref=14001;//Gas specific heat at 155 bar and 618K
38 _href=4.175e5;//water enthalpy at 155 bars and 618K
39 _cpref=8950;//water specific heat at 155 bar and 618K
41 //saturation data for water at 155 bars and 618K
42 _hsatv=2.6e6;//Gas enthalpy at saturation at 155 bars
43 _hsatl=1.63e6;//water enthalpy at saturation at 155 bars
44 _rhosatv=101.9;//Gas density at saturation at 155 bars
45 _rhosatl=594.4;//water density at saturation at 155 bars
47 _Ndim=vitesseTransport.size();
48 _vitesseTransport=Vector(_Ndim);
49 for(int i=0;i<_Ndim;i++)
50 _vitesseTransport[i]=vitesseTransport[i];
54 _transportMatrixSet=false;
57 void TransportEquation::initialize()
60 throw CdmathException("TransportEquation::initialize() set initial data first");
62 cout<<"Initialising the transport of a fluid enthalpy"<<endl;
63 /**************** Field creation *********************/
65 //post processing fields used only for saving results
66 _TT=Field ("Temperature", CELLS, _mesh, 1);
67 _Alpha=Field ("Void fraction", CELLS, _mesh, 1);
68 _Rho=Field ("Mixture density", CELLS, _mesh, 1);
69 //Construction des champs de post-traitement
70 VecCreate(PETSC_COMM_SELF, &_Hn);
71 VecSetSizes(_Hn,PETSC_DECIDE,_Nmailles);
72 VecSetFromOptions(_Hn);
73 for(int i =0; i<_Nmailles;i++){
74 _TT(i)=temperature(_VV(i));
75 _Alpha(i)=voidFraction(_VV(i));
76 _Rho(i)=density(_Alpha(i));
77 VecSetValue(_Hn,i,_VV(i), INSERT_VALUES);
79 if(!_heatPowerFieldSet){
80 _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
81 for(int i =0; i<_Nmailles; i++)
82 _heatPowerField(i) = _heatSource;
84 if(!_rodTemperatureFieldSet){
85 _rodTemperatureField=Field("Rod temperature",CELLS,_mesh,1);
86 for(int i =0; i<_Nmailles; i++)
87 _rodTemperatureField(i) = _rodTemperature;
90 //creation de la matrice
91 MatCreateSeqAIJ(PETSC_COMM_SELF, _Nmailles, _Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
92 VecDuplicate(_Hn, &_Hk);
93 VecDuplicate(_Hn, &_Hkm1);
94 VecDuplicate(_Hn, &_deltaH);
95 VecDuplicate(_Hn, &_b);//RHS of the linear system: _b=Hn/dt + _b0 + puisance
96 VecDuplicate(_Hn, &_b0);//part of the RHS that comes from the boundary conditions
99 KSPCreate(PETSC_COMM_SELF, &_ksp);
100 KSPSetType(_ksp, _ksptype);
101 // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
102 KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
103 KSPGetPC(_ksp, &_pc);
104 PCSetType(_pc, _pctype);
106 _initializedMemory=true;
107 save();//save initial data
110 double TransportEquation::computeTimeStep(bool & stop){
111 if(!_transportMatrixSet)
112 _dt_transport=computeTransportMatrix();
113 _dt_src=computeRHS();
116 return min(_dt_transport,_dt_src);
118 double TransportEquation::computeTransportMatrix(){
119 long nbFaces = _mesh.getNumberOfFaces();
123 double inv_dxi, inv_dxj;
124 Vector normale(_Ndim);
127 std::vector< int > idCells;
130 for (int j=0; j<nbFaces;j++){
131 Fj = _mesh.getFace(j);
133 // compute the normal vector corresponding to face j : from idCells[0] to idCells[1]
134 idCells = Fj.getCellsId();
135 Cell1 = _mesh.getCell(idCells[0]);
138 for(int l=0; l<Cell1.getNumberOfFaces(); l++){
139 if (j == Cell1.getFacesId()[l]){
140 for (int idim = 0; idim < _Ndim; ++idim)
141 normale[idim] = Cell1.getNormalVector(l,idim);
145 }else{ // _Ndim = 1 : assume that this is normal mesh : the face index increases in positive direction
146 if (Fj.getNumberOfCells()<2) {
149 else if (j==nbFaces-1)
152 throw CdmathException("TransportEquation::ComputeTimeStep(): computation of normal vector failed");
153 } else if(Fj.getNumberOfCells()==2){
154 if (idCells[0] < idCells[1])
160 //Compute velocity at the face Fj
161 un=normale*_vitesseTransport;
165 // compute 1/dxi = volume of Ci/area of Fj
167 inv_dxi = Fj.getMeasure()/Cell1.getMeasure();
169 inv_dxi = 1/Cell1.getMeasure();
171 // If Fj is on the boundary
172 if (Fj.getNumberOfCells()==1) {
173 if(_verbose && _nbTimeStep%_freqSave ==0)
175 cout << "face numero " << j << " cellule frontiere " << idCells[0] << " ; vecteur normal=(";
176 for(int p=0; p<_Ndim; p++)
177 cout << normale[p] << ",";
180 nameOfGroup = Fj.getGroupName();
182 if (_limitField[nameOfGroup].bcType==NeumannTransport){
183 MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
185 else if(_limitField[nameOfGroup].bcType==InletTransport){
187 MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
190 hk=_limitField[nameOfGroup].h;
191 VecSetValue(_b0,idm,-inv_dxi*un*hk, ADD_VALUES);
195 cout<<"!!!!!!!!!!!!!!! Error TransportEquation::computeTransportMatrix() !!!!!!!!!!"<<endl;
196 cout<<"!!!!!!!!! Boundary condition not treated for boundary named "<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<" !!!!!!!!!!!!!! "<<endl;
197 cout<<"Accepted boundary conditions are NeumannTransport "<<NeumannTransport<< " and InletTransport "<< InletTransport <<endl;
198 throw CdmathException("Boundary condition not accepted");
200 // if Fj is inside the domain
201 } else if (Fj.getNumberOfCells()==2 ){
202 if(_verbose && _nbTimeStep%_freqSave ==0)
204 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
205 cout << " ; vecteur normal=(";
206 for(int p=0; p<_Ndim; p++)
207 cout << normale[p] << ",";
210 Cell2 = _mesh.getCell(idCells[1]);
213 inv_dxj = Fj.getMeasure()/Cell2.getMeasure();
215 inv_dxj = 1/Cell2.getMeasure();
218 MatSetValue(_A,idm,idm,inv_dxi*un, ADD_VALUES);
219 MatSetValue(_A,idn,idm,-inv_dxj*un, ADD_VALUES);
222 MatSetValue(_A,idm,idn,inv_dxi*un, ADD_VALUES);
223 MatSetValue(_A,idn,idn,-inv_dxj*un, ADD_VALUES);
227 throw CdmathException("TransportEquation::ComputeTimeStep(): incompatible number of cells around a face");
229 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
230 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
231 VecAssemblyBegin(_b0);
233 _transportMatrixSet=true;
234 if(abs(_maxvp)<_precision)
235 throw CdmathException("TransportEquation::computeTransportMatrix(): maximum eigenvalue for time step is zero");
237 return _cfl*_minl/_maxvp;
239 double TransportEquation::computeRHS(){
240 double rhomin=INFINITY;
242 VecAssemblyBegin(_b);
244 cout<<"second membre of transport problem"<<endl;
245 for (int i=0; i<_Nmailles;i++){
246 VecSetValue(_b,i,_heatTransfertCoeff*(_rodTemperatureField(i)-_TT(i))/_Rho(i),ADD_VALUES);
247 VecSetValue(_b,i,_heatPowerField(i)/_Rho(i),ADD_VALUES);
249 cout<<_heatPowerField(i)/_Rho(i)<<endl;
255 VecView(_b, PETSC_VIEWER_STDOUT_WORLD);
257 return rhomin*_cpref/_heatTransfertCoeff;
259 void TransportEquation::updatePrimitives()
262 for(int i=0; i<_Nmailles; i++)
264 VecGetValues(_Hk, 1, &i, &hi);
266 _TT(i)=temperature(hi);
267 _Alpha(i)=voidFraction(hi);
268 _Rho(i)=density(_Alpha(i));
272 bool TransportEquation::initTimeStep(double dt){
274 if(_verbose && _nbTimeStep%_freqSave ==0)
275 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
276 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
277 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
279 if(_timeScheme == Implicit)
285 void TransportEquation::abortTimeStep(){
286 VecAXPY(_b, -1/_dt, _Hn);
287 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
288 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
290 if(_timeScheme == Implicit)
295 bool TransportEquation::iterateTimeStep(bool &converged)
299 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
301 computeTimeStep(stop);
307 VecAXPY(_b, 1/_dt, _Hn);
310 cout << "Vecteur Hn: " << endl;
311 VecView(_Hn, PETSC_VIEWER_STDOUT_WORLD);
313 cout<<"Vecteur _b "<<endl;
314 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
315 cout << "Matrice A "<<endl;
316 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
319 if(_timeScheme == Explicit)
321 MatMult(_A, _Hn, _Hk);
324 cout << "Nouveau vecteur Hk: " << endl;
325 VecView(_Hk, PETSC_VIEWER_STDOUT_WORLD);
328 VecAXPY(_Hk, -1, _b);
331 cout << "Nouveau vecteur Hk-b: " << endl;
332 VecView(_Hk, PETSC_VIEWER_STDOUT_WORLD);
338 cout << "Nouveau vecteur dt*(Hk-b): " << endl;
339 VecView(_Hk, PETSC_VIEWER_STDOUT_WORLD);
346 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
347 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
349 #if PETSC_VERSION_GREATER_3_5
350 KSPSetOperators(_ksp, _A, _A);
352 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
356 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
358 KSPSolve(_ksp, _b, _Hk);
361 KSPGetIterationNumber(_ksp, &_PetscIts);
362 if( _MaxIterLinearSolver < _PetscIts)
363 _MaxIterLinearSolver = _PetscIts;
364 if(_PetscIts>=_maxPetscIts)
366 cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
372 VecCopy(_Hk, _deltaH);//ici on a deltaH=Hk
373 VecAXPY(_deltaH, -1, _Hkm1);//On obtient deltaH=Hk-Hkm1
377 for(int i=0; i<_Nmailles; i++)
379 VecGetValues(_deltaH, 1, &i, &dhi);
380 VecGetValues(_Hk, 1, &i, &hi);
381 if(_erreur_rel < fabs(dhi/hi))
382 _erreur_rel = fabs(dhi/hi);
386 converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
396 void TransportEquation::validateTimeStep()
398 VecCopy(_Hk, _deltaH);//ici Hk=Hnp1 donc on a deltaH=Hnp1
399 VecAXPY(_deltaH, -1, _Hn);//On obtient deltaH=Hnp1-Hn
404 for(int i=0; i<_Nmailles; i++)
406 VecGetValues(_deltaH, 1, &i, &dhi);
407 VecGetValues(_Hk, 1, &i, &hi);
408 if(_erreur_rel < fabs(dhi/hi))
409 _erreur_rel = fabs(dhi/hi);
411 _isStationary =(_erreur_rel <_precision);
416 if(_nbTimeStep%_freqSave ==0)
418 cout <<"Valeur propre locale max: " << _maxvp << endl;
419 //Find minimum and maximum void fractions
420 double alphamin=INFINITY;
421 double alphamax=-INFINITY;
422 for(int i=0; i<_Nmailles; i++)
424 if(_Alpha(i)>alphamax)
426 if(_Alpha(i)<alphamin)
429 cout<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<endl;
437 void TransportEquation::terminate(){
441 VecDestroy(&_deltaH);
447 void TransportEquation::save(){
448 string resultFile(_path+"/TransportEquation_");///Results
449 resultFile+=_fileName;
451 _VV.setTime(_time,_nbTimeStep);
452 _TT.setTime(_time,_nbTimeStep);
453 _Alpha.setTime(_time,_nbTimeStep);
454 _Rho.setTime(_time,_nbTimeStep);
456 // create mesh and component info
457 if (_nbTimeStep ==0 || _restartWithNewFileName){
458 if (_restartWithNewFileName)
459 _restartWithNewFileName=false;
460 string suppress ="rm -rf "+resultFile+"_*";
461 system(suppress.c_str());//Nettoyage des précédents calculs identiques
466 _VV.writeVTK(resultFile+"Enthalpy");
467 _TT.writeVTK(resultFile+"Temperature");
468 _Alpha.writeVTK(resultFile+"GasFraction");
469 _Rho.writeVTK(resultFile+"MixtureDensity");
472 _VV.writeMED(resultFile+"Enthalpy");
473 _TT.writeMED(resultFile+"Temperature");
474 _Alpha.writeMED(resultFile+"GasFraction");
475 _Rho.writeMED(resultFile+"MixtureDensity");
478 _VV.writeCSV(resultFile+"Enthalpy");
479 _TT.writeCSV(resultFile+"Temperature");
480 _Alpha.writeCSV(resultFile+"GasFraction");
481 _Rho.writeCSV(resultFile+"MixtureDensity");
485 // do not create mesh
490 _VV.writeVTK(resultFile+"Enthalpy",false);
491 _TT.writeVTK(resultFile+"Temperature",false);
492 _Alpha.writeVTK(resultFile+"GasFraction",false);
493 _Rho.writeVTK(resultFile+"MixtureDensity",false);
496 _VV.writeMED(resultFile+"Enthalpy",false);
497 _TT.writeMED(resultFile+"Temperature",false);
498 _Alpha.writeMED(resultFile+"GasFraction",false);
499 _Rho.writeMED(resultFile+"MixtureDensity",false);
502 _VV.writeCSV(resultFile+"Enthalpy");
503 _TT.writeCSV(resultFile+"Temperature");
504 _Alpha.writeCSV(resultFile+"GasFraction");
505 _Rho.writeCSV(resultFile+"MixtureDensity");
511 vector<string> TransportEquation::getOutputFieldsNames()
513 vector<string> result(2);
515 result[0]="Enthalpy";
516 result[1]="FluidTemperature";
521 Field& TransportEquation::getOutputField(const string& nameField )
523 if(nameField=="FluidTemperature" || nameField=="FLUIDTEMPERATURE" )
524 return getFluidTemperatureField();
525 else if(nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPY" )
526 return getEnthalpyField();
529 cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
530 throw CdmathException("TransportEquation::getOutputField error : Unknown Field name");