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() const
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*deltaU
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);
683 if(_nbTimeStep+1<_cfl)
684 return (_nbTimeStep+1)*_minl/_maxvp;
687 return _cfl*_minl/_maxvp;
690 void ProblemFluid::computeNewtonVariation()
694 cout<<"Vecteur courant Uk "<<endl;
695 VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
698 if(_timeScheme == Explicit)
700 VecCopy(_b,_newtonVariation);
701 VecScale(_newtonVariation, _dt);
702 if(_verbose && _nbTimeStep%_freqSave ==0)
704 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
705 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
713 cout << "Matrice du système linéaire avant contribution delta t" << endl;
714 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
716 cout << "Second membre du système linéaire avant contribution delta t" << endl;
717 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
720 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
721 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
723 VecAXPY(_b, 1/_dt, _old);
724 VecAXPY(_b, -1/_dt, _conservativeVars);
727 #if PETSC_VERSION_GREATER_3_5
728 KSPSetOperators(_ksp, _A, _A);
730 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
735 cout << "Matrice du système linéaire après contribution delta t" << endl;
736 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
738 cout << "Second membre du système linéaire après contribution delta t" << endl;
739 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
744 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
747 KSPSolve(_ksp, _b, _newtonVariation);
751 computeScaling(_maxvp);
753 VecAssemblyBegin(_vecScaling);
754 VecAssemblyBegin(_invVecScaling);
755 for(int imaille = 0; imaille<_Nmailles; imaille++)
758 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
759 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
761 VecAssemblyEnd(_vecScaling);
762 VecAssemblyEnd(_invVecScaling);
766 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
767 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
769 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
770 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
773 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
776 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
777 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
780 VecCopy(_b,_bScaling);
781 VecPointwiseMult(_b,_vecScaling,_bScaling);
784 cout << "Produit du second membre par le preconditionneur bloc diagonal a gauche" << endl;
785 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
789 KSPSolve(_ksp,_b, _bScaling);
790 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
794 cout << "solution du systeme lineaire local:" << endl;
795 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
801 void ProblemFluid::validateTimeStep()
805 cout<<" Vecteur Un"<<endl;
806 VecView(_old, PETSC_VIEWER_STDOUT_WORLD);
807 cout<<" Vecteur Un+1"<<endl;
808 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_WORLD);
810 VecAXPY(_old, -1, _conservativeVars);//old contient old-courant
812 //Calcul de la variation Un+1-Un
817 for(int j=1; j<=_Nmailles; j++)
819 for(int k=0; k<_nVar; k++)
822 VecGetValues(_old, 1, &I, &dx);
823 VecGetValues(_conservativeVars, 1, &I, &x);
824 if (fabs(x)< _precision)
826 if(_erreur_rel < fabs(dx))
827 _erreur_rel = fabs(dx);
829 else if(_erreur_rel < fabs(dx/x))
830 _erreur_rel = fabs(dx/x);
834 _isStationary =_erreur_rel <_precision;
836 VecCopy(_conservativeVars, _old);
838 if(_verbose && _nbTimeStep%_freqSave ==0){
839 if(!_usePrimitiveVarsInNewton)
841 cout <<"Valeur propre locale max: " << _maxvp << endl;
844 if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
845 //Find minimum and maximum void fractions
846 double alphamin=1e30;
847 double alphamax=-1e30;
848 double T, Tmax=-1e30;
851 for(int j=0; j<_Nmailles; j++)
853 VecGetValues(_primitiveVars, 1, &I, &x);//extract void fraction
859 VecGetValues(_primitiveVars, 1, &J, &T);//extract void fraction
864 cout<<"Alpha min = " << alphamin << ", Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
866 *_runLogFile<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
871 if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
875 void ProblemFluid::abortTimeStep(){
879 void ProblemFluid::addConvectionToSecondMember
881 const int &j, bool isBord, string groupname
884 if(_verbose && _nbTimeStep%_freqSave ==0)
885 cout<<"ProblemFluid::addConvectionToSecondMember start"<<endl;
887 //extraction des valeurs
888 for(int k=0; k<_nVar; k++)
889 _idm[k] = _nVar*i + k;
890 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
891 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
894 for(int k=0; k<_nVar; k++)
895 _idn[k] = _nVar*j + k;
896 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
897 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
900 for(int k=0; k<_nVar; k++)
902 VecGetValues(_Uext, _nVar, _idn, _Uj);
903 consToPrim(_Uj, _Vj,_porosityj);
908 if(_verbose && _nbTimeStep%_freqSave ==0)
910 cout << "addConvectionToSecondMember : état i= " << i << ", _Ui=" << endl;
911 for(int q=0; q<_nVar; q++)
912 cout << _Ui[q] << endl;
914 cout << "addConvectionToSecondMember : état j= " << j << ", _Uj=" << endl;
915 for(int q=0; q<_nVar; q++)
916 cout << _Uj[q] << endl;
919 if(_nonLinearFormulation!=reducedRoe){//VFRoe, Roe or VFFC
920 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar),Fij(_nVar);
921 for(int i1=0;i1<_nVar;i1++)
928 Vector normale(_Ndim);
929 for(int i1=0;i1<_Ndim;i1++)
930 normale(i1)=_vec_normal[i1];
932 Matrix signAroe(_nVar);
933 for(int i1=0;i1<_nVar;i1++)
934 for(int i2=0;i2<_nVar;i2++)
935 signAroe(i1,i2)=_signAroe[i1*_nVar+i2];
937 Matrix absAroe(_nVar);
938 for(int i1=0;i1<_nVar;i1++)
939 for(int i2=0;i2<_nVar;i2++)
940 absAroe(i1,i2)=_absAroe[i1*_nVar+i2];
942 if(_nonLinearFormulation==VFRoe)//VFRoe
944 Vector Uij(_nVar),Vij(_nVar);
945 double porosityij=sqrt(_porosityi*_porosityj);
947 Uij=(Ui+Uj)/2+signAroe*(Ui-Uj)/2;
949 for(int i1=0;i1<_nVar;i1++)
952 consToPrim(_Uij, _Vij,porosityij);
954 applyVFRoeLowMachCorrections(isBord, groupname);
956 for(int i1=0;i1<_nVar;i1++)
962 Fij=convectionFlux(Uij,Vij,normale,porosityij);
964 if(_verbose && _nbTimeStep%_freqSave ==0)
966 cout<<"Etat interfacial conservatif "<<i<<", "<<j<< endl;
968 cout<<"Etat interfacial primitif "<<i<<", "<<j<< endl;
974 if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
975 Fij=staggeredVFFCFlux();
978 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
979 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
980 if(_nonLinearFormulation==VFFC)//VFFC
981 Fij=(Fi+Fj)/2+signAroe*(Fi-Fj)/2;
982 else if(_nonLinearFormulation==Roe)//Roe
983 Fij=(Fi+Fj)/2+absAroe*(Ui-Uj)/2;
985 if(_verbose && _nbTimeStep%_freqSave ==0)
987 cout<<"Flux cellule "<<i<<" = "<< endl;
989 cout<<"Flux cellule "<<j<<" = "<< endl;
994 for(int i1=0;i1<_nVar;i1++)
995 _phi[i1]=-Fij(i1)*_inv_dxi;
996 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
997 if(_verbose && _nbTimeStep%_freqSave ==0)
999 cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1000 cout<<"Flux interfacial "<<i<<", "<<j<< endl;
1002 cout << "Contribution convection à " << i << ", -Fij(i1)*_inv_dxi= "<<endl;
1003 for(int q=0; q<_nVar; q++)
1004 cout << _phi[q] << endl;
1008 for(int i1=0;i1<_nVar;i1++)
1009 _phi[i1]*=-_inv_dxj/_inv_dxi;
1010 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1011 if(_verbose && _nbTimeStep%_freqSave ==0)
1013 cout << "Contribution convection à " << j << ", Fij(i1)*_inv_dxj= "<<endl;
1014 for(int q=0; q<_nVar; q++)
1015 cout << _phi[q] << endl;
1019 }else //_nonLinearFormulation==reducedRoe)
1021 for(int k=0; k<_nVar; k++)
1022 _temp[k]=(_Ui[k] - _Uj[k])*_inv_dxi;//(Ui-Uj)*_inv_dxi
1024 Poly.matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
1025 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1027 if(_verbose && _nbTimeStep%_freqSave ==0)
1029 cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1030 cout << "(Ui - Uj)*_inv_dxi= "<<endl;;
1031 for(int q=0; q<_nVar; q++)
1032 cout << _temp[q] << endl;
1034 cout << "Contribution convection à " << i << ", A^-*(Ui - Uj)*_inv_dxi= "<<endl;
1035 for(int q=0; q<_nVar; q++)
1036 cout << _phi[q] << endl;
1042 for(int k=0; k<_nVar; k++)
1043 _temp[k]*=_inv_dxj/_inv_dxi;//(Ui-Uj)*_inv_dxj
1044 Poly.matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
1045 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1047 if(_verbose && _nbTimeStep%_freqSave ==0)
1049 cout << "Contribution convection à " << j << ", A^+*(Ui - Uj)*_inv_dxi= "<<endl;
1050 for(int q=0; q<_nVar; q++)
1051 cout << _phi[q] << endl;
1056 if(_verbose && _nbTimeStep%_freqSave ==0)
1058 cout<<"ProblemFluid::addConvectionToSecondMember end : matrices de décentrement cellules i= " << i << ", et j= " << j<< "):"<<endl;
1059 displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
1060 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
1061 displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
1062 displayMatrix(_signAroe, _nVar,"Signe de la matrice de Roe");
1066 void ProblemFluid::addSourceTermToSecondMember
1067 ( const int i, int nbVoisinsi,
1068 const int j, int nbVoisinsj,
1069 bool isBord, int ij, double mesureFace)//To do : generalise to unstructured meshes
1071 if(_verbose && _nbTimeStep%_freqSave ==0)
1072 cout<<"ProblemFluid::addSourceTerm cell i= "<<i<< " cell j= "<< j<< " isbord "<<isBord<<endl;
1075 for(int k=1; k<_nVar;k++)
1076 _idm[k] = _idm[k-1] + 1;
1077 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1078 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1079 sourceVector(_Si,_Ui,_Vi,i);
1081 if (_verbose && _nbTimeStep%_freqSave ==0)
1083 cout << "Terme source S(Ui), i= " << i<<endl;
1084 for(int q=0; q<_nVar; q++)
1085 cout << _Si[q] << endl;
1089 for(int k=0; k<_nVar; k++)
1090 _idn[k] = _nVar*j + k;
1091 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1092 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1093 sourceVector(_Sj,_Uj,_Vj,j);
1096 for(int k=0; k<_nVar; k++)
1098 VecGetValues(_Uext, _nVar, _idn, _Uj);
1099 consToPrim(_Uj, _Vj,_porosityj);
1100 sourceVector(_Sj,_Uj,_Vj,i);
1102 if (_verbose && _nbTimeStep%_freqSave ==0)
1105 cout << "Terme source S(Uj), j= " << j<<endl;
1107 cout << "Terme source S(Uj), cellule fantôme "<<endl;
1108 for(int q=0; q<_nVar; q++)
1109 cout << _Sj[q] << endl;
1113 //Compute pressure loss vector
1115 if(_pressureLossFieldSet){
1116 K=_pressureLossField(ij);
1117 pressureLossVector(_pressureLossVector, K,_Ui, _Vi, _Uj, _Vj);
1121 for(int k=0; k<_nVar;k++)
1122 _pressureLossVector[k]=0;
1124 if (_verbose && _nbTimeStep%_freqSave ==0)
1126 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", K="<<K<<endl;
1127 for(int k=0; k<_nVar;k++)
1128 cout<< _pressureLossVector[k]<<", ";
1131 //Contribution of the porosityField gradient:
1133 porosityGradientSourceVector();
1135 for(int k=0; k<_nVar;k++)
1136 _porosityGradientSourceVector[k]=0;
1139 if (_verbose && _nbTimeStep%_freqSave ==0)
1142 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<", dxj= "<<1/_inv_dxj<<endl;
1144 cout<<"interface frontière i= "<<i<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<endl;
1145 cout<<"Gradient de porosite à l'interface"<<endl;
1146 for(int k=0; k<_nVar;k++)
1147 cout<< _porosityGradientSourceVector[k]<<", ";
1152 if(_wellBalancedCorrection){
1153 for(int k=0; k<_nVar;k++)
1154 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1155 Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1156 for(int k=0; k<_nVar;k++){
1157 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1158 _Sj[k]=(_phi[k]+_l[k])*mesureFace/_perimeters(j);///nbVoisinsj;
1160 if (_verbose && _nbTimeStep%_freqSave ==0)
1162 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant (après décentrement) de la face (i,j), j="<<j<<endl;
1163 for(int q=0; q<_nVar; q++)
1164 cout << _Si[q] << endl;
1165 cout << "Contribution au terme source Sj de la cellule j= " << j<<" venant (après décentrement) de la face (i,j), i="<<i<<endl;
1166 for(int q=0; q<_nVar; q++)
1167 cout << _Sj[q] << endl;
1169 cout<<"ratio surface sur volume i = "<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1170 cout<<"ratio surface sur volume j = "<<mesureFace/_perimeters(j)<<" perimeter = "<< _perimeters(j)<<endl;
1175 for(int k=0; k<_nVar;k++){
1176 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i)
1177 _Sj[k]=_Sj[k]/nbVoisinsj+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(j)
1179 if (_verbose && _nbTimeStep%_freqSave ==0)
1181 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,j), j="<<j<<endl;
1182 for(int q=0; q<_nVar; q++)
1183 cout << _Si[q] << endl;
1184 cout << "Contribution au terme source Sj de la cellule j = " << j<<" venant de la face (i,j), i="<<i <<endl;
1185 for(int q=0; q<_nVar; q++)
1186 cout << _Sj[q] << endl;
1191 VecSetValuesBlocked(_b, 1, _idn, _Sj, ADD_VALUES);
1193 if(_wellBalancedCorrection){
1194 for(int k=0; k<_nVar;k++)
1195 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1196 Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1197 for(int k=0; k<_nVar;k++)
1198 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1199 if (_verbose && _nbTimeStep%_freqSave ==0)
1201 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant (après décentrement) de la face (i,bord)"<<endl;
1202 for(int q=0; q<_nVar; q++)
1203 cout << _Si[q] << endl;
1204 cout<<"ratio surface sur volume i ="<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1210 for(int k=0; k<_nVar;k++)
1211 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i);//
1212 if (_verbose && _nbTimeStep%_freqSave ==0)
1214 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,bord) "<<endl;
1215 for(int q=0; q<_nVar; q++)
1216 cout << _Si[q] << endl;
1222 VecSetValuesBlocked(_b, 1, _idm, _Si, ADD_VALUES);
1224 if(_verbose && _nbTimeStep%_freqSave ==0 && _wellBalancedCorrection)
1225 displayMatrix( _signAroe,_nVar,"Signe matrice de Roe");
1228 void ProblemFluid::updatePrimitives()
1230 VecAssemblyBegin(_primitiveVars);
1231 for(int i=1; i<=_Nmailles; i++)
1233 _idm[0] = _nVar*( (i-1));
1234 for(int k=1; k<_nVar; k++)
1235 _idm[k] = _idm[k-1] + 1;
1237 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1238 consToPrim(_Ui,_Vi,_porosityField(i-1));
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 if(_nbPhases==2 && _Psat>-1e30){//Cas simulation flashing
1254 VecGetValues(_primitiveVars, 1, _idm+1, &pressure);
1255 _dp_over_dt(i-1)=(_Vi[1]-pressure)/_dt;//pn+1-pn
1258 VecSetValuesBlocked(_primitiveVars, 1, _idm, _Vi, INSERT_VALUES);
1260 VecAssemblyEnd(_primitiveVars);
1264 cout << "Nouvelles variables primitives : " << endl;
1265 VecView(_primitiveVars, PETSC_VIEWER_STDOUT_WORLD);
1270 void ProblemFluid::updateConservatives()
1272 VecAssemblyBegin(_conservativeVars);
1273 for(int i=1; i<=_Nmailles; i++)
1275 _idm[0] = _nVar*( (i-1));
1276 for(int k=1; k<_nVar; k++)
1277 _idm[k] = _idm[k-1] + 1;
1279 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1281 primToCons(_Vi,0,_Ui,0);
1283 VecSetValuesBlocked(_conservativeVars, 1, _idm, _Ui, INSERT_VALUES);
1285 if(_verbose && _nbTimeStep%_freqSave ==0)
1287 cout << "ProblemFluid::updateConservatives() cell " << i-1 << endl;
1289 for(int q=0; q<_nVar; q++)
1290 cout << _Vi[q] << "\t";
1293 for(int q=0; q<_nVar; q++)
1294 cout << _Ui[q] << "\t";
1298 VecAssemblyEnd(_conservativeVars);
1302 cout << "Nouvelles variables conservatives : " << endl;
1303 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_WORLD);
1308 vector< complex<double> > ProblemFluid::getRacines(vector< double > pol_car){
1309 int degre_polynom = pol_car.size()-1;
1311 vector< complex< double > > valeurs_propres (degre_polynom);
1312 vector< double > zeror(degre_polynom);
1313 vector< double > zeroi(degre_polynom);
1314 for(int j=0; j<(degre_polynom+1)/2; j++){//coefficients in order of decreasing powers for rpoly
1316 pol_car[j] =pol_car[degre_polynom-j];
1317 pol_car[degre_polynom-j]=tmp;
1320 //Calcul des racines du polynome
1321 roots_polynoms roots;
1322 int size_vp = roots.rpoly(&pol_car[0],degre_polynom,&zeror[0],&zeroi[0]);
1324 //On ordonne les valeurs propres
1325 if(zeror[1]<zeror[0])
1335 if(size_vp == degre_polynom) {
1336 for(int ct =0; ct<degre_polynom; ct++)
1337 valeurs_propres[ct] = complex<double> (zeror[ct],zeroi[ct]); //vp non triviales
1340 cout << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1341 *_runLogFile << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
1342 cout<<"Coefficients polynome caracteristique: "<<endl;
1343 for(int ct =0; ct<degre_polynom+1; ct++)
1344 cout<<pol_car[ct]<< " , " <<endl;
1346 *_runLogFile<<"getRacines computation failed"<<endl;
1347 _runLogFile->close();
1348 throw CdmathException("getRacines computation failed");
1351 return valeurs_propres;
1354 void ProblemFluid::AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1357 int nbVp_dist=valeurs_propres_dist.size();
1358 vector< complex< double > > y (nbVp_dist,0);
1359 for( int i=0 ; i<nbVp_dist ; i++)
1360 y[i] = Poly.abs_generalise(valeurs_propres_dist[i]);
1361 Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
1362 if(_verbose && _nbTimeStep%_freqSave ==0)
1364 cout<< endl<<"ProblemFluid::AbsMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1365 for(int ct =0; ct<nbVp_dist; ct++)
1366 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
1368 cout<<"ProblemFluid::AbsMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1369 for(int ct =0; ct<nbVp_dist; ct++)
1370 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
1372 displayMatrix(_absAroe,_nVar,"Valeur absolue de la matrice de Roe");
1376 void ProblemFluid::SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1378 int nbVp_dist=valeurs_propres_dist.size();
1379 vector< complex< double > > y (nbVp_dist,0);
1380 for( int i=0 ; i<nbVp_dist ; i++)
1382 if(valeurs_propres_dist[i].real()>0)
1384 else if(valeurs_propres_dist[i].real()<0)
1390 Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
1391 if(_verbose && _nbTimeStep%_freqSave ==0)
1393 cout<< endl<<"ProblemFluid::SigneMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
1394 for(int ct =0; ct<nbVp_dist; ct++)
1395 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
1397 cout<<"ProblemFluid::SigneMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
1398 for(int ct =0; ct<nbVp_dist; ct++)
1399 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
1401 displayMatrix(_signAroe,_nVar,"Signe matrice de Roe");
1404 void ProblemFluid::InvMatriceRoe(vector< complex<double> > valeurs_propres_dist)
1406 int nbVp_dist=valeurs_propres_dist.size();
1407 vector< complex< double > > y (nbVp_dist,0);
1409 for( int i=0 ; i<nbVp_dist ; i++)
1411 if(Poly.module(valeurs_propres_dist[i])>_precision)
1412 y[i] = 1./valeurs_propres_dist[i];
1414 y[i] = 1./_precision;
1416 Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
1417 if(_verbose && _nbTimeStep%_freqSave ==0)
1419 cout<< endl<<"ProblemFluid::InvMatriceRoe : Valeurs propres :" << nbVp_dist<<endl;
1420 for(int ct =0; ct<nbVp_dist; ct++)
1421 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
1423 cout<<"ProblemFluid::InvMatriceRoe : Valeurs à atteindre pour l'interpolation"<<endl;
1424 for(int ct =0; ct<nbVp_dist; ct++)
1425 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
1427 displayMatrix(_invAroe,_nVar,"Inverse matrice de Roe");
1431 Field ProblemFluid::getConservativeField() const
1433 if(!_initializedMemory)
1435 _runLogFile->close();
1436 throw CdmathException("ProblemFluid::getConservativeField call initialize() first");
1441 Field ProblemFluid::getPrimitiveField() const
1443 if(!_initializedMemory)
1445 _runLogFile->close();
1446 throw CdmathException("ProblemFluid::getPrimitiveField call initialize() first");
1451 void ProblemFluid::terminate(){
1454 delete[] _Diffusion;
1455 delete[] _GravityImplicitationMatrix;
1456 delete[] _AroeMinus;
1461 delete[] _AroeImplicit;
1462 delete[] _AroeMinusImplicit;
1463 delete[] _AroePlusImplicit;
1464 delete[] _absAroeImplicit;
1469 delete[] _primToConsJacoMat;
1476 delete[] _externalStates;
1479 delete[] _vec_normal;
1484 if(_nonLinearFormulation==VFRoe){
1490 delete[] _pressureLossVector;
1491 delete[] _porosityGradientSourceVector;
1494 delete[] _blockDiag;
1495 delete[] _invBlockDiag;
1497 VecDestroy(&_vecScaling);
1498 VecDestroy(&_invVecScaling);
1499 VecDestroy(&_bScaling);
1502 VecDestroy(&_conservativeVars);
1503 VecDestroy(&_newtonVariation);
1505 VecDestroy(&_primitiveVars);
1508 VecDestroy(&_Uextdiff);
1512 for(int i=0;i<_nbPhases;i++)