4 * Created on: Sep 15, 2014
8 #include "SinglePhase.hxx"
12 SinglePhase::SinglePhase(phaseType fluid, pressureEstimate pEstimate, int dim, bool useDellacherieEOS){
16 _dragCoeffs=vector<double>(1,0);
18 _useDellacherieEOS=useDellacherieEOS;
21 if(pEstimate==around1bar300K){//EOS at 1 bar and 300K
25 cout<<"Fluid is air around 1 bar and 300 K (27°C)"<<endl;
26 *_runLogFile<<"Fluid is air around 1 bar and 300 K (27°C)"<<endl;
27 _fluides[0] = new StiffenedGas(1.4,743,_Tref,2.22e5); //ideal gas law for nitrogen at pressure 1 bar and temperature 27°C, e=2.22e5, c_v=743
30 cout<<"Fluid is water around 1 bar and 300 K (27°C)"<<endl;
31 *_runLogFile<<"Fluid is water around 1 bar and 300 K (27°C)"<<endl;
32 _fluides[0] = new StiffenedGas(996,_Pref,_Tref,1.12e5,1501,4130); //stiffened gas law for water at pressure 1 bar and temperature 27°C, e=1.12e5, c_v=4130
35 else{//EOS at 155 bars and 618K
39 cout<<"Fluid is Gas around saturation point 155 bars and 618 K (345°C)"<<endl;
40 *_runLogFile<<"Fluid is Gas around saturation point 155 bars and 618 K (345°C)"<<endl;
41 _fluides[0] = new StiffenedGas(102,_Pref,_Tref,2.44e6, 433,3633); //stiffened gas law for Gas at pressure 155 bar and temperature 345°C
43 else{//To do : change to normal regime: 155 bars and 573K
44 cout<<"Fluid is water around saturation point 155 bars and 618 K (345°C)"<<endl;
45 *_runLogFile<<"Fluid is water around saturation point 155 bars and 618 K (345°C)"<<endl;
46 if(_useDellacherieEOS)
47 _fluides[0]= new StiffenedGasDellacherie(2.35,1e9,-1.167e6,1816); //stiffened gas law for water from S. Dellacherie
49 _fluides[0]= new StiffenedGas(594.,_Pref,_Tref,1.6e6, 621.,3100.); //stiffened gas law for water at pressure 155 bar, and temperature 345°C
53 void SinglePhase::initialize(){
54 cout<<"Initialising the Navier-Stokes model"<<endl;
55 *_runLogFile<<"Initialising the Navier-Stokes model"<<endl;
57 _Uroe = new double[_nVar];
58 _gravite = vector<double>(_nVar,0);//Not to be confused with _GravityField3d (size _Ndim). _gravite (size _Nvar) is usefull for dealing with source term and implicitation of gravity vector
59 for(int i=0; i<_Ndim; i++)
60 _gravite[i+1]=_GravityField3d[i];
62 _GravityImplicitationMatrix = new PetscScalar[_nVar*_nVar];
63 if(_saveVelocity || _saveAllFields)
64 _Vitesse=Field("Velocity",CELLS,_mesh,3);//Forcement en dimension 3 pour le posttraitement des lignes de courant
68 _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
69 _Pressure=Field("Pressure",CELLS,_mesh,1);
70 _Density=Field("Density",CELLS,_mesh,1);
71 _Temperature=Field("Temperature",CELLS,_mesh,1);
72 _MachNumber=Field("MachNumber",CELLS,_mesh,1);
73 _VitesseX=Field("Velocity x",CELLS,_mesh,1);
76 _VitesseY=Field("Velocity y",CELLS,_mesh,1);
78 _VitesseZ=Field("Velocity z",CELLS,_mesh,1);
82 if(_entropicCorrection)
83 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
85 ProblemFluid::initialize();
88 bool SinglePhase::iterateTimeStep(bool &converged)
90 if(_timeScheme == Explicit || !_usePrimitiveVarsInNewton)
91 return ProblemFluid::iterateTimeStep(converged);
96 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
98 computeTimeStep(stop);
100 if(stop){//Le compute time step ne s'est pas bien passé
101 cout<<"ComputeTimeStep failed"<<endl;
106 computeNewtonVariation();
108 //converged=convergence des iterations
109 if(_timeScheme == Explicit)
111 else{//Implicit scheme
113 KSPGetIterationNumber(_ksp, &_PetscIts);
114 if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
115 _MaxIterLinearSolver = _PetscIts;
116 if(_PetscIts>=_maxPetscIts)//solving the linear system failed
118 cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
119 *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
123 else{//solving the linear system succeeded
124 //Calcul de la variation relative Uk+1-Uk
128 for(int j=1; j<=_Nmailles; j++)
130 for(int k=0; k<_nVar; k++)
133 VecGetValues(_newtonVariation, 1, &I, &dx);
134 VecGetValues(_primitiveVars, 1, &I, &x);
135 if (fabs(x)*fabs(x)< _precision)
137 if(_erreur_rel < fabs(dx))
138 _erreur_rel = fabs(dx);
140 else if(_erreur_rel < fabs(dx/x))
141 _erreur_rel = fabs(dx/x);
145 converged = _erreur_rel <= _precision_Newton;
148 double relaxation=1;//Vk+1=Vk+relaxation*deltaV
150 VecAXPY(_primitiveVars, relaxation, _newtonVariation);
152 //mise a jour du champ primitif
153 updateConservatives();
155 if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
157 cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
158 *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
164 cout<<"Vecteur Vkp1-Vk "<<endl;
165 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
166 cout << "Nouvel etat courant Vk de l'iteration Newton: " << endl;
167 VecView(_primitiveVars, PETSC_VIEWER_STDOUT_SELF);
170 if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
171 if(_minm1<-_precision || _minm2<-_precision)
173 cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
174 *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
178 cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
179 *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
191 void SinglePhase::computeNewtonVariation()
193 if(!_usePrimitiveVarsInNewton)
194 ProblemFluid::computeNewtonVariation();
199 cout<<"Vecteur courant Vk "<<endl;
200 VecView(_primitiveVars,PETSC_VIEWER_STDOUT_SELF);
202 cout << "Matrice du système linéaire avant contribution delta t" << endl;
203 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
205 cout << "Second membre du système linéaire avant contribution delta t" << endl;
206 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
209 if(_timeScheme == Explicit)
211 VecCopy(_b,_newtonVariation);
212 VecScale(_newtonVariation, _dt);
213 if(_verbose && _nbTimeStep%_freqSave ==0)
215 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
216 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
222 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
224 VecAXPY(_b, 1/_dt, _old);
225 VecAXPY(_b, -1/_dt, _conservativeVars);
227 for(int imaille = 0; imaille<_Nmailles; imaille++){
228 _idm[0] = _nVar*imaille;
229 for(int k=1; k<_nVar; k++)
230 _idm[k] = _idm[k-1] + 1;
231 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
232 primToConsJacobianMatrix(_Vi);
233 for(int k=0; k<_nVar*_nVar; k++)
234 _primToConsJacoMat[k]*=1/_dt;
235 MatSetValuesBlocked(_A, 1, &imaille, 1, &imaille, _primToConsJacoMat, ADD_VALUES);
237 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
239 #if PETSC_VERSION_GREATER_3_5
240 KSPSetOperators(_ksp, _A, _A);
242 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
247 cout << "Matrice du système linéaire" << endl;
248 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
250 cout << "Second membre du système linéaire" << endl;
251 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
256 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
259 KSPSolve(_ksp, _b, _newtonVariation);
263 computeScaling(_maxvp);
265 VecAssemblyBegin(_vecScaling);
266 VecAssemblyBegin(_invVecScaling);
267 for(int imaille = 0; imaille<_Nmailles; imaille++)
270 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
271 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
273 VecAssemblyEnd(_vecScaling);
274 VecAssemblyEnd(_invVecScaling);
278 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
279 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
281 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
282 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
285 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
288 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
289 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
292 VecCopy(_b,_bScaling);
293 VecPointwiseMult(_b,_vecScaling,_bScaling);
296 cout << "Produit du second membre par le preconditionneur bloc diagonal a gauche" << endl;
297 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
301 KSPSolve(_ksp,_b, _bScaling);
302 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
306 cout << "solution du systeme lineaire local:" << endl;
307 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
313 void SinglePhase::convectionState( const long &i, const long &j, const bool &IsBord){
314 //First conservative state then further down we will compute interface (Roe) state and then compute primitive state
316 for(int k=1; k<_nVar; k++)
317 _idm[k] = _idm[k-1] + 1;
318 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
321 for(int k=1; k<_nVar; k++)
322 _idm[k] = _idm[k-1] + 1;
324 VecGetValues(_Uext, _nVar, _idm, _Uj);
326 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
327 if(_verbose && _nbTimeStep%_freqSave ==0)
329 cout<<"Convection Left state cell " << i<< ": "<<endl;
330 for(int k =0; k<_nVar; k++)
332 cout<<"Convection Right state cell " << j<< ": "<<endl;
333 for(int k =0; k<_nVar; k++)
336 if(_Ui[0]<0||_Uj[0]<0)
338 cout<<"!!!!!!!!!!!!!!!!!!!!!!!!densite negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
339 *_runLogFile<<"!!!!!!!!!!!!!!!!!!!!!!!!densite negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
340 _runLogFile->close();
341 throw CdmathException("densite negative, arret de calcul");
344 //Computation of Roe state
345 PetscScalar ri, rj, xi, xj, pi, pj;
347 ri = sqrt(_Ui[0]);//racine carre de phi_i rho_i
349 _Uroe[0] = ri*rj; //moyenne geometrique des densites
350 if(_verbose && _nbTimeStep%_freqSave ==0)
351 cout << "Densité moyenne Roe gauche " << i << ": " << ri*ri << ", droite " << j << ": " << rj*rj << "->" << _Uroe[0] << endl;
352 for(int k=0;k<_Ndim;k++){
355 _Uroe[1+k] = (xi/ri + xj/rj)/(ri + rj);
356 //"moyenne" des vitesses
357 if(_verbose && _nbTimeStep%_freqSave ==0)
358 cout << "Vitesse de Roe composante "<< k<<" gauche " << i << ": " << xi/(ri*ri) << ", droite " << j << ": " << xj/(rj*rj) << "->" << _Uroe[k+1] << endl;
360 // H = (rho E + p)/rho
361 xi = _Ui[_nVar-1];//phi rho E
363 Ii = i*_nVar; // correct Kieu
364 VecGetValues(_primitiveVars, 1, &Ii, &pi);// _primitiveVars pour p
368 for(int k=1;k<=_Ndim;k++)
369 q_2 += _Uj[k]*_Uj[k];
370 q_2 /= _Uj[0]; //phi rho u²
371 pj = _fluides[0]->getPressure((_Uj[(_Ndim+2)-1]-q_2/2)/_porosityj,_Uj[0]/_porosityj);
375 Ii = j*_nVar; // correct Kieu
376 VecGetValues(_primitiveVars, 1, &Ii, &pj);
378 xi = (xi + pi)/(ri*ri);
379 xj = (xj + pj)/(rj*rj);
380 _Uroe[_nVar-1] = (ri*xi + rj*xj)/(ri + rj);
381 //on se donne l enthalpie ici
382 if(_verbose && _nbTimeStep%_freqSave ==0)
383 cout << "Enthalpie totale de Roe H gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[_nVar-1] << endl;
385 if(_verbose && _nbTimeStep%_freqSave ==0)
387 cout<<"Convection interfacial state"<<endl;
388 for(int k=0;k<_nVar;k++)
389 cout<< _Uroe[k]<<" , "<<endl;
392 //Extraction of primitive states
393 _idm[0] = _nVar*i; // Kieu
394 for(int k=1; k<_nVar; k++)
395 _idm[k] = _idm[k-1] + 1;
397 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
398 if (_verbose && _nbTimeStep%_freqSave ==0)
400 cout << "Convection state: variables primitives maille " << i<<endl;
401 for(int q=0; q<_nVar; q++)
402 cout << _Vi[q] << endl;
407 for(int k=0; k<_nVar; k++)
408 _idn[k] = _nVar*j + k;
410 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
414 for(int k=0; k<_nVar; k++)
417 VecGetValues(_Uextdiff, _nVar, _idn, _phi);
418 consToPrim(_phi,_Vj,1);
421 if (_verbose && _nbTimeStep%_freqSave ==0)
423 cout << "Convection state: variables primitives maille " <<j <<endl;
424 for(int q=0; q<_nVar; q++)
425 cout << _Vj[q] << endl;
430 void SinglePhase::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
431 //sortie: matrices et etat de diffusion (rho, q, T)
433 for(int k=1; k<_nVar; k++)
434 _idm[k] = _idm[k-1] + 1;
436 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
437 for(int k=0; k<_nVar; k++)
441 VecGetValues(_Uextdiff, _nVar, _idn, _Uj);
443 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
445 if(_verbose && _nbTimeStep%_freqSave ==0)
447 cout << "SinglePhase::diffusionStateAndMatrices cellule gauche" << i << endl;
449 for(int q=0; q<_nVar; q++)
450 cout << _Ui[q] << "\t";
452 cout << "SinglePhase::diffusionStateAndMatrices cellule droite" << j << endl;
454 for(int q=0; q<_nVar; q++)
455 cout << _Uj[q] << "\t";
459 for(int k=0; k<_nVar; k++)
460 _Udiff[k] = (_Ui[k]/_porosityi+_Uj[k]/_porosityj)/2;
462 if(_verbose && _nbTimeStep%_freqSave ==0)
464 cout << "SinglePhase::diffusionStateAndMatrices conservative diffusion state" << endl;
466 for(int q=0; q<_nVar; q++)
467 cout << _Udiff[q] << "\t";
469 cout << "porosite gauche= "<<_porosityi<< ", porosite droite= "<<_porosityj<<endl;
471 consToPrim(_Udiff,_phi,1);
472 _Udiff[_nVar-1]=_phi[_nVar-1];
474 if(_timeScheme==Implicit)
477 for (int i = 0; i<_Ndim;i++)
478 q_2+=_Udiff[i+1]*_Udiff[i+1];
479 double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
480 double lambda = _fluides[0]->getConductivity(_Udiff[_nVar-1]);
481 double Cv= _fluides[0]->constante("Cv");
482 for(int i=0; i<_nVar*_nVar;i++)
484 for(int i=1;i<(_nVar-1);i++)
486 _Diffusion[i*_nVar] = mu*_Udiff[i]/(_Udiff[0]*_Udiff[0]);
487 _Diffusion[i*_nVar+i] = -mu/_Udiff[0];
489 int i = (_nVar-1)*_nVar;
490 _Diffusion[i]=lambda*(_Udiff[_nVar-1]/_Udiff[0]-q_2/(2*Cv*_Udiff[0]*_Udiff[0]*_Udiff[0]));
491 for(int k=1;k<(_nVar-1);k++)
493 _Diffusion[i+k]= lambda*_Udiff[k]/(_Udiff[0]*_Udiff[0]*Cv);
495 _Diffusion[_nVar*_nVar-1]=-lambda/(_Udiff[0]*Cv);
498 void SinglePhase::setBoundaryState(string nameOfGroup, const int &j,double *normale){
500 for(int k=1; k<_nVar; k++)
501 _idm[k] = _idm[k-1] + 1;
503 double porosityj=_porosityField(j);
505 if(_verbose && _nbTimeStep%_freqSave ==0)
507 cout << "setBoundaryState for group "<< nameOfGroup << ", inner cell j= "<<j<< " face unit normal vector "<<endl;
508 for(int k=0; k<_Ndim; k++){
509 cout<<normale[k]<<", ";
514 if (_limitField[nameOfGroup].bcType==Wall){
515 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne conservatif
516 double q_n=0;//q_n=quantité de mouvement normale à la face frontière;
517 for(int k=0; k<_Ndim; k++)
518 q_n+=_externalStates[(k+1)]*normale[k];
520 //Pour la convection, inversion du sens de la vitesse normale
521 for(int k=0; k<_Ndim; k++)
522 _externalStates[(k+1)]-= 2*q_n*normale[k];
525 for(int k=1; k<_nVar; k++)
526 _idm[k] = _idm[k-1] + 1;
528 VecAssemblyBegin(_Uext);
529 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
530 VecAssemblyEnd(_Uext);
532 //Pour la diffusion, paroi à vitesse et temperature imposees
534 for(int k=1; k<_nVar; k++)
535 _idm[k] = _idm[k-1] + 1;
536 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);//L'état fantome contient à présent les variables primitives internes
537 double pression=_externalStates[0];
538 double T=_limitField[nameOfGroup].T;
539 double rho=_fluides[0]->getDensity(pression,T);
541 _externalStates[0]=porosityj*rho;
542 _externalStates[1]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
544 v2 +=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
547 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
548 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
551 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
552 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
555 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
557 for(int k=1; k<_nVar; k++)
558 _idm[k] = _idm[k-1] + 1;
559 VecAssemblyBegin(_Uextdiff);
560 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
561 VecAssemblyEnd(_Uextdiff);
563 else if (_limitField[nameOfGroup].bcType==Neumann){
564 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On prend l'état fantôme égal à l'état interne (conditions limites de Neumann)
567 for(int k=1; k<_nVar; k++)
568 _idm[k] = _idm[k-1] + 1;
570 VecAssemblyBegin(_Uext);
571 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
572 VecAssemblyEnd(_Uext);
574 VecAssemblyBegin(_Uextdiff);
575 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
576 VecAssemblyEnd(_Uextdiff);
578 else if (_limitField[nameOfGroup].bcType==Inlet){
579 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne (conditions limites de Neumann)
580 double q_int_n=0;//q_int_n=composante normale de la quantité de mouvement à la face frontière;
581 for(int k=0; k<_Ndim; k++)
582 q_int_n+=_externalStates[(k+1)]*normale[k];//On calcule la vitesse normale sortante
584 double q_ext_n=_limitField[nameOfGroup].v_x[0]*normale[0];
587 q_ext_n+=_limitField[nameOfGroup].v_y[0]*normale[1];
589 q_ext_n+=_limitField[nameOfGroup].v_z[0]*normale[2];
592 if(q_int_n+q_ext_n<=0){//Interfacial velocity goes inward
593 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);//On met à jour l'état fantome avec les variables primitives internes
594 double pression=_externalStates[0];
595 double T=_limitField[nameOfGroup].T;
596 double rho=_fluides[0]->getDensity(pression,T);
598 _externalStates[0]=porosityj*rho;//Composante fantome de masse
599 _externalStates[1]=_externalStates[0]*(_limitField[nameOfGroup].v_x[0]);//Composante fantome de qdm x
601 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
604 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
605 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];//Composante fantome de qdm y
608 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];//Composante fantome de qdm z
609 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
612 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);//Composante fantome de l'nrj
614 else if(_nbTimeStep%_freqSave ==0)
615 cout<< "Warning : fluid possibly going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
618 for(int k=1; k<_nVar; k++)
619 _idm[k] = _idm[k-1] + 1;
620 VecAssemblyBegin(_Uext);
621 VecAssemblyBegin(_Uextdiff);
622 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
623 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
624 VecAssemblyEnd(_Uext);
625 VecAssemblyEnd(_Uextdiff);
627 else if (_limitField[nameOfGroup].bcType==InletRotationVelocity){
628 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
629 double u_int_n=0;//u_int_n=composante normale de la vitesse intérieure à la face frontière;
630 for(int k=0; k<_Ndim; k++)
631 u_int_n+=_externalStates[(k+1)]*normale[k];
638 omega[0]=_limitField[nameOfGroup].v_x[0];
639 omega[1]=_limitField[nameOfGroup].v_y[0];
642 Normale[0]=normale[0];
643 Normale[1]=normale[1];
647 omega[2]=_limitField[nameOfGroup].v_z[0];
648 Normale[2]=normale[2];
651 Vector tangent_vel=omega%Normale;
652 u_ext_n=-0.01*tangent_vel.norm();
653 //Changing external state velocity
654 for(int k=0; k<_Ndim; k++)
655 _externalStates[(k+1)]=u_ext_n*normale[k] + tangent_vel[k];
658 if(u_ext_n + u_int_n <= 0){
659 double pression=_externalStates[0];
660 double T=_limitField[nameOfGroup].T;
661 double rho=_fluides[0]->getDensity(pression,T);
664 v2 +=_externalStates[1]*_externalStates[1];//v_x*v_x
665 _externalStates[0]=porosityj*rho;
666 _externalStates[1]*=_externalStates[0];
669 v2 +=_externalStates[2]*_externalStates[2];//+v_y*v_y
670 _externalStates[2]*=_externalStates[0];
673 v2 +=_externalStates[3]*_externalStates[3];//+v_z*v_z
674 _externalStates[3]*=_externalStates[0];
677 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
679 else if(_nbTimeStep%_freqSave ==0)
682 * cout<< "Warning : fluid going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
684 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On définit l'état fantôme avec l'état interne
685 if(_nbTimeStep%_freqSave ==0)
686 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Wall boundary condition."<<endl;
688 //Changing external state momentum
689 for(int k=0; k<_Ndim; k++)
690 _externalStates[(k+1)]-=2*_externalStates[0]*u_int_n*normale[k];
694 for(int k=1; k<_nVar; k++)
695 _idm[k] = _idm[k-1] + 1;
696 VecAssemblyBegin(_Uext);
697 VecAssemblyBegin(_Uextdiff);
698 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
699 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
700 VecAssemblyEnd(_Uext);
701 VecAssemblyEnd(_Uextdiff);
703 else if (_limitField[nameOfGroup].bcType==InletPressure){
704 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
706 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
707 Cell Cj=_mesh.getCell(j);
708 double hydroPress=Cj.x()*_GravityField3d[0];
710 hydroPress+=Cj.y()*_GravityField3d[1];
712 hydroPress+=Cj.z()*_GravityField3d[2];
714 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
716 //Building the primitive external state
717 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
718 double u_n=0;//u_n=vitesse normale à la face frontière;
719 for(int k=0; k<_Ndim; k++)
720 u_n+=_externalStates[(k+1)]*normale[k];
724 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress,_limitField[nameOfGroup].T);
726 //Contribution from the tangential velocity
730 omega[0]=_limitField[nameOfGroup].v_x[0];
731 omega[1]=_limitField[nameOfGroup].v_y[0];
734 Normale[0]=normale[0];
735 Normale[1]=normale[1];
739 omega[2]=_limitField[nameOfGroup].v_z[0];
740 Normale[2]=normale[2];
743 Vector tangent_vel=omega%Normale;
745 //Changing external state velocity
746 for(int k=0; k<_Ndim; k++)
747 _externalStates[(k+1)]=u_n*normale[k] + abs(u_n)*tangent_vel[k];
752 if(_nbTimeStep%_freqSave ==0)
753 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Neumann boundary condition for velocity and temperature (only pressure value is imposed as in outlet BC)."<<endl;
754 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
756 if(_nbTimeStep%_freqSave ==0)
757 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Wall boundary condition."<<endl;
758 _externalStates[0]=porosityj*_fluides[0]->getDensity(_externalStates[0]+hydroPress, _externalStates[_nVar-1]);
759 //Changing external state velocity
760 for(int k=0; k<_Ndim; k++)
761 _externalStates[(k+1)]-=2*u_n*normale[k];
765 for(int k=0; k<_Ndim; k++)
767 v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
768 _externalStates[(k+1)]*=_externalStates[0] ;//qdm component
770 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);//nrj component
774 for(int k=1; k<_nVar; k++)
775 _idm[k] = _idm[k-1] + 1;
776 VecAssemblyBegin(_Uext);
777 VecAssemblyBegin(_Uextdiff);
778 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
779 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
780 VecAssemblyEnd(_Uext);
781 VecAssemblyEnd(_Uextdiff);
783 else if (_limitField[nameOfGroup].bcType==Outlet){
784 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne conservatif
785 double q_n=0;//q_n=quantité de mouvement normale à la face frontière;
786 for(int k=0; k<_Ndim; k++)
787 q_n+=_externalStates[(k+1)]*normale[k];
789 if(q_n < -_precision && _nbTimeStep%_freqSave ==0)
791 cout<< "Warning : fluid going in through outlet boundary "<<nameOfGroup<<" with flow rate "<< q_n<<endl;
792 cout<< "Applying Neumann boundary condition for velocity and temperature"<<endl;
794 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
795 Cell Cj=_mesh.getCell(j);
796 double hydroPress=Cj.x()*_GravityField3d[0];
798 hydroPress+=Cj.y()*_GravityField3d[1];
800 hydroPress+=Cj.z()*_GravityField3d[2];
802 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
804 if(_verbose && _nbTimeStep%_freqSave ==0)
806 cout<<"Cond lim outlet densite= "<<_externalStates[0]<<" gravite= "<<_GravityField3d[0]<<" Cj.x()= "<<Cj.x()<<endl;
807 cout<<"Cond lim outlet pression ref= "<<_limitField[nameOfGroup].p<<" pression hydro= "<<hydroPress<<" total= "<<_limitField[nameOfGroup].p+hydroPress<<endl;
809 //Building the external state
810 _idm[0] = _nVar*j;// Kieu
811 for(int k=1; k<_nVar; k++)
812 _idm[k] = _idm[k-1] + 1;
813 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
815 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
817 for(int k=0; k<_Ndim; k++)
819 v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
820 _externalStates[(k+1)]*=_externalStates[0] ;
822 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);
824 for(int k=1; k<_nVar; k++)
825 _idm[k] = _idm[k-1] + 1;
826 VecAssemblyBegin(_Uext);
827 VecAssemblyBegin(_Uextdiff);
828 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
829 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
830 VecAssemblyEnd(_Uext);
831 VecAssemblyEnd(_Uextdiff);
833 cout<<"Boundary condition not set for boundary named "<<nameOfGroup<< " _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
834 cout<<"Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
835 *_runLogFile<<"Boundary condition not set for boundary named. Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
836 _runLogFile->close();
837 throw CdmathException("Unknown boundary condition");
841 void SinglePhase::convectionMatrices()
843 //entree: URoe = rho, u, H
844 //sortie: matrices Roe+ et Roe-
846 if(_verbose && _nbTimeStep%_freqSave ==0)
847 cout<<"SinglePhase::convectionMatrices()"<<endl;
849 double u_n=0, u_2=0;//vitesse normale et carré du module
851 for(int i=0;i<_Ndim;i++)
853 u_2 += _Uroe[1+i]*_Uroe[1+i];
854 u_n += _Uroe[1+i]*_vec_normal[i];
857 vector<complex<double> > vp_dist(3,0);
859 if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
861 staggeredVFFCMatricesConservativeVariables(u_n);//Computation of classical upwinding matrices
862 if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//For use in implicit matrix
863 staggeredVFFCMatricesPrimitiveVariables(u_n);
867 Vector vitesse(_Ndim);
868 for(int idim=0;idim<_Ndim;idim++)
869 vitesse[idim]=_Uroe[1+idim];
872 /***********Calcul des valeurs propres ********/
874 c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
875 k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
876 K = u_2*k/2; //g-1/2 *|u|²
878 vp_dist[0]=u_n-c;vp_dist[1]=u_n;vp_dist[2]=u_n+c;
880 _maxvploc=fabs(u_n)+c;
884 if(_verbose && _nbTimeStep%_freqSave ==0)
885 cout<<"SinglePhase::convectionMatrices Eigenvalues "<<u_n-c<<" , "<<u_n<<" , "<<u_n+c<<endl;
887 RoeMatrixConservativeVariables( u_n, H,vitesse,k,K);
889 /******** Construction des matrices de decentrement ********/
890 if( _spaceScheme ==centered){
891 if(_entropicCorrection)
893 *_runLogFile<<"SinglePhase::convectionMatrices: entropy scheme not available for centered scheme"<<endl;
894 _runLogFile->close();
895 throw CdmathException("SinglePhase::convectionMatrices: entropy scheme not available for centered scheme");
898 for(int i=0; i<_nVar*_nVar;i++)
901 else if(_spaceScheme == upwind || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach){
902 if(_entropicCorrection)
903 entropicShift(_vec_normal);
905 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
907 vector< complex< double > > y (3,0);
908 for( int i=0 ; i<3 ; i++)
909 y[i] = Polynoms::abs_generalise(vp_dist[i])+_entropicShift[i];
910 Polynoms::abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
912 if( _spaceScheme ==pressureCorrection)
913 for( int i=0 ; i<_Ndim ; i++)
914 for( int j=0 ; j<_Ndim ; j++)
915 _absAroe[(1+i)*_nVar+1+j]-=(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
916 else if( _spaceScheme ==lowMach){
917 double M=sqrt(u_2)/c;
918 for( int i=0 ; i<_Ndim ; i++)
919 for( int j=0 ; j<_Ndim ; j++)
920 _absAroe[(1+i)*_nVar+1+j]-=(1-M)*(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
923 else if( _spaceScheme ==staggered ){
924 if(_entropicCorrection)//To do: study entropic correction for staggered
926 *_runLogFile<<"SinglePhase::convectionMatrices: entropy scheme not available for staggered scheme"<<endl;
927 _runLogFile->close();
928 throw CdmathException("SinglePhase::convectionMatrices: entropy scheme not available for staggered scheme");
931 staggeredRoeUpwindingMatrixConservativeVariables( u_n, H, vitesse, k, K);
935 *_runLogFile<<"SinglePhase::convectionMatrices: scheme not treated"<<endl;
936 _runLogFile->close();
937 throw CdmathException("SinglePhase::convectionMatrices: scheme not treated");
940 for(int i=0; i<_nVar*_nVar;i++)
942 _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
943 _AroePlus[i] = (_Aroe[i]+_absAroe[i])/2;
945 if(_timeScheme==Implicit)
947 if(_usePrimitiveVarsInNewton)//Implicitation using primitive variables
949 _Vij[0]=_fluides[0]->getPressureFromEnthalpy(_Uroe[_nVar-1]-u_2/2, _Uroe[0]);//pressure
950 _Vij[_nVar-1]=_fluides[0]->getTemperatureFromPressure( _Vij[0], _Uroe[0]);//Temperature
951 for(int idim=0;idim<_Ndim; idim++)
952 _Vij[1+idim]=_Uroe[1+idim];
953 primToConsJacobianMatrix(_Vij);
954 Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
955 Polynoms::matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
958 for(int i=0; i<_nVar*_nVar;i++)
960 _AroeMinusImplicit[i] = _AroeMinus[i];
961 _AroePlusImplicit[i] = _AroePlus[i];
964 if(_verbose && _nbTimeStep%_freqSave ==0)
966 displayMatrix(_Aroe, _nVar,"Matrice de Roe");
967 displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
968 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
969 displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
973 if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
975 displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
976 displayMatrix(_AroePlusImplicit, _nVar,"Matrice _AroePlusImplicit");
979 /*********Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source*****/
980 if(_entropicCorrection)
982 InvMatriceRoe( vp_dist);
983 Polynoms::matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
985 else if (_spaceScheme==upwind || (_spaceScheme ==pressureCorrection ) || (_spaceScheme ==lowMach ))//upwind sans entropic
986 SigneMatriceRoe( vp_dist);
987 else if (_spaceScheme==centered)//centre sans entropic
988 for(int i=0; i<_nVar*_nVar;i++)
990 else if( _spaceScheme ==staggered )//à tester
999 for(int i=0; i<_nVar*_nVar;i++)
1001 _signAroe[0] = signu;
1002 for(int i=1; i<_nVar-1;i++)
1003 _signAroe[i*_nVar+i] = -signu;
1004 _signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
1008 *_runLogFile<<"SinglePhase::convectionMatrices: well balanced option not treated"<<endl;
1009 _runLogFile->close();
1010 throw CdmathException("SinglePhase::convectionMatrices: well balanced option not treated");
1013 void SinglePhase::computeScaling(double maxvp)
1017 for(int q=1; q<_nVar-1; q++)
1019 _blockDiag[q]=1./maxvp;//
1020 _invBlockDiag[q]= maxvp;//1.;//
1022 _blockDiag[_nVar - 1]=(_fluides[0]->constante("gamma")-1)/(maxvp*maxvp);//1
1023 _invBlockDiag[_nVar - 1]= 1./_blockDiag[_nVar - 1] ;// 1.;//
1026 void SinglePhase::addDiffusionToSecondMember
1032 double lambda = _fluides[0]->getConductivity(_Udiff[_nVar-1]);
1033 double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
1036 lambda=max(lambda,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
1038 if(lambda==0 && mu ==0 && _heatTransfertCoeff==0)
1042 for(int k=1; k<_nVar; k++)
1043 _idm[k] = _idm[k-1] + 1;
1046 VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
1050 //on n'a pas de contribution sur la masse
1052 //contribution visqueuse sur la quantite de mouvement
1053 for(int k=1; k<_nVar-1; k++)
1054 _phi[k] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu*(_porosityj*_Vj[k] - _porosityi*_Vi[k]);
1055 //contribution visqueuse sur l'energie
1056 _phi[_nVar-1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*lambda*(_porosityj*_Vj[_nVar-1] - _porosityi*_Vi[_nVar-1]);
1059 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1061 if(_verbose && _nbTimeStep%_freqSave ==0)
1063 cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
1064 for(int q=0; q<_nVar; q++)
1065 cout << _phi[q] << endl;
1071 //On change de signe pour l'autre contribution
1072 for(int k=0; k<_nVar; k++)
1073 _phi[k] *= -_inv_dxj/_inv_dxi;
1076 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1077 if(_verbose && _nbTimeStep%_freqSave ==0)
1079 cout << "Contribution diffusion au 2nd membre pour la maille " << j << ": "<<endl;
1080 for(int q=0; q<_nVar; q++)
1081 cout << _phi[q] << endl;
1086 if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
1088 cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
1089 for(int i=0; i<_nVar; i++)
1091 for(int j=0; j<_nVar; j++)
1092 cout << _Diffusion[i*_nVar+j]<<", ";
1099 void SinglePhase::sourceVector(PetscScalar * Si, PetscScalar * Ui, PetscScalar * Vi, int i)
1101 double phirho=Ui[0], T=Vi[_nVar-1];
1103 for(int k=0; k<_Ndim; k++)
1104 norm_u+=Vi[1+k]*Vi[1+k];
1105 norm_u=sqrt(norm_u);
1107 Si[0]=_heatPowerField(i)/_latentHeat;
1110 for(int k=1; k<_nVar-1; k++)
1111 Si[k] =(_gravite[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*phirho;
1113 Si[_nVar-1]=_heatPowerField(i);
1115 for(int k=0; k<_Ndim; k++)
1116 Si[_nVar-1] +=(_GravityField3d[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*Vi[1+k]*phirho;
1118 if(_timeScheme==Implicit)
1120 for(int k=0; k<_nVar*_nVar;k++)
1121 _GravityImplicitationMatrix[k] = 0;
1122 if(!_usePrimitiveVarsInNewton)
1123 for(int k=0; k<_nVar;k++)
1124 _GravityImplicitationMatrix[k*_nVar]=-_gravite[k];
1127 double pression=Vi[0];
1128 getDensityDerivatives( pression, T, norm_u*norm_u);
1129 for(int k=0; k<_nVar;k++)
1131 _GravityImplicitationMatrix[k*_nVar+0] =-_gravite[k]*_drho_sur_dp;
1132 _GravityImplicitationMatrix[k*_nVar+_nVar-1]=-_gravite[k]*_drho_sur_dT;
1136 if(_verbose && _nbTimeStep%_freqSave ==0)
1138 cout<<"SinglePhase::sourceVector"<<endl;
1140 for(int k=0;k<_nVar;k++)
1144 for(int k=0;k<_nVar;k++)
1148 for(int k=0;k<_nVar;k++)
1151 if(_timeScheme==Implicit)
1152 displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
1156 void SinglePhase::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
1158 double norm_u=0, u_n=0, rho;
1159 for(int i=0;i<_Ndim;i++)
1160 u_n += _Uroe[1+i]*_vec_normal[i];
1164 for(int i=0;i<_Ndim;i++)
1165 norm_u += Vi[1+i]*Vi[1+i];
1166 norm_u=sqrt(norm_u);
1168 for(int i=0;i<_Ndim;i++)
1169 pressureLoss[1+i]=-1/2*K*rho*norm_u*Vi[1+i];
1172 for(int i=0;i<_Ndim;i++)
1173 norm_u += Vj[1+i]*Vj[1+i];
1174 norm_u=sqrt(norm_u);
1176 for(int i=0;i<_Ndim;i++)
1177 pressureLoss[1+i]=-1/2*K*rho*norm_u*Vj[1+i];
1179 pressureLoss[_nVar-1]=-1/2*K*rho*norm_u*norm_u*norm_u;
1181 if(_verbose && _nbTimeStep%_freqSave ==0)
1183 cout<<"SinglePhase::pressureLossVector K= "<<K<<endl;
1185 for(int k=0;k<_nVar;k++)
1189 for(int k=0;k<_nVar;k++)
1193 for(int k=0;k<_nVar;k++)
1197 for(int k=0;k<_nVar;k++)
1200 cout<<"pressureLoss="<<endl;
1201 for(int k=0;k<_nVar;k++)
1202 cout<<pressureLoss[k]<<", ";
1207 void SinglePhase::porosityGradientSourceVector()
1209 double u_ni=0, u_nj=0, rhoi,rhoj, pi=_Vi[0], pj=_Vj[0],pij;
1210 for(int i=0;i<_Ndim;i++) {
1211 u_ni += _Vi[1+i]*_vec_normal[i];
1212 u_nj += _Vj[1+i]*_vec_normal[i];
1214 _porosityGradientSourceVector[0]=0;
1215 rhoj=_Uj[0]/_porosityj;
1216 rhoi=_Ui[0]/_porosityi;
1217 pij=(pi+pj)/2+rhoi*rhoj/2/(rhoi+rhoj)*(u_ni-u_nj)*(u_ni-u_nj);
1218 for(int i=0;i<_Ndim;i++)
1219 _porosityGradientSourceVector[1+i]=pij*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
1220 _porosityGradientSourceVector[_nVar-1]=0;
1223 void SinglePhase::jacobian(const int &j, string nameOfGroup,double * normale)
1225 if(_verbose && _nbTimeStep%_freqSave ==0)
1226 cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
1229 for(k=0; k<_nVar*_nVar;k++)
1230 _Jcb[k] = 0;//No implicitation at this stage
1233 for(k=1; k<_nVar; k++)
1234 _idm[k] = _idm[k-1] + 1;
1235 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);
1236 double q_n=0;//quantité de mouvement normale à la paroi
1237 for(k=0; k<_Ndim; k++)
1238 q_n+=_externalStates[(k+1)]*normale[k];
1240 // loop of boundary types
1241 if (_limitField[nameOfGroup].bcType==Wall)
1243 for(k=0; k<_nVar;k++)
1244 _Jcb[k*_nVar + k] = 1;
1245 for(k=1; k<_nVar-1;k++)
1246 for(int l=1; l<_nVar-1;l++)
1247 _Jcb[k*_nVar + l] -= 2*normale[k-1]*normale[l-1];
1249 else if (_limitField[nameOfGroup].bcType==Inlet)
1253 double v[_Ndim], ve[_Ndim], v2, ve2;
1256 for(k=1; k<_nVar; k++)
1257 _idm[k] = _idm[k-1] + 1;
1258 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1259 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1261 ve[0] = _limitField[nameOfGroup].v_x[0];
1266 ve[1] = _limitField[nameOfGroup].v_y[0];
1272 ve[2] = _limitField[nameOfGroup].v_z[0];
1277 double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,_Uj[0]);
1278 double total_energy=internal_energy+ve2/2;
1281 _Jcb[0]=v2/(2*internal_energy);
1282 for(k=0; k<_Ndim;k++)
1283 _Jcb[1+k]=-v[k]/internal_energy;
1284 _Jcb[_nVar-1]=1/internal_energy;
1286 for(int l =1;l<1+_Ndim;l++){
1287 _Jcb[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
1288 for(k=0; k<_Ndim;k++)
1289 _Jcb[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
1290 _Jcb[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
1293 _Jcb[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
1294 for(k=0; k<_Ndim;k++)
1295 _Jcb[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
1296 _Jcb[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
1299 for(k=0;k<_nVar;k++)
1304 _Jcb[_nVar]=_limitField[nameOfGroup].v_x[0];//Kieu
1305 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
1308 _Jcb[2*_nVar]= _limitField[nameOfGroup].v_y[0];
1309 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
1311 _Jcb[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1312 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
1315 _Jcb[(_nVar-1)*_nVar]=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2;
1318 else if (_limitField[nameOfGroup].bcType==InletPressure && q_n<0){
1320 double v[_Ndim], v2=0;
1322 for(k=1; k<_nVar; k++)
1323 _idm[k] = _idm[k-1] + 1;
1324 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1326 for(k=0; k<_Ndim;k++){
1331 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _limitField[nameOfGroup].T);
1332 double rho_int = _externalStates[0];
1333 double density_ratio=rho_ext/rho_int;
1335 for(int l =1;l<1+_Ndim;l++){
1336 _Jcb[l*_nVar]=-density_ratio*v[l-1];
1337 _Jcb[l*_nVar+l]=density_ratio;
1340 _Jcb[(_nVar-1)*_nVar]=-v2*density_ratio;
1341 for(k=0; k<_Ndim;k++)
1342 _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k];
1344 // not wall, not inlet, not inletPressure
1345 else if(_limitField[nameOfGroup].bcType==Outlet || (_limitField[nameOfGroup].bcType==InletPressure && q_n>=0))
1348 double v[_Ndim], v2=0;
1350 for(k=1; k<_nVar; k++)
1351 _idm[k] = _idm[k-1] + 1;
1352 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1354 for(k=0; k<_Ndim;k++){
1359 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _externalStates[_nVar-1]);
1360 double rho_int = _externalStates[0];
1361 double density_ratio=rho_ext/rho_int;
1362 double internal_energy=_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho_int);
1363 double total_energy=internal_energy+v2/2;
1366 _Jcb[0]=density_ratio*(1-v2/(2*internal_energy));
1367 for(k=0; k<_Ndim;k++)
1368 _Jcb[1+k]=density_ratio*v[k]/internal_energy;
1369 _Jcb[_nVar-1]=-density_ratio/internal_energy;
1371 for(int l =1;l<1+_Ndim;l++){
1372 _Jcb[l*_nVar]=density_ratio*v2*v[l-1]/(2*internal_energy);
1373 for(k=0; k<_Ndim;k++)
1374 _Jcb[l*_nVar+1+k]=density_ratio*v[k]*v[l-1]/internal_energy;
1375 _Jcb[l*_nVar+1+k]-=density_ratio;
1376 _Jcb[l*_nVar+_nVar-1]=-density_ratio*v[l-1]/internal_energy;
1379 _Jcb[(_nVar-1)*_nVar]=density_ratio*v2*total_energy/(2*internal_energy);
1380 for(k=0; k<_Ndim;k++)
1381 _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k]*total_energy/internal_energy;
1382 _Jcb[(_nVar-1)*_nVar+_nVar-1]=density_ratio*(1-total_energy/internal_energy);
1386 double cd = 1,cn=0,p0, gamma;
1387 _idm[0] = j*_nVar;// Kieu
1388 for(k=1; k<_nVar;k++)
1389 _idm[k] = _idm[k-1] + 1;
1390 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1391 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1393 // compute the common numerator and common denominator
1394 p0=_fluides[0]->constante("p0");
1395 gamma =_fluides[0]->constante("gamma");
1396 cn =_limitField[nameOfGroup].p +p0;
1397 cd = _phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0;
1401 for(k=1; k<_nVar-1;k++)
1402 v2+=_externalStates[k]*_externalStates[k];
1404 _JcbDiff[0] = cn*(_phi[_nVar-1] -v2 -p0)/cd;
1405 for(k=1; k<_nVar-1;k++)
1406 _JcbDiff[k]=cn*_phi[k]/cd;
1407 _JcbDiff[_nVar-1]= -cn*_phi[0]/cd;
1409 for(idim=0; idim<_Ndim;idim++)
1412 _JcbDiff[(1+idim)*_nVar]=-(v2*cn*_phi[idim+1])/(2*cd);
1413 //colonnes intermediaire
1414 for(jdim=0; jdim<_Ndim;jdim++)
1416 _JcbDiff[(1+idim)*_nVar + jdim + 1] =_externalStates[idim+1]*_phi[jdim+1];
1417 _JcbDiff[(1+idim)*_nVar + jdim + 1]*=cn/cd;
1419 //matrice identite*cn*(rhoe- p0)
1420 _JcbDiff[(1+idim)*_nVar + idim + 1] +=( cn*(_phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0))/cd;
1423 _JcbDiff[(1+idim)*_nVar + _nVar-1]=-_phi[idim+1]*cn/cd;
1426 _JcbDiff[_nVar*(_nVar-1)] = -(v2*_phi[_nVar -1]*cn)/(2*cd);
1427 for(int idim=0; idim<_Ndim;idim++)
1428 _JcbDiff[_nVar*(_nVar-1)+idim+1]=_externalStates[idim+1]*_phi[_nVar -1]*cn/cd;
1429 _JcbDiff[_nVar*_nVar -1] = -(v2/2+p0)*cn/cd;
1432 else if (_limitField[nameOfGroup].bcType!=Neumann)// not wall, not inlet, not outlet
1434 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1435 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1436 _runLogFile->close();
1437 throw CdmathException("SinglePhase::jacobian: This boundary condition is not treated");
1441 //calcule la jacobienne pour une CL de diffusion
1442 void SinglePhase::jacobianDiff(const int &j, string nameOfGroup)
1444 if(_verbose && _nbTimeStep%_freqSave ==0)
1445 cout<<"Jacobienne condition limite diffusion bord "<< nameOfGroup<<endl;
1448 for(k=0; k<_nVar*_nVar;k++)
1451 if (_limitField[nameOfGroup].bcType==Wall){
1452 double v[_Ndim], ve[_Ndim], v2, ve2;
1455 for(k=1; k<_nVar; k++)
1456 _idm[k] = _idm[k-1] + 1;
1457 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1458 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1460 ve[0] = _limitField[nameOfGroup].v_x[0];
1465 ve[1] = _limitField[nameOfGroup].v_y[0];
1471 ve[2] = _limitField[nameOfGroup].v_z[0];
1477 double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho);
1478 double total_energy=internal_energy+ve2/2;
1481 _JcbDiff[0]=v2/(2*internal_energy);
1482 for(k=0; k<_Ndim;k++)
1483 _JcbDiff[1+k]=-v[k]/internal_energy;
1484 _JcbDiff[_nVar-1]=1/internal_energy;
1486 for(int l =1;l<1+_Ndim;l++){
1487 _JcbDiff[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
1488 for(k=0; k<_Ndim;k++)
1489 _JcbDiff[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
1490 _JcbDiff[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
1493 _JcbDiff[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
1494 for(k=0; k<_Ndim;k++)
1495 _JcbDiff[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
1496 _JcbDiff[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
1498 else if (_limitField[nameOfGroup].bcType==Outlet || _limitField[nameOfGroup].bcType==Neumann
1499 ||_limitField[nameOfGroup].bcType==Inlet || _limitField[nameOfGroup].bcType==InletPressure)
1501 for(k=0;k<_nVar;k++)
1502 _JcbDiff[k*_nVar+k]=1;
1505 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1506 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1507 _runLogFile->close();
1508 throw CdmathException("SinglePhase::jacobianDiff: This boundary condition is not recognised");
1512 void SinglePhase::primToCons(const double *P, const int &i, double *W, const int &j){
1513 //cout<<"SinglePhase::primToCons i="<<i<<", j="<<j<<", P[i*(_Ndim+2)]="<<P[i*(_Ndim+2)]<<", P[i*(_Ndim+2)+(_Ndim+1)]="<<P[i*(_Ndim+2)+(_Ndim+1)]<<endl;
1515 double rho=_fluides[0]->getDensity(P[i*(_Ndim+2)], P[i*(_Ndim+2)+(_Ndim+1)]);
1516 W[j*(_Ndim+2)] = _porosityField(j)*rho;//phi*rho
1517 for(int k=0; k<_Ndim; k++)
1518 W[j*(_Ndim+2)+(k+1)] = W[j*(_Ndim+2)]*P[i*(_Ndim+2)+(k+1)];//phi*rho*u
1520 W[j*(_Ndim+2)+(_Ndim+1)] = W[j*(_Ndim+2)]*_fluides[0]->getInternalEnergy(P[i*(_Ndim+2)+ (_Ndim+1)],rho);//rho*e
1521 for(int k=0; k<_Ndim; k++)
1522 W[j*(_Ndim+2)+(_Ndim+1)] += W[j*(_Ndim+2)]*P[i*(_Ndim+2)+(k+1)]*P[i*(_Ndim+2)+(k+1)]*0.5;//phi*rho*e+0.5*phi*rho*u^2
1525 void SinglePhase::primToConsJacobianMatrix(double *V)
1527 double pression=V[0];
1528 double temperature=V[_nVar-1];
1529 double vitesse[_Ndim];
1530 for(int idim=0;idim<_Ndim;idim++)
1531 vitesse[idim]=V[1+idim];
1533 for(int idim=0;idim<_Ndim;idim++)
1534 v2+=vitesse[idim]*vitesse[idim];
1536 double rho=_fluides[0]->getDensity(pression,temperature);
1537 double gamma=_fluides[0]->constante("gamma");
1538 double Pinf=_fluides[0]->constante("p0");
1539 double q=_fluides[0]->constante("q");
1540 double cv=_fluides[0]->constante("cv");
1542 for(int k=0;k<_nVar*_nVar; k++)
1543 _primToConsJacoMat[k]=0;
1545 if( !_useDellacherieEOS)
1547 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1548 double e=fluide0->getInternalEnergy(temperature);
1551 _primToConsJacoMat[0]=1/((gamma-1)*(e-q));
1552 _primToConsJacoMat[_nVar-1]=-rho*cv/(e-q);
1554 for(int idim=0;idim<_Ndim;idim++)
1556 _primToConsJacoMat[_nVar+idim*_nVar]=vitesse[idim]/((gamma-1)*(e-q));
1557 _primToConsJacoMat[_nVar+idim*_nVar+1+idim]=rho;
1558 _primToConsJacoMat[_nVar+idim*_nVar+_nVar-1]=-rho*vitesse[idim]*cv/(e-q);
1560 _primToConsJacoMat[(_nVar-1)*_nVar]=E/((gamma-1)*(e-q));
1561 for(int idim=0;idim<_Ndim;idim++)
1562 _primToConsJacoMat[(_nVar-1)*_nVar+1+idim]=rho*vitesse[idim];
1563 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*cv*(1-E/(e-q));
1565 else if( _useDellacherieEOS)
1567 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1568 double h=fluide0->getEnthalpy(temperature);
1570 double cp=_fluides[0]->constante("cp");
1572 _primToConsJacoMat[0]=gamma/((gamma-1)*(h-q));
1573 _primToConsJacoMat[_nVar-1]=-rho*cp/(h-q);
1575 for(int idim=0;idim<_Ndim;idim++)
1577 _primToConsJacoMat[_nVar+idim*_nVar]=gamma*vitesse[idim]/((gamma-1)*(h-q));
1578 _primToConsJacoMat[_nVar+idim*_nVar+1+idim]=rho;
1579 _primToConsJacoMat[_nVar+idim*_nVar+_nVar-1]=-rho*vitesse[idim]*cp/(h-q);
1581 _primToConsJacoMat[(_nVar-1)*_nVar]=gamma*H/((gamma-1)*(h-q))-1;
1582 for(int idim=0;idim<_Ndim;idim++)
1583 _primToConsJacoMat[(_nVar-1)*_nVar+1+idim]=rho*vitesse[idim];
1584 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*cp*(1-H/(h-q));
1588 *_runLogFile<<"SinglePhase::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie"<<endl;
1589 _runLogFile->close();
1590 throw CdmathException("SinglePhase::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
1593 if(_verbose && _nbTimeStep%_freqSave ==0)
1595 cout<<" SinglePhase::primToConsJacobianMatrix" << endl;
1596 displayVector(_Vi,_nVar," _Vi " );
1597 cout<<" Jacobienne primToCons: " << endl;
1598 displayMatrix(_primToConsJacoMat,_nVar," Jacobienne primToCons: ");
1602 void SinglePhase::consToPrim(const double *Wcons, double* Wprim,double porosity)
1605 for(int k=1;k<=_Ndim;k++)
1606 q_2 += Wcons[k]*Wcons[k];
1607 q_2 /= Wcons[0]; //phi rho u²
1608 double rhoe=(Wcons[(_Ndim+2)-1]-q_2/2)/porosity;
1609 double rho=Wcons[0]/porosity;
1610 Wprim[0] = _fluides[0]->getPressure(rhoe,rho);//pressure p
1612 cout << "pressure = "<< Wprim[0] << " < 0 " << endl;
1613 *_runLogFile<< "pressure = "<< Wprim[0] << " < 0 " << endl;
1614 _runLogFile->close();
1615 throw CdmathException("SinglePhase::consToPrim: negative pressure");
1617 for(int k=1;k<=_Ndim;k++)
1618 Wprim[k] = Wcons[k]/Wcons[0];//velocity u
1619 Wprim[(_Ndim+2)-1] = _fluides[0]->getTemperatureFromPressure(Wprim[0],Wcons[0]/porosity);
1621 if(_verbose && _nbTimeStep%_freqSave ==0)
1623 cout<<"ConsToPrim Vecteur conservatif"<<endl;
1624 for(int k=0;k<_nVar;k++)
1625 cout<<Wcons[k]<<endl;
1626 cout<<"ConsToPrim Vecteur primitif"<<endl;
1627 for(int k=0;k<_nVar;k++)
1628 cout<<Wprim[k]<<endl;
1632 void SinglePhase::entropicShift(double* n)//TO do: make sure _Vi and _Vj are well set
1634 /*Left and right values */
1635 double ul_n = 0, ul_2=0, ur_n=0, ur_2 = 0; //valeurs de vitesse normale et |u|² a droite et a gauche
1636 for(int i=0;i<_Ndim;i++)
1638 ul_n += _Vi[1+i]*n[i];
1639 ul_2 += _Vi[1+i]*_Vi[1+i];
1640 ur_n += _Vj[1+i]*n[i];
1641 ur_2 += _Vj[1+i]*_Vj[1+i];
1644 double cl = _fluides[0]->vitesseSonEnthalpie(_Vi[_Ndim+1]-ul_2/2);//vitesse du son a l'interface
1645 double cr = _fluides[0]->vitesseSonEnthalpie(_Vj[_Ndim+1]-ur_2/2);//vitesse du son a l'interface
1647 _entropicShift[0]=fabs(ul_n-cl - (ur_n-cr));
1648 _entropicShift[1]=fabs(ul_n - ur_n);
1649 _entropicShift[2]=fabs(ul_n+cl - (ur_n+cr));
1651 if(_verbose && _nbTimeStep%_freqSave ==0)
1653 cout<<"Entropic shift left eigenvalues: "<<endl;
1654 cout<<"("<< ul_n-cl<< ", " <<ul_n<<", "<<ul_n+cl << ")";
1656 cout<<"Entropic shift right eigenvalues: "<<endl;
1657 cout<<"("<< ur_n-cr<< ", " <<ur_n<<", "<<ur_n+cr << ")";
1659 cout<<"eigenvalue jumps "<<endl;
1660 cout<< _entropicShift[0] << ", " << _entropicShift[1] << ", "<< _entropicShift[2] <<endl;
1664 Vector SinglePhase::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
1665 if(_verbose && _nbTimeStep%_freqSave ==0)
1667 cout<<"SinglePhase::convectionFlux start"<<endl;
1668 cout<<"Ucons"<<endl;
1670 cout<<"Vprim"<<endl;
1674 double phirho=U(0);//phi rho
1675 Vector phiq(_Ndim);//phi rho u
1676 for(int i=0;i<_Ndim;i++)
1679 double pression=V(0);
1680 Vector vitesse(_Ndim);
1681 for(int i=0;i<_Ndim;i++)
1683 double Temperature= V(1+_Ndim);
1685 double vitessen=vitesse*normale;
1686 double rho=phirho/porosity;
1687 double e_int=_fluides[0]->getInternalEnergy(Temperature,rho);
1690 F(0)=phirho*vitessen;
1691 for(int i=0;i<_Ndim;i++)
1692 F(1+i)=phirho*vitessen*vitesse(i)+pression*normale(i)*porosity;
1693 F(1+_Ndim)=phirho*(e_int+0.5*vitesse*vitesse+pression/rho)*vitessen;
1695 if(_verbose && _nbTimeStep%_freqSave ==0)
1697 cout<<"SinglePhase::convectionFlux end"<<endl;
1698 cout<<"Flux F(U,V)"<<endl;
1705 void SinglePhase::RoeMatrixConservativeVariables(double u_n, double H,Vector velocity, double k, double K)
1707 /******** Construction de la matrice de Roe *********/
1708 //premiere ligne (masse)
1710 for(int idim=0; idim<_Ndim;idim++)
1711 _Aroe[1+idim]=_vec_normal[idim];
1714 //lignes intermadiaires(qdm)
1715 for(int idim=0; idim<_Ndim;idim++)
1718 _Aroe[(1+idim)*_nVar]=K*_vec_normal[idim] - u_n*_Uroe[1+idim];
1719 //colonnes intermediaires
1720 for(int jdim=0; jdim<_Ndim;jdim++)
1721 _Aroe[(1+idim)*_nVar + jdim + 1] = _Uroe[1+idim]*_vec_normal[jdim]-k*_vec_normal[idim]*_Uroe[1+jdim];
1723 _Aroe[(1+idim)*_nVar + idim + 1] += u_n;
1725 _Aroe[(1+idim)*_nVar + _nVar-1]=k*_vec_normal[idim];
1728 //derniere ligne (energie)
1729 _Aroe[_nVar*(_nVar-1)] = (K - H)*u_n;
1730 for(int idim=0; idim<_Ndim;idim++)
1731 _Aroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] - k*u_n*_Uroe[idim+1];
1732 _Aroe[_nVar*_nVar -1] = (1 + k)*u_n;
1734 void SinglePhase::convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector vitesse)
1736 //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1737 //On remplit la matrice de Roe en variables primitives : F(V_L)-F(V_R)=Aroe (V_L-V_R)
1738 //EOS is more involved with primitive variables
1739 // call to getDensityDerivatives(double concentration, double pression, double temperature,double v2) needed
1740 _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1741 for(int i=0;i<_Ndim;i++)
1742 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1743 _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1744 for(int i=0;i<_Ndim;i++)
1746 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]+_vec_normal[i];
1747 for(int j=0;j<_Ndim;j++)
1748 _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1749 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1750 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1752 _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1753 for(int i=0;i<_Ndim;i++)
1754 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1755 _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1757 void SinglePhase::staggeredRoeUpwindingMatrixConservativeVariables( double u_n, double H,Vector velocity, double k, double K)
1759 //Calcul de décentrement de type décalé pour formulation de Roe
1760 if(fabs(u_n)>_precision)
1762 //premiere ligne (masse)
1764 for(int idim=0; idim<_Ndim;idim++)
1765 _absAroe[1+idim]=_vec_normal[idim];
1766 _absAroe[_nVar-1]=0;
1768 //lignes intermadiaires(qdm)
1769 for(int idim=0; idim<_Ndim;idim++)
1772 _absAroe[(1+idim)*_nVar]=-K*_vec_normal[idim] - u_n*_Uroe[1+idim];
1773 //colonnes intermediaires
1774 for(int jdim=0; jdim<_Ndim;jdim++)
1775 _absAroe[(1+idim)*_nVar + jdim + 1] = _Uroe[1+idim]*_vec_normal[jdim]+k*_vec_normal[idim]*_Uroe[1+jdim];
1777 _absAroe[(1+idim)*_nVar + idim + 1] += u_n;
1779 _absAroe[(1+idim)*_nVar + _nVar-1]=-k*_vec_normal[idim];
1782 //derniere ligne (energie)
1783 _absAroe[_nVar*(_nVar-1)] = (-K - H)*u_n;
1784 for(int idim=0; idim<_Ndim;idim++)
1785 _absAroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] + k*u_n*_Uroe[idim+1];
1786 _absAroe[_nVar*_nVar -1] = (1 - k)*u_n;
1794 for(int i=0; i<_nVar*_nVar;i++)
1795 _absAroe[i] *= signu;
1797 else//umn=0 ->centered scheme
1799 for(int i=0; i<_nVar*_nVar;i++)
1803 void SinglePhase::staggeredRoeUpwindingMatrixPrimitiveVariables(double rho, double u_n,double H, Vector vitesse)
1805 //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1806 //Calcul de décentrement de type décalé pour formulation Roe
1807 _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1808 for(int i=0;i<_Ndim;i++)
1809 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1810 _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1811 for(int i=0;i<_Ndim;i++)
1813 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]-_vec_normal[i];
1814 for(int j=0;j<_Ndim;j++)
1815 _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1816 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1817 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1819 _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1820 for(int i=0;i<_Ndim;i++)
1821 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1822 _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1825 Vector SinglePhase::staggeredVFFCFlux()
1827 if(_verbose && _nbTimeStep%_freqSave ==0)
1828 cout<<"SinglePhase::staggeredVFFCFlux() start"<<endl;
1830 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1832 *_runLogFile<< "SinglePhase::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding, pressure = "<< endl;
1833 _runLogFile->close();
1834 throw CdmathException("SinglePhase::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding");
1836 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
1840 double uijn=0, phiqn=0, uin=0, ujn=0;
1841 for(int idim=0;idim<_Ndim;idim++)
1843 uijn+=_vec_normal[idim]*_Uroe[1+idim];//URoe = rho, u, H
1844 uin +=_vec_normal[idim]*_Ui[1+idim];
1845 ujn +=_vec_normal[idim]*_Uj[1+idim];
1848 if( (uin>0 && ujn >0) || (uin>=0 && ujn <=0 && uijn>0) ) // formerly (uijn>_precision)
1850 for(int idim=0;idim<_Ndim;idim++)
1851 phiqn+=_vec_normal[idim]*_Ui[1+idim];//phi rho u n
1853 for(int idim=0;idim<_Ndim;idim++)
1854 Fij(1+idim)=phiqn*_Vi[1+idim]+_Vj[0]*_vec_normal[idim]*_porosityj;
1855 Fij(_nVar-1)=phiqn/_Ui[0]*(_Ui[_nVar-1]+_Vj[0]*sqrt(_porosityj/_porosityi));
1857 else if( (uin<0 && ujn <0) || (uin>=0 && ujn <=0 && uijn<0) ) // formerly (uijn<-_precision)
1859 for(int idim=0;idim<_Ndim;idim++)
1860 phiqn+=_vec_normal[idim]*_Uj[1+idim];//phi rho u n
1862 for(int idim=0;idim<_Ndim;idim++)
1863 Fij(1+idim)=phiqn*_Vj[1+idim]+_Vi[0]*_vec_normal[idim]*_porosityi;
1864 Fij(_nVar-1)=phiqn/_Uj[0]*(_Uj[_nVar-1]+_Vi[0]*sqrt(_porosityi/_porosityj));
1866 else//case (uin<=0 && ujn >=0) or (uin>=0 && ujn <=0 && uijn==0), apply centered scheme
1868 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar);
1869 Vector normale(_Ndim);
1870 for(int i1=0;i1<_Ndim;i1++)
1871 normale(i1)=_vec_normal[i1];
1872 for(int i1=0;i1<_nVar;i1++)
1879 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
1880 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
1881 Fij=(Fi+Fj)/2;//+_maxvploc*(Ui-Uj)/2;
1883 if(_verbose && _nbTimeStep%_freqSave ==0)
1885 cout<<"SinglePhase::staggeredVFFCFlux() endf uijn="<<uijn<<endl;
1892 void SinglePhase::staggeredVFFCMatricesConservativeVariables(double un)//vitesse normale de Roe en entrée
1894 if(_verbose && _nbTimeStep%_freqSave ==0)
1895 cout<<"SinglePhase::staggeredVFFCMatrices()"<<endl;
1897 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1899 *_runLogFile<< "SinglePhase::staggeredVFFCMatrices: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding"<< endl;
1900 _runLogFile->close();
1901 throw CdmathException("SinglePhase::staggeredVFFCMatrices: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding");
1903 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
1905 double ui_n=0, ui_2=0, uj_n=0, uj_2=0;//vitesse normale et carré du module
1906 double H;//enthalpie totale (expression particulière)
1907 consToPrim(_Ui,_Vi,_porosityi);
1908 consToPrim(_Uj,_Vj,_porosityj);
1909 for(int i=0;i<_Ndim;i++)
1911 ui_2 += _Vi[1+i]*_Vi[1+i];
1912 ui_n += _Vi[1+i]*_vec_normal[i];
1913 uj_2 += _Vj[1+i]*_Vj[1+i];
1914 uj_n += _Vj[1+i]*_vec_normal[i];
1917 double rhoi,pj, Ei, rhoj;
1918 double cj, Kj, kj;//dérivées de la pression
1919 rhoi=_Ui[0]/_porosityi;
1920 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
1922 rhoj=_Uj[0]/_porosityj;
1924 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
1925 kj = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
1926 Kj = uj_2*kj/2; //g-1/2 *|u|²
1929 double ci, Ki, ki;//dérivées de la pression
1930 Ej= _Uj[_Ndim+1]/rhoj;
1933 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
1934 ki = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
1935 Ki = ui_2*ki/2; //g-1/2 *|u|²
1939 /***********Calcul des valeurs propres ********/
1940 vector<complex<double> > vp_dist(3,0);
1941 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
1943 _maxvploc=fabs(ui_n)+cj;
1944 if(_maxvploc>_maxvp)
1947 if(_verbose && _nbTimeStep%_freqSave ==0)
1948 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
1950 /******** Construction de la matrice A^+ *********/
1951 //premiere ligne (masse)
1953 for(int idim=0; idim<_Ndim;idim++)
1954 _AroePlus[1+idim]=_vec_normal[idim];
1955 _AroePlus[_nVar-1]=0;
1957 //lignes intermadiaires(qdm)
1958 for(int idim=0; idim<_Ndim;idim++)
1961 _AroePlus[(1+idim)*_nVar]=- ui_n*_Vi[1+idim];
1962 //colonnes intermediaires
1963 for(int jdim=0; jdim<_Ndim;jdim++)
1964 _AroePlus[(1+idim)*_nVar + jdim + 1] = _Vi[1+idim]*_vec_normal[jdim];
1966 _AroePlus[(1+idim)*_nVar + idim + 1] += ui_n;
1968 _AroePlus[(1+idim)*_nVar + _nVar-1]=0;
1971 //derniere ligne (energie)
1972 _AroePlus[_nVar*(_nVar-1)] = - H*ui_n;
1973 for(int idim=0; idim<_Ndim;idim++)
1974 _AroePlus[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] ;
1975 _AroePlus[_nVar*_nVar -1] = ui_n;
1977 /******** Construction de la matrice A^- *********/
1978 //premiere ligne (masse)
1980 for(int idim=0; idim<_Ndim;idim++)
1981 _AroeMinus[1+idim]=0;
1982 _AroeMinus[_nVar-1]=0;
1984 //lignes intermadiaires(qdm)
1985 for(int idim=0; idim<_Ndim;idim++)
1988 _AroeMinus[(1+idim)*_nVar]=Kj*_vec_normal[idim] ;
1989 //colonnes intermediaires
1990 for(int jdim=0; jdim<_Ndim;jdim++)
1991 _AroeMinus[(1+idim)*_nVar + jdim + 1] = -kj*_vec_normal[idim]*_Vj[1+jdim];
1993 _AroeMinus[(1+idim)*_nVar + idim + 1] += 0;
1995 _AroeMinus[(1+idim)*_nVar + _nVar-1]=kj*_vec_normal[idim];
1998 //derniere ligne (energie)
1999 _AroeMinus[_nVar*(_nVar-1)] = Kj *ui_n;
2000 for(int idim=0; idim<_Ndim;idim++)
2001 _AroeMinus[_nVar*(_nVar-1)+idim+1]= - kj*uj_n*_Vi[idim+1];
2002 _AroeMinus[_nVar*_nVar -1] = kj*ui_n;
2004 else if(un<-_precision)
2006 /***********Calcul des valeurs propres ********/
2007 vector<complex<double> > vp_dist(3,0);
2008 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2010 _maxvploc=fabs(uj_n)+ci;
2011 if(_maxvploc>_maxvp)
2014 if(_verbose && _nbTimeStep%_freqSave ==0)
2015 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2017 /******** Construction de la matrice A^+ *********/
2018 //premiere ligne (masse)
2020 for(int idim=0; idim<_Ndim;idim++)
2021 _AroePlus[1+idim]=0;
2022 _AroePlus[_nVar-1]=0;
2024 //lignes intermadiaires(qdm)
2025 for(int idim=0; idim<_Ndim;idim++)
2028 _AroePlus[(1+idim)*_nVar]=Ki*_vec_normal[idim] ;
2029 //colonnes intermediaires
2030 for(int jdim=0; jdim<_Ndim;jdim++)
2031 _AroePlus[(1+idim)*_nVar + jdim + 1] = -ki*_vec_normal[idim]*_Vi[1+jdim];
2033 _AroePlus[(1+idim)*_nVar + idim + 1] += 0;
2035 _AroePlus[(1+idim)*_nVar + _nVar-1]=ki*_vec_normal[idim];
2038 //derniere ligne (energie)
2039 _AroePlus[_nVar*(_nVar-1)] = Ki *uj_n;
2040 for(int idim=0; idim<_Ndim;idim++)
2041 _AroePlus[_nVar*(_nVar-1)+idim+1]= - ki*ui_n*_Vj[idim+1];
2042 _AroePlus[_nVar*_nVar -1] = ki*uj_n;
2044 /******** Construction de la matrice A^- *********/
2045 //premiere ligne (masse)
2047 for(int idim=0; idim<_Ndim;idim++)
2048 _AroeMinus[1+idim]=_vec_normal[idim];
2049 _AroeMinus[_nVar-1]=0;
2051 //lignes intermadiaires(qdm)
2052 for(int idim=0; idim<_Ndim;idim++)
2055 _AroeMinus[(1+idim)*_nVar]= - uj_n*_Vj[1+idim];
2056 //colonnes intermediaires
2057 for(int jdim=0; jdim<_Ndim;jdim++)
2058 _AroeMinus[(1+idim)*_nVar + jdim + 1] = _Vj[1+idim]*_vec_normal[jdim];
2060 _AroeMinus[(1+idim)*_nVar + idim + 1] += uj_n;
2062 _AroeMinus[(1+idim)*_nVar + _nVar-1]=0;
2065 //derniere ligne (energie)
2066 _AroeMinus[_nVar*(_nVar-1)] = - H*uj_n;
2067 for(int idim=0; idim<_Ndim;idim++)
2068 _AroeMinus[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] ;
2069 _AroeMinus[_nVar*_nVar -1] = uj_n;
2071 else//case nil velocity on the interface, apply centered scheme
2073 double u_n=0, u_2=0;//vitesse normale et carré du module
2074 for(int i=0;i<_Ndim;i++)
2076 u_2 += _Uroe[1+i]*_Uroe[1+i];
2077 u_n += _Uroe[1+i]*_vec_normal[i];
2079 Vector vitesse(_Ndim);
2080 for(int idim=0;idim<_Ndim;idim++)
2081 vitesse[idim]=_Uroe[1+idim];
2084 /***********Calcul des valeurs propres ********/
2086 c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
2087 k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
2088 K = u_2*k/2; //g-1/2 *|u|²
2090 _maxvploc=fabs(u_n)+c;
2091 if(_maxvploc>_maxvp)
2094 /******** Construction de la matrice A^+ *********/
2095 //premiere ligne (masse)
2097 for(int idim=0; idim<_Ndim;idim++)
2098 _AroePlus[1+idim]=0;
2099 _AroePlus[_nVar-1]=0;
2101 //lignes intermadiaires(qdm)
2102 for(int idim=0; idim<_Ndim;idim++)
2105 _AroePlus[(1+idim)*_nVar]=- ui_n*_Vi[1+idim];
2106 //colonnes intermediaires
2107 for(int jdim=0; jdim<_Ndim;jdim++)
2108 _AroePlus[(1+idim)*_nVar + jdim + 1] = _Vi[1+idim]*_vec_normal[jdim]-0.5*ki*_vec_normal[idim]*_Vi[1+jdim];
2110 _AroePlus[(1+idim)*_nVar + idim + 1] += 0.5*ui_n;
2112 _AroePlus[(1+idim)*_nVar + _nVar-1]=0.5*ki*_vec_normal[idim];
2115 //derniere ligne (energie)
2116 _AroePlus[_nVar*(_nVar-1)] = 0;
2117 for(int idim=0; idim<_Ndim;idim++)
2118 _AroePlus[_nVar*(_nVar-1)+idim+1]=0 ;
2119 _AroePlus[_nVar*_nVar -1] = 0;
2121 /******** Construction de la matrice A^- *********/
2122 //premiere ligne (masse)
2124 for(int idim=0; idim<_Ndim;idim++)
2125 _AroeMinus[1+idim]=0;
2126 _AroeMinus[_nVar-1]=0;
2128 //lignes intermadiaires(qdm)
2129 for(int idim=0; idim<_Ndim;idim++)
2132 _AroeMinus[(1+idim)*_nVar]=Kj*_vec_normal[idim] - uj_n*_Vj[1+idim];
2133 //colonnes intermediaires
2134 for(int jdim=0; jdim<_Ndim;jdim++)
2135 _AroeMinus[(1+idim)*_nVar + jdim + 1] = -0.5*kj*_vec_normal[idim]*_Vj[1+jdim];
2137 _AroeMinus[(1+idim)*_nVar + idim + 1] += 0.5*uj_n;
2139 _AroeMinus[(1+idim)*_nVar + _nVar-1]=0.5*kj*_vec_normal[idim];
2142 //derniere ligne (energie)
2143 _AroeMinus[_nVar*(_nVar-1)] = 0;
2144 for(int idim=0; idim<_Ndim;idim++)
2145 _AroeMinus[_nVar*(_nVar-1)+idim+1]= 0;
2146 _AroeMinus[_nVar*_nVar -1] = 0;
2149 if(_timeScheme==Implicit)
2150 for(int i=0; i<_nVar*_nVar;i++)
2152 _AroeMinusImplicit[i] = _AroeMinus[i];
2153 _AroePlusImplicit[i] = _AroePlus[i];
2156 /******** Construction de la matrices Aroe *********/
2158 //premiere ligne (masse)
2160 for(int idim=0; idim<_Ndim;idim++)
2161 _Aroe[1+idim]=_vec_normal[idim];
2164 //lignes intermadiaires(qdm)
2165 for(int idim=0; idim<_Ndim;idim++)
2168 _Aroe[(1+idim)*_nVar]=Ki*_vec_normal[idim] - uj_n*_Uj[1+idim];
2169 //colonnes intermediaires
2170 for(int jdim=0; jdim<_Ndim;jdim++)
2171 _Aroe[(1+idim)*_nVar + jdim + 1] = _Uj[1+idim]*_vec_normal[jdim]-ki*_vec_normal[idim]*_Ui[1+jdim];
2173 _Aroe[(1+idim)*_nVar + idim + 1] += uj_n;
2175 _Aroe[(1+idim)*_nVar + _nVar-1]=ki*_vec_normal[idim];
2178 //derniere ligne (energie)
2179 _Aroe[_nVar*(_nVar-1)] = (Ki - H)*uj_n;
2180 for(int idim=0; idim<_Ndim;idim++)
2181 _Aroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] - ki*ui_n*_Uj[idim+1];
2182 _Aroe[_nVar*_nVar -1] = (1 + ki)*uj_n;
2186 void SinglePhase::staggeredVFFCMatricesPrimitiveVariables(double un)//vitesse normale de Roe en entrée
2188 if(_verbose && _nbTimeStep%_freqSave ==0)
2189 cout<<"SinglePhase::staggeredVFFCMatricesPrimitiveVariables()"<<endl;
2191 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2193 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding" << endl;
2194 _runLogFile->close();
2195 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding");
2197 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
2199 double ui_n=0., ui_2=0., uj_n=0., uj_2=0.;//vitesse normale et carré du module
2200 double H;//enthalpie totale (expression particulière)
2201 consToPrim(_Ui,_Vi,_porosityi);
2202 consToPrim(_Uj,_Vj,_porosityj);
2204 for(int i=0;i<_Ndim;i++)
2206 ui_2 += _Vi[1+i]*_Vi[1+i];
2207 ui_n += _Vi[1+i]*_vec_normal[i];
2208 uj_2 += _Vj[1+i]*_Vj[1+i];
2209 uj_n += _Vj[1+i]*_vec_normal[i];
2212 if(_verbose && _nbTimeStep%_freqSave ==0){
2213 cout <<"SinglePhase::staggeredVFFCMatricesPrimitiveVariables " << endl;
2214 cout<<"Vecteur primitif _Vi" << endl;
2215 for(int i=0;i<_nVar;i++)
2218 cout<<"Vecteur primitif _Vj" << endl;
2219 for(int i=0;i<_nVar;i++)
2224 double gamma=_fluides[0]->constante("gamma");
2225 double q=_fluides[0]->constante("q");
2227 if(fabs(un)>_precision)//non zero velocity on the interface
2229 if( !_useDellacherieEOS)
2231 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2232 double cv=fluide0->constante("cv");
2236 double rhoi,rhoj,pj, Ei, ei;
2237 double cj;//vitesse du son pour calcul valeurs propres
2238 rhoi=_Ui[0]/_porosityi;
2239 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
2242 rhoj=_Uj[0]/_porosityj;
2243 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
2245 /***********Calcul des valeurs propres ********/
2246 vector<complex<double> > vp_dist(3,0);
2247 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
2249 _maxvploc=fabs(ui_n)+cj;
2250 if(_maxvploc>_maxvp)
2253 if(_verbose && _nbTimeStep%_freqSave ==0)
2254 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
2256 /******** Construction de la matrice A^+ *********/
2257 //premiere ligne (masse)
2258 _AroePlusImplicit[0]=ui_n/((gamma-1)*(ei-q));
2259 for(int idim=0; idim<_Ndim;idim++)
2260 _AroePlusImplicit[1+idim]=rhoi*_vec_normal[idim];
2261 _AroePlusImplicit[_nVar-1]=-rhoi*ui_n*cv/(ei-q);
2263 //lignes intermadiaires(qdm)
2264 for(int idim=0; idim<_Ndim;idim++)
2267 _AroePlusImplicit[(1+idim)*_nVar]=ui_n/((gamma-1)*(ei-q))*_Vi[1+idim];
2268 //colonnes intermediaires
2269 for(int jdim=0; jdim<_Ndim;jdim++)
2270 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = rhoi*_Vi[1+idim]*_vec_normal[jdim];
2272 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += rhoi*ui_n;
2274 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoi*ui_n*cv/(ei-q)*_Vi[1+idim];
2277 //derniere ligne (energie)
2278 _AroePlusImplicit[_nVar*(_nVar-1)] = Ei*ui_n/((gamma-1)*(ei-q));
2279 for(int idim=0; idim<_Ndim;idim++)
2280 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoi*Ei+pj)*_vec_normal[idim]+rhoi*ui_n*_Vi[1+idim];
2281 _AroePlusImplicit[_nVar*_nVar -1] = rhoi*ui_n*(1-Ei/(ei-q))*cv;
2283 /******** Construction de la matrice A^- *********/
2284 //premiere ligne (masse)
2285 _AroeMinusImplicit[0]=0;
2286 for(int idim=0; idim<_Ndim;idim++)
2287 _AroeMinusImplicit[1+idim]=0;
2288 _AroeMinusImplicit[_nVar-1]=0;
2290 //lignes intermadiaires(qdm)
2291 for(int idim=0; idim<_Ndim;idim++)
2294 _AroeMinusImplicit[(1+idim)*_nVar]=_vec_normal[idim] ;
2295 //colonnes intermediaires
2296 for(int jdim=0; jdim<_Ndim;jdim++)
2297 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2299 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2301 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2304 //derniere ligne (energie)
2305 _AroeMinusImplicit[_nVar*(_nVar-1)] = ui_n;
2306 for(int idim=0; idim<_Ndim;idim++)
2307 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2308 _AroeMinusImplicit[_nVar*_nVar -1] = 0;
2310 else if(un<-_precision)
2312 double pi, Ej, rhoi, rhoj, ej;
2313 double ci;//vitesse du son pour calcul valeurs propres
2314 rhoj=_Uj[0]/_porosityj;
2315 Ej= _Uj[_Ndim+1]/rhoj;
2318 rhoi=_Ui[0]/_porosityi;
2319 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
2321 /***********Calcul des valeurs propres ********/
2322 vector<complex<double> > vp_dist(3,0);
2323 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2325 _maxvploc=fabs(uj_n)+ci;
2326 if(_maxvploc>_maxvp)
2329 if(_verbose && _nbTimeStep%_freqSave ==0)
2330 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2332 /******** Construction de la matrice A^+ *********/
2333 //premiere ligne (masse)
2334 _AroePlusImplicit[0]=0;
2335 for(int idim=0; idim<_Ndim;idim++)
2336 _AroePlusImplicit[1+idim]=0;
2337 _AroePlusImplicit[_nVar-1]=0;
2339 //lignes intermadiaires(qdm)
2340 for(int idim=0; idim<_Ndim;idim++)
2343 _AroePlusImplicit[(1+idim)*_nVar]=0;
2344 //colonnes intermediaires
2345 for(int jdim=0; jdim<_Ndim;jdim++)
2346 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2348 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2350 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2353 //derniere ligne (energie)
2354 _AroePlusImplicit[_nVar*(_nVar-1)] = uj_n;
2355 for(int idim=0; idim<_Ndim;idim++)
2356 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2357 _AroePlusImplicit[_nVar*_nVar -1] = 0;
2359 /******** Construction de la matrice A^- *********/
2360 //premiere ligne (masse)
2361 _AroeMinusImplicit[0]=uj_n/((gamma-1)*(ej-q));
2362 for(int idim=0; idim<_Ndim;idim++)
2363 _AroeMinusImplicit[1+idim]=rhoj*_vec_normal[idim];
2364 _AroeMinusImplicit[_nVar-1]=-rhoj*uj_n*cv/(ej-q);
2366 //lignes intermadiaires(qdm)
2367 for(int idim=0; idim<_Ndim;idim++)
2370 _AroeMinusImplicit[(1+idim)*_nVar]= uj_n/((gamma-1)*(ej-q))*_Vj[1+idim];
2371 //colonnes intermediaires
2372 for(int jdim=0; jdim<_Ndim;jdim++)
2373 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = rhoj*_Vj[1+idim]*_vec_normal[jdim];
2375 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += rhoj*uj_n;
2377 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoj*uj_n*cv/(ej-q)*_Vj[1+idim];
2380 //derniere ligne (energie)
2381 _AroeMinusImplicit[_nVar*(_nVar-1)] = Ej*uj_n/((gamma-1)*(ej-q));
2382 for(int idim=0; idim<_Ndim;idim++)
2383 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoj*Ej+pi)*_vec_normal[idim]+rhoj*uj_n*_Vj[1+idim];
2384 _AroeMinusImplicit[_nVar*_nVar -1] = rhoj*uj_n*(1-Ej/(ej-q))*cv;
2388 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero" << endl;
2389 _runLogFile->close();
2390 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero");
2393 else if(_useDellacherieEOS )
2395 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2396 double cp=fluide0->constante("cp");
2400 double rhoi,rhoj,pj, Ei, hi, Hi;
2401 double cj;//vitesse du son pour calcul valeurs propres
2402 rhoi=_Ui[0]/_porosityi;
2403 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
2407 rhoj=_Uj[0]/_porosityj;
2408 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
2410 /***********Calcul des valeurs propres ********/
2411 vector<complex<double> > vp_dist(3,0);
2412 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
2414 _maxvploc=fabs(ui_n)+cj;
2415 if(_maxvploc>_maxvp)
2418 if(_verbose && _nbTimeStep%_freqSave ==0)
2419 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
2421 /******** Construction de la matrice A^+ *********/
2422 //premiere ligne (masse)
2423 _AroePlusImplicit[0]=ui_n*gamma/((gamma-1)*(hi-q));
2424 for(int idim=0; idim<_Ndim;idim++)
2425 _AroePlusImplicit[1+idim]=rhoi*_vec_normal[idim];
2426 _AroePlusImplicit[_nVar-1]=-rhoi*ui_n*cp/(hi-q);
2428 //lignes intermadiaires(qdm)
2429 for(int idim=0; idim<_Ndim;idim++)
2432 _AroePlusImplicit[(1+idim)*_nVar]=ui_n*gamma/((gamma-1)*(hi-q))*_Vi[1+idim];
2433 //colonnes intermediaires
2434 for(int jdim=0; jdim<_Ndim;jdim++)
2435 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = rhoi*_Vi[1+idim]*_vec_normal[jdim];
2437 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += rhoi*ui_n;
2439 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoi*ui_n*cp/(hi-q)*_Vi[1+idim];
2442 //derniere ligne (energie)
2443 _AroePlusImplicit[_nVar*(_nVar-1)] = ui_n*(Hi*gamma/((gamma-1)*(hi-q))-1);
2444 for(int idim=0; idim<_Ndim;idim++)
2445 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoi*Ei+pj)*_vec_normal[idim]+rhoi*ui_n*_Vi[1+idim];
2446 _AroePlusImplicit[_nVar*_nVar -1] = rhoi*ui_n*(1-Hi/(hi-q))*cp;
2448 /******** Construction de la matrice A^- *********/
2449 //premiere ligne (masse)
2450 _AroeMinusImplicit[0]=0;
2451 for(int idim=0; idim<_Ndim;idim++)
2452 _AroeMinusImplicit[1+idim]=0;
2453 _AroeMinusImplicit[_nVar-1]=0;
2455 //lignes intermadiaires(qdm)
2456 for(int idim=0; idim<_Ndim;idim++)
2459 _AroeMinusImplicit[(1+idim)*_nVar]=_vec_normal[idim] ;
2460 //colonnes intermediaires
2461 for(int jdim=0; jdim<_Ndim;jdim++)
2462 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2464 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2466 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2469 //derniere ligne (energie)
2470 _AroeMinusImplicit[_nVar*(_nVar-1)] = ui_n;
2471 for(int idim=0; idim<_Ndim;idim++)
2472 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2473 _AroeMinusImplicit[_nVar*_nVar -1] = 0;
2475 else if(un<-_precision)
2477 double pi, Ej, rhoi,rhoj, Hj, hj;
2478 double ci;//vitesse du son pour calcul valeurs propres
2479 rhoj=_Uj[0]/_porosityj;
2480 Ej= _Uj[_Ndim+1]/rhoj;
2484 rhoi=_Ui[0]/_porosityi;
2485 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
2487 /***********Calcul des valeurs propres ********/
2488 vector<complex<double> > vp_dist(3,0);
2489 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2491 _maxvploc=fabs(uj_n)+ci;
2492 if(_maxvploc>_maxvp)
2495 if(_verbose && _nbTimeStep%_freqSave ==0)
2496 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2498 /******** Construction de la matrice A^+ *********/
2499 //premiere ligne (masse)
2500 _AroePlusImplicit[0]=0;
2501 for(int idim=0; idim<_Ndim;idim++)
2502 _AroePlusImplicit[1+idim]=0;
2503 _AroePlusImplicit[_nVar-1]=0;
2505 //lignes intermadiaires(qdm)
2506 for(int idim=0; idim<_Ndim;idim++)
2509 _AroePlusImplicit[(1+idim)*_nVar]=0;
2510 //colonnes intermediaires
2511 for(int jdim=0; jdim<_Ndim;jdim++)
2512 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2514 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2516 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2519 //derniere ligne (energie)
2520 _AroePlusImplicit[_nVar*(_nVar-1)] = uj_n;
2521 for(int idim=0; idim<_Ndim;idim++)
2522 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2523 _AroePlusImplicit[_nVar*_nVar -1] = 0;
2525 /******** Construction de la matrice A^- *********/
2526 //premiere ligne (masse)
2527 _AroeMinusImplicit[0]=uj_n*gamma/((gamma-1)*(hj-q));
2528 for(int idim=0; idim<_Ndim;idim++)
2529 _AroeMinusImplicit[1+idim]=rhoj*_vec_normal[idim];
2530 _AroeMinusImplicit[_nVar-1]=-rhoj*uj_n*cp/(hj-q);
2532 //lignes intermadiaires(qdm)
2533 for(int idim=0; idim<_Ndim;idim++)
2536 _AroeMinusImplicit[(1+idim)*_nVar]= uj_n*gamma/((gamma-1)*(hj-q))*_Vj[1+idim];
2537 //colonnes intermediaires
2538 for(int jdim=0; jdim<_Ndim;jdim++)
2539 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = rhoj*_Vj[1+idim]*_vec_normal[jdim];
2541 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += rhoj*uj_n;
2543 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoj*uj_n*cp/(hj-q)*_Vj[1+idim];
2546 //derniere ligne (energie)
2547 _AroeMinusImplicit[_nVar*(_nVar-1)] = uj_n*(Hj*gamma/((gamma-1)*(hj-q))-1);
2548 for(int idim=0; idim<_Ndim;idim++)
2549 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoj*Ej+pi)*_vec_normal[idim]+rhoj*uj_n*_Vj[1+idim];
2550 _AroeMinusImplicit[_nVar*_nVar -1] = rhoj*uj_n*(1-Hj/(hj-q))*cp;
2554 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero" << endl;
2555 _runLogFile->close();
2556 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero");
2561 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
2562 _runLogFile->close();
2563 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
2566 else//case nil velocity on the interface, apply centered scheme
2568 primToConsJacobianMatrix(_Vj);
2569 Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
2570 primToConsJacobianMatrix(_Vi);
2571 Polynoms::matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
2575 void SinglePhase::applyVFRoeLowMachCorrections(bool isBord, string groupname)
2577 if(_nonLinearFormulation!=VFRoe)
2579 *_runLogFile<< "SinglePhase::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation" << endl;
2580 _runLogFile->close();
2581 throw CdmathException("SinglePhase::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation");
2583 else//_nonLinearFormulation==VFRoe
2585 if(_spaceScheme==lowMach){
2587 for(int i=0;i<_Ndim;i++)
2588 u_2 += _Uroe[1+i]*_Uroe[1+i];
2589 double c = _maxvploc;//vitesse du son a l'interface
2590 double M=sqrt(u_2)/c;//Mach number
2591 _Vij[0]=M*_Vij[0]+(1-M)*(_Vi[0]+_Vj[0])/2;
2593 else if(_spaceScheme==pressureCorrection)
2594 {//order 1 : no correction, oarder 2 : correction everywhere, order 3 : correction only inside, orders 4 and 5 : special correction at boundaries
2595 if(_pressureCorrectionOrder==2 || (!isBord && _pressureCorrectionOrder==3) || (!isBord && _pressureCorrectionOrder==4) || (!isBord && _pressureCorrectionOrder==5) )
2597 double norm_uij=0, uij_n=0, ui_n=0, uj_n=0;
2598 for(int i=0;i<_Ndim;i++)
2600 norm_uij += _Uroe[1+i]*_Uroe[1+i];
2601 uij_n += _Uroe[1+i]*_vec_normal[i];
2602 ui_n += _Vi[1+i]*_vec_normal[i];
2603 uj_n += _Vj[1+i]*_vec_normal[i];
2605 norm_uij=sqrt(norm_uij);
2606 if(norm_uij>_precision)//avoid division by zero
2607 _Vij[0]=(_Vi[0]+_Vj[0])/2 + uij_n/norm_uij*(_Vj[0]-_Vi[0])/4 - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
2609 _Vij[0]=(_Vi[0]+_Vj[0])/2 - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
2611 else if(_pressureCorrectionOrder==4 && isBord)
2613 else if(_pressureCorrectionOrder==5 && isBord)
2615 double g_n=0;//scalar product of gravity and normal vector
2616 for(int i=0;i<_Ndim;i++)
2617 g_n += _GravityField3d[i]*_vec_normal[i];
2618 _Vij[0]=_Vi[0]- _Ui[0]*g_n/_inv_dxi/2;
2621 else if(_spaceScheme==staggered)
2624 for(int i=0;i<_Ndim;i++)
2625 uij_n += _Uroe[1+i]*_vec_normal[i];
2627 if(uij_n>_precision){
2629 for(int i=0;i<_Ndim;i++)
2631 _Vij[_nVar-1]=_Vi[_nVar-1];
2633 else if(uij_n<-_precision){
2635 for(int i=0;i<_Ndim;i++)
2637 _Vij[_nVar-1]=_Vj[_nVar-1];
2640 _Vij[0]=(_Vi[0]+_Vi[0])/2;
2641 for(int i=0;i<_Ndim;i++)
2642 _Vij[1+i]=(_Vj[1+i]+_Vj[1+i])/2;
2643 _Vij[_nVar-1]=(_Vj[_nVar-1]+_Vj[_nVar-1])/2;
2646 primToCons(_Vij,0,_Uij,0);
2650 void SinglePhase::testConservation()
2652 double SUM, DELTA, x;
2654 for(int i=0; i<_nVar; i++)
2658 cout << "Masse totale (kg): ";
2662 cout << "Energie totale (J): ";
2664 cout << "Quantite de mouvement totale (kg.m.s^-1): ";
2670 for(int j=0; j<_Nmailles; j++)
2672 if(!_usePrimitiveVarsInNewton)
2673 VecGetValues(_conservativeVars, 1, &I, &x);//on recupere la valeur du champ
2675 VecGetValues(_primitiveVars, 1, &I, &x);//on recupere la valeur du champ
2676 SUM += x*_mesh.getCell(j).getMeasure();
2677 VecGetValues(_newtonVariation, 1, &I, &x);//on recupere la variation du champ
2678 DELTA += x*_mesh.getCell(j).getMeasure();
2681 if(fabs(SUM)>_precision)
2682 cout << SUM << ", variation relative: " << fabs(DELTA /SUM) << endl;
2684 cout << " a une somme quasi nulle, variation absolue: " << fabs(DELTA) << endl;
2688 void SinglePhase::getDensityDerivatives( double pressure, double temperature, double v2)
2690 double rho=_fluides[0]->getDensity(pressure,temperature);
2691 double gamma=_fluides[0]->constante("gamma");
2692 double q=_fluides[0]->constante("q");
2694 if( !_useDellacherieEOS)
2696 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2697 double e = fluide0->getInternalEnergy(temperature);
2698 double cv=fluide0->constante("cv");
2701 _drho_sur_dp=1/((gamma-1)*(e-q));
2702 _drho_sur_dT=-rho*cv/(e-q);
2703 _drhoE_sur_dp=E/((gamma-1)*(e-q));
2704 _drhoE_sur_dT=rho*cv*(1-E/(e-q));
2706 else if(_useDellacherieEOS )
2708 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2709 double h=fluide0->getEnthalpy(temperature);
2711 double cp=fluide0->constante("cp");
2713 _drho_sur_dp=gamma/((gamma-1)*(h-q));
2714 _drho_sur_dT=-rho*cp/(h-q);
2715 _drhoE_sur_dp=gamma*H/((gamma-1)*(h-q))-1;
2716 _drhoE_sur_dT=rho*cp*(1-H/(h-q));
2720 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
2721 _runLogFile->close();
2722 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
2725 if(_verbose && _nbTimeStep%_freqSave ==0)
2727 cout<<"_drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
2728 cout<<"_drhoE_sur_dp= "<<_drhoE_sur_dp<<", _drhoE_sur_dT= "<<_drhoE_sur_dT<<endl;
2731 void SinglePhase::save(){
2732 string prim(_path+"/SinglePhasePrim_");///Results
2733 string cons(_path+"/SinglePhaseCons_");
2734 string allFields(_path+"/");
2737 allFields+=_fileName;
2740 for (long i = 0; i < _Nmailles; i++){
2741 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
2742 for (int j = 0; j < _nVar; j++){
2744 VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
2747 if(_saveConservativeField){
2748 for (long i = 0; i < _Nmailles; i++){
2749 // j = 0 : density; j = _nVar - 1 : energy j = 1,..,_nVar-2: momentum
2750 for (int j = 0; j < _nVar; j++){
2752 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
2755 _UU.setTime(_time,_nbTimeStep);
2757 _VV.setTime(_time,_nbTimeStep);
2759 // create mesh and component info
2760 if (_nbTimeStep ==0 || _restartWithNewFileName){
2761 string prim_suppress ="rm -rf "+prim+"_*";
2762 string cons_suppress ="rm -rf "+cons+"_*";
2764 system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
2765 system(cons_suppress.c_str());//Nettoyage des précédents calculs identiques
2767 if(_saveConservativeField){
2768 _UU.setInfoOnComponent(0,"Density_(kg/m^3)");
2769 _UU.setInfoOnComponent(1,"Momentum_x");// (kg/m^2/s)
2771 _UU.setInfoOnComponent(2,"Momentum_y");// (kg/m^2/s)
2773 _UU.setInfoOnComponent(3,"Momentum_z");// (kg/m^2/s)
2775 _UU.setInfoOnComponent(_nVar-1,"Energy_(J/m^3)");
2790 _VV.setInfoOnComponent(0,"Pressure_(Pa)");
2791 _VV.setInfoOnComponent(1,"Velocity_x_(m/s)");
2793 _VV.setInfoOnComponent(2,"Velocity_y_(m/s)");
2795 _VV.setInfoOnComponent(3,"Velocity_z_(m/s)");
2796 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
2812 // do not create mesh
2817 _VV.writeVTK(prim,false);
2820 _VV.writeMED(prim,false);
2826 if(_saveConservativeField){
2830 _UU.writeVTK(cons,false);
2833 _UU.writeMED(cons,false);
2841 if(_saveVelocity || _saveAllFields){
2842 for (long i = 0; i < _Nmailles; i++){
2843 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
2844 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
2846 int Ii = i*_nVar +1+j;
2847 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
2849 for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
2852 _Vitesse.setTime(_time,_nbTimeStep);
2853 if (_nbTimeStep ==0 || _restartWithNewFileName){
2854 _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
2855 _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
2856 _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
2861 _Vitesse.writeVTK(prim+"_Velocity");
2864 _Vitesse.writeMED(prim+"_Velocity");
2867 _Vitesse.writeCSV(prim+"_Velocity");
2875 _Vitesse.writeVTK(prim+"_Velocity",false);
2878 _Vitesse.writeMED(prim+"_Velocity",false);
2881 _Vitesse.writeCSV(prim+"_Velocity");
2889 double p,T,rho, h, vx,vy,vz,v2;
2891 for (long i = 0; i < _Nmailles; i++){
2893 VecGetValues(_conservativeVars,1,&Ii,&rho);
2895 VecGetValues(_primitiveVars,1,&Ii,&p);
2896 Ii = i*_nVar +_nVar-1;
2897 VecGetValues(_primitiveVars,1,&Ii,&T);
2899 VecGetValues(_primitiveVars,1,&Ii,&vx);
2903 VecGetValues(_primitiveVars,1,&Ii,&vy);
2906 VecGetValues(_primitiveVars,1,&Ii,&vz);
2910 h = _fluides[0]->getEnthalpy(T,rho);
2928 _MachNumber(i)=sqrt(v2)/_fluides[0]->vitesseSonEnthalpie(h);
2930 _Enthalpy.setTime(_time,_nbTimeStep);
2931 _Density.setTime(_time,_nbTimeStep);
2932 _Pressure.setTime(_time,_nbTimeStep);
2933 _Temperature.setTime(_time,_nbTimeStep);
2934 _MachNumber.setTime(_time,_nbTimeStep);
2935 _VitesseX.setTime(_time,_nbTimeStep);
2938 _VitesseY.setTime(_time,_nbTimeStep);
2940 _VitesseZ.setTime(_time,_nbTimeStep);
2942 if (_nbTimeStep ==0 || _restartWithNewFileName){
2946 _Enthalpy.writeVTK(allFields+"_Enthalpy");
2947 _Density.writeVTK(allFields+"_Density");
2948 _Pressure.writeVTK(allFields+"_Pressure");
2949 _Temperature.writeVTK(allFields+"_Temperature");
2950 _MachNumber.writeVTK(allFields+"_MachNumber");
2951 _VitesseX.writeVTK(allFields+"_VelocityX");
2954 _VitesseY.writeVTK(allFields+"_VelocityY");
2956 _VitesseZ.writeVTK(allFields+"_VelocityZ");
2960 _Enthalpy.writeMED(allFields+"_Enthalpy");
2961 _Density.writeMED(allFields+"_Density");
2962 _Pressure.writeMED(allFields+"_Pressure");
2963 _Temperature.writeMED(allFields+"_Temperature");
2964 _MachNumber.writeMED(allFields+"_MachNumber");
2965 _VitesseX.writeMED(allFields+"_VelocityX");
2968 _VitesseY.writeMED(allFields+"_VelocityY");
2970 _VitesseZ.writeMED(allFields+"_VelocityZ");
2974 _Enthalpy.writeCSV(allFields+"_Enthalpy");
2975 _Density.writeCSV(allFields+"_Density");
2976 _Pressure.writeCSV(allFields+"_Pressure");
2977 _Temperature.writeCSV(allFields+"_Temperature");
2978 _MachNumber.writeCSV(allFields+"_MachNumber");
2979 _VitesseX.writeCSV(allFields+"_VelocityX");
2982 _VitesseY.writeCSV(allFields+"_VelocityY");
2984 _VitesseZ.writeCSV(allFields+"_VelocityZ");
2993 _Enthalpy.writeVTK(allFields+"_Enthalpy",false);
2994 _Density.writeVTK(allFields+"_Density",false);
2995 _Pressure.writeVTK(allFields+"_Pressure",false);
2996 _Temperature.writeVTK(allFields+"_Temperature",false);
2997 _MachNumber.writeVTK(allFields+"_MachNumber",false);
2998 _VitesseX.writeVTK(allFields+"_VelocityX",false);
3001 _VitesseY.writeVTK(allFields+"_VelocityY",false);
3003 _VitesseZ.writeVTK(allFields+"_VelocityZ",false);
3007 _Enthalpy.writeMED(allFields+"_Enthalpy",false);
3008 _Density.writeMED(allFields+"_Density",false);
3009 _Pressure.writeMED(allFields+"_Pressure",false);
3010 _Temperature.writeMED(allFields+"_Temperature",false);
3011 _MachNumber.writeMED(allFields+"_MachNumber",false);
3012 _VitesseX.writeMED(allFields+"_VelocityX",false);
3015 _VitesseY.writeMED(allFields+"_VelocityY",false);
3017 _VitesseZ.writeMED(allFields+"_VelocityZ",false);
3021 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3022 _Density.writeCSV(allFields+"_Density");
3023 _Pressure.writeCSV(allFields+"_Pressure");
3024 _Temperature.writeCSV(allFields+"_Temperature");
3025 _MachNumber.writeCSV(allFields+"_MachNumber");
3026 _VitesseX.writeCSV(allFields+"_VelocityX");
3029 _VitesseY.writeCSV(allFields+"_VelocityY");
3031 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3056 if(_saveConservativeField){
3071 if(_saveVelocity || _saveAllFields){
3075 _Vitesse.writeVTK(prim+"_Velocity");
3078 _Vitesse.writeMED(prim+"_Velocity");
3081 _Vitesse.writeCSV(prim+"_Velocity");
3087 if (_restartWithNewFileName)
3088 _restartWithNewFileName=false;
3091 Field& SinglePhase::getPressureField()
3095 _Pressure=Field("Pressure",CELLS,_mesh,1);
3097 for (long i = 0; i < _Nmailles; i++){
3099 VecGetValues(_primitiveVars,1,&Ii,&_Pressure(i));
3101 _Pressure.setTime(_time,_nbTimeStep);
3106 Field& SinglePhase::getTemperatureField()
3110 _Temperature=Field("Temperature",CELLS,_mesh,1);
3112 for (long i = 0; i < _Nmailles; i++){
3113 Ii = i*_nVar +_nVar-1;
3114 VecGetValues(_primitiveVars,1,&Ii,&_Temperature(i));
3116 _Temperature.setTime(_time,_nbTimeStep);
3118 return _Temperature;
3121 Field& SinglePhase::getVelocityField()
3123 if(!_saveAllFields )
3125 _Vitesse=Field("Vitesse",CELLS,_mesh,3);
3127 for (long i = 0; i < _Nmailles; i++)
3129 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
3131 int Ii = i*_nVar +1+j;
3132 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
3134 for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
3137 _Vitesse.setTime(_time,_nbTimeStep);
3138 _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
3139 _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
3140 _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
3146 Field& SinglePhase::getVelocityXField()
3148 if(!_saveAllFields )
3150 _VitesseX=Field("Velocity X",CELLS,_mesh,1);
3152 for (long i = 0; i < _Nmailles; i++)
3154 int Ii = i*_nVar +1;
3155 VecGetValues(_primitiveVars,1,&Ii,&_VitesseX(i));
3157 _VitesseX.setTime(_time,_nbTimeStep);
3158 _VitesseX.setInfoOnComponent(0,"Velocity_x_(m/s)");
3164 Field& SinglePhase::getDensityField()
3166 if(!_saveAllFields )
3168 _Density=Field("Density",CELLS,_mesh,1);
3170 for (long i = 0; i < _Nmailles; i++){
3172 VecGetValues(_conservativeVars,1,&Ii,&_Density(i));
3174 _Density.setTime(_time,_nbTimeStep);
3179 Field& SinglePhase::getMomentumField()//not yet managed by parameter _saveAllFields
3181 _Momentum=Field("Momentum",CELLS,_mesh,_Ndim);
3183 for (long i = 0; i < _Nmailles; i++)
3184 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de qdm
3186 int Ii = i*_nVar +1+j;
3187 VecGetValues(_conservativeVars,1,&Ii,&_Momentum(i,j));
3189 _Momentum.setTime(_time,_nbTimeStep);
3194 Field& SinglePhase::getTotalEnergyField()//not yet managed by parameter _saveAllFields
3196 _TotalEnergy=Field("TotalEnergy",CELLS,_mesh,1);
3198 for (long i = 0; i < _Nmailles; i++){
3199 Ii = i*_nVar +_nVar-1;
3200 VecGetValues(_conservativeVars,1,&Ii,&_TotalEnergy(i));
3202 _TotalEnergy.setTime(_time,_nbTimeStep);
3204 return _TotalEnergy;
3207 Field& SinglePhase::getEnthalpyField()
3209 if(!_saveAllFields )
3211 _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
3214 for (long i = 0; i < _Nmailles; i++){
3216 VecGetValues(_primitiveVars,1,&Ii,&p);
3217 Ii = i*_nVar +_nVar-1;
3218 VecGetValues(_primitiveVars,1,&Ii,&T);
3220 rho=_fluides[0]->getDensity(p,T);
3221 _Enthalpy(i)=_fluides[0]->getEnthalpy(T,rho);
3223 _Enthalpy.setTime(_time,_nbTimeStep);
3229 vector<string> SinglePhase::getOutputFieldsNames()
3231 vector<string> result(8);
3233 result[0]="Pressure";
3234 result[1]="Velocity";
3235 result[2]="Temperature";
3236 result[3]="Density";
3237 result[4]="Momentum";
3238 result[5]="TotalEnergy";
3239 result[6]="Enthalpy";
3240 result[7]="VelocityX";
3245 Field& SinglePhase::getOutputField(const string& nameField )
3247 if(nameField=="pressure" || nameField=="Pressure" || nameField=="PRESSURE" || nameField=="PRESSION" || nameField=="Pression" || nameField=="pression" )
3248 return getPressureField();
3249 else if(nameField=="velocity" || nameField=="Velocity" || nameField=="VELOCITY" || nameField=="Vitesse" || nameField=="VITESSE" || nameField=="vitesse" )
3250 return getVelocityField();
3251 else if(nameField=="velocityX" || nameField=="VelocityX" || nameField=="VELOCITYX" || nameField=="VitesseX" || nameField=="VITESSEX" || nameField=="vitesseX" )
3252 return getVelocityXField();
3253 else if(nameField=="temperature" || nameField=="Temperature" || nameField=="TEMPERATURE" || nameField=="temperature" )
3254 return getTemperatureField();
3255 else if(nameField=="density" || nameField=="Density" || nameField=="DENSITY" || nameField=="Densite" || nameField=="DENSITE" || nameField=="densite" )
3256 return getDensityField();
3257 else if(nameField=="momentum" || nameField=="Momentum" || nameField=="MOMENTUM" || nameField=="Qdm" || nameField=="QDM" || nameField=="qdm" )
3258 return getMomentumField();
3259 else if(nameField=="enthalpy" || nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPIE" || nameField=="enthalpie" )
3260 return getEnthalpyField();
3261 else if(nameField=="totalenergy" || nameField=="TotalEnergy" || nameField=="TOTALENERGY" || nameField=="ENERGIETOTALE" || nameField=="EnergieTotale" || nameField=="energietotale" )
3262 return getTotalEnergyField();
3265 cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
3266 throw CdmathException("SinglePhase::getOutputField error : Unknown Field name");