1 #include "ProblemFluid.hxx"
7 ProblemFluid::ProblemFluid(void){
12 _porosityFieldSet=false;
13 _pressureLossFieldSet=false;
14 _sectionFieldSet=false;
15 _GravityField3d=vector<double>(3,0);
16 _gravityReferencePoint=vector<double>(3,0);
17 _Uroe=NULL;_Udiff=NULL;_temp=NULL;_l=NULL;_vec_normal=NULL;;
20 _saveConservativeField=false;
21 _usePrimitiveVarsInNewton=false;
22 _saveInterfacialField=false;
23 _err_press_max=0; _part_imag_max=0; _nbMaillesNeg=0; _nbVpCplx=0;_minm1=1e30;_minm2=1e30;//Pour affichage paramètres diphasiques dans IterateTimeStep
25 _entropicCorrection=false;
26 _pressureCorrectionOrder=2;
27 _nonLinearFormulation=Roe;
32 SpaceScheme ProblemFluid::getSpaceScheme()
36 void ProblemFluid::setNumericalScheme(SpaceScheme spaceScheme, TimeScheme timeScheme)
38 if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
39 _restartWithNewTimeScheme=true;
40 _timeScheme = timeScheme;
41 _spaceScheme = spaceScheme;
44 void ProblemFluid::initialize()
47 *_runLogFile<<"ProblemFluid::initialize() set initial data first"<<endl;
49 throw CdmathException("ProblemFluid::initialize() set initial data first");
51 cout << "Number of Phases = " << _nbPhases << " mesh dimension = "<<_Ndim<<" number of variables = "<<_nVar<<endl;
52 *_runLogFile << "Number of Phases = " << _nbPhases << " spaceDim= "<<_Ndim<<" number of variables= "<<_nVar<<endl;
54 /********* local arrays ****************/
55 _AroePlus = new PetscScalar[_nVar*_nVar];
56 _Diffusion = new PetscScalar[_nVar*_nVar];
57 _AroeMinus = new PetscScalar[_nVar*_nVar];
58 _Aroe = new PetscScalar[_nVar*_nVar];
59 _absAroe = new PetscScalar[_nVar*_nVar];
60 _signAroe = new PetscScalar[_nVar*_nVar];
61 _invAroe = new PetscScalar[_nVar*_nVar];
62 _AroeImplicit = new PetscScalar[_nVar*_nVar];
63 _AroeMinusImplicit = new PetscScalar[_nVar*_nVar];
64 _AroePlusImplicit = new PetscScalar[_nVar*_nVar];
65 _absAroeImplicit = new PetscScalar[_nVar*_nVar];
66 _primToConsJacoMat = new PetscScalar[_nVar*_nVar];
67 _phi = new PetscScalar[_nVar];
68 _Jcb = new PetscScalar[_nVar*_nVar];
69 _JcbDiff = new PetscScalar[_nVar*_nVar];
70 _a = new PetscScalar[_nVar*_nVar];
71 _externalStates = new PetscScalar[_nVar];
72 _Ui = new PetscScalar[_nVar];
73 _Uj = new PetscScalar[_nVar];
74 _Vi = new PetscScalar[_nVar];
75 _Vj = new PetscScalar[_nVar];
76 _Uij = new PetscScalar[_nVar];//used for VFRoe scheme
77 _Vij = new PetscScalar[_nVar];//used for VFRoe scheme
78 _Si = new PetscScalar[_nVar];
79 _Sj = new PetscScalar[_nVar];
80 _pressureLossVector= new PetscScalar[_nVar];
81 _porosityGradientSourceVector= new PetscScalar[_nVar];
83 _l = new double[_nVar];
84 _r = new double[_nVar];
85 _Udiff = new double[_nVar];
86 _temp = new double[_nVar*_nVar];
88 _idm = new int[_nVar];
89 _idn = new int[_nVar];
90 _vec_normal = new double[_Ndim];
92 for(int k=0;k<_nVar*_nVar;k++)
97 for(int k=0; k<_nVar; k++){
102 /**************** Field creation *********************/
103 if(!_heatPowerFieldSet){
104 _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
105 for(int i =0; i<_Nmailles; i++)
106 _heatPowerField(i) = _heatSource;
109 _dp_over_dt=Field("dP/dt",CELLS,_mesh,1);
110 if(!_porosityFieldSet){
111 _porosityField=Field("Porosity field",CELLS,_mesh,1);
112 for(int i =0; i<_Nmailles; i++)
113 _porosityField(i) = 1;
114 _porosityFieldSet=true;
117 //conservative field used only for saving results
118 _UU=Field ("Conservative vec", CELLS, _mesh, _nVar);
119 if(_saveInterfacialField && _nonLinearFormulation==VFRoe)
121 _UUstar=Field ("Interfacial U", CELLS, _mesh, _nVar);
122 _VVstar=Field ("Interfacial V", CELLS, _mesh, _nVar);
125 //Construction des champs primitifs et conservatifs initiaux comme avant dans ParaFlow
126 double * initialFieldPrim = new double[_nVar*_Nmailles];
127 for(int i =0; i<_Nmailles;i++)
128 for(int j =0; j<_nVar;j++)
129 initialFieldPrim[i*_nVar+j]=_VV(i,j);
130 double *initialFieldCons = new double[_nVar*_Nmailles];
131 for(int i=0; i<_Nmailles; i++)
132 primToCons(initialFieldPrim, i, initialFieldCons, i);
133 for(int i =0; i<_Nmailles;i++)
134 for(int j =0; j<_nVar;j++)
135 _UU(i,j)=initialFieldCons[i*_nVar+j];
137 /**********Petsc structures: ****************/
139 //creation de la matrice
140 if(_timeScheme == Implicit)
141 MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
143 //creation des vecteurs
144 VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uext);
145 VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Vext);
146 VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uextdiff);
147 // VecCreateSeq(PETSC_COMM_SELF, _nVar*_Nmailles, &_conservativeVars);
148 VecCreate(PETSC_COMM_SELF, &_conservativeVars);
149 VecSetSizes(_conservativeVars,PETSC_DECIDE,_nVar*_Nmailles);
150 VecSetBlockSize(_conservativeVars,_nVar);
151 VecSetFromOptions(_conservativeVars);
152 VecDuplicate(_conservativeVars, &_old);
153 VecDuplicate(_conservativeVars, &_newtonVariation);
154 VecDuplicate(_conservativeVars, &_b);
155 VecDuplicate(_conservativeVars, &_primitiveVars);
159 VecDuplicate(_conservativeVars, &_vecScaling);
160 VecDuplicate(_conservativeVars, &_invVecScaling);
161 VecDuplicate(_conservativeVars, &_bScaling);
163 _blockDiag = new PetscScalar[_nVar];
164 _invBlockDiag = new PetscScalar[_nVar];
168 int *indices = new int[_Nmailles];
169 for(int i=0; i<_Nmailles; i++)
171 VecSetValuesBlocked(_conservativeVars, _Nmailles, indices, initialFieldCons, INSERT_VALUES);
172 VecAssemblyBegin(_conservativeVars);
173 VecAssemblyEnd(_conservativeVars);
174 VecCopy(_conservativeVars, _old);
175 VecAssemblyBegin(_old);
176 VecAssemblyEnd(_old);
177 VecSetValuesBlocked(_primitiveVars, _Nmailles, indices, initialFieldPrim, INSERT_VALUES);
178 VecAssemblyBegin(_primitiveVars);
179 VecAssemblyEnd(_primitiveVars);
182 cout << "Variables primitives initiales : " << endl;
183 VecView(_primitiveVars, PETSC_VIEWER_STDOUT_WORLD);
185 cout<<"Variables conservatives initiales : "<<endl;
186 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_SELF);
189 delete[] initialFieldPrim;
190 delete[] initialFieldCons;
194 KSPCreate(PETSC_COMM_SELF, &_ksp);
195 KSPSetType(_ksp, _ksptype);
196 // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
197 KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
198 KSPGetPC(_ksp, &_pc);
199 PCSetType(_pc, _pctype);
201 _initializedMemory=true;
202 save();//save initial data
205 bool ProblemFluid::initTimeStep(double dt){
207 return _dt>0;//No need to call MatShift as the linear system matrix is filled at each Newton iteration (unlike linear problem)
210 bool ProblemFluid::iterateTimeStep(bool &converged)
214 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
216 computeTimeStep(stop);//This compute timestep is just to update the linear system. The time step was imposed befor starting the Newton iterations
218 if(stop){//Le compute time step ne s'est pas bien passé
219 cout<<"ComputeTimeStep failed"<<endl;
224 computeNewtonVariation();
226 //converged=convergence des iterations
227 if(_timeScheme == Explicit)
229 else{//Implicit scheme
231 KSPGetIterationNumber(_ksp, &_PetscIts);
232 if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
233 _MaxIterLinearSolver = _PetscIts;
235 KSPConvergedReason reason;
236 KSPGetConvergedReason(_ksp,&reason);
238 if(reason<0)//solving the linear system failed
240 cout<<"Systeme lineaire : pas de convergence de Petsc. Raison PETSC numéro "<<reason<<endl;
241 cout<<"Nombre d'itérations effectuées "<< _PetscIts<<" nombre maximal Itérations autorisées "<<_maxPetscIts<<endl;
242 *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Raison PETSC numéro "<<reason<<endl;
243 *_runLogFile<<"Nombre d'itérations effectuées "<< _PetscIts<<" nombre maximal Itérations autorisées "<<_maxPetscIts<<endl;
247 else{//solving the linear system succeeded
248 //Calcul de la variation relative Uk+1-Uk
252 for(int j=1; j<=_Nmailles; j++)
254 for(int k=0; k<_nVar; k++)
257 VecGetValues(_newtonVariation, 1, &I, &dx);
258 VecGetValues(_conservativeVars, 1, &I, &x);
259 if (fabs(x)*fabs(x)< _precision)
261 if(_erreur_rel < fabs(dx))
262 _erreur_rel = fabs(dx);
264 else if(_erreur_rel < fabs(dx/x))
265 _erreur_rel = fabs(dx/x);
269 converged = _erreur_rel <= _precision_Newton;
272 double relaxation=1;//Uk+1=Uk+relaxation*daltaU
274 VecAXPY(_conservativeVars, relaxation, _newtonVariation);
276 //mise a jour du champ primitif
279 if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
281 cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
282 *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
288 cout<<"Vecteur Ukp1-Uk "<<endl;
289 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
290 cout << "Nouvel etat courant Uk de l'iteration Newton: " << endl;
291 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_SELF);
294 if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
295 if(_minm1<-_precision || _minm2<-_precision)
297 cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
298 *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
302 cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
303 *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
315 double ProblemFluid::computeTimeStep(bool & stop){
318 if(_restartWithNewTimeScheme)//This is a change of time scheme during a simulation
320 if(_timeScheme == Implicit)
321 MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
324 _restartWithNewTimeScheme=false;
326 if(_timeScheme == Implicit)
329 VecAssemblyBegin(_b);
330 VecAssemblyBegin(_conservativeVars);
331 std::vector< int > idCells;
332 PetscInt idm, idn, size = 1;
334 long nbFaces = _mesh.getNumberOfFaces();
339 for (int j=0; j<nbFaces;j++){
340 Fj = _mesh.getFace(j);
341 _isBoundary=Fj.isBorder();
342 idCells = Fj.getCellsId();
344 // If Fj is on the boundary
347 for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
349 // compute the normal vector corresponding to face j : from Ctemp1 outward
350 Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
352 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
353 {//we look for l the index of the face Fj for the cell Ctemp1
354 if (j == Ctemp1.getFacesId()[l])
356 for (int idim = 0; idim < _Ndim; ++idim)
357 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
361 }else{ // _Ndim = 1, build normal vector (bug cdmath)
362 if(!_sectionFieldSet)
364 if (Fj.x()<Ctemp1.x())
377 if(_verbose && _nbTimeStep%_freqSave ==0)
379 cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
380 for(int p=0; p<_Ndim; p++)
381 cout << _vec_normal[p] << ",";
384 nameOfGroup = Fj.getGroupName();
385 _porosityi=_porosityField(idCells[k]);
386 _porosityj=_porosityi;
387 setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
388 convectionState(idCells[k],0,true);
389 convectionMatrices();
390 diffusionStateAndMatrices(idCells[k], 0, true);
393 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
395 _inv_dxi = 1/Ctemp1.getMeasure();
397 addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
398 addDiffusionToSecondMember(idCells[k],-1,true);
399 addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
401 if(_timeScheme == Implicit){
402 for(int l=0; l<_nVar*_nVar;l++){
403 _AroeMinusImplicit[l] *= _inv_dxi;
404 _Diffusion[l] *=_inv_dxi*_inv_dxi;
407 jacobian(idCells[k],nameOfGroup,_vec_normal);
408 jacobianDiff(idCells[k],nameOfGroup);
409 if(_verbose && _nbTimeStep%_freqSave ==0){
410 cout << "Matrice Jacobienne CL convection:" << endl;
411 for(int p=0; p<_nVar; p++){
412 for(int q=0; q<_nVar; q++)
413 cout << _Jcb[p*_nVar+q] << "\t";
417 cout << "Matrice Jacobienne CL diffusion:" << endl;
418 for(int p=0; p<_nVar; p++){
419 for(int q=0; q<_nVar; q++)
420 cout << _JcbDiff[p*_nVar+q] << "\t";
427 //calcul et insertion de A^-*Jcb
428 Poly.matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
429 MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
432 displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
435 for(int k=0; k<_nVar*_nVar;k++){
436 _AroeMinusImplicit[k] *= -1;
438 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
440 displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
442 //calcul et insertion de D*JcbDiff
443 Poly.matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
444 MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
445 for(int k=0; k<_nVar*_nVar;k++)
447 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
450 } else if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
451 // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
452 Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
453 Ctemp2 = _mesh.getCell(idCells[1]);
455 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
456 if (j == Ctemp1.getFacesId()[l]){
457 for (int idim = 0; idim < _Ndim; ++idim)
458 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
462 }else{ // _Ndim = 1, build normal vector (bug cdmath)
463 if(!_sectionFieldSet)
465 if (Fj.x()<Ctemp1.x())
472 if(idCells[0]>idCells[1])
478 if(_verbose && _nbTimeStep%_freqSave ==0)
480 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
481 cout<<" Normal vector= ";
482 for (int idim = 0; idim < _Ndim; ++idim)
483 cout<<_vec_normal[idim]<<", ";
486 _porosityi=_porosityField(idCells[0]);
487 _porosityj=_porosityField(idCells[1]);
488 convectionState(idCells[0],idCells[1],false);
489 convectionMatrices();
490 diffusionStateAndMatrices(idCells[0], idCells[1], false);
492 // compute 1/dxi and 1/dxj
495 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
496 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
500 _inv_dxi = 1/Ctemp1.getMeasure();
501 _inv_dxj = 1/Ctemp2.getMeasure();
504 addConvectionToSecondMember( idCells[0],idCells[1], false);
505 addDiffusionToSecondMember( idCells[0],idCells[1], false);
506 addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
508 if(_timeScheme == Implicit){
509 for(int k=0; k<_nVar*_nVar;k++)
511 _AroeMinusImplicit[k] *= _inv_dxi;
512 _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);
516 //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNb<<endl;
517 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
518 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
521 displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
522 displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
524 for(int k=0;k<_nVar*_nVar;k++){
525 _AroeMinusImplicit[k] *= -1;
528 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
529 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
531 displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
532 displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
534 for(int k=0; k<_nVar*_nVar;k++)
536 _AroePlusImplicit[k] *= _inv_dxj;
537 _Diffusion[k] *=_inv_dxj/_inv_dxi;
539 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
540 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
542 displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
544 for(int k=0;k<_nVar*_nVar;k++){
545 _AroePlusImplicit[k] *= -1;
548 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
549 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
552 displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
555 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
556 if(_verbose && _nbTimeStep%_freqSave ==0)
557 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
558 *_runLogFile<<"Warning: treatment of a junction node"<<endl;
560 if(!_sectionFieldSet)
562 _runLogFile->close();
563 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
565 int largestSectionCellIndex=0;
566 for(int i=1;i<Fj.getNumberOfCells();i++){
567 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
568 largestSectionCellIndex=i;
570 idm = idCells[largestSectionCellIndex];
571 Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
572 _porosityi=_porosityField(idm);
574 if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
578 if(_verbose && _nbTimeStep%_freqSave ==0)
580 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
581 cout << " ; vecteur normal=(";
582 for(int p=0; p<_Ndim; p++)
583 cout << _vec_normal[p] << ",";
586 for(int i=0;i<Fj.getNumberOfCells();i++){
587 if(i != largestSectionCellIndex){
589 Ctemp2 = _mesh.getCell(idn);
590 _porosityj=_porosityField(idn);
591 convectionState(idm,idn,false);
592 convectionMatrices();
593 diffusionStateAndMatrices(idm, idn,false);
595 if(_verbose && _nbTimeStep%_freqSave ==0)
596 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
598 _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
599 _inv_dxj = 1/Ctemp2.getMeasure();
601 addConvectionToSecondMember(idm,idn, false);
602 _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
603 addDiffusionToSecondMember(idm,idn, false);
604 _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
605 addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
607 if(_timeScheme == Implicit){
608 for(int k=0; k<_nVar*_nVar;k++)
610 _AroeMinusImplicit[k] *= _inv_dxi;
611 _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
613 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
614 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
617 displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
618 displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
620 for(int k=0;k<_nVar*_nVar;k++){
621 _AroeMinusImplicit[k] *= -1;
624 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
625 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
627 displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
628 displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
630 for(int k=0; k<_nVar*_nVar;k++)
632 _AroePlusImplicit[k] *= _inv_dxj;
633 _Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
635 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
636 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
638 displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
640 for(int k=0;k<_nVar*_nVar;k++){
641 _AroePlusImplicit[k] *= -1;
644 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
645 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
648 displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
655 cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
656 _runLogFile->close();
657 throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
661 VecAssemblyEnd(_conservativeVars);
664 if(_timeScheme == Implicit){
665 for(int imaille = 0; imaille<_Nmailles; imaille++)
666 MatSetValuesBlocked(_A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
668 if(_verbose && _nbTimeStep%_freqSave ==0)
669 displayMatrix(_GravityImplicitationMatrix,_nVar,"Gravity matrix:");
671 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
672 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
673 if(_verbose && _nbTimeStep%_freqSave ==0){
674 cout << "Fin ProblemFluid.cxx : matrice implicite"<<endl;
675 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
684 if(_nbTimeStep+1<_cfl)
685 return (_nbTimeStep+1)*_minl/_maxvp;
688 return _cfl*_minl/_maxvp;
691 void ProblemFluid::computeNewtonVariation()
695 cout<<"Vecteur courant Uk "<<endl;
696 VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
699 if(_timeScheme == Explicit)
701 VecCopy(_b,_newtonVariation);
702 VecScale(_newtonVariation, _dt);
703 if(_verbose && _nbTimeStep%_freqSave ==0)
705 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
706 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
714 cout << "Matrice du système linéaire avant contribution delta t" << endl;
715 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
717 cout << "Second membre du système linéaire avant contribution delta t" << endl;
718 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
721 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
722 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
724 VecAXPY(_b, 1/_dt, _old);
725 VecAXPY(_b, -1/_dt, _conservativeVars);
728 #if PETSC_VERSION_GREATER_3_5
729 KSPSetOperators(_ksp, _A, _A);
731 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
736 cout << "Matrice du système linéaire après contribution delta t" << endl;
737 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
739 cout << "Second membre du système linéaire après contribution delta t" << endl;
740 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
745 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
748 KSPSolve(_ksp, _b, _newtonVariation);
752 computeScaling(_maxvp);
754 VecAssemblyBegin(_vecScaling);
755 VecAssemblyBegin(_invVecScaling);
756 for(int imaille = 0; imaille<_Nmailles; imaille++)
759 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
760 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
762 VecAssemblyEnd(_vecScaling);
763 VecAssemblyEnd(_invVecScaling);
767 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
768 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
770 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
771 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
774 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
777 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
778 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
781 VecCopy(_b,_bScaling);
782 VecPointwiseMult(_b,_vecScaling,_bScaling);
785 cout << "Produit du second membre par le preconditionneur bloc diagonal a gauche" << endl;
786 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
790 KSPSolve(_ksp,_b, _bScaling);
791 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
795 cout << "solution du systeme lineaire local:" << endl;
796 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
802 void ProblemFluid::validateTimeStep()
806 cout<<" Vecteur Un"<<endl;
807 VecView(_old, PETSC_VIEWER_STDOUT_WORLD);
808 cout<<" Vecteur Un+1"<<endl;
809 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_WORLD);
811 VecAXPY(_old, -1, _conservativeVars);//old contient old-courant
813 //Calcul de la variation Un+1-Un
818 for(int j=1; j<=_Nmailles; j++)
820 for(int k=0; k<_nVar; k++)
823 VecGetValues(_old, 1, &I, &dx);
824 VecGetValues(_conservativeVars, 1, &I, &x);
825 if (fabs(x)< _precision)
827 if(_erreur_rel < fabs(dx))
828 _erreur_rel = fabs(dx);
830 else if(_erreur_rel < fabs(dx/x))
831 _erreur_rel = fabs(dx/x);
835 _isStationary =_erreur_rel <_precision;
837 VecCopy(_conservativeVars, _old);
839 if(_verbose && _nbTimeStep%_freqSave ==0){
840 if(!_usePrimitiveVarsInNewton)
842 cout <<"Valeur propre locale max: " << _maxvp << endl;
845 if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
846 //Find minimum and maximum void fractions
847 double alphamin=1e30;
848 double alphamax=-1e30;
849 double T, Tmax=-1e30;
852 for(int j=0; j<_Nmailles; j++)
854 VecGetValues(_primitiveVars, 1, &I, &x);//extract void fraction
860 VecGetValues(_primitiveVars, 1, &J, &T);//extract void fraction
865 cout<<"Alpha min = " << alphamin << ", Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
867 *_runLogFile<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
872 if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
876 void ProblemFluid::abortTimeStep(){
880 void ProblemFluid::addConvectionToSecondMember
882 const int &j, bool isBord, string groupname
885 if(_verbose && _nbTimeStep%_freqSave ==0)
886 cout<<"ProblemFluid::addConvectionToSecondMember start"<<endl;
888 //extraction des valeurs
889 for(int k=0; k<_nVar; k++)
890 _idm[k] = _nVar*i + k;
891 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
892 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
895 for(int k=0; k<_nVar; k++)
896 _idn[k] = _nVar*j + k;
897 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
898 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
901 for(int k=0; k<_nVar; k++)
903 VecGetValues(_Uext, _nVar, _idn, _Uj);
904 consToPrim(_Uj, _Vj,_porosityj);
909 if(_verbose && _nbTimeStep%_freqSave ==0)
911 cout << "addConvectionToSecondMember : état i= " << i << ", _Ui=" << endl;
912 for(int q=0; q<_nVar; q++)
913 cout << _Ui[q] << endl;
915 cout << "addConvectionToSecondMember : état j= " << j << ", _Uj=" << endl;
916 for(int q=0; q<_nVar; q++)
917 cout << _Uj[q] << endl;
920 if(_nonLinearFormulation!=reducedRoe){//VFRoe, Roe or VFFC
921 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar),Fij(_nVar);
922 for(int i1=0;i1<_nVar;i1++)
929 Vector normale(_Ndim);
930 for(int i1=0;i1<_Ndim;i1++)
931 normale(i1)=_vec_normal[i1];
933 Matrix signAroe(_nVar);
934 for(int i1=0;i1<_nVar;i1++)
935 for(int i2=0;i2<_nVar;i2++)
936 signAroe(i1,i2)=_signAroe[i1*_nVar+i2];
938 Matrix absAroe(_nVar);
939 for(int i1=0;i1<_nVar;i1++)
940 for(int i2=0;i2<_nVar;i2++)
941 absAroe(i1,i2)=_absAroe[i1*_nVar+i2];
943 if(_nonLinearFormulation==VFRoe)//VFRoe
945 Vector Uij(_nVar),Vij(_nVar);
946 double porosityij=sqrt(_porosityi*_porosityj);
948 Uij=(Ui+Uj)/2+signAroe*(Ui-Uj)/2;
950 for(int i1=0;i1<_nVar;i1++)
953 consToPrim(_Uij, _Vij,porosityij);
955 applyVFRoeLowMachCorrections(isBord, groupname);
957 for(int i1=0;i1<_nVar;i1++)
963 Fij=convectionFlux(Uij,Vij,normale,porosityij);
965 if(_verbose && _nbTimeStep%_freqSave ==0)
967 cout<<"Etat interfacial conservatif "<<i<<", "<<j<< endl;
969 cout<<"Etat interfacial primitif "<<i<<", "<<j<< endl;
975 if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
976 Fij=staggeredVFFCFlux();
979 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
980 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
981 if(_nonLinearFormulation==VFFC)//VFFC
982 Fij=(Fi+Fj)/2+signAroe*(Fi-Fj)/2;
983 else if(_nonLinearFormulation==Roe)//Roe
984 Fij=(Fi+Fj)/2+absAroe*(Ui-Uj)/2;
986 if(_verbose && _nbTimeStep%_freqSave ==0)
988 cout<<"Flux cellule "<<i<<" = "<< endl;
990 cout<<"Flux cellule "<<j<<" = "<< endl;
995 for(int i1=0;i1<_nVar;i1++)
996 _phi[i1]=-Fij(i1)*_inv_dxi;
997 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
998 if(_verbose && _nbTimeStep%_freqSave ==0)
1000 cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1001 cout<<"Flux interfacial "<<i<<", "<<j<< endl;
1003 cout << "Contribution convection à " << i << ", -Fij(i1)*_inv_dxi= "<<endl;
1004 for(int q=0; q<_nVar; q++)
1005 cout << _phi[q] << endl;
1009 for(int i1=0;i1<_nVar;i1++)
1010 _phi[i1]*=-_inv_dxj/_inv_dxi;
1011 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1012 if(_verbose && _nbTimeStep%_freqSave ==0)
1014 cout << "Contribution convection à " << j << ", Fij(i1)*_inv_dxj= "<<endl;
1015 for(int q=0; q<_nVar; q++)
1016 cout << _phi[q] << endl;
1020 }else //_nonLinearFormulation==reducedRoe)
1022 for(int k=0; k<_nVar; k++)
1023 _temp[k]=(_Ui[k] - _Uj[k])*_inv_dxi;//(Ui-Uj)*_inv_dxi
1025 Poly.matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
1026 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1028 if(_verbose && _nbTimeStep%_freqSave ==0)
1030 cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1031 cout << "(Ui - Uj)*_inv_dxi= "<<endl;;
1032 for(int q=0; q<_nVar; q++)
1033 cout << _temp[q] << endl;
1035 cout << "Contribution convection à " << i << ", A^-*(Ui - Uj)*_inv_dxi= "<<endl;
1036 for(int q=0; q<_nVar; q++)
1037 cout << _phi[q] << endl;
1043 for(int k=0; k<_nVar; k++)
1044 _temp[k]*=_inv_dxj/_inv_dxi;//(Ui-Uj)*_inv_dxj
1045 Poly.matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
1046 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1048 if(_verbose && _nbTimeStep%_freqSave ==0)
1050 cout << "Contribution convection à " << j << ", A^+*(Ui - Uj)*_inv_dxi= "<<endl;
1051 for(int q=0; q<_nVar; q++)
1052 cout << _phi[q] << endl;
1057 if(_verbose && _nbTimeStep%_freqSave ==0)
1059 cout<<"ProblemFluid::addConvectionToSecondMember end : matrices de décentrement cellules i= " << i << ", et j= " << j<< "):"<<endl;
1060 displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
1061 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
1062 displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
1063 displayMatrix(_signAroe, _nVar,"Signe de la matrice de Roe");
1067 void ProblemFluid::addSourceTermToSecondMember
1068 ( const int i, int nbVoisinsi,
1069 const int j, int nbVoisinsj,
1070 bool isBord, int ij, double mesureFace)//To do : generalise to unstructured meshes
1072 if(_verbose && _nbTimeStep%_freqSave ==0)
1073 cout<<"ProblemFluid::addSourceTerm cell i= "<<i<< " cell j= "<< j<< " isbord "<<isBord<<endl;
1076 for(int k=1; k<_nVar;k++)
1077 _idm[k] = _idm[k-1] + 1;
1078 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1079 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1080 sourceVector(_Si,_Ui,_Vi,i);
1082 if (_verbose && _nbTimeStep%_freqSave ==0)
1084 cout << "Terme source S(Ui), i= " << i<<endl;
1085 for(int q=0; q<_nVar; q++)
1086 cout << _Si[q] << endl;
1090 for(int k=0; k<_nVar; k++)
1091 _idn[k] = _nVar*j + k;
1092 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1093 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1094 sourceVector(_Sj,_Uj,_Vj,j);
1097 for(int k=0; k<_nVar; k++)
1099 VecGetValues(_Uext, _nVar, _idn, _Uj);
1100 consToPrim(_Uj, _Vj,_porosityj);
1101 sourceVector(_Sj,_Uj,_Vj,i);
1103 if (_verbose && _nbTimeStep%_freqSave ==0)
1106 cout << "Terme source S(Uj), j= " << j<<endl;
1108 cout << "Terme source S(Uj), cellule fantôme "<<endl;
1109 for(int q=0; q<_nVar; q++)
1110 cout << _Sj[q] << endl;
1114 //Compute pressure loss vector
1116 if(_pressureLossFieldSet){
1117 K=_pressureLossField(ij);
1118 pressureLossVector(_pressureLossVector, K,_Ui, _Vi, _Uj, _Vj);
1122 for(int k=0; k<_nVar;k++)
1123 _pressureLossVector[k]=0;
1125 if (_verbose && _nbTimeStep%_freqSave ==0)
1127 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", K="<<K<<endl;
1128 for(int k=0; k<_nVar;k++)
1129 cout<< _pressureLossVector[k]<<", ";
1132 //Contribution of the porosityField gradient:
1134 porosityGradientSourceVector();
1136 for(int k=0; k<_nVar;k++)
1137 _porosityGradientSourceVector[k]=0;
1140 if (_verbose && _nbTimeStep%_freqSave ==0)
1143 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<", dxj= "<<1/_inv_dxj<<endl;
1145 cout<<"interface frontière i= "<<i<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<endl;
1146 cout<<"Gradient de porosite à l'interface"<<endl;
1147 for(int k=0; k<_nVar;k++)
1148 cout<< _porosityGradientSourceVector[k]<<", ";
1153 if(_wellBalancedCorrection){
1154 for(int k=0; k<_nVar;k++)
1155 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1156 Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1157 for(int k=0; k<_nVar;k++){
1158 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1159 _Sj[k]=(_phi[k]+_l[k])*mesureFace/_perimeters(j);///nbVoisinsj;
1161 if (_verbose && _nbTimeStep%_freqSave ==0)
1163 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant (après décentrement) de la face (i,j), j="<<j<<endl;
1164 for(int q=0; q<_nVar; q++)
1165 cout << _Si[q] << endl;
1166 cout << "Contribution au terme source Sj de la cellule j= " << j<<" venant (après décentrement) de la face (i,j), i="<<i<<endl;
1167 for(int q=0; q<_nVar; q++)
1168 cout << _Sj[q] << endl;
1170 cout<<"ratio surface sur volume i = "<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1171 cout<<"ratio surface sur volume j = "<<mesureFace/_perimeters(j)<<" perimeter = "<< _perimeters(j)<<endl;
1176 for(int k=0; k<_nVar;k++){
1177 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i)
1178 _Sj[k]=_Sj[k]/nbVoisinsj+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(j)
1180 if (_verbose && _nbTimeStep%_freqSave ==0)
1182 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,j), j="<<j<<endl;
1183 for(int q=0; q<_nVar; q++)
1184 cout << _Si[q] << endl;
1185 cout << "Contribution au terme source Sj de la cellule j = " << j<<" venant de la face (i,j), i="<<i <<endl;
1186 for(int q=0; q<_nVar; q++)
1187 cout << _Sj[q] << endl;
1192 VecSetValuesBlocked(_b, 1, _idn, _Sj, ADD_VALUES);
1194 if(_wellBalancedCorrection){
1195 for(int k=0; k<_nVar;k++)
1196 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1197 Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1198 for(int k=0; k<_nVar;k++)
1199 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1200 if (_verbose && _nbTimeStep%_freqSave ==0)
1202 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant (après décentrement) de la face (i,bord)"<<endl;
1203 for(int q=0; q<_nVar; q++)
1204 cout << _Si[q] << endl;
1205 cout<<"ratio surface sur volume i ="<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1211 for(int k=0; k<_nVar;k++)
1212 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i);//
1213 if (_verbose && _nbTimeStep%_freqSave ==0)
1215 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,bord) "<<endl;
1216 for(int q=0; q<_nVar; q++)
1217 cout << _Si[q] << endl;
1223 VecSetValuesBlocked(_b, 1, _idm, _Si, ADD_VALUES);
1225 if(_verbose && _nbTimeStep%_freqSave ==0 && _wellBalancedCorrection)
1226 displayMatrix( _signAroe,_nVar,"Signe matrice de Roe");
1229 void ProblemFluid::updatePrimitives()
1231 VecAssemblyBegin(_primitiveVars);
1232 for(int i=1; i<=_Nmailles; i++)
1234 _idm[0] = _nVar*( (i-1));
1235 for(int k=1; k<_nVar; k++)
1236 _idm[k] = _idm[k-1] + 1;
1238 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1239 if(_verbose && _nbTimeStep%_freqSave ==0)
1241 cout << "ProblemFluid::updatePrimitives() cell " << i-1 << endl;
1243 for(int q=0; q<_nVar; q++)
1244 cout << _Ui[q] << "\t";
1247 for(int q=0; q<_nVar; q++)
1248 cout << _Vi[q] << "\t";
1252 consToPrim(_Ui,_Vi,_porosityField(i-1));
1253 if(_nbPhases==2 && _Psat>-1e30){//Cas simulation flashing
1255 VecGetValues(_primitiveVars, 1, _idm+1, &pressure);
1256 _dp_over_dt(i-1)=(_Vi[1]-pressure)/_dt;//pn+1-pn
1259 VecSetValuesBlocked(_primitiveVars, 1, _idm, _Vi, INSERT_VALUES);
1261 VecAssemblyEnd(_primitiveVars);
1265 cout << "Nouvelles variables primitives : " << endl;
1266 VecView(_primitiveVars, PETSC_VIEWER_STDOUT_WORLD);
1271 void ProblemFluid::updateConservatives()
1273 VecAssemblyBegin(_conservativeVars);
1274 for(int i=1; i<=_Nmailles; i++)
1276 _idm[0] = _nVar*( (i-1));
1277 for(int k=1; k<_nVar; k++)
1278 _idm[k] = _idm[k-1] + 1;
1280 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1282 primToCons(_Vi,0,_Ui,0);
1284 VecSetValuesBlocked(_conservativeVars, 1, _idm, _Ui, INSERT_VALUES);
1286 if(_verbose && _nbTimeStep%_freqSave ==0)
1288 cout << "ProblemFluid::updateConservatives() cell " << i-1 << endl;
1290 for(int q=0; q<_nVar; q++)
1291 cout << _Vi[q] << "\t";
1294 for(int q=0; q<_nVar; q++)
1295 cout << _Ui[q] << "\t";
1299 VecAssemblyEnd(_conservativeVars);
1303 cout << "Nouvelles variables conservatives : " << endl;
1304 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_WORLD);
1309 vector< complex<double> > ProblemFluid::getRacines(vector< double > pol_car){
1310 int degre_polynom = pol_car.size()-1;
1312 vector< complex< double > > valeurs_propres (degre_polynom);
1313 vector< double > zeror(degre_polynom);
1314 vector< double > zeroi(degre_polynom);
1315 for(int j=0; j<(degre_polynom+1)/2; j++){//coefficients in order of decreasing powers for rpoly
1317 pol_car[j] =pol_car[degre_polynom-j];
1318 pol_car[degre_polynom-j]=tmp;
1321 //Calcul des racines du polynome
1322 roots_polynoms roots;
1323 int size_vp = roots.rpoly(&pol_car[0],degre_polynom,&zeror[0],&zeroi[0]);
1325 //On ordonne les valeurs propres
1326 if(zeror[1]<zeror[0])
1336 if(size_vp == degre_polynom) {
1337 for(int ct =0; ct<degre_polynom; ct++)
1338 valeurs_propres[ct] = complex<double> (zeror[ct],zeroi[ct]); //vp non triviales
1341 cout << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1342 *_runLogFile << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1343 cout<<"Coefficients polynome caracteristique: "<<endl;
1344 for(int ct =0; ct<degre_polynom+1; ct++)
1345 cout<<pol_car[ct]<< " , " <<endl;
1347 *_runLogFile<<"getRacines computation failed"<<endl;
1348 _runLogFile->close();
1349 throw CdmathException("getRacines computation failed");
1352 return valeurs_propres;
1355 void ProblemFluid::AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1358 int nbVp_dist=valeurs_propres_dist.size();
1359 vector< complex< double > > y (nbVp_dist,0);
1360 for( int i=0 ; i<nbVp_dist ; i++)
1361 y[i] = Poly.abs_generalise(valeurs_propres_dist[i]);
1362 Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
1363 if(_verbose && _nbTimeStep%_freqSave ==0)
1365 cout<< endl<<"ProblemFluid::AbsMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1366 for(int ct =0; ct<nbVp_dist; ct++)
1367 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
1369 cout<<"ProblemFluid::AbsMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1370 for(int ct =0; ct<nbVp_dist; ct++)
1371 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
1373 displayMatrix(_absAroe,_nVar,"Valeur absolue de la matrice de Roe");
1377 void ProblemFluid::SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1379 int nbVp_dist=valeurs_propres_dist.size();
1380 vector< complex< double > > y (nbVp_dist,0);
1381 for( int i=0 ; i<nbVp_dist ; i++)
1383 if(valeurs_propres_dist[i].real()>0)
1385 else if(valeurs_propres_dist[i].real()<0)
1391 Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
1392 if(_verbose && _nbTimeStep%_freqSave ==0)
1394 cout<< endl<<"ProblemFluid::SigneMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1395 for(int ct =0; ct<nbVp_dist; ct++)
1396 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
1398 cout<<"ProblemFluid::SigneMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1399 for(int ct =0; ct<nbVp_dist; ct++)
1400 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
1402 displayMatrix(_signAroe,_nVar,"Signe matrice de Roe");
1405 void ProblemFluid::InvMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1407 int nbVp_dist=valeurs_propres_dist.size();
1408 vector< complex< double > > y (nbVp_dist,0);
1410 for( int i=0 ; i<nbVp_dist ; i++)
1412 if(Poly.module(valeurs_propres_dist[i])>_precision)
1413 y[i] = 1./valeurs_propres_dist[i];
1415 y[i] = 1./_precision;
1417 Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
1418 if(_verbose && _nbTimeStep%_freqSave ==0)
1420 cout<< endl<<"ProblemFluid::InvMatriceRoe : Valeurs propres :" << nbVp_dist<<endl;
1421 for(int ct =0; ct<nbVp_dist; ct++)
1422 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
1424 cout<<"ProblemFluid::InvMatriceRoe : Valeurs à atteindre pour l'interpolation"<<endl;
1425 for(int ct =0; ct<nbVp_dist; ct++)
1426 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
1428 displayMatrix(_invAroe,_nVar,"Inverse matrice de Roe");
1432 void ProblemFluid::terminate(){
1435 delete[] _Diffusion;
1436 delete[] _GravityImplicitationMatrix;
1437 delete[] _AroeMinus;
1442 delete[] _AroeImplicit;
1443 delete[] _AroeMinusImplicit;
1444 delete[] _AroePlusImplicit;
1445 delete[] _absAroeImplicit;
1450 delete[] _primToConsJacoMat;
1457 delete[] _externalStates;
1460 delete[] _vec_normal;
1465 if(_nonLinearFormulation==VFRoe){
1471 delete[] _pressureLossVector;
1472 delete[] _porosityGradientSourceVector;
1475 delete[] _blockDiag;
1476 delete[] _invBlockDiag;
1478 VecDestroy(&_vecScaling);
1479 VecDestroy(&_invVecScaling);
1480 VecDestroy(&_bScaling);
1483 VecDestroy(&_conservativeVars);
1484 VecDestroy(&_newtonVariation);
1486 VecDestroy(&_primitiveVars);
1489 VecDestroy(&_Uextdiff);
1493 for(int i=0;i<_nbPhases;i++)