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 _VitesseX=Field("Velocity x",CELLS,_mesh,1);
75 _VitesseY=Field("Velocity y",CELLS,_mesh,1);
77 _VitesseZ=Field("Velocity z",CELLS,_mesh,1);
81 if(_entropicCorrection)
82 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
84 ProblemFluid::initialize();
87 bool SinglePhase::iterateTimeStep(bool &converged)
89 if(_timeScheme == Explicit || !_usePrimitiveVarsInNewton)
90 return ProblemFluid::iterateTimeStep(converged);
95 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
97 computeTimeStep(stop);
99 if(stop){//Le compute time step ne s'est pas bien passé
100 cout<<"ComputeTimeStep failed"<<endl;
105 computeNewtonVariation();
107 //converged=convergence des iterations
108 if(_timeScheme == Explicit)
110 else{//Implicit scheme
112 KSPGetIterationNumber(_ksp, &_PetscIts);
113 if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
114 _MaxIterLinearSolver = _PetscIts;
115 if(_PetscIts>=_maxPetscIts)//solving the linear system failed
117 cout<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
118 *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Itérations maximales "<<_maxPetscIts<<" atteintes"<<endl;
122 else{//solving the linear system succeeded
123 //Calcul de la variation relative Uk+1-Uk
127 for(int j=1; j<=_Nmailles; j++)
129 for(int k=0; k<_nVar; k++)
132 VecGetValues(_newtonVariation, 1, &I, &dx);
133 VecGetValues(_primitiveVars, 1, &I, &x);
134 if (fabs(x)*fabs(x)< _precision)
136 if(_erreur_rel < fabs(dx))
137 _erreur_rel = fabs(dx);
139 else if(_erreur_rel < fabs(dx/x))
140 _erreur_rel = fabs(dx/x);
144 converged = _erreur_rel <= _precision_Newton;
147 double relaxation=1;//Vk+1=Vk+relaxation*deltaV
149 VecAXPY(_primitiveVars, relaxation, _newtonVariation);
151 //mise a jour du champ primitif
152 updateConservatives();
154 if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
156 cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
157 *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
163 cout<<"Vecteur Vkp1-Vk "<<endl;
164 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
165 cout << "Nouvel etat courant Vk de l'iteration Newton: " << endl;
166 VecView(_primitiveVars, PETSC_VIEWER_STDOUT_SELF);
169 if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
170 if(_minm1<-_precision || _minm2<-_precision)
172 cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
173 *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
177 cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
178 *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
190 void SinglePhase::computeNewtonVariation()
192 if(!_usePrimitiveVarsInNewton)
193 ProblemFluid::computeNewtonVariation();
198 cout<<"Vecteur courant Vk "<<endl;
199 VecView(_primitiveVars,PETSC_VIEWER_STDOUT_SELF);
201 cout << "Matrice du système linéaire avant contribution delta t" << endl;
202 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
204 cout << "Second membre du système linéaire avant contribution delta t" << endl;
205 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
208 if(_timeScheme == Explicit)
210 VecCopy(_b,_newtonVariation);
211 VecScale(_newtonVariation, _dt);
212 if(_verbose && _nbTimeStep%_freqSave ==0)
214 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
215 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
221 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
223 VecAXPY(_b, 1/_dt, _old);
224 VecAXPY(_b, -1/_dt, _conservativeVars);
226 for(int imaille = 0; imaille<_Nmailles; imaille++){
227 _idm[0] = _nVar*imaille;
228 for(int k=1; k<_nVar; k++)
229 _idm[k] = _idm[k-1] + 1;
230 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
231 primToConsJacobianMatrix(_Vi);
232 for(int k=0; k<_nVar*_nVar; k++)
233 _primToConsJacoMat[k]*=1/_dt;
234 MatSetValuesBlocked(_A, 1, &imaille, 1, &imaille, _primToConsJacoMat, ADD_VALUES);
236 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
238 #if PETSC_VERSION_GREATER_3_5
239 KSPSetOperators(_ksp, _A, _A);
241 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
246 cout << "Matrice du système linéaire" << endl;
247 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
249 cout << "Second membre du système linéaire" << endl;
250 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
255 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
258 KSPSolve(_ksp, _b, _newtonVariation);
262 computeScaling(_maxvp);
264 VecAssemblyBegin(_vecScaling);
265 VecAssemblyBegin(_invVecScaling);
266 for(int imaille = 0; imaille<_Nmailles; imaille++)
269 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
270 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
272 VecAssemblyEnd(_vecScaling);
273 VecAssemblyEnd(_invVecScaling);
277 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
278 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
280 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
281 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
284 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
287 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
288 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
291 VecCopy(_b,_bScaling);
292 VecPointwiseMult(_b,_vecScaling,_bScaling);
295 cout << "Produit du second membre par le preconditionneur bloc diagonal a gauche" << endl;
296 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
300 KSPSolve(_ksp,_b, _bScaling);
301 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
305 cout << "solution du systeme lineaire local:" << endl;
306 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
312 void SinglePhase::convectionState( const long &i, const long &j, const bool &IsBord){
315 for(int k=1; k<_nVar; k++)
316 _idm[k] = _idm[k-1] + 1;
317 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
320 for(int k=1; k<_nVar; k++)
321 _idm[k] = _idm[k-1] + 1;
323 VecGetValues(_Uext, _nVar, _idm, _Uj);
325 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
326 if(_verbose && _nbTimeStep%_freqSave ==0)
328 cout<<"Convection Left state cell " << i<< ": "<<endl;
329 for(int k =0; k<_nVar; k++)
331 cout<<"Convection Right state cell " << j<< ": "<<endl;
332 for(int k =0; k<_nVar; k++)
335 if(_Ui[0]<0||_Uj[0]<0)
337 cout<<"!!!!!!!!!!!!!!!!!!!!!!!!densite negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
338 *_runLogFile<<"!!!!!!!!!!!!!!!!!!!!!!!!densite negative, arret de calcul!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl;
339 _runLogFile->close();
340 throw CdmathException("densite negative, arret de calcul");
342 PetscScalar ri, rj, xi, xj, pi, pj;
344 ri = sqrt(_Ui[0]);//racine carre de phi_i rho_i
346 _Uroe[0] = ri*rj; //moyenne geometrique des densites
347 if(_verbose && _nbTimeStep%_freqSave ==0)
348 cout << "Densité moyenne Roe gauche " << i << ": " << ri*ri << ", droite " << j << ": " << rj*rj << "->" << _Uroe[0] << endl;
349 for(int k=0;k<_Ndim;k++){
352 _Uroe[1+k] = (xi/ri + xj/rj)/(ri + rj);
353 //"moyenne" des vitesses
354 if(_verbose && _nbTimeStep%_freqSave ==0)
355 cout << "Vitesse de Roe composante "<< k<<" gauche " << i << ": " << xi/(ri*ri) << ", droite " << j << ": " << xj/(rj*rj) << "->" << _Uroe[k+1] << endl;
357 // H = (rho E + p)/rho
358 xi = _Ui[_nVar-1];//phi rho E
360 Ii = i*_nVar; // correct Kieu
361 VecGetValues(_primitiveVars, 1, &Ii, &pi);// _primitiveVars pour p
365 for(int k=1;k<=_Ndim;k++)
366 q_2 += _Uj[k]*_Uj[k];
367 q_2 /= _Uj[0]; //phi rho u²
368 pj = _fluides[0]->getPressure((_Uj[(_Ndim+2)-1]-q_2/2)/_porosityj,_Uj[0]/_porosityj);
372 Ii = j*_nVar; // correct Kieu
373 VecGetValues(_primitiveVars, 1, &Ii, &pj);
375 xi = (xi + pi)/(ri*ri);
376 xj = (xj + pj)/(rj*rj);
377 _Uroe[_nVar-1] = (ri*xi + rj*xj)/(ri + rj);
378 //on se donne l enthalpie ici
379 if(_verbose && _nbTimeStep%_freqSave ==0)
380 cout << "Enthalpie totale de Roe H gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[_nVar-1] << endl;
382 if(_verbose && _nbTimeStep%_freqSave ==0)
384 cout<<"Convection interfacial state"<<endl;
385 for(int k=0;k<_nVar;k++)
386 cout<< _Uroe[k]<<" , "<<endl;
390 void SinglePhase::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
391 //sortie: matrices et etat de diffusion (rho, q, T)
393 for(int k=1; k<_nVar; k++)
394 _idm[k] = _idm[k-1] + 1;
396 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
398 for(int k=1; k<_nVar; k++)
399 _idm[k] = _idm[k-1] + 1;
402 VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
404 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
406 if(_verbose && _nbTimeStep%_freqSave ==0)
408 cout << "SinglePhase::diffusionStateAndMatrices cellule gauche" << i << endl;
410 for(int q=0; q<_nVar; q++)
411 cout << _Ui[q] << "\t";
413 cout << "SinglePhase::diffusionStateAndMatrices cellule droite" << j << endl;
415 for(int q=0; q<_nVar; q++)
416 cout << _Uj[q] << "\t";
420 for(int k=0; k<_nVar; k++)
421 _Udiff[k] = (_Ui[k]/_porosityi+_Uj[k]/_porosityj)/2;
423 if(_verbose && _nbTimeStep%_freqSave ==0)
425 cout << "SinglePhase::diffusionStateAndMatrices conservative diffusion state" << endl;
427 for(int q=0; q<_nVar; q++)
428 cout << _Udiff[q] << "\t";
430 cout << "porosite gauche= "<<_porosityi<< ", porosite droite= "<<_porosityj<<endl;
432 consToPrim(_Udiff,_phi,1);
433 _Udiff[_nVar-1]=_phi[_nVar-1];
435 if(_timeScheme==Implicit)
438 for (int i = 0; i<_Ndim;i++)
439 q_2+=_Udiff[i+1]*_Udiff[i+1];
440 double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
441 double lambda = _fluides[0]->getConductivity(_Udiff[_nVar-1]);
442 double Cv= _fluides[0]->constante("Cv");
443 for(int i=0; i<_nVar*_nVar;i++)
445 for(int i=1;i<(_nVar-1);i++)
447 _Diffusion[i*_nVar] = mu*_Udiff[i]/(_Udiff[0]*_Udiff[0]);
448 _Diffusion[i*_nVar+i] = -mu/_Udiff[0];
450 int i = (_nVar-1)*_nVar;
451 _Diffusion[i]=lambda*(_Udiff[_nVar-1]/_Udiff[0]-q_2/(2*Cv*_Udiff[0]*_Udiff[0]*_Udiff[0]));
452 for(int k=1;k<(_nVar-1);k++)
454 _Diffusion[i+k]= lambda*_Udiff[k]/(_Udiff[0]*_Udiff[0]*Cv);
456 _Diffusion[_nVar*_nVar-1]=-lambda/(_Udiff[0]*Cv);
460 void SinglePhase::setBoundaryState(string nameOfGroup, const int &j,double *normale){
462 for(int k=1; k<_nVar; k++)
463 _idm[k] = _idm[k-1] + 1;
465 double porosityj=_porosityField(j);
467 if(_verbose && _nbTimeStep%_freqSave ==0)
469 cout << "setBoundaryState for group "<< nameOfGroup << ", inner cell j= "<<j<< " face unit normal vector "<<endl;
470 for(int k=0; k<_Ndim; k++){
471 cout<<normale[k]<<", ";
476 if (_limitField[nameOfGroup].bcType==Wall){
477 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne conservatif
478 double q_n=0;//q_n=quantité de mouvement normale à la face frontière;
479 for(int k=0; k<_Ndim; k++)
480 q_n+=_externalStates[(k+1)]*normale[k];
482 //Pour la convection, inversion du sens de la vitesse normale
483 for(int k=0; k<_Ndim; k++)
484 _externalStates[(k+1)]-= 2*q_n*normale[k];
487 for(int k=1; k<_nVar; k++)
488 _idm[k] = _idm[k-1] + 1;
490 VecAssemblyBegin(_Uext);
491 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
492 VecAssemblyEnd(_Uext);
494 //Pour la diffusion, paroi à vitesse et temperature imposees
496 for(int k=1; k<_nVar; k++)
497 _idm[k] = _idm[k-1] + 1;
498 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);//L'état fantome contient à présent les variables primitives internes
499 double pression=_externalStates[0];
500 double T=_limitField[nameOfGroup].T;
501 double rho=_fluides[0]->getDensity(pression,T);
503 _externalStates[0]=porosityj*rho;
504 _externalStates[1]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
506 v2 +=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
509 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
510 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
513 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
514 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
517 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
519 for(int k=1; k<_nVar; k++)
520 _idm[k] = _idm[k-1] + 1;
521 VecAssemblyBegin(_Uextdiff);
522 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
523 VecAssemblyEnd(_Uextdiff);
525 else if (_limitField[nameOfGroup].bcType==Neumann){
526 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On prend l'état fantôme égal à l'état interne (conditions limites de Neumann)
529 for(int k=1; k<_nVar; k++)
530 _idm[k] = _idm[k-1] + 1;
532 VecAssemblyBegin(_Uext);
533 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
534 VecAssemblyEnd(_Uext);
536 VecAssemblyBegin(_Uextdiff);
537 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
538 VecAssemblyEnd(_Uextdiff);
540 else if (_limitField[nameOfGroup].bcType==Inlet){
541 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne (conditions limites de Neumann)
542 double q_int_n=0;//q_int_n=composante normale de la quantité de mouvement à la face frontière;
543 for(int k=0; k<_Ndim; k++)
544 q_int_n+=_externalStates[(k+1)]*normale[k];//On calcule la vitesse normale sortante
546 double q_ext_n=_limitField[nameOfGroup].v_x[0]*normale[0];
549 q_ext_n+=_limitField[nameOfGroup].v_y[0]*normale[1];
551 q_ext_n+=_limitField[nameOfGroup].v_z[0]*normale[2];
554 if(q_int_n+q_ext_n<=0){//Interfacial velocity goes inward
555 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);//On met à jour l'état fantome avec les variables primitives internes
556 double pression=_externalStates[0];
557 double T=_limitField[nameOfGroup].T;
558 double rho=_fluides[0]->getDensity(pression,T);
560 _externalStates[0]=porosityj*rho;//Composante fantome de masse
561 _externalStates[1]=_externalStates[0]*(_limitField[nameOfGroup].v_x[0]);//Composante fantome de qdm x
563 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
566 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
567 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];//Composante fantome de qdm y
570 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];//Composante fantome de qdm z
571 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
574 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);//Composante fantome de l'nrj
576 else if(_nbTimeStep%_freqSave ==0)
577 cout<< "Warning : fluid possibly going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
580 for(int k=1; k<_nVar; k++)
581 _idm[k] = _idm[k-1] + 1;
582 VecAssemblyBegin(_Uext);
583 VecAssemblyBegin(_Uextdiff);
584 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
585 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
586 VecAssemblyEnd(_Uext);
587 VecAssemblyEnd(_Uextdiff);
589 else if (_limitField[nameOfGroup].bcType==InletRotationVelocity){
590 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
591 double u_int_n=0;//u_int_n=composante normale de la vitesse intérieure à la face frontière;
592 for(int k=0; k<_Ndim; k++)
593 u_int_n+=_externalStates[(k+1)]*normale[k];
600 omega[0]=_limitField[nameOfGroup].v_x[0];
601 omega[1]=_limitField[nameOfGroup].v_y[0];
604 Normale[0]=normale[0];
605 Normale[1]=normale[1];
609 omega[2]=_limitField[nameOfGroup].v_z[0];
610 Normale[2]=normale[2];
613 Vector tangent_vel=omega%Normale;
614 u_ext_n=-0.01*tangent_vel.norm();
615 //Changing external state velocity
616 for(int k=0; k<_Ndim; k++)
617 _externalStates[(k+1)]=u_ext_n*normale[k] + tangent_vel[k];
620 if(u_ext_n + u_int_n <= 0){
621 double pression=_externalStates[0];
622 double T=_limitField[nameOfGroup].T;
623 double rho=_fluides[0]->getDensity(pression,T);
626 v2 +=_externalStates[1]*_externalStates[1];//v_x*v_x
627 _externalStates[0]=porosityj*rho;
628 _externalStates[1]*=_externalStates[0];
631 v2 +=_externalStates[2]*_externalStates[2];//+v_y*v_y
632 _externalStates[2]*=_externalStates[0];
635 v2 +=_externalStates[3]*_externalStates[3];//+v_z*v_z
636 _externalStates[3]*=_externalStates[0];
639 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
641 else if(_nbTimeStep%_freqSave ==0)
644 * cout<< "Warning : fluid going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
646 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On définit l'état fantôme avec l'état interne
647 if(_nbTimeStep%_freqSave ==0)
648 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Wall boundary condition."<<endl;
650 //Changing external state momentum
651 for(int k=0; k<_Ndim; k++)
652 _externalStates[(k+1)]-=2*_externalStates[0]*u_int_n*normale[k];
656 for(int k=1; k<_nVar; k++)
657 _idm[k] = _idm[k-1] + 1;
658 VecAssemblyBegin(_Uext);
659 VecAssemblyBegin(_Uextdiff);
660 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
661 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
662 VecAssemblyEnd(_Uext);
663 VecAssemblyEnd(_Uextdiff);
665 else if (_limitField[nameOfGroup].bcType==InletPressure){
666 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
668 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
669 Cell Cj=_mesh.getCell(j);
670 double hydroPress=Cj.x()*_GravityField3d[0];
672 hydroPress+=Cj.y()*_GravityField3d[1];
674 hydroPress+=Cj.z()*_GravityField3d[2];
676 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
678 //Building the primitive external state
679 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
680 double u_n=0;//u_n=vitesse normale à la face frontière;
681 for(int k=0; k<_Ndim; k++)
682 u_n+=_externalStates[(k+1)]*normale[k];
686 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress,_limitField[nameOfGroup].T);
688 //Contribution from the tangential velocity
692 omega[0]=_limitField[nameOfGroup].v_x[0];
693 omega[1]=_limitField[nameOfGroup].v_y[0];
696 Normale[0]=normale[0];
697 Normale[1]=normale[1];
701 omega[2]=_limitField[nameOfGroup].v_z[0];
702 Normale[2]=normale[2];
705 Vector tangent_vel=omega%Normale;
707 //Changing external state velocity
708 for(int k=0; k<_Ndim; k++)
709 _externalStates[(k+1)]=u_n*normale[k] + abs(u_n)*tangent_vel[k];
714 if(_nbTimeStep%_freqSave ==0)
715 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;
716 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
718 if(_nbTimeStep%_freqSave ==0)
719 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Wall boundary condition."<<endl;
720 _externalStates[0]=porosityj*_fluides[0]->getDensity(_externalStates[0]+hydroPress, _externalStates[_nVar-1]);
721 //Changing external state velocity
722 for(int k=0; k<_Ndim; k++)
723 _externalStates[(k+1)]-=2*u_n*normale[k];
727 for(int k=0; k<_Ndim; k++)
729 v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
730 _externalStates[(k+1)]*=_externalStates[0] ;//qdm component
732 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);//nrj component
736 for(int k=1; k<_nVar; k++)
737 _idm[k] = _idm[k-1] + 1;
738 VecAssemblyBegin(_Uext);
739 VecAssemblyBegin(_Uextdiff);
740 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
741 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
742 VecAssemblyEnd(_Uext);
743 VecAssemblyEnd(_Uextdiff);
745 else if (_limitField[nameOfGroup].bcType==Outlet){
746 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne conservatif
747 double q_n=0;//q_n=quantité de mouvement normale à la face frontière;
748 for(int k=0; k<_Ndim; k++)
749 q_n+=_externalStates[(k+1)]*normale[k];
751 if(q_n < -_precision && _nbTimeStep%_freqSave ==0)
753 cout<< "Warning : fluid going in through outlet boundary "<<nameOfGroup<<" with flow rate "<< q_n<<endl;
754 cout<< "Applying Neumann boundary condition for velocity and temperature"<<endl;
756 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
757 Cell Cj=_mesh.getCell(j);
758 double hydroPress=Cj.x()*_GravityField3d[0];
760 hydroPress+=Cj.y()*_GravityField3d[1];
762 hydroPress+=Cj.z()*_GravityField3d[2];
764 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
766 if(_verbose && _nbTimeStep%_freqSave ==0)
768 cout<<"Cond lim outlet densite= "<<_externalStates[0]<<" gravite= "<<_GravityField3d[0]<<" Cj.x()= "<<Cj.x()<<endl;
769 cout<<"Cond lim outlet pression ref= "<<_limitField[nameOfGroup].p<<" pression hydro= "<<hydroPress<<" total= "<<_limitField[nameOfGroup].p+hydroPress<<endl;
771 //Building the external state
772 _idm[0] = _nVar*j;// Kieu
773 for(int k=1; k<_nVar; k++)
774 _idm[k] = _idm[k-1] + 1;
775 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
777 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
779 for(int k=0; k<_Ndim; k++)
781 v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
782 _externalStates[(k+1)]*=_externalStates[0] ;
784 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);
786 for(int k=1; k<_nVar; k++)
787 _idm[k] = _idm[k-1] + 1;
788 VecAssemblyBegin(_Uext);
789 VecAssemblyBegin(_Uextdiff);
790 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
791 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
792 VecAssemblyEnd(_Uext);
793 VecAssemblyEnd(_Uextdiff);
795 cout<<"Boundary condition not set for boundary named "<<nameOfGroup<< " _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
796 cout<<"Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
797 *_runLogFile<<"Boundary condition not set for boundary named. Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
798 _runLogFile->close();
799 throw CdmathException("Unknown boundary condition");
803 void SinglePhase::convectionMatrices()
805 //entree: URoe = rho, u, H
806 //sortie: matrices Roe+ et Roe-
808 if(_verbose && _nbTimeStep%_freqSave ==0)
809 cout<<"SinglePhase::convectionMatrices()"<<endl;
811 double u_n=0, u_2=0;//vitesse normale et carré du module
813 for(int i=0;i<_Ndim;i++)
815 u_2 += _Uroe[1+i]*_Uroe[1+i];
816 u_n += _Uroe[1+i]*_vec_normal[i];
819 vector<complex<double> > vp_dist(3,0);
821 if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
823 staggeredVFFCMatricesConservativeVariables(u_n);//Computation of classical upwinding matrices
824 if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//For use in implicit matrix
825 staggeredVFFCMatricesPrimitiveVariables(u_n);
829 Vector vitesse(_Ndim);
830 for(int idim=0;idim<_Ndim;idim++)
831 vitesse[idim]=_Uroe[1+idim];
834 /***********Calcul des valeurs propres ********/
836 c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
837 k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
838 K = u_2*k/2; //g-1/2 *|u|²
840 vp_dist[0]=u_n-c;vp_dist[1]=u_n;vp_dist[2]=u_n+c;
842 _maxvploc=fabs(u_n)+c;
846 if(_verbose && _nbTimeStep%_freqSave ==0)
847 cout<<"SinglePhase::convectionMatrices Eigenvalues "<<u_n-c<<" , "<<u_n<<" , "<<u_n+c<<endl;
849 RoeMatrixConservativeVariables( u_n, H,vitesse,k,K);
851 /******** Construction des matrices de decentrement ********/
852 if( _spaceScheme ==centered){
853 if(_entropicCorrection)
855 *_runLogFile<<"SinglePhase::convectionMatrices: entropy scheme not available for centered scheme"<<endl;
856 _runLogFile->close();
857 throw CdmathException("SinglePhase::convectionMatrices: entropy scheme not available for centered scheme");
860 for(int i=0; i<_nVar*_nVar;i++)
863 else if(_spaceScheme == upwind || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach){
864 if(_entropicCorrection)
865 entropicShift(_vec_normal);
867 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
869 vector< complex< double > > y (3,0);
871 for( int i=0 ; i<3 ; i++)
872 y[i] = Poly.abs_generalise(vp_dist[i])+_entropicShift[i];
873 Poly.abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
875 if( _spaceScheme ==pressureCorrection)
876 for( int i=0 ; i<_Ndim ; i++)
877 for( int j=0 ; j<_Ndim ; j++)
878 _absAroe[(1+i)*_nVar+1+j]-=(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
879 else if( _spaceScheme ==lowMach){
880 double M=sqrt(u_2)/c;
881 for( int i=0 ; i<_Ndim ; i++)
882 for( int j=0 ; j<_Ndim ; j++)
883 _absAroe[(1+i)*_nVar+1+j]-=(1-M)*(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
886 else if( _spaceScheme ==staggered ){
887 if(_entropicCorrection)//To do: study entropic correction for staggered
889 *_runLogFile<<"SinglePhase::convectionMatrices: entropy scheme not available for staggered scheme"<<endl;
890 _runLogFile->close();
891 throw CdmathException("SinglePhase::convectionMatrices: entropy scheme not available for staggered scheme");
894 staggeredRoeUpwindingMatrixConservativeVariables( u_n, H, vitesse, k, K);
898 *_runLogFile<<"SinglePhase::convectionMatrices: scheme not treated"<<endl;
899 _runLogFile->close();
900 throw CdmathException("SinglePhase::convectionMatrices: scheme not treated");
903 for(int i=0; i<_nVar*_nVar;i++)
905 _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
906 _AroePlus[i] = (_Aroe[i]+_absAroe[i])/2;
908 if(_timeScheme==Implicit)
910 if(_usePrimitiveVarsInNewton)//Implicitation using primitive variables
912 _Vij[0]=_fluides[0]->getPressureFromEnthalpy(_Uroe[_nVar-1]-u_2/2, _Uroe[0]);//pressure
913 _Vij[_nVar-1]=_fluides[0]->getTemperatureFromPressure( _Vij[0], _Uroe[0]);//Temperature
914 for(int idim=0;idim<_Ndim; idim++)
915 _Vij[1+idim]=_Uroe[1+idim];
916 primToConsJacobianMatrix(_Vij);
918 Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
919 Poly.matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
922 for(int i=0; i<_nVar*_nVar;i++)
924 _AroeMinusImplicit[i] = _AroeMinus[i];
925 _AroePlusImplicit[i] = _AroePlus[i];
928 if(_verbose && _nbTimeStep%_freqSave ==0)
930 displayMatrix(_Aroe, _nVar,"Matrice de Roe");
931 displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
932 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
933 displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
937 if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
939 displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
940 displayMatrix(_AroePlusImplicit, _nVar,"Matrice _AroePlusImplicit");
943 /*********Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source*****/
944 if(_entropicCorrection)
946 InvMatriceRoe( vp_dist);
948 Poly.matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
950 else if (_spaceScheme==upwind || (_spaceScheme ==pressureCorrection ) || (_spaceScheme ==lowMach ))//upwind sans entropic
951 SigneMatriceRoe( vp_dist);
952 else if (_spaceScheme==centered)//centre sans entropic
953 for(int i=0; i<_nVar*_nVar;i++)
955 else if( _spaceScheme ==staggered )//à tester
964 for(int i=0; i<_nVar*_nVar;i++)
966 _signAroe[0] = signu;
967 for(int i=1; i<_nVar-1;i++)
968 _signAroe[i*_nVar+i] = -signu;
969 _signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
973 *_runLogFile<<"SinglePhase::convectionMatrices: well balanced option not treated"<<endl;
974 _runLogFile->close();
975 throw CdmathException("SinglePhase::convectionMatrices: well balanced option not treated");
978 void SinglePhase::computeScaling(double maxvp)
982 for(int q=1; q<_nVar-1; q++)
984 _blockDiag[q]=1./maxvp;//
985 _invBlockDiag[q]= maxvp;//1.;//
987 _blockDiag[_nVar - 1]=(_fluides[0]->constante("gamma")-1)/(maxvp*maxvp);//1
988 _invBlockDiag[_nVar - 1]= 1./_blockDiag[_nVar - 1] ;// 1.;//
991 void SinglePhase::addDiffusionToSecondMember
997 double lambda=_fluides[0]->getConductivity(_Udiff[_nVar-1]);
998 double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
1000 if(lambda==0 && mu ==0 && _heatTransfertCoeff==0)
1003 //extraction des valeurs
1004 _idm[0] = _nVar*i; // Kieu
1005 for(int k=1; k<_nVar; k++)
1006 _idm[k] = _idm[k-1] + 1;
1008 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1009 if (_verbose && _nbTimeStep%_freqSave ==0)
1011 cout << "Calcul diffusion: variables primitives maille " << i<<endl;
1012 for(int q=0; q<_nVar; q++)
1013 cout << _Vi[q] << endl;
1018 for(int k=0; k<_nVar; k++)
1019 _idn[k] = _nVar*j + k;
1021 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1025 lambda=max(lambda,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
1026 for(int k=0; k<_nVar; k++)
1029 VecGetValues(_Uextdiff, _nVar, _idn, _phi);
1030 consToPrim(_phi,_Vj,1);
1033 if (_verbose && _nbTimeStep%_freqSave ==0)
1035 cout << "Calcul diffusion: variables primitives maille " <<j <<endl;
1036 for(int q=0; q<_nVar; q++)
1037 cout << _Vj[q] << endl;
1040 //on n'a pas de contribution sur la masse
1042 //contribution visqueuse sur la quantite de mouvement
1043 for(int k=1; k<_nVar-1; k++)
1044 _phi[k] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu*(_porosityj*_Vj[k] - _porosityi*_Vi[k]);
1046 //contribution visqueuse sur l'energie
1047 _phi[_nVar-1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*lambda*(_porosityj*_Vj[_nVar-1] - _porosityi*_Vi[_nVar-1]);
1050 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1052 if(_verbose && _nbTimeStep%_freqSave ==0)
1054 cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
1055 for(int q=0; q<_nVar; q++)
1056 cout << _phi[q] << endl;
1062 //On change de signe pour l'autre contribution
1063 for(int k=0; k<_nVar; k++)
1064 _phi[k] *= -_inv_dxj/_inv_dxi;
1067 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1068 if(_verbose && _nbTimeStep%_freqSave ==0)
1070 cout << "Contribution diffusion au 2nd membre pour la maille " << j << ": "<<endl;
1071 for(int q=0; q<_nVar; q++)
1072 cout << _phi[q] << endl;
1077 if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
1079 cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
1080 for(int i=0; i<_nVar; i++)
1082 for(int j=0; j<_nVar; j++)
1083 cout << _Diffusion[i*_nVar+j]<<", ";
1090 void SinglePhase::sourceVector(PetscScalar * Si, PetscScalar * Ui, PetscScalar * Vi, int i)
1092 double phirho=Ui[0], T=Vi[_nVar-1];
1094 for(int k=0; k<_Ndim; k++)
1095 norm_u+=Vi[1+k]*Vi[1+k];
1096 norm_u=sqrt(norm_u);
1098 Si[0]=_heatPowerField(i)/_latentHeat;
1101 for(int k=1; k<_nVar-1; k++)
1102 Si[k] =(_gravite[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*phirho;
1104 Si[_nVar-1]=_heatPowerField(i);
1106 for(int k=0; k<_Ndim; k++)
1107 Si[_nVar-1] +=(_GravityField3d[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*Vi[1+k]*phirho;
1109 if(_timeScheme==Implicit)
1111 for(int k=0; k<_nVar*_nVar;k++)
1112 _GravityImplicitationMatrix[k] = 0;
1113 if(!_usePrimitiveVarsInNewton)
1114 for(int k=0; k<_nVar;k++)
1115 _GravityImplicitationMatrix[k*_nVar]=-_gravite[k];
1118 double pression=Vi[0];
1119 getDensityDerivatives( pression, T, norm_u*norm_u);
1120 for(int k=0; k<_nVar;k++)
1122 _GravityImplicitationMatrix[k*_nVar+0] =-_gravite[k]*_drho_sur_dp;
1123 _GravityImplicitationMatrix[k*_nVar+_nVar-1]=-_gravite[k]*_drho_sur_dT;
1127 if(_verbose && _nbTimeStep%_freqSave ==0)
1129 cout<<"SinglePhase::sourceVector"<<endl;
1131 for(int k=0;k<_nVar;k++)
1135 for(int k=0;k<_nVar;k++)
1139 for(int k=0;k<_nVar;k++)
1142 if(_timeScheme==Implicit)
1143 displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
1147 void SinglePhase::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
1149 double norm_u=0, u_n=0, rho;
1150 for(int i=0;i<_Ndim;i++)
1151 u_n += _Uroe[1+i]*_vec_normal[i];
1155 for(int i=0;i<_Ndim;i++)
1156 norm_u += Vi[1+i]*Vi[1+i];
1157 norm_u=sqrt(norm_u);
1159 for(int i=0;i<_Ndim;i++)
1160 pressureLoss[1+i]=-1/2*K*rho*norm_u*Vi[1+i];
1163 for(int i=0;i<_Ndim;i++)
1164 norm_u += Vj[1+i]*Vj[1+i];
1165 norm_u=sqrt(norm_u);
1167 for(int i=0;i<_Ndim;i++)
1168 pressureLoss[1+i]=-1/2*K*rho*norm_u*Vj[1+i];
1170 pressureLoss[_nVar-1]=-1/2*K*rho*norm_u*norm_u*norm_u;
1172 if(_verbose && _nbTimeStep%_freqSave ==0)
1174 cout<<"SinglePhase::pressureLossVector K= "<<K<<endl;
1176 for(int k=0;k<_nVar;k++)
1180 for(int k=0;k<_nVar;k++)
1184 for(int k=0;k<_nVar;k++)
1188 for(int k=0;k<_nVar;k++)
1191 cout<<"pressureLoss="<<endl;
1192 for(int k=0;k<_nVar;k++)
1193 cout<<pressureLoss[k]<<", ";
1198 void SinglePhase::porosityGradientSourceVector()
1200 double u_ni=0, u_nj=0, rhoi,rhoj, pi=_Vi[0], pj=_Vj[0],pij;
1201 for(int i=0;i<_Ndim;i++) {
1202 u_ni += _Vi[1+i]*_vec_normal[i];
1203 u_nj += _Vj[1+i]*_vec_normal[i];
1205 _porosityGradientSourceVector[0]=0;
1206 rhoj=_Uj[0]/_porosityj;
1207 rhoi=_Ui[0]/_porosityi;
1208 pij=(pi+pj)/2+rhoi*rhoj/2/(rhoi+rhoj)*(u_ni-u_nj)*(u_ni-u_nj);
1209 for(int i=0;i<_Ndim;i++)
1210 _porosityGradientSourceVector[1+i]=pij*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
1211 _porosityGradientSourceVector[_nVar-1]=0;
1214 void SinglePhase::jacobian(const int &j, string nameOfGroup,double * normale)
1216 if(_verbose && _nbTimeStep%_freqSave ==0)
1217 cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
1220 for(k=0; k<_nVar*_nVar;k++)
1221 _Jcb[k] = 0;//No implicitation at this stage
1224 for(k=1; k<_nVar; k++)
1225 _idm[k] = _idm[k-1] + 1;
1226 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);
1227 double q_n=0;//quantité de mouvement normale à la paroi
1228 for(k=0; k<_Ndim; k++)
1229 q_n+=_externalStates[(k+1)]*normale[k];
1231 // loop of boundary types
1232 if (_limitField[nameOfGroup].bcType==Wall)
1234 for(k=0; k<_nVar;k++)
1235 _Jcb[k*_nVar + k] = 1;
1236 for(k=1; k<_nVar-1;k++)
1237 for(int l=1; l<_nVar-1;l++)
1238 _Jcb[k*_nVar + l] -= 2*normale[k-1]*normale[l-1];
1240 else if (_limitField[nameOfGroup].bcType==Inlet)
1244 double v[_Ndim], ve[_Ndim], v2, ve2;
1247 for(k=1; k<_nVar; k++)
1248 _idm[k] = _idm[k-1] + 1;
1249 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1250 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1252 ve[0] = _limitField[nameOfGroup].v_x[0];
1257 ve[1] = _limitField[nameOfGroup].v_y[0];
1263 ve[2] = _limitField[nameOfGroup].v_z[0];
1268 double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,_Uj[0]);
1269 double total_energy=internal_energy+ve2/2;
1272 _Jcb[0]=v2/(2*internal_energy);
1273 for(k=0; k<_Ndim;k++)
1274 _Jcb[1+k]=-v[k]/internal_energy;
1275 _Jcb[_nVar-1]=1/internal_energy;
1277 for(int l =1;l<1+_Ndim;l++){
1278 _Jcb[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
1279 for(k=0; k<_Ndim;k++)
1280 _Jcb[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
1281 _Jcb[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
1284 _Jcb[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
1285 for(k=0; k<_Ndim;k++)
1286 _Jcb[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
1287 _Jcb[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
1290 for(k=0;k<_nVar;k++)
1295 _Jcb[_nVar]=_limitField[nameOfGroup].v_x[0];//Kieu
1296 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
1299 _Jcb[2*_nVar]= _limitField[nameOfGroup].v_y[0];
1300 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
1302 _Jcb[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1303 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
1306 _Jcb[(_nVar-1)*_nVar]=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2;
1309 else if (_limitField[nameOfGroup].bcType==InletPressure && q_n<0){
1311 double v[_Ndim], v2=0;
1313 for(k=1; k<_nVar; k++)
1314 _idm[k] = _idm[k-1] + 1;
1315 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1317 for(k=0; k<_Ndim;k++){
1322 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _limitField[nameOfGroup].T);
1323 double rho_int = _externalStates[0];
1324 double density_ratio=rho_ext/rho_int;
1326 for(int l =1;l<1+_Ndim;l++){
1327 _Jcb[l*_nVar]=-density_ratio*v[l-1];
1328 _Jcb[l*_nVar+l]=density_ratio;
1331 _Jcb[(_nVar-1)*_nVar]=-v2*density_ratio;
1332 for(k=0; k<_Ndim;k++)
1333 _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k];
1335 // not wall, not inlet, not inletPressure
1336 else if(_limitField[nameOfGroup].bcType==Outlet || (_limitField[nameOfGroup].bcType==InletPressure && q_n>=0))
1339 double v[_Ndim], v2=0;
1341 for(k=1; k<_nVar; k++)
1342 _idm[k] = _idm[k-1] + 1;
1343 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1345 for(k=0; k<_Ndim;k++){
1350 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _externalStates[_nVar-1]);
1351 double rho_int = _externalStates[0];
1352 double density_ratio=rho_ext/rho_int;
1353 double internal_energy=_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho_int);
1354 double total_energy=internal_energy+v2/2;
1357 _Jcb[0]=density_ratio*(1-v2/(2*internal_energy));
1358 for(k=0; k<_Ndim;k++)
1359 _Jcb[1+k]=density_ratio*v[k]/internal_energy;
1360 _Jcb[_nVar-1]=-density_ratio/internal_energy;
1362 for(int l =1;l<1+_Ndim;l++){
1363 _Jcb[l*_nVar]=density_ratio*v2*v[l-1]/(2*internal_energy);
1364 for(k=0; k<_Ndim;k++)
1365 _Jcb[l*_nVar+1+k]=density_ratio*v[k]*v[l-1]/internal_energy;
1366 _Jcb[l*_nVar+1+k]-=density_ratio;
1367 _Jcb[l*_nVar+_nVar-1]=-density_ratio*v[l-1]/internal_energy;
1370 _Jcb[(_nVar-1)*_nVar]=density_ratio*v2*total_energy/(2*internal_energy);
1371 for(k=0; k<_Ndim;k++)
1372 _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k]*total_energy/internal_energy;
1373 _Jcb[(_nVar-1)*_nVar+_nVar-1]=density_ratio*(1-total_energy/internal_energy);
1377 double cd = 1,cn=0,p0, gamma;
1378 _idm[0] = j*_nVar;// Kieu
1379 for(k=1; k<_nVar;k++)
1380 _idm[k] = _idm[k-1] + 1;
1381 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1382 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1384 // compute the common numerator and common denominator
1385 p0=_fluides[0]->constante("p0");
1386 gamma =_fluides[0]->constante("gamma");
1387 cn =_limitField[nameOfGroup].p +p0;
1388 cd = _phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0;
1392 for(k=1; k<_nVar-1;k++)
1393 v2+=_externalStates[k]*_externalStates[k];
1395 _JcbDiff[0] = cn*(_phi[_nVar-1] -v2 -p0)/cd;
1396 for(k=1; k<_nVar-1;k++)
1397 _JcbDiff[k]=cn*_phi[k]/cd;
1398 _JcbDiff[_nVar-1]= -cn*_phi[0]/cd;
1400 for(idim=0; idim<_Ndim;idim++)
1403 _JcbDiff[(1+idim)*_nVar]=-(v2*cn*_phi[idim+1])/(2*cd);
1404 //colonnes intermediaire
1405 for(jdim=0; jdim<_Ndim;jdim++)
1407 _JcbDiff[(1+idim)*_nVar + jdim + 1] =_externalStates[idim+1]*_phi[jdim+1];
1408 _JcbDiff[(1+idim)*_nVar + jdim + 1]*=cn/cd;
1410 //matrice identite*cn*(rhoe- p0)
1411 _JcbDiff[(1+idim)*_nVar + idim + 1] +=( cn*(_phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0))/cd;
1414 _JcbDiff[(1+idim)*_nVar + _nVar-1]=-_phi[idim+1]*cn/cd;
1417 _JcbDiff[_nVar*(_nVar-1)] = -(v2*_phi[_nVar -1]*cn)/(2*cd);
1418 for(int idim=0; idim<_Ndim;idim++)
1419 _JcbDiff[_nVar*(_nVar-1)+idim+1]=_externalStates[idim+1]*_phi[_nVar -1]*cn/cd;
1420 _JcbDiff[_nVar*_nVar -1] = -(v2/2+p0)*cn/cd;
1423 else if (_limitField[nameOfGroup].bcType!=Neumann)// not wall, not inlet, not outlet
1425 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1426 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1427 _runLogFile->close();
1428 throw CdmathException("SinglePhase::jacobian: This boundary condition is not treated");
1432 //calcule la jacobienne pour une CL de diffusion
1433 void SinglePhase::jacobianDiff(const int &j, string nameOfGroup)
1435 if(_verbose && _nbTimeStep%_freqSave ==0)
1436 cout<<"Jacobienne condition limite diffusion bord "<< nameOfGroup<<endl;
1439 for(k=0; k<_nVar*_nVar;k++)
1442 if (_limitField[nameOfGroup].bcType==Wall){
1443 double v[_Ndim], ve[_Ndim], v2, ve2;
1446 for(k=1; k<_nVar; k++)
1447 _idm[k] = _idm[k-1] + 1;
1448 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1449 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1451 ve[0] = _limitField[nameOfGroup].v_x[0];
1456 ve[1] = _limitField[nameOfGroup].v_y[0];
1462 ve[2] = _limitField[nameOfGroup].v_z[0];
1468 double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho);
1469 double total_energy=internal_energy+ve2/2;
1472 _JcbDiff[0]=v2/(2*internal_energy);
1473 for(k=0; k<_Ndim;k++)
1474 _JcbDiff[1+k]=-v[k]/internal_energy;
1475 _JcbDiff[_nVar-1]=1/internal_energy;
1477 for(int l =1;l<1+_Ndim;l++){
1478 _JcbDiff[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
1479 for(k=0; k<_Ndim;k++)
1480 _JcbDiff[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
1481 _JcbDiff[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
1484 _JcbDiff[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
1485 for(k=0; k<_Ndim;k++)
1486 _JcbDiff[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
1487 _JcbDiff[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
1489 else if (_limitField[nameOfGroup].bcType==Outlet || _limitField[nameOfGroup].bcType==Neumann
1490 ||_limitField[nameOfGroup].bcType==Inlet || _limitField[nameOfGroup].bcType==InletPressure)
1492 for(k=0;k<_nVar;k++)
1493 _JcbDiff[k*_nVar+k]=1;
1496 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1497 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1498 _runLogFile->close();
1499 throw CdmathException("SinglePhase::jacobianDiff: This boundary condition is not recognised");
1503 void SinglePhase::primToCons(const double *P, const int &i, double *W, const int &j){
1504 //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;
1506 double rho=_fluides[0]->getDensity(P[i*(_Ndim+2)], P[i*(_Ndim+2)+(_Ndim+1)]);
1507 W[j*(_Ndim+2)] = _porosityField(j)*rho;//phi*rho
1508 for(int k=0; k<_Ndim; k++)
1509 W[j*(_Ndim+2)+(k+1)] = W[j*(_Ndim+2)]*P[i*(_Ndim+2)+(k+1)];//phi*rho*u
1511 W[j*(_Ndim+2)+(_Ndim+1)] = W[j*(_Ndim+2)]*_fluides[0]->getInternalEnergy(P[i*(_Ndim+2)+ (_Ndim+1)],rho);//rho*e
1512 for(int k=0; k<_Ndim; k++)
1513 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
1516 void SinglePhase::primToConsJacobianMatrix(double *V)
1518 double pression=V[0];
1519 double temperature=V[_nVar-1];
1520 double vitesse[_Ndim];
1521 for(int idim=0;idim<_Ndim;idim++)
1522 vitesse[idim]=V[1+idim];
1524 for(int idim=0;idim<_Ndim;idim++)
1525 v2+=vitesse[idim]*vitesse[idim];
1527 double rho=_fluides[0]->getDensity(pression,temperature);
1528 double gamma=_fluides[0]->constante("gamma");
1529 double Pinf=_fluides[0]->constante("p0");
1530 double q=_fluides[0]->constante("q");
1531 double cv=_fluides[0]->constante("cv");
1533 for(int k=0;k<_nVar*_nVar; k++)
1534 _primToConsJacoMat[k]=0;
1536 if( !_useDellacherieEOS)
1538 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1539 double e=fluide0->getInternalEnergy(temperature);
1542 _primToConsJacoMat[0]=1/((gamma-1)*(e-q));
1543 _primToConsJacoMat[_nVar-1]=-rho*cv/(e-q);
1545 for(int idim=0;idim<_Ndim;idim++)
1547 _primToConsJacoMat[_nVar+idim*_nVar]=vitesse[idim]/((gamma-1)*(e-q));
1548 _primToConsJacoMat[_nVar+idim*_nVar+1+idim]=rho;
1549 _primToConsJacoMat[_nVar+idim*_nVar+_nVar-1]=-rho*vitesse[idim]*cv/(e-q);
1551 _primToConsJacoMat[(_nVar-1)*_nVar]=E/((gamma-1)*(e-q));
1552 for(int idim=0;idim<_Ndim;idim++)
1553 _primToConsJacoMat[(_nVar-1)*_nVar+1+idim]=rho*vitesse[idim];
1554 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*cv*(1-E/(e-q));
1556 else if( _useDellacherieEOS)
1558 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1559 double h=fluide0->getEnthalpy(temperature);
1561 double cp=_fluides[0]->constante("cp");
1563 _primToConsJacoMat[0]=gamma/((gamma-1)*(h-q));
1564 _primToConsJacoMat[_nVar-1]=-rho*cp/(h-q);
1566 for(int idim=0;idim<_Ndim;idim++)
1568 _primToConsJacoMat[_nVar+idim*_nVar]=gamma*vitesse[idim]/((gamma-1)*(h-q));
1569 _primToConsJacoMat[_nVar+idim*_nVar+1+idim]=rho;
1570 _primToConsJacoMat[_nVar+idim*_nVar+_nVar-1]=-rho*vitesse[idim]*cp/(h-q);
1572 _primToConsJacoMat[(_nVar-1)*_nVar]=gamma*H/((gamma-1)*(h-q))-1;
1573 for(int idim=0;idim<_Ndim;idim++)
1574 _primToConsJacoMat[(_nVar-1)*_nVar+1+idim]=rho*vitesse[idim];
1575 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*cp*(1-H/(h-q));
1579 *_runLogFile<<"SinglePhase::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie"<<endl;
1580 _runLogFile->close();
1581 throw CdmathException("SinglePhase::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
1584 if(_verbose && _nbTimeStep%_freqSave ==0)
1586 cout<<" SinglePhase::primToConsJacobianMatrix" << endl;
1587 displayVector(_Vi,_nVar," _Vi " );
1588 cout<<" Jacobienne primToCons: " << endl;
1589 displayMatrix(_primToConsJacoMat,_nVar," Jacobienne primToCons: ");
1593 void SinglePhase::consToPrim(const double *Wcons, double* Wprim,double porosity)
1596 for(int k=1;k<=_Ndim;k++)
1597 q_2 += Wcons[k]*Wcons[k];
1598 q_2 /= Wcons[0]; //phi rho u²
1599 double rhoe=(Wcons[(_Ndim+2)-1]-q_2/2)/porosity;
1600 double rho=Wcons[0]/porosity;
1601 Wprim[0] = _fluides[0]->getPressure(rhoe,rho);//pressure p
1603 cout << "pressure = "<< Wprim[0] << " < 0 " << endl;
1604 *_runLogFile<< "pressure = "<< Wprim[0] << " < 0 " << endl;
1605 _runLogFile->close();
1606 throw CdmathException("SinglePhase::consToPrim: negative pressure");
1608 for(int k=1;k<=_Ndim;k++)
1609 Wprim[k] = Wcons[k]/Wcons[0];//velocity u
1610 Wprim[(_Ndim+2)-1] = _fluides[0]->getTemperatureFromPressure(Wprim[0],Wcons[0]/porosity);
1612 if(_verbose && _nbTimeStep%_freqSave ==0)
1614 cout<<"ConsToPrim Vecteur conservatif"<<endl;
1615 for(int k=0;k<_nVar;k++)
1616 cout<<Wcons[k]<<endl;
1617 cout<<"ConsToPrim Vecteur primitif"<<endl;
1618 for(int k=0;k<_nVar;k++)
1619 cout<<Wprim[k]<<endl;
1623 void SinglePhase::entropicShift(double* n)//TO do: make sure _Vi and _Vj are well set
1625 /*Left and right values */
1626 double ul_n = 0, ul_2=0, ur_n=0, ur_2 = 0; //valeurs de vitesse normale et |u|² a droite et a gauche
1627 for(int i=0;i<_Ndim;i++)
1629 ul_n += _Vi[1+i]*n[i];
1630 ul_2 += _Vi[1+i]*_Vi[1+i];
1631 ur_n += _Vj[1+i]*n[i];
1632 ur_2 += _Vj[1+i]*_Vj[1+i];
1636 double cl = _fluides[0]->vitesseSonEnthalpie(_Vi[_Ndim+1]-ul_2/2);//vitesse du son a l'interface
1637 double cr = _fluides[0]->vitesseSonEnthalpie(_Vj[_Ndim+1]-ur_2/2);//vitesse du son a l'interface
1639 _entropicShift[0]=fabs(ul_n-cl - (ur_n-cr));
1640 _entropicShift[1]=fabs(ul_n - ur_n);
1641 _entropicShift[2]=fabs(ul_n+cl - (ur_n+cr));
1643 if(_verbose && _nbTimeStep%_freqSave ==0)
1645 cout<<"Entropic shift left eigenvalues: "<<endl;
1646 cout<<"("<< ul_n-cl<< ", " <<ul_n<<", "<<ul_n+cl << ")";
1648 cout<<"Entropic shift right eigenvalues: "<<endl;
1649 cout<<"("<< ur_n-cr<< ", " <<ur_n<<", "<<ur_n+cr << ")";
1651 cout<<"eigenvalue jumps "<<endl;
1652 cout<< _entropicShift[0] << ", " << _entropicShift[1] << ", "<< _entropicShift[2] <<endl;
1656 Vector SinglePhase::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
1657 if(_verbose && _nbTimeStep%_freqSave ==0)
1659 cout<<"SinglePhase::convectionFlux start"<<endl;
1660 cout<<"Ucons"<<endl;
1662 cout<<"Vprim"<<endl;
1666 double phirho=U(0);//phi rho
1667 Vector phiq(_Ndim);//phi rho u
1668 for(int i=0;i<_Ndim;i++)
1671 double pression=V(0);
1672 Vector vitesse(_Ndim);
1673 for(int i=0;i<_Ndim;i++)
1675 double Temperature= V(1+_Ndim);
1677 double vitessen=vitesse*normale;
1678 double rho=phirho/porosity;
1679 double e_int=_fluides[0]->getInternalEnergy(Temperature,rho);
1682 F(0)=phirho*vitessen;
1683 for(int i=0;i<_Ndim;i++)
1684 F(1+i)=phirho*vitessen*vitesse(i)+pression*normale(i)*porosity;
1685 F(1+_Ndim)=phirho*(e_int+0.5*vitesse*vitesse+pression/rho)*vitessen;
1687 if(_verbose && _nbTimeStep%_freqSave ==0)
1689 cout<<"SinglePhase::convectionFlux end"<<endl;
1690 cout<<"Flux F(U,V)"<<endl;
1697 void SinglePhase::RoeMatrixConservativeVariables(double u_n, double H,Vector velocity, double k, double K)
1699 /******** Construction de la matrice de Roe *********/
1700 //premiere ligne (masse)
1702 for(int idim=0; idim<_Ndim;idim++)
1703 _Aroe[1+idim]=_vec_normal[idim];
1706 //lignes intermadiaires(qdm)
1707 for(int idim=0; idim<_Ndim;idim++)
1710 _Aroe[(1+idim)*_nVar]=K*_vec_normal[idim] - u_n*_Uroe[1+idim];
1711 //colonnes intermediaires
1712 for(int jdim=0; jdim<_Ndim;jdim++)
1713 _Aroe[(1+idim)*_nVar + jdim + 1] = _Uroe[1+idim]*_vec_normal[jdim]-k*_vec_normal[idim]*_Uroe[1+jdim];
1715 _Aroe[(1+idim)*_nVar + idim + 1] += u_n;
1717 _Aroe[(1+idim)*_nVar + _nVar-1]=k*_vec_normal[idim];
1720 //derniere ligne (energie)
1721 _Aroe[_nVar*(_nVar-1)] = (K - H)*u_n;
1722 for(int idim=0; idim<_Ndim;idim++)
1723 _Aroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] - k*u_n*_Uroe[idim+1];
1724 _Aroe[_nVar*_nVar -1] = (1 + k)*u_n;
1726 void SinglePhase::convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector vitesse)
1728 //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1729 //On remplit la matrice de Roe en variables primitives : F(V_L)-F(V_R)=Aroe (V_L-V_R)
1730 //EOS is more involved with primitive variables
1731 // call to getDensityDerivatives(double concentration, double pression, double temperature,double v2) needed
1732 _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1733 for(int i=0;i<_Ndim;i++)
1734 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1735 _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1736 for(int i=0;i<_Ndim;i++)
1738 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]+_vec_normal[i];
1739 for(int j=0;j<_Ndim;j++)
1740 _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1741 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1742 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1744 _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1745 for(int i=0;i<_Ndim;i++)
1746 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1747 _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1749 void SinglePhase::staggeredRoeUpwindingMatrixConservativeVariables( double u_n, double H,Vector velocity, double k, double K)
1751 //Calcul de décentrement de type décalé pour formulation de Roe
1752 if(fabs(u_n)>_precision)
1754 //premiere ligne (masse)
1756 for(int idim=0; idim<_Ndim;idim++)
1757 _absAroe[1+idim]=_vec_normal[idim];
1758 _absAroe[_nVar-1]=0;
1760 //lignes intermadiaires(qdm)
1761 for(int idim=0; idim<_Ndim;idim++)
1764 _absAroe[(1+idim)*_nVar]=-K*_vec_normal[idim] - u_n*_Uroe[1+idim];
1765 //colonnes intermediaires
1766 for(int jdim=0; jdim<_Ndim;jdim++)
1767 _absAroe[(1+idim)*_nVar + jdim + 1] = _Uroe[1+idim]*_vec_normal[jdim]+k*_vec_normal[idim]*_Uroe[1+jdim];
1769 _absAroe[(1+idim)*_nVar + idim + 1] += u_n;
1771 _absAroe[(1+idim)*_nVar + _nVar-1]=-k*_vec_normal[idim];
1774 //derniere ligne (energie)
1775 _absAroe[_nVar*(_nVar-1)] = (-K - H)*u_n;
1776 for(int idim=0; idim<_Ndim;idim++)
1777 _absAroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] + k*u_n*_Uroe[idim+1];
1778 _absAroe[_nVar*_nVar -1] = (1 - k)*u_n;
1786 for(int i=0; i<_nVar*_nVar;i++)
1787 _absAroe[i] *= signu;
1789 else//umn=0 ->centered scheme
1791 for(int i=0; i<_nVar*_nVar;i++)
1795 void SinglePhase::staggeredRoeUpwindingMatrixPrimitiveVariables(double rho, double u_n,double H, Vector vitesse)
1797 //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1798 //Calcul de décentrement de type décalé pour formulation Roe
1799 _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1800 for(int i=0;i<_Ndim;i++)
1801 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1802 _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1803 for(int i=0;i<_Ndim;i++)
1805 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]-_vec_normal[i];
1806 for(int j=0;j<_Ndim;j++)
1807 _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1808 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1809 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1811 _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1812 for(int i=0;i<_Ndim;i++)
1813 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1814 _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1817 Vector SinglePhase::staggeredVFFCFlux()
1819 if(_verbose && _nbTimeStep%_freqSave ==0)
1820 cout<<"SinglePhase::staggeredVFFCFlux() start"<<endl;
1822 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1824 *_runLogFile<< "SinglePhase::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding, pressure = "<< endl;
1825 _runLogFile->close();
1826 throw CdmathException("SinglePhase::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding");
1828 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
1832 double uijn=0, phiqn=0, uin=0, ujn=0;
1833 for(int idim=0;idim<_Ndim;idim++)
1835 uijn+=_vec_normal[idim]*_Uroe[1+idim];//URoe = rho, u, H
1836 uin +=_vec_normal[idim]*_Ui[1+idim];
1837 ujn +=_vec_normal[idim]*_Uj[1+idim];
1840 if( (uin>0 && ujn >0) || (uin>=0 && ujn <=0 && uijn>0) ) // formerly (uijn>_precision)
1842 for(int idim=0;idim<_Ndim;idim++)
1843 phiqn+=_vec_normal[idim]*_Ui[1+idim];//phi rho u n
1845 for(int idim=0;idim<_Ndim;idim++)
1846 Fij(1+idim)=phiqn*_Vi[1+idim]+_Vj[0]*_vec_normal[idim]*_porosityj;
1847 Fij(_nVar-1)=phiqn/_Ui[0]*(_Ui[_nVar-1]+_Vj[0]*sqrt(_porosityj/_porosityi));
1849 else if( (uin<0 && ujn <0) || (uin>=0 && ujn <=0 && uijn<0) ) // formerly (uijn<-_precision)
1851 for(int idim=0;idim<_Ndim;idim++)
1852 phiqn+=_vec_normal[idim]*_Uj[1+idim];//phi rho u n
1854 for(int idim=0;idim<_Ndim;idim++)
1855 Fij(1+idim)=phiqn*_Vj[1+idim]+_Vi[0]*_vec_normal[idim]*_porosityi;
1856 Fij(_nVar-1)=phiqn/_Uj[0]*(_Uj[_nVar-1]+_Vi[0]*sqrt(_porosityi/_porosityj));
1858 else//case (uin<=0 && ujn >=0) or (uin>=0 && ujn <=0 && uijn==0), apply centered scheme
1860 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar);
1861 Vector normale(_Ndim);
1862 for(int i1=0;i1<_Ndim;i1++)
1863 normale(i1)=_vec_normal[i1];
1864 for(int i1=0;i1<_nVar;i1++)
1871 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
1872 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
1873 Fij=(Fi+Fj)/2;//+_maxvploc*(Ui-Uj)/2;
1875 if(_verbose && _nbTimeStep%_freqSave ==0)
1877 cout<<"SinglePhase::staggeredVFFCFlux() endf uijn="<<uijn<<endl;
1884 void SinglePhase::staggeredVFFCMatricesConservativeVariables(double un)//vitesse normale de Roe en entrée
1886 if(_verbose && _nbTimeStep%_freqSave ==0)
1887 cout<<"SinglePhase::staggeredVFFCMatrices()"<<endl;
1889 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1891 *_runLogFile<< "SinglePhase::staggeredVFFCMatrices: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding"<< endl;
1892 _runLogFile->close();
1893 throw CdmathException("SinglePhase::staggeredVFFCMatrices: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding");
1895 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
1897 double ui_n=0, ui_2=0, uj_n=0, uj_2=0;//vitesse normale et carré du module
1898 double H;//enthalpie totale (expression particulière)
1899 consToPrim(_Ui,_Vi,_porosityi);
1900 consToPrim(_Uj,_Vj,_porosityj);
1901 for(int i=0;i<_Ndim;i++)
1903 ui_2 += _Vi[1+i]*_Vi[1+i];
1904 ui_n += _Vi[1+i]*_vec_normal[i];
1905 uj_2 += _Vj[1+i]*_Vj[1+i];
1906 uj_n += _Vj[1+i]*_vec_normal[i];
1909 double rhoi,pj, Ei, rhoj;
1910 double cj, Kj, kj;//dérivées de la pression
1911 rhoi=_Ui[0]/_porosityi;
1912 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
1914 rhoj=_Uj[0]/_porosityj;
1916 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
1917 kj = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
1918 Kj = uj_2*kj/2; //g-1/2 *|u|²
1921 double ci, Ki, ki;//dérivées de la pression
1922 Ej= _Uj[_Ndim+1]/rhoj;
1925 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
1926 ki = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
1927 Ki = ui_2*ki/2; //g-1/2 *|u|²
1931 /***********Calcul des valeurs propres ********/
1932 vector<complex<double> > vp_dist(3,0);
1933 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
1935 _maxvploc=fabs(ui_n)+cj;
1936 if(_maxvploc>_maxvp)
1939 if(_verbose && _nbTimeStep%_freqSave ==0)
1940 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
1942 /******** Construction de la matrice A^+ *********/
1943 //premiere ligne (masse)
1945 for(int idim=0; idim<_Ndim;idim++)
1946 _AroePlus[1+idim]=_vec_normal[idim];
1947 _AroePlus[_nVar-1]=0;
1949 //lignes intermadiaires(qdm)
1950 for(int idim=0; idim<_Ndim;idim++)
1953 _AroePlus[(1+idim)*_nVar]=- ui_n*_Vi[1+idim];
1954 //colonnes intermediaires
1955 for(int jdim=0; jdim<_Ndim;jdim++)
1956 _AroePlus[(1+idim)*_nVar + jdim + 1] = _Vi[1+idim]*_vec_normal[jdim];
1958 _AroePlus[(1+idim)*_nVar + idim + 1] += ui_n;
1960 _AroePlus[(1+idim)*_nVar + _nVar-1]=0;
1963 //derniere ligne (energie)
1964 _AroePlus[_nVar*(_nVar-1)] = - H*ui_n;
1965 for(int idim=0; idim<_Ndim;idim++)
1966 _AroePlus[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] ;
1967 _AroePlus[_nVar*_nVar -1] = ui_n;
1969 /******** Construction de la matrice A^- *********/
1970 //premiere ligne (masse)
1972 for(int idim=0; idim<_Ndim;idim++)
1973 _AroeMinus[1+idim]=0;
1974 _AroeMinus[_nVar-1]=0;
1976 //lignes intermadiaires(qdm)
1977 for(int idim=0; idim<_Ndim;idim++)
1980 _AroeMinus[(1+idim)*_nVar]=Kj*_vec_normal[idim] ;
1981 //colonnes intermediaires
1982 for(int jdim=0; jdim<_Ndim;jdim++)
1983 _AroeMinus[(1+idim)*_nVar + jdim + 1] = -kj*_vec_normal[idim]*_Vj[1+jdim];
1985 _AroeMinus[(1+idim)*_nVar + idim + 1] += 0;
1987 _AroeMinus[(1+idim)*_nVar + _nVar-1]=kj*_vec_normal[idim];
1990 //derniere ligne (energie)
1991 _AroeMinus[_nVar*(_nVar-1)] = Kj *ui_n;
1992 for(int idim=0; idim<_Ndim;idim++)
1993 _AroeMinus[_nVar*(_nVar-1)+idim+1]= - kj*uj_n*_Vi[idim+1];
1994 _AroeMinus[_nVar*_nVar -1] = kj*ui_n;
1996 else if(un<-_precision)
1998 /***********Calcul des valeurs propres ********/
1999 vector<complex<double> > vp_dist(3,0);
2000 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2002 _maxvploc=fabs(uj_n)+ci;
2003 if(_maxvploc>_maxvp)
2006 if(_verbose && _nbTimeStep%_freqSave ==0)
2007 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2009 /******** Construction de la matrice A^+ *********/
2010 //premiere ligne (masse)
2012 for(int idim=0; idim<_Ndim;idim++)
2013 _AroePlus[1+idim]=0;
2014 _AroePlus[_nVar-1]=0;
2016 //lignes intermadiaires(qdm)
2017 for(int idim=0; idim<_Ndim;idim++)
2020 _AroePlus[(1+idim)*_nVar]=Ki*_vec_normal[idim] ;
2021 //colonnes intermediaires
2022 for(int jdim=0; jdim<_Ndim;jdim++)
2023 _AroePlus[(1+idim)*_nVar + jdim + 1] = -ki*_vec_normal[idim]*_Vi[1+jdim];
2025 _AroePlus[(1+idim)*_nVar + idim + 1] += 0;
2027 _AroePlus[(1+idim)*_nVar + _nVar-1]=ki*_vec_normal[idim];
2030 //derniere ligne (energie)
2031 _AroePlus[_nVar*(_nVar-1)] = Ki *uj_n;
2032 for(int idim=0; idim<_Ndim;idim++)
2033 _AroePlus[_nVar*(_nVar-1)+idim+1]= - ki*ui_n*_Vj[idim+1];
2034 _AroePlus[_nVar*_nVar -1] = ki*uj_n;
2036 /******** Construction de la matrice A^- *********/
2037 //premiere ligne (masse)
2039 for(int idim=0; idim<_Ndim;idim++)
2040 _AroeMinus[1+idim]=_vec_normal[idim];
2041 _AroeMinus[_nVar-1]=0;
2043 //lignes intermadiaires(qdm)
2044 for(int idim=0; idim<_Ndim;idim++)
2047 _AroeMinus[(1+idim)*_nVar]= - uj_n*_Vj[1+idim];
2048 //colonnes intermediaires
2049 for(int jdim=0; jdim<_Ndim;jdim++)
2050 _AroeMinus[(1+idim)*_nVar + jdim + 1] = _Vj[1+idim]*_vec_normal[jdim];
2052 _AroeMinus[(1+idim)*_nVar + idim + 1] += uj_n;
2054 _AroeMinus[(1+idim)*_nVar + _nVar-1]=0;
2057 //derniere ligne (energie)
2058 _AroeMinus[_nVar*(_nVar-1)] = - H*uj_n;
2059 for(int idim=0; idim<_Ndim;idim++)
2060 _AroeMinus[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] ;
2061 _AroeMinus[_nVar*_nVar -1] = uj_n;
2063 else//case nil velocity on the interface, apply centered scheme
2065 double u_n=0, u_2=0;//vitesse normale et carré du module
2066 for(int i=0;i<_Ndim;i++)
2068 u_2 += _Uroe[1+i]*_Uroe[1+i];
2069 u_n += _Uroe[1+i]*_vec_normal[i];
2071 Vector vitesse(_Ndim);
2072 for(int idim=0;idim<_Ndim;idim++)
2073 vitesse[idim]=_Uroe[1+idim];
2076 /***********Calcul des valeurs propres ********/
2078 c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
2079 k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
2080 K = u_2*k/2; //g-1/2 *|u|²
2082 _maxvploc=fabs(u_n)+c;
2083 if(_maxvploc>_maxvp)
2086 /******** Construction de la matrice A^+ *********/
2087 //premiere ligne (masse)
2089 for(int idim=0; idim<_Ndim;idim++)
2090 _AroePlus[1+idim]=0;
2091 _AroePlus[_nVar-1]=0;
2093 //lignes intermadiaires(qdm)
2094 for(int idim=0; idim<_Ndim;idim++)
2097 _AroePlus[(1+idim)*_nVar]=- ui_n*_Vi[1+idim];
2098 //colonnes intermediaires
2099 for(int jdim=0; jdim<_Ndim;jdim++)
2100 _AroePlus[(1+idim)*_nVar + jdim + 1] = _Vi[1+idim]*_vec_normal[jdim]-0.5*ki*_vec_normal[idim]*_Vi[1+jdim];
2102 _AroePlus[(1+idim)*_nVar + idim + 1] += 0.5*ui_n;
2104 _AroePlus[(1+idim)*_nVar + _nVar-1]=0.5*ki*_vec_normal[idim];
2107 //derniere ligne (energie)
2108 _AroePlus[_nVar*(_nVar-1)] = 0;
2109 for(int idim=0; idim<_Ndim;idim++)
2110 _AroePlus[_nVar*(_nVar-1)+idim+1]=0 ;
2111 _AroePlus[_nVar*_nVar -1] = 0;
2113 /******** Construction de la matrice A^- *********/
2114 //premiere ligne (masse)
2116 for(int idim=0; idim<_Ndim;idim++)
2117 _AroeMinus[1+idim]=0;
2118 _AroeMinus[_nVar-1]=0;
2120 //lignes intermadiaires(qdm)
2121 for(int idim=0; idim<_Ndim;idim++)
2124 _AroeMinus[(1+idim)*_nVar]=Kj*_vec_normal[idim] - uj_n*_Vj[1+idim];
2125 //colonnes intermediaires
2126 for(int jdim=0; jdim<_Ndim;jdim++)
2127 _AroeMinus[(1+idim)*_nVar + jdim + 1] = -0.5*kj*_vec_normal[idim]*_Vj[1+jdim];
2129 _AroeMinus[(1+idim)*_nVar + idim + 1] += 0.5*uj_n;
2131 _AroeMinus[(1+idim)*_nVar + _nVar-1]=0.5*kj*_vec_normal[idim];
2134 //derniere ligne (energie)
2135 _AroeMinus[_nVar*(_nVar-1)] = 0;
2136 for(int idim=0; idim<_Ndim;idim++)
2137 _AroeMinus[_nVar*(_nVar-1)+idim+1]= 0;
2138 _AroeMinus[_nVar*_nVar -1] = 0;
2141 if(_timeScheme==Implicit)
2142 for(int i=0; i<_nVar*_nVar;i++)
2144 _AroeMinusImplicit[i] = _AroeMinus[i];
2145 _AroePlusImplicit[i] = _AroePlus[i];
2148 /******** Construction de la matrices Aroe *********/
2150 //premiere ligne (masse)
2152 for(int idim=0; idim<_Ndim;idim++)
2153 _Aroe[1+idim]=_vec_normal[idim];
2156 //lignes intermadiaires(qdm)
2157 for(int idim=0; idim<_Ndim;idim++)
2160 _Aroe[(1+idim)*_nVar]=Ki*_vec_normal[idim] - uj_n*_Uj[1+idim];
2161 //colonnes intermediaires
2162 for(int jdim=0; jdim<_Ndim;jdim++)
2163 _Aroe[(1+idim)*_nVar + jdim + 1] = _Uj[1+idim]*_vec_normal[jdim]-ki*_vec_normal[idim]*_Ui[1+jdim];
2165 _Aroe[(1+idim)*_nVar + idim + 1] += uj_n;
2167 _Aroe[(1+idim)*_nVar + _nVar-1]=ki*_vec_normal[idim];
2170 //derniere ligne (energie)
2171 _Aroe[_nVar*(_nVar-1)] = (Ki - H)*uj_n;
2172 for(int idim=0; idim<_Ndim;idim++)
2173 _Aroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] - ki*ui_n*_Uj[idim+1];
2174 _Aroe[_nVar*_nVar -1] = (1 + ki)*uj_n;
2178 void SinglePhase::staggeredVFFCMatricesPrimitiveVariables(double un)//vitesse normale de Roe en entrée
2180 if(_verbose && _nbTimeStep%_freqSave ==0)
2181 cout<<"SinglePhase::staggeredVFFCMatricesPrimitiveVariables()"<<endl;
2183 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2185 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding" << endl;
2186 _runLogFile->close();
2187 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding");
2189 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
2191 double ui_n=0., ui_2=0., uj_n=0., uj_2=0.;//vitesse normale et carré du module
2192 double H;//enthalpie totale (expression particulière)
2193 consToPrim(_Ui,_Vi,_porosityi);
2194 consToPrim(_Uj,_Vj,_porosityj);
2196 for(int i=0;i<_Ndim;i++)
2198 ui_2 += _Vi[1+i]*_Vi[1+i];
2199 ui_n += _Vi[1+i]*_vec_normal[i];
2200 uj_2 += _Vj[1+i]*_Vj[1+i];
2201 uj_n += _Vj[1+i]*_vec_normal[i];
2204 if(_verbose && _nbTimeStep%_freqSave ==0){
2205 cout <<"SinglePhase::staggeredVFFCMatricesPrimitiveVariables " << endl;
2206 cout<<"Vecteur primitif _Vi" << endl;
2207 for(int i=0;i<_nVar;i++)
2210 cout<<"Vecteur primitif _Vj" << endl;
2211 for(int i=0;i<_nVar;i++)
2216 double gamma=_fluides[0]->constante("gamma");
2217 double q=_fluides[0]->constante("q");
2219 if(fabs(un)>_precision)//non zero velocity on the interface
2221 if( !_useDellacherieEOS)
2223 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2224 double cv=fluide0->constante("cv");
2228 double rhoi,rhoj,pj, Ei, ei;
2229 double cj;//vitesse du son pour calcul valeurs propres
2230 rhoi=_Ui[0]/_porosityi;
2231 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
2234 rhoj=_Uj[0]/_porosityj;
2235 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
2237 /***********Calcul des valeurs propres ********/
2238 vector<complex<double> > vp_dist(3,0);
2239 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
2241 _maxvploc=fabs(ui_n)+cj;
2242 if(_maxvploc>_maxvp)
2245 if(_verbose && _nbTimeStep%_freqSave ==0)
2246 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
2248 /******** Construction de la matrice A^+ *********/
2249 //premiere ligne (masse)
2250 _AroePlusImplicit[0]=ui_n/((gamma-1)*(ei-q));
2251 for(int idim=0; idim<_Ndim;idim++)
2252 _AroePlusImplicit[1+idim]=rhoi*_vec_normal[idim];
2253 _AroePlusImplicit[_nVar-1]=-rhoi*ui_n*cv/(ei-q);
2255 //lignes intermadiaires(qdm)
2256 for(int idim=0; idim<_Ndim;idim++)
2259 _AroePlusImplicit[(1+idim)*_nVar]=ui_n/((gamma-1)*(ei-q))*_Vi[1+idim];
2260 //colonnes intermediaires
2261 for(int jdim=0; jdim<_Ndim;jdim++)
2262 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = rhoi*_Vi[1+idim]*_vec_normal[jdim];
2264 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += rhoi*ui_n;
2266 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoi*ui_n*cv/(ei-q)*_Vi[1+idim];
2269 //derniere ligne (energie)
2270 _AroePlusImplicit[_nVar*(_nVar-1)] = Ei*ui_n/((gamma-1)*(ei-q));
2271 for(int idim=0; idim<_Ndim;idim++)
2272 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoi*Ei+pj)*_vec_normal[idim]+rhoi*ui_n*_Vi[1+idim];
2273 _AroePlusImplicit[_nVar*_nVar -1] = rhoi*ui_n*(1-Ei/(ei-q))*cv;
2275 /******** Construction de la matrice A^- *********/
2276 //premiere ligne (masse)
2277 _AroeMinusImplicit[0]=0;
2278 for(int idim=0; idim<_Ndim;idim++)
2279 _AroeMinusImplicit[1+idim]=0;
2280 _AroeMinusImplicit[_nVar-1]=0;
2282 //lignes intermadiaires(qdm)
2283 for(int idim=0; idim<_Ndim;idim++)
2286 _AroeMinusImplicit[(1+idim)*_nVar]=_vec_normal[idim] ;
2287 //colonnes intermediaires
2288 for(int jdim=0; jdim<_Ndim;jdim++)
2289 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2291 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2293 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2296 //derniere ligne (energie)
2297 _AroeMinusImplicit[_nVar*(_nVar-1)] = ui_n;
2298 for(int idim=0; idim<_Ndim;idim++)
2299 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2300 _AroeMinusImplicit[_nVar*_nVar -1] = 0;
2302 else if(un<-_precision)
2304 double pi, Ej, rhoi, rhoj, ej;
2305 double ci;//vitesse du son pour calcul valeurs propres
2306 rhoj=_Uj[0]/_porosityj;
2307 Ej= _Uj[_Ndim+1]/rhoj;
2310 rhoi=_Ui[0]/_porosityi;
2311 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
2313 /***********Calcul des valeurs propres ********/
2314 vector<complex<double> > vp_dist(3,0);
2315 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2317 _maxvploc=fabs(uj_n)+ci;
2318 if(_maxvploc>_maxvp)
2321 if(_verbose && _nbTimeStep%_freqSave ==0)
2322 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2324 /******** Construction de la matrice A^+ *********/
2325 //premiere ligne (masse)
2326 _AroePlusImplicit[0]=0;
2327 for(int idim=0; idim<_Ndim;idim++)
2328 _AroePlusImplicit[1+idim]=0;
2329 _AroePlusImplicit[_nVar-1]=0;
2331 //lignes intermadiaires(qdm)
2332 for(int idim=0; idim<_Ndim;idim++)
2335 _AroePlusImplicit[(1+idim)*_nVar]=0;
2336 //colonnes intermediaires
2337 for(int jdim=0; jdim<_Ndim;jdim++)
2338 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2340 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2342 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2345 //derniere ligne (energie)
2346 _AroePlusImplicit[_nVar*(_nVar-1)] = uj_n;
2347 for(int idim=0; idim<_Ndim;idim++)
2348 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2349 _AroePlusImplicit[_nVar*_nVar -1] = 0;
2351 /******** Construction de la matrice A^- *********/
2352 //premiere ligne (masse)
2353 _AroeMinusImplicit[0]=uj_n/((gamma-1)*(ej-q));
2354 for(int idim=0; idim<_Ndim;idim++)
2355 _AroeMinusImplicit[1+idim]=rhoj*_vec_normal[idim];
2356 _AroeMinusImplicit[_nVar-1]=-rhoj*uj_n*cv/(ej-q);
2358 //lignes intermadiaires(qdm)
2359 for(int idim=0; idim<_Ndim;idim++)
2362 _AroeMinusImplicit[(1+idim)*_nVar]= uj_n/((gamma-1)*(ej-q))*_Vj[1+idim];
2363 //colonnes intermediaires
2364 for(int jdim=0; jdim<_Ndim;jdim++)
2365 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = rhoj*_Vj[1+idim]*_vec_normal[jdim];
2367 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += rhoj*uj_n;
2369 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoj*uj_n*cv/(ej-q)*_Vj[1+idim];
2372 //derniere ligne (energie)
2373 _AroeMinusImplicit[_nVar*(_nVar-1)] = Ej*uj_n/((gamma-1)*(ej-q));
2374 for(int idim=0; idim<_Ndim;idim++)
2375 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoj*Ej+pi)*_vec_normal[idim]+rhoj*uj_n*_Vj[1+idim];
2376 _AroeMinusImplicit[_nVar*_nVar -1] = rhoj*uj_n*(1-Ej/(ej-q))*cv;
2380 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero" << endl;
2381 _runLogFile->close();
2382 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero");
2385 else if(_useDellacherieEOS )
2387 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2388 double cp=fluide0->constante("cp");
2392 double rhoi,rhoj,pj, Ei, hi, Hi;
2393 double cj;//vitesse du son pour calcul valeurs propres
2394 rhoi=_Ui[0]/_porosityi;
2395 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
2399 rhoj=_Uj[0]/_porosityj;
2400 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
2402 /***********Calcul des valeurs propres ********/
2403 vector<complex<double> > vp_dist(3,0);
2404 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
2406 _maxvploc=fabs(ui_n)+cj;
2407 if(_maxvploc>_maxvp)
2410 if(_verbose && _nbTimeStep%_freqSave ==0)
2411 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
2413 /******** Construction de la matrice A^+ *********/
2414 //premiere ligne (masse)
2415 _AroePlusImplicit[0]=ui_n*gamma/((gamma-1)*(hi-q));
2416 for(int idim=0; idim<_Ndim;idim++)
2417 _AroePlusImplicit[1+idim]=rhoi*_vec_normal[idim];
2418 _AroePlusImplicit[_nVar-1]=-rhoi*ui_n*cp/(hi-q);
2420 //lignes intermadiaires(qdm)
2421 for(int idim=0; idim<_Ndim;idim++)
2424 _AroePlusImplicit[(1+idim)*_nVar]=ui_n*gamma/((gamma-1)*(hi-q))*_Vi[1+idim];
2425 //colonnes intermediaires
2426 for(int jdim=0; jdim<_Ndim;jdim++)
2427 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = rhoi*_Vi[1+idim]*_vec_normal[jdim];
2429 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += rhoi*ui_n;
2431 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoi*ui_n*cp/(hi-q)*_Vi[1+idim];
2434 //derniere ligne (energie)
2435 _AroePlusImplicit[_nVar*(_nVar-1)] = ui_n*(Hi*gamma/((gamma-1)*(hi-q))-1);
2436 for(int idim=0; idim<_Ndim;idim++)
2437 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoi*Ei+pj)*_vec_normal[idim]+rhoi*ui_n*_Vi[1+idim];
2438 _AroePlusImplicit[_nVar*_nVar -1] = rhoi*ui_n*(1-Hi/(hi-q))*cp;
2440 /******** Construction de la matrice A^- *********/
2441 //premiere ligne (masse)
2442 _AroeMinusImplicit[0]=0;
2443 for(int idim=0; idim<_Ndim;idim++)
2444 _AroeMinusImplicit[1+idim]=0;
2445 _AroeMinusImplicit[_nVar-1]=0;
2447 //lignes intermadiaires(qdm)
2448 for(int idim=0; idim<_Ndim;idim++)
2451 _AroeMinusImplicit[(1+idim)*_nVar]=_vec_normal[idim] ;
2452 //colonnes intermediaires
2453 for(int jdim=0; jdim<_Ndim;jdim++)
2454 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2456 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2458 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2461 //derniere ligne (energie)
2462 _AroeMinusImplicit[_nVar*(_nVar-1)] = ui_n;
2463 for(int idim=0; idim<_Ndim;idim++)
2464 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2465 _AroeMinusImplicit[_nVar*_nVar -1] = 0;
2467 else if(un<-_precision)
2469 double pi, Ej, rhoi,rhoj, Hj, hj;
2470 double ci;//vitesse du son pour calcul valeurs propres
2471 rhoj=_Uj[0]/_porosityj;
2472 Ej= _Uj[_Ndim+1]/rhoj;
2476 rhoi=_Ui[0]/_porosityi;
2477 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
2479 /***********Calcul des valeurs propres ********/
2480 vector<complex<double> > vp_dist(3,0);
2481 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2483 _maxvploc=fabs(uj_n)+ci;
2484 if(_maxvploc>_maxvp)
2487 if(_verbose && _nbTimeStep%_freqSave ==0)
2488 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2490 /******** Construction de la matrice A^+ *********/
2491 //premiere ligne (masse)
2492 _AroePlusImplicit[0]=0;
2493 for(int idim=0; idim<_Ndim;idim++)
2494 _AroePlusImplicit[1+idim]=0;
2495 _AroePlusImplicit[_nVar-1]=0;
2497 //lignes intermadiaires(qdm)
2498 for(int idim=0; idim<_Ndim;idim++)
2501 _AroePlusImplicit[(1+idim)*_nVar]=0;
2502 //colonnes intermediaires
2503 for(int jdim=0; jdim<_Ndim;jdim++)
2504 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2506 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2508 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2511 //derniere ligne (energie)
2512 _AroePlusImplicit[_nVar*(_nVar-1)] = uj_n;
2513 for(int idim=0; idim<_Ndim;idim++)
2514 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2515 _AroePlusImplicit[_nVar*_nVar -1] = 0;
2517 /******** Construction de la matrice A^- *********/
2518 //premiere ligne (masse)
2519 _AroeMinusImplicit[0]=uj_n*gamma/((gamma-1)*(hj-q));
2520 for(int idim=0; idim<_Ndim;idim++)
2521 _AroeMinusImplicit[1+idim]=rhoj*_vec_normal[idim];
2522 _AroeMinusImplicit[_nVar-1]=-rhoj*uj_n*cp/(hj-q);
2524 //lignes intermadiaires(qdm)
2525 for(int idim=0; idim<_Ndim;idim++)
2528 _AroeMinusImplicit[(1+idim)*_nVar]= uj_n*gamma/((gamma-1)*(hj-q))*_Vj[1+idim];
2529 //colonnes intermediaires
2530 for(int jdim=0; jdim<_Ndim;jdim++)
2531 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = rhoj*_Vj[1+idim]*_vec_normal[jdim];
2533 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += rhoj*uj_n;
2535 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoj*uj_n*cp/(hj-q)*_Vj[1+idim];
2538 //derniere ligne (energie)
2539 _AroeMinusImplicit[_nVar*(_nVar-1)] = uj_n*(Hj*gamma/((gamma-1)*(hj-q))-1);
2540 for(int idim=0; idim<_Ndim;idim++)
2541 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoj*Ej+pi)*_vec_normal[idim]+rhoj*uj_n*_Vj[1+idim];
2542 _AroeMinusImplicit[_nVar*_nVar -1] = rhoj*uj_n*(1-Hj/(hj-q))*cp;
2546 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero" << endl;
2547 _runLogFile->close();
2548 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero");
2553 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
2554 _runLogFile->close();
2555 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
2558 else//case nil velocity on the interface, apply centered scheme
2561 primToConsJacobianMatrix(_Vj);
2562 Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
2563 primToConsJacobianMatrix(_Vi);
2564 Poly.matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
2568 void SinglePhase::applyVFRoeLowMachCorrections(bool isBord, string groupname)
2570 if(_nonLinearFormulation!=VFRoe)
2572 *_runLogFile<< "SinglePhase::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation" << endl;
2573 _runLogFile->close();
2574 throw CdmathException("SinglePhase::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation");
2576 else//_nonLinearFormulation==VFRoe
2578 if(_spaceScheme==lowMach){
2580 for(int i=0;i<_Ndim;i++)
2581 u_2 += _Uroe[1+i]*_Uroe[1+i];
2582 double c = _maxvploc;//vitesse du son a l'interface
2583 double M=sqrt(u_2)/c;//Mach number
2584 _Vij[0]=M*_Vij[0]+(1-M)*(_Vi[0]+_Vj[0])/2;
2586 else if(_spaceScheme==pressureCorrection)
2587 {//order 1 : no correction, oarder 2 : correction everywhere, order 3 : correction only inside, orders 4 and 5 : special correction at boundaries
2588 if(_pressureCorrectionOrder==2 || (!isBord && _pressureCorrectionOrder==3) || (!isBord && _pressureCorrectionOrder==4) || (!isBord && _pressureCorrectionOrder==5) )
2590 double norm_uij=0, uij_n=0, ui_n=0, uj_n=0;
2591 for(int i=0;i<_Ndim;i++)
2593 norm_uij += _Uroe[1+i]*_Uroe[1+i];
2594 uij_n += _Uroe[1+i]*_vec_normal[i];
2595 ui_n += _Vi[1+i]*_vec_normal[i];
2596 uj_n += _Vj[1+i]*_vec_normal[i];
2598 norm_uij=sqrt(norm_uij);
2599 if(norm_uij>_precision)//avoid division by zero
2600 _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;
2602 _Vij[0]=(_Vi[0]+_Vj[0])/2 - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
2604 else if(_pressureCorrectionOrder==4 && isBord)
2606 else if(_pressureCorrectionOrder==5 && isBord)
2608 double g_n=0;//scalar product of gravity and normal vector
2609 for(int i=0;i<_Ndim;i++)
2610 g_n += _GravityField3d[i]*_vec_normal[i];
2611 _Vij[0]=_Vi[0]- _Ui[0]*g_n/_inv_dxi/2;
2614 else if(_spaceScheme==staggered)
2617 for(int i=0;i<_Ndim;i++)
2618 uij_n += _Uroe[1+i]*_vec_normal[i];
2620 if(uij_n>_precision){
2622 for(int i=0;i<_Ndim;i++)
2624 _Vij[_nVar-1]=_Vi[_nVar-1];
2626 else if(uij_n<-_precision){
2628 for(int i=0;i<_Ndim;i++)
2630 _Vij[_nVar-1]=_Vj[_nVar-1];
2633 _Vij[0]=(_Vi[0]+_Vi[0])/2;
2634 for(int i=0;i<_Ndim;i++)
2635 _Vij[1+i]=(_Vj[1+i]+_Vj[1+i])/2;
2636 _Vij[_nVar-1]=(_Vj[_nVar-1]+_Vj[_nVar-1])/2;
2639 primToCons(_Vij,0,_Uij,0);
2643 void SinglePhase::testConservation()
2645 double SUM, DELTA, x;
2647 for(int i=0; i<_nVar; i++)
2651 cout << "Masse totale (kg): ";
2655 cout << "Energie totale (J): ";
2657 cout << "Quantite de mouvement totale (kg.m.s^-1): ";
2663 for(int j=0; j<_Nmailles; j++)
2665 if(!_usePrimitiveVarsInNewton)
2666 VecGetValues(_conservativeVars, 1, &I, &x);//on recupere la valeur du champ
2668 VecGetValues(_primitiveVars, 1, &I, &x);//on recupere la valeur du champ
2669 SUM += x*_mesh.getCell(j).getMeasure();
2670 VecGetValues(_newtonVariation, 1, &I, &x);//on recupere la variation du champ
2671 DELTA += x*_mesh.getCell(j).getMeasure();
2674 if(fabs(SUM)>_precision)
2675 cout << SUM << ", variation relative: " << fabs(DELTA /SUM) << endl;
2677 cout << " a une somme quasi nulle, variation absolue: " << fabs(DELTA) << endl;
2681 void SinglePhase::getDensityDerivatives( double pressure, double temperature, double v2)
2683 double rho=_fluides[0]->getDensity(pressure,temperature);
2684 double gamma=_fluides[0]->constante("gamma");
2685 double q=_fluides[0]->constante("q");
2687 if( !_useDellacherieEOS)
2689 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2690 double e = fluide0->getInternalEnergy(temperature);
2691 double cv=fluide0->constante("cv");
2694 _drho_sur_dp=1/((gamma-1)*(e-q));
2695 _drho_sur_dT=-rho*cv/(e-q);
2696 _drhoE_sur_dp=E/((gamma-1)*(e-q));
2697 _drhoE_sur_dT=rho*cv*(1-E/(e-q));
2699 else if(_useDellacherieEOS )
2701 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2702 double h=fluide0->getEnthalpy(temperature);
2704 double cp=fluide0->constante("cp");
2706 _drho_sur_dp=gamma/((gamma-1)*(h-q));
2707 _drho_sur_dT=-rho*cp/(h-q);
2708 _drhoE_sur_dp=gamma*H/((gamma-1)*(h-q))-1;
2709 _drhoE_sur_dT=rho*cp*(1-H/(h-q));
2713 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
2714 _runLogFile->close();
2715 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
2718 if(_verbose && _nbTimeStep%_freqSave ==0)
2720 cout<<"_drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
2721 cout<<"_drhoE_sur_dp= "<<_drhoE_sur_dp<<", _drhoE_sur_dT= "<<_drhoE_sur_dT<<endl;
2724 void SinglePhase::save(){
2725 string prim(_path+"/SinglePhasePrim_");///Results
2726 string cons(_path+"/SinglePhaseCons_");
2727 string allFields(_path+"/");
2730 allFields+=_fileName;
2733 for (long i = 0; i < _Nmailles; i++){
2734 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
2735 for (int j = 0; j < _nVar; j++){
2737 VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
2740 if(_saveConservativeField){
2741 for (long i = 0; i < _Nmailles; i++){
2742 // j = 0 : density; j = _nVar - 1 : energy j = 1,..,_nVar-2: momentum
2743 for (int j = 0; j < _nVar; j++){
2745 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
2748 _UU.setTime(_time,_nbTimeStep);
2750 _VV.setTime(_time,_nbTimeStep);
2752 // create mesh and component info
2753 if (_nbTimeStep ==0 || _restartWithNewFileName){
2754 string prim_suppress ="rm -rf "+prim+"_*";
2755 string cons_suppress ="rm -rf "+cons+"_*";
2757 system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
2758 system(cons_suppress.c_str());//Nettoyage des précédents calculs identiques
2760 if(_saveConservativeField){
2761 _UU.setInfoOnComponent(0,"Density_(kg/m^3)");
2762 _UU.setInfoOnComponent(1,"Momentum_x");// (kg/m^2/s)
2764 _UU.setInfoOnComponent(2,"Momentum_y");// (kg/m^2/s)
2766 _UU.setInfoOnComponent(3,"Momentum_z");// (kg/m^2/s)
2768 _UU.setInfoOnComponent(_nVar-1,"Energy_(J/m^3)");
2783 _VV.setInfoOnComponent(0,"Pressure_(Pa)");
2784 _VV.setInfoOnComponent(1,"Velocity_x_(m/s)");
2786 _VV.setInfoOnComponent(2,"Velocity_y_(m/s)");
2788 _VV.setInfoOnComponent(3,"Velocity_z_(m/s)");
2789 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
2805 // do not create mesh
2810 _VV.writeVTK(prim,false);
2813 _VV.writeMED(prim,false);
2819 if(_saveConservativeField){
2823 _UU.writeVTK(cons,false);
2826 _UU.writeMED(cons,false);
2834 if(_saveVelocity || _saveAllFields){
2835 for (long i = 0; i < _Nmailles; i++){
2836 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
2837 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
2839 int Ii = i*_nVar +1+j;
2840 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
2842 for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
2845 _Vitesse.setTime(_time,_nbTimeStep);
2846 if (_nbTimeStep ==0 || _restartWithNewFileName){
2847 _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
2848 _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
2849 _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
2854 _Vitesse.writeVTK(prim+"_Velocity");
2857 _Vitesse.writeMED(prim+"_Velocity");
2860 _Vitesse.writeCSV(prim+"_Velocity");
2868 _Vitesse.writeVTK(prim+"_Velocity",false);
2871 _Vitesse.writeMED(prim+"_Velocity",false);
2874 _Vitesse.writeCSV(prim+"_Velocity");
2882 double p,T,rho, h, vx,vy,vz;
2884 for (long i = 0; i < _Nmailles; i++){
2886 VecGetValues(_conservativeVars,1,&Ii,&rho);
2888 VecGetValues(_primitiveVars,1,&Ii,&p);
2889 Ii = i*_nVar +_nVar-1;
2890 VecGetValues(_primitiveVars,1,&Ii,&T);
2892 VecGetValues(_primitiveVars,1,&Ii,&vx);
2896 VecGetValues(_primitiveVars,1,&Ii,&vy);
2899 VecGetValues(_primitiveVars,1,&Ii,&vz);
2903 h = _fluides[0]->getEnthalpy(T,rho);
2917 _Enthalpy.setTime(_time,_nbTimeStep);
2918 _Density.setTime(_time,_nbTimeStep);
2919 _Pressure.setTime(_time,_nbTimeStep);
2920 _Temperature.setTime(_time,_nbTimeStep);
2921 _VitesseX.setTime(_time,_nbTimeStep);
2924 _VitesseY.setTime(_time,_nbTimeStep);
2926 _VitesseZ.setTime(_time,_nbTimeStep);
2928 if (_nbTimeStep ==0 || _restartWithNewFileName){
2932 _Enthalpy.writeVTK(allFields+"_Enthalpy");
2933 _Density.writeVTK(allFields+"_Density");
2934 _Pressure.writeVTK(allFields+"_Pressure");
2935 _Temperature.writeVTK(allFields+"_Temperature");
2936 _VitesseX.writeVTK(allFields+"_VelocityX");
2939 _VitesseY.writeVTK(allFields+"_VelocityY");
2941 _VitesseZ.writeVTK(allFields+"_VelocityZ");
2945 _Enthalpy.writeMED(allFields+"_Enthalpy");
2946 _Density.writeMED(allFields+"_Density");
2947 _Pressure.writeMED(allFields+"_Pressure");
2948 _Temperature.writeMED(allFields+"_Temperature");
2949 _VitesseX.writeMED(allFields+"_VelocityX");
2952 _VitesseY.writeMED(allFields+"_VelocityY");
2954 _VitesseZ.writeMED(allFields+"_VelocityZ");
2958 _Enthalpy.writeCSV(allFields+"_Enthalpy");
2959 _Density.writeCSV(allFields+"_Density");
2960 _Pressure.writeCSV(allFields+"_Pressure");
2961 _Temperature.writeCSV(allFields+"_Temperature");
2962 _VitesseX.writeCSV(allFields+"_VelocityX");
2965 _VitesseY.writeCSV(allFields+"_VelocityY");
2967 _VitesseZ.writeCSV(allFields+"_VelocityZ");
2976 _Enthalpy.writeVTK(allFields+"_Enthalpy",false);
2977 _Density.writeVTK(allFields+"_Density",false);
2978 _Pressure.writeVTK(allFields+"_Pressure",false);
2979 _Temperature.writeVTK(allFields+"_Temperature",false);
2980 _VitesseX.writeVTK(allFields+"_VelocityX",false);
2983 _VitesseY.writeVTK(allFields+"_VelocityY",false);
2985 _VitesseZ.writeVTK(allFields+"_VelocityZ",false);
2989 _Enthalpy.writeMED(allFields+"_Enthalpy",false);
2990 _Density.writeMED(allFields+"_Density",false);
2991 _Pressure.writeMED(allFields+"_Pressure",false);
2992 _Temperature.writeMED(allFields+"_Temperature",false);
2993 _VitesseX.writeMED(allFields+"_VelocityX",false);
2996 _VitesseY.writeMED(allFields+"_VelocityY",false);
2998 _VitesseZ.writeMED(allFields+"_VelocityZ",false);
3002 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3003 _Density.writeCSV(allFields+"_Density");
3004 _Pressure.writeCSV(allFields+"_Pressure");
3005 _Temperature.writeCSV(allFields+"_Temperature");
3006 _VitesseX.writeCSV(allFields+"_VelocityX");
3009 _VitesseY.writeCSV(allFields+"_VelocityY");
3011 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3036 if(_saveConservativeField){
3051 if(_saveVelocity || _saveAllFields){
3055 _Vitesse.writeVTK(prim+"_Velocity");
3058 _Vitesse.writeMED(prim+"_Velocity");
3061 _Vitesse.writeCSV(prim+"_Velocity");
3067 if (_restartWithNewFileName)
3068 _restartWithNewFileName=false;
3071 Field& SinglePhase::getPressureField()
3075 _Pressure=Field("Pressure",CELLS,_mesh,1);
3077 for (long i = 0; i < _Nmailles; i++){
3079 VecGetValues(_primitiveVars,1,&Ii,&_Pressure(i));
3081 _Pressure.setTime(_time,_nbTimeStep);
3086 Field& SinglePhase::getTemperatureField()
3090 _Temperature=Field("Temperature",CELLS,_mesh,1);
3092 for (long i = 0; i < _Nmailles; i++){
3093 Ii = i*_nVar +_nVar-1;
3094 VecGetValues(_primitiveVars,1,&Ii,&_Temperature(i));
3096 _Temperature.setTime(_time,_nbTimeStep);
3098 return _Temperature;
3101 Field& SinglePhase::getVelocityField()
3103 if(!_saveAllFields )
3105 _Vitesse=Field("Vitesse",CELLS,_mesh,3);
3107 for (long i = 0; i < _Nmailles; i++)
3109 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
3111 int Ii = i*_nVar +1+j;
3112 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
3114 for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
3117 _Vitesse.setTime(_time,_nbTimeStep);
3118 _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
3119 _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
3120 _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
3126 Field& SinglePhase::getVelocityXField()
3128 if(!_saveAllFields )
3130 _VitesseX=Field("Velocity X",CELLS,_mesh,1);
3132 for (long i = 0; i < _Nmailles; i++)
3134 int Ii = i*_nVar +1;
3135 VecGetValues(_primitiveVars,1,&Ii,&_VitesseX(i));
3137 _VitesseX.setTime(_time,_nbTimeStep);
3138 _VitesseX.setInfoOnComponent(0,"Velocity_x_(m/s)");
3144 Field& SinglePhase::getDensityField()
3146 if(!_saveAllFields )
3148 _Density=Field("Density",CELLS,_mesh,1);
3150 for (long i = 0; i < _Nmailles; i++){
3152 VecGetValues(_conservativeVars,1,&Ii,&_Density(i));
3154 _Density.setTime(_time,_nbTimeStep);
3159 Field& SinglePhase::getMomentumField()//not yet managed by parameter _saveAllFields
3161 _Momentum=Field("Momentum",CELLS,_mesh,_Ndim);
3163 for (long i = 0; i < _Nmailles; i++)
3164 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de qdm
3166 int Ii = i*_nVar +1+j;
3167 VecGetValues(_conservativeVars,1,&Ii,&_Momentum(i,j));
3169 _Momentum.setTime(_time,_nbTimeStep);
3174 Field& SinglePhase::getTotalEnergyField()//not yet managed by parameter _saveAllFields
3176 _TotalEnergy=Field("TotalEnergy",CELLS,_mesh,1);
3178 for (long i = 0; i < _Nmailles; i++){
3179 Ii = i*_nVar +_nVar-1;
3180 VecGetValues(_conservativeVars,1,&Ii,&_TotalEnergy(i));
3182 _TotalEnergy.setTime(_time,_nbTimeStep);
3184 return _TotalEnergy;
3187 Field& SinglePhase::getEnthalpyField()
3189 if(!_saveAllFields )
3191 _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
3194 for (long i = 0; i < _Nmailles; i++){
3196 VecGetValues(_primitiveVars,1,&Ii,&p);
3197 Ii = i*_nVar +_nVar-1;
3198 VecGetValues(_primitiveVars,1,&Ii,&T);
3200 rho=_fluides[0]->getDensity(p,T);
3201 _Enthalpy(i)=_fluides[0]->getEnthalpy(T,rho);
3203 _Enthalpy.setTime(_time,_nbTimeStep);
3209 vector<string> SinglePhase::getOutputFieldsNames()
3211 vector<string> result(8);
3213 result[0]="Pressure";
3214 result[1]="Velocity";
3215 result[2]="Temperature";
3216 result[3]="Density";
3217 result[4]="Momentum";
3218 result[5]="TotalEnergy";
3219 result[6]="Enthalpy";
3220 result[7]="VelocityX";
3225 Field& SinglePhase::getOutputField(const string& nameField )
3227 if(nameField=="pressure" || nameField=="Pressure" || nameField=="PRESSURE" || nameField=="PRESSION" || nameField=="Pression" || nameField=="pression" )
3228 return getPressureField();
3229 else if(nameField=="velocity" || nameField=="Velocity" || nameField=="VELOCITY" || nameField=="Vitesse" || nameField=="VITESSE" || nameField=="vitesse" )
3230 return getVelocityField();
3231 else if(nameField=="velocityX" || nameField=="VelocityX" || nameField=="VELOCITYX" || nameField=="VitesseX" || nameField=="VITESSEX" || nameField=="vitesseX" )
3232 return getVelocityXField();
3233 else if(nameField=="temperature" || nameField=="Temperature" || nameField=="TEMPERATURE" || nameField=="temperature" )
3234 return getTemperatureField();
3235 else if(nameField=="density" || nameField=="Density" || nameField=="DENSITY" || nameField=="Densite" || nameField=="DENSITE" || nameField=="densite" )
3236 return getDensityField();
3237 else if(nameField=="momentum" || nameField=="Momentum" || nameField=="MOMENTUM" || nameField=="Qdm" || nameField=="QDM" || nameField=="qdm" )
3238 return getMomentumField();
3239 else if(nameField=="enthalpy" || nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPIE" || nameField=="enthalpie" )
3240 return getEnthalpyField();
3241 else if(nameField=="totalenergy" || nameField=="TotalEnergy" || nameField=="TOTALENERGY" || nameField=="ENERGIETOTALE" || nameField=="EnergieTotale" || nameField=="energietotale" )
3242 return getTotalEnergyField();
3245 cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
3246 throw CdmathException("SinglePhase::getOutputField error : Unknown Field name");