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){
313 //First conservative state then further down we will compute interface (Roe) state and then compute primitive state
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");
343 //Computation of Roe state
344 PetscScalar ri, rj, xi, xj, pi, pj;
346 ri = sqrt(_Ui[0]);//racine carre de phi_i rho_i
348 _Uroe[0] = ri*rj; //moyenne geometrique des densites
349 if(_verbose && _nbTimeStep%_freqSave ==0)
350 cout << "Densité moyenne Roe gauche " << i << ": " << ri*ri << ", droite " << j << ": " << rj*rj << "->" << _Uroe[0] << endl;
351 for(int k=0;k<_Ndim;k++){
354 _Uroe[1+k] = (xi/ri + xj/rj)/(ri + rj);
355 //"moyenne" des vitesses
356 if(_verbose && _nbTimeStep%_freqSave ==0)
357 cout << "Vitesse de Roe composante "<< k<<" gauche " << i << ": " << xi/(ri*ri) << ", droite " << j << ": " << xj/(rj*rj) << "->" << _Uroe[k+1] << endl;
359 // H = (rho E + p)/rho
360 xi = _Ui[_nVar-1];//phi rho E
362 Ii = i*_nVar; // correct Kieu
363 VecGetValues(_primitiveVars, 1, &Ii, &pi);// _primitiveVars pour p
367 for(int k=1;k<=_Ndim;k++)
368 q_2 += _Uj[k]*_Uj[k];
369 q_2 /= _Uj[0]; //phi rho u²
370 pj = _fluides[0]->getPressure((_Uj[(_Ndim+2)-1]-q_2/2)/_porosityj,_Uj[0]/_porosityj);
374 Ii = j*_nVar; // correct Kieu
375 VecGetValues(_primitiveVars, 1, &Ii, &pj);
377 xi = (xi + pi)/(ri*ri);
378 xj = (xj + pj)/(rj*rj);
379 _Uroe[_nVar-1] = (ri*xi + rj*xj)/(ri + rj);
380 //on se donne l enthalpie ici
381 if(_verbose && _nbTimeStep%_freqSave ==0)
382 cout << "Enthalpie totale de Roe H gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[_nVar-1] << endl;
384 if(_verbose && _nbTimeStep%_freqSave ==0)
386 cout<<"Convection interfacial state"<<endl;
387 for(int k=0;k<_nVar;k++)
388 cout<< _Uroe[k]<<" , "<<endl;
391 //Extraction of primitive states
392 _idm[0] = _nVar*i; // Kieu
393 for(int k=1; k<_nVar; k++)
394 _idm[k] = _idm[k-1] + 1;
396 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
397 if (_verbose && _nbTimeStep%_freqSave ==0)
399 cout << "Convection state: variables primitives maille " << i<<endl;
400 for(int q=0; q<_nVar; q++)
401 cout << _Vi[q] << endl;
406 for(int k=0; k<_nVar; k++)
407 _idn[k] = _nVar*j + k;
409 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
413 for(int k=0; k<_nVar; k++)
416 VecGetValues(_Uextdiff, _nVar, _idn, _phi);
417 consToPrim(_phi,_Vj,1);
420 if (_verbose && _nbTimeStep%_freqSave ==0)
422 cout << "Convection state: variables primitives maille " <<j <<endl;
423 for(int q=0; q<_nVar; q++)
424 cout << _Vj[q] << endl;
429 void SinglePhase::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
430 //sortie: matrices et etat de diffusion (rho, q, T)
432 for(int k=1; k<_nVar; k++)
433 _idm[k] = _idm[k-1] + 1;
435 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
436 for(int k=0; k<_nVar; k++)
440 VecGetValues(_Uextdiff, _nVar, _idn, _Uj);
442 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
444 if(_verbose && _nbTimeStep%_freqSave ==0)
446 cout << "SinglePhase::diffusionStateAndMatrices cellule gauche" << i << endl;
448 for(int q=0; q<_nVar; q++)
449 cout << _Ui[q] << "\t";
451 cout << "SinglePhase::diffusionStateAndMatrices cellule droite" << j << endl;
453 for(int q=0; q<_nVar; q++)
454 cout << _Uj[q] << "\t";
458 for(int k=0; k<_nVar; k++)
459 _Udiff[k] = (_Ui[k]/_porosityi+_Uj[k]/_porosityj)/2;
461 if(_verbose && _nbTimeStep%_freqSave ==0)
463 cout << "SinglePhase::diffusionStateAndMatrices conservative diffusion state" << endl;
465 for(int q=0; q<_nVar; q++)
466 cout << _Udiff[q] << "\t";
468 cout << "porosite gauche= "<<_porosityi<< ", porosite droite= "<<_porosityj<<endl;
470 consToPrim(_Udiff,_phi,1);
471 _Udiff[_nVar-1]=_phi[_nVar-1];
473 if(_timeScheme==Implicit)
476 for (int i = 0; i<_Ndim;i++)
477 q_2+=_Udiff[i+1]*_Udiff[i+1];
478 double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
479 double lambda = _fluides[0]->getConductivity(_Udiff[_nVar-1]);
480 double Cv= _fluides[0]->constante("Cv");
481 for(int i=0; i<_nVar*_nVar;i++)
483 for(int i=1;i<(_nVar-1);i++)
485 _Diffusion[i*_nVar] = mu*_Udiff[i]/(_Udiff[0]*_Udiff[0]);
486 _Diffusion[i*_nVar+i] = -mu/_Udiff[0];
488 int i = (_nVar-1)*_nVar;
489 _Diffusion[i]=lambda*(_Udiff[_nVar-1]/_Udiff[0]-q_2/(2*Cv*_Udiff[0]*_Udiff[0]*_Udiff[0]));
490 for(int k=1;k<(_nVar-1);k++)
492 _Diffusion[i+k]= lambda*_Udiff[k]/(_Udiff[0]*_Udiff[0]*Cv);
494 _Diffusion[_nVar*_nVar-1]=-lambda/(_Udiff[0]*Cv);
497 void SinglePhase::setBoundaryState(string nameOfGroup, const int &j,double *normale){
499 for(int k=1; k<_nVar; k++)
500 _idm[k] = _idm[k-1] + 1;
502 double porosityj=_porosityField(j);
504 if(_verbose && _nbTimeStep%_freqSave ==0)
506 cout << "setBoundaryState for group "<< nameOfGroup << ", inner cell j= "<<j<< " face unit normal vector "<<endl;
507 for(int k=0; k<_Ndim; k++){
508 cout<<normale[k]<<", ";
513 if (_limitField[nameOfGroup].bcType==Wall){
514 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne conservatif
515 double q_n=0;//q_n=quantité de mouvement normale à la face frontière;
516 for(int k=0; k<_Ndim; k++)
517 q_n+=_externalStates[(k+1)]*normale[k];
519 //Pour la convection, inversion du sens de la vitesse normale
520 for(int k=0; k<_Ndim; k++)
521 _externalStates[(k+1)]-= 2*q_n*normale[k];
524 for(int k=1; k<_nVar; k++)
525 _idm[k] = _idm[k-1] + 1;
527 VecAssemblyBegin(_Uext);
528 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
529 VecAssemblyEnd(_Uext);
531 //Pour la diffusion, paroi à vitesse et temperature imposees
533 for(int k=1; k<_nVar; k++)
534 _idm[k] = _idm[k-1] + 1;
535 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);//L'état fantome contient à présent les variables primitives internes
536 double pression=_externalStates[0];
537 double T=_limitField[nameOfGroup].T;
538 double rho=_fluides[0]->getDensity(pression,T);
540 _externalStates[0]=porosityj*rho;
541 _externalStates[1]=_externalStates[0]*_limitField[nameOfGroup].v_x[0];
543 v2 +=_limitField[nameOfGroup].v_x[0]*_limitField[nameOfGroup].v_x[0];
546 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
547 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];
550 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];
551 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
554 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
556 for(int k=1; k<_nVar; k++)
557 _idm[k] = _idm[k-1] + 1;
558 VecAssemblyBegin(_Uextdiff);
559 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
560 VecAssemblyEnd(_Uextdiff);
562 else if (_limitField[nameOfGroup].bcType==Neumann){
563 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On prend l'état fantôme égal à l'état interne (conditions limites de Neumann)
566 for(int k=1; k<_nVar; k++)
567 _idm[k] = _idm[k-1] + 1;
569 VecAssemblyBegin(_Uext);
570 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
571 VecAssemblyEnd(_Uext);
573 VecAssemblyBegin(_Uextdiff);
574 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
575 VecAssemblyEnd(_Uextdiff);
577 else if (_limitField[nameOfGroup].bcType==Inlet){
578 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne (conditions limites de Neumann)
579 double q_int_n=0;//q_int_n=composante normale de la quantité de mouvement à la face frontière;
580 for(int k=0; k<_Ndim; k++)
581 q_int_n+=_externalStates[(k+1)]*normale[k];//On calcule la vitesse normale sortante
583 double q_ext_n=_limitField[nameOfGroup].v_x[0]*normale[0];
586 q_ext_n+=_limitField[nameOfGroup].v_y[0]*normale[1];
588 q_ext_n+=_limitField[nameOfGroup].v_z[0]*normale[2];
591 if(q_int_n+q_ext_n<=0){//Interfacial velocity goes inward
592 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);//On met à jour l'état fantome avec les variables primitives internes
593 double pression=_externalStates[0];
594 double T=_limitField[nameOfGroup].T;
595 double rho=_fluides[0]->getDensity(pression,T);
597 _externalStates[0]=porosityj*rho;//Composante fantome de masse
598 _externalStates[1]=_externalStates[0]*(_limitField[nameOfGroup].v_x[0]);//Composante fantome de qdm x
600 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
603 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
604 _externalStates[2]=_externalStates[0]*_limitField[nameOfGroup].v_y[0];//Composante fantome de qdm y
607 _externalStates[3]=_externalStates[0]*_limitField[nameOfGroup].v_z[0];//Composante fantome de qdm z
608 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
611 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);//Composante fantome de l'nrj
613 else if(_nbTimeStep%_freqSave ==0)
614 cout<< "Warning : fluid possibly going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
617 for(int k=1; k<_nVar; k++)
618 _idm[k] = _idm[k-1] + 1;
619 VecAssemblyBegin(_Uext);
620 VecAssemblyBegin(_Uextdiff);
621 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
622 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
623 VecAssemblyEnd(_Uext);
624 VecAssemblyEnd(_Uextdiff);
626 else if (_limitField[nameOfGroup].bcType==InletRotationVelocity){
627 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
628 double u_int_n=0;//u_int_n=composante normale de la vitesse intérieure à la face frontière;
629 for(int k=0; k<_Ndim; k++)
630 u_int_n+=_externalStates[(k+1)]*normale[k];
637 omega[0]=_limitField[nameOfGroup].v_x[0];
638 omega[1]=_limitField[nameOfGroup].v_y[0];
641 Normale[0]=normale[0];
642 Normale[1]=normale[1];
646 omega[2]=_limitField[nameOfGroup].v_z[0];
647 Normale[2]=normale[2];
650 Vector tangent_vel=omega%Normale;
651 u_ext_n=-0.01*tangent_vel.norm();
652 //Changing external state velocity
653 for(int k=0; k<_Ndim; k++)
654 _externalStates[(k+1)]=u_ext_n*normale[k] + tangent_vel[k];
657 if(u_ext_n + u_int_n <= 0){
658 double pression=_externalStates[0];
659 double T=_limitField[nameOfGroup].T;
660 double rho=_fluides[0]->getDensity(pression,T);
663 v2 +=_externalStates[1]*_externalStates[1];//v_x*v_x
664 _externalStates[0]=porosityj*rho;
665 _externalStates[1]*=_externalStates[0];
668 v2 +=_externalStates[2]*_externalStates[2];//+v_y*v_y
669 _externalStates[2]*=_externalStates[0];
672 v2 +=_externalStates[3]*_externalStates[3];//+v_z*v_z
673 _externalStates[3]*=_externalStates[0];
676 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2);
678 else if(_nbTimeStep%_freqSave ==0)
681 * cout<< "Warning : fluid going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
683 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On définit l'état fantôme avec l'état interne
684 if(_nbTimeStep%_freqSave ==0)
685 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Wall boundary condition."<<endl;
687 //Changing external state momentum
688 for(int k=0; k<_Ndim; k++)
689 _externalStates[(k+1)]-=2*_externalStates[0]*u_int_n*normale[k];
693 for(int k=1; k<_nVar; k++)
694 _idm[k] = _idm[k-1] + 1;
695 VecAssemblyBegin(_Uext);
696 VecAssemblyBegin(_Uextdiff);
697 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
698 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
699 VecAssemblyEnd(_Uext);
700 VecAssemblyEnd(_Uextdiff);
702 else if (_limitField[nameOfGroup].bcType==InletPressure){
703 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
705 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
706 Cell Cj=_mesh.getCell(j);
707 double hydroPress=Cj.x()*_GravityField3d[0];
709 hydroPress+=Cj.y()*_GravityField3d[1];
711 hydroPress+=Cj.z()*_GravityField3d[2];
713 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
715 //Building the primitive external state
716 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
717 double u_n=0;//u_n=vitesse normale à la face frontière;
718 for(int k=0; k<_Ndim; k++)
719 u_n+=_externalStates[(k+1)]*normale[k];
723 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress,_limitField[nameOfGroup].T);
725 //Contribution from the tangential velocity
729 omega[0]=_limitField[nameOfGroup].v_x[0];
730 omega[1]=_limitField[nameOfGroup].v_y[0];
733 Normale[0]=normale[0];
734 Normale[1]=normale[1];
738 omega[2]=_limitField[nameOfGroup].v_z[0];
739 Normale[2]=normale[2];
742 Vector tangent_vel=omega%Normale;
744 //Changing external state velocity
745 for(int k=0; k<_Ndim; k++)
746 _externalStates[(k+1)]=u_n*normale[k] + abs(u_n)*tangent_vel[k];
751 if(_nbTimeStep%_freqSave ==0)
752 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;
753 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
755 if(_nbTimeStep%_freqSave ==0)
756 cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Wall boundary condition."<<endl;
757 _externalStates[0]=porosityj*_fluides[0]->getDensity(_externalStates[0]+hydroPress, _externalStates[_nVar-1]);
758 //Changing external state velocity
759 for(int k=0; k<_Ndim; k++)
760 _externalStates[(k+1)]-=2*u_n*normale[k];
764 for(int k=0; k<_Ndim; k++)
766 v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
767 _externalStates[(k+1)]*=_externalStates[0] ;//qdm component
769 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);//nrj component
773 for(int k=1; k<_nVar; k++)
774 _idm[k] = _idm[k-1] + 1;
775 VecAssemblyBegin(_Uext);
776 VecAssemblyBegin(_Uextdiff);
777 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
778 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
779 VecAssemblyEnd(_Uext);
780 VecAssemblyEnd(_Uextdiff);
782 else if (_limitField[nameOfGroup].bcType==Outlet){
783 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne conservatif
784 double q_n=0;//q_n=quantité de mouvement normale à la face frontière;
785 for(int k=0; k<_Ndim; k++)
786 q_n+=_externalStates[(k+1)]*normale[k];
788 if(q_n < -_precision && _nbTimeStep%_freqSave ==0)
790 cout<< "Warning : fluid going in through outlet boundary "<<nameOfGroup<<" with flow rate "<< q_n<<endl;
791 cout<< "Applying Neumann boundary condition for velocity and temperature"<<endl;
793 //Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
794 Cell Cj=_mesh.getCell(j);
795 double hydroPress=Cj.x()*_GravityField3d[0];
797 hydroPress+=Cj.y()*_GravityField3d[1];
799 hydroPress+=Cj.z()*_GravityField3d[2];
801 hydroPress*=_externalStates[0]/porosityj;//multiplication by rho the total density
803 if(_verbose && _nbTimeStep%_freqSave ==0)
805 cout<<"Cond lim outlet densite= "<<_externalStates[0]<<" gravite= "<<_GravityField3d[0]<<" Cj.x()= "<<Cj.x()<<endl;
806 cout<<"Cond lim outlet pression ref= "<<_limitField[nameOfGroup].p<<" pression hydro= "<<hydroPress<<" total= "<<_limitField[nameOfGroup].p+hydroPress<<endl;
808 //Building the external state
809 _idm[0] = _nVar*j;// Kieu
810 for(int k=1; k<_nVar; k++)
811 _idm[k] = _idm[k-1] + 1;
812 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
814 _externalStates[0]=porosityj*_fluides[0]->getDensity(_limitField[nameOfGroup].p+hydroPress, _externalStates[_nVar-1]);
816 for(int k=0; k<_Ndim; k++)
818 v2+=_externalStates[(k+1)]*_externalStates[(k+1)];
819 _externalStates[(k+1)]*=_externalStates[0] ;
821 _externalStates[_nVar-1] = _externalStates[0]*(_fluides[0]->getInternalEnergy( _externalStates[_nVar-1],_externalStates[0]) + v2/2);
823 for(int k=1; k<_nVar; k++)
824 _idm[k] = _idm[k-1] + 1;
825 VecAssemblyBegin(_Uext);
826 VecAssemblyBegin(_Uextdiff);
827 VecSetValues(_Uext, _nVar, _idm, _externalStates, INSERT_VALUES);
828 VecSetValues(_Uextdiff, _nVar, _idm, _externalStates, INSERT_VALUES);
829 VecAssemblyEnd(_Uext);
830 VecAssemblyEnd(_Uextdiff);
832 cout<<"Boundary condition not set for boundary named "<<nameOfGroup<< " _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
833 cout<<"Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
834 *_runLogFile<<"Boundary condition not set for boundary named. Accepted boundary condition are Neumann, Wall, Inlet, and Outlet"<<endl;
835 _runLogFile->close();
836 throw CdmathException("Unknown boundary condition");
840 void SinglePhase::convectionMatrices()
842 //entree: URoe = rho, u, H
843 //sortie: matrices Roe+ et Roe-
845 if(_verbose && _nbTimeStep%_freqSave ==0)
846 cout<<"SinglePhase::convectionMatrices()"<<endl;
848 double u_n=0, u_2=0;//vitesse normale et carré du module
850 for(int i=0;i<_Ndim;i++)
852 u_2 += _Uroe[1+i]*_Uroe[1+i];
853 u_n += _Uroe[1+i]*_vec_normal[i];
856 vector<complex<double> > vp_dist(3,0);
858 if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
860 staggeredVFFCMatricesConservativeVariables(u_n);//Computation of classical upwinding matrices
861 if(_timeScheme==Implicit && _usePrimitiveVarsInNewton)//For use in implicit matrix
862 staggeredVFFCMatricesPrimitiveVariables(u_n);
866 Vector vitesse(_Ndim);
867 for(int idim=0;idim<_Ndim;idim++)
868 vitesse[idim]=_Uroe[1+idim];
871 /***********Calcul des valeurs propres ********/
873 c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
874 k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
875 K = u_2*k/2; //g-1/2 *|u|²
877 vp_dist[0]=u_n-c;vp_dist[1]=u_n;vp_dist[2]=u_n+c;
879 _maxvploc=fabs(u_n)+c;
883 if(_verbose && _nbTimeStep%_freqSave ==0)
884 cout<<"SinglePhase::convectionMatrices Eigenvalues "<<u_n-c<<" , "<<u_n<<" , "<<u_n+c<<endl;
886 RoeMatrixConservativeVariables( u_n, H,vitesse,k,K);
888 /******** Construction des matrices de decentrement ********/
889 if( _spaceScheme ==centered){
890 if(_entropicCorrection)
892 *_runLogFile<<"SinglePhase::convectionMatrices: entropy scheme not available for centered scheme"<<endl;
893 _runLogFile->close();
894 throw CdmathException("SinglePhase::convectionMatrices: entropy scheme not available for centered scheme");
897 for(int i=0; i<_nVar*_nVar;i++)
900 else if(_spaceScheme == upwind || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach){
901 if(_entropicCorrection)
902 entropicShift(_vec_normal);
904 _entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
906 vector< complex< double > > y (3,0);
908 for( int i=0 ; i<3 ; i++)
909 y[i] = Poly.abs_generalise(vp_dist[i])+_entropicShift[i];
910 Poly.abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
912 if( _spaceScheme ==pressureCorrection)
913 for( int i=0 ; i<_Ndim ; i++)
914 for( int j=0 ; j<_Ndim ; j++)
915 _absAroe[(1+i)*_nVar+1+j]-=(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
916 else if( _spaceScheme ==lowMach){
917 double M=sqrt(u_2)/c;
918 for( int i=0 ; i<_Ndim ; i++)
919 for( int j=0 ; j<_Ndim ; j++)
920 _absAroe[(1+i)*_nVar+1+j]-=(1-M)*(vp_dist[2].real()-vp_dist[0].real())/2*_vec_normal[i]*_vec_normal[j];
923 else if( _spaceScheme ==staggered ){
924 if(_entropicCorrection)//To do: study entropic correction for staggered
926 *_runLogFile<<"SinglePhase::convectionMatrices: entropy scheme not available for staggered scheme"<<endl;
927 _runLogFile->close();
928 throw CdmathException("SinglePhase::convectionMatrices: entropy scheme not available for staggered scheme");
931 staggeredRoeUpwindingMatrixConservativeVariables( u_n, H, vitesse, k, K);
935 *_runLogFile<<"SinglePhase::convectionMatrices: scheme not treated"<<endl;
936 _runLogFile->close();
937 throw CdmathException("SinglePhase::convectionMatrices: scheme not treated");
940 for(int i=0; i<_nVar*_nVar;i++)
942 _AroeMinus[i] = (_Aroe[i]-_absAroe[i])/2;
943 _AroePlus[i] = (_Aroe[i]+_absAroe[i])/2;
945 if(_timeScheme==Implicit)
947 if(_usePrimitiveVarsInNewton)//Implicitation using primitive variables
949 _Vij[0]=_fluides[0]->getPressureFromEnthalpy(_Uroe[_nVar-1]-u_2/2, _Uroe[0]);//pressure
950 _Vij[_nVar-1]=_fluides[0]->getTemperatureFromPressure( _Vij[0], _Uroe[0]);//Temperature
951 for(int idim=0;idim<_Ndim; idim++)
952 _Vij[1+idim]=_Uroe[1+idim];
953 primToConsJacobianMatrix(_Vij);
955 Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
956 Poly.matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
959 for(int i=0; i<_nVar*_nVar;i++)
961 _AroeMinusImplicit[i] = _AroeMinus[i];
962 _AroePlusImplicit[i] = _AroePlus[i];
965 if(_verbose && _nbTimeStep%_freqSave ==0)
967 displayMatrix(_Aroe, _nVar,"Matrice de Roe");
968 displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
969 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
970 displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
974 if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
976 displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
977 displayMatrix(_AroePlusImplicit, _nVar,"Matrice _AroePlusImplicit");
980 /*********Calcul de la matrice signe pour VFFC, VFRoe et décentrement des termes source*****/
981 if(_entropicCorrection)
983 InvMatriceRoe( vp_dist);
985 Poly.matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
987 else if (_spaceScheme==upwind || (_spaceScheme ==pressureCorrection ) || (_spaceScheme ==lowMach ))//upwind sans entropic
988 SigneMatriceRoe( vp_dist);
989 else if (_spaceScheme==centered)//centre sans entropic
990 for(int i=0; i<_nVar*_nVar;i++)
992 else if( _spaceScheme ==staggered )//à tester
1001 for(int i=0; i<_nVar*_nVar;i++)
1003 _signAroe[0] = signu;
1004 for(int i=1; i<_nVar-1;i++)
1005 _signAroe[i*_nVar+i] = -signu;
1006 _signAroe[_nVar*(_nVar-1)+_nVar-1] = signu;
1010 *_runLogFile<<"SinglePhase::convectionMatrices: well balanced option not treated"<<endl;
1011 _runLogFile->close();
1012 throw CdmathException("SinglePhase::convectionMatrices: well balanced option not treated");
1015 void SinglePhase::computeScaling(double maxvp)
1019 for(int q=1; q<_nVar-1; q++)
1021 _blockDiag[q]=1./maxvp;//
1022 _invBlockDiag[q]= maxvp;//1.;//
1024 _blockDiag[_nVar - 1]=(_fluides[0]->constante("gamma")-1)/(maxvp*maxvp);//1
1025 _invBlockDiag[_nVar - 1]= 1./_blockDiag[_nVar - 1] ;// 1.;//
1028 void SinglePhase::addDiffusionToSecondMember
1034 double lambda = _fluides[0]->getConductivity(_Udiff[_nVar-1]);
1035 double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
1038 lambda=max(lambda,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
1040 if(lambda==0 && mu ==0 && _heatTransfertCoeff==0)
1044 for(int k=1; k<_nVar; k++)
1045 _idm[k] = _idm[k-1] + 1;
1048 VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
1052 //on n'a pas de contribution sur la masse
1054 //contribution visqueuse sur la quantite de mouvement
1055 for(int k=1; k<_nVar-1; k++)
1056 _phi[k] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu*(_porosityj*_Vj[k] - _porosityi*_Vi[k]);
1057 //contribution visqueuse sur l'energie
1058 _phi[_nVar-1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*lambda*(_porosityj*_Vj[_nVar-1] - _porosityi*_Vi[_nVar-1]);
1061 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1063 if(_verbose && _nbTimeStep%_freqSave ==0)
1065 cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
1066 for(int q=0; q<_nVar; q++)
1067 cout << _phi[q] << endl;
1073 //On change de signe pour l'autre contribution
1074 for(int k=0; k<_nVar; k++)
1075 _phi[k] *= -_inv_dxj/_inv_dxi;
1078 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1079 if(_verbose && _nbTimeStep%_freqSave ==0)
1081 cout << "Contribution diffusion au 2nd membre pour la maille " << j << ": "<<endl;
1082 for(int q=0; q<_nVar; q++)
1083 cout << _phi[q] << endl;
1088 if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
1090 cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
1091 for(int i=0; i<_nVar; i++)
1093 for(int j=0; j<_nVar; j++)
1094 cout << _Diffusion[i*_nVar+j]<<", ";
1101 void SinglePhase::sourceVector(PetscScalar * Si, PetscScalar * Ui, PetscScalar * Vi, int i)
1103 double phirho=Ui[0], T=Vi[_nVar-1];
1105 for(int k=0; k<_Ndim; k++)
1106 norm_u+=Vi[1+k]*Vi[1+k];
1107 norm_u=sqrt(norm_u);
1109 Si[0]=_heatPowerField(i)/_latentHeat;
1112 for(int k=1; k<_nVar-1; k++)
1113 Si[k] =(_gravite[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*phirho;
1115 Si[_nVar-1]=_heatPowerField(i);
1117 for(int k=0; k<_Ndim; k++)
1118 Si[_nVar-1] +=(_GravityField3d[k]-_dragCoeffs[0]*norm_u*Vi[1+k])*Vi[1+k]*phirho;
1120 if(_timeScheme==Implicit)
1122 for(int k=0; k<_nVar*_nVar;k++)
1123 _GravityImplicitationMatrix[k] = 0;
1124 if(!_usePrimitiveVarsInNewton)
1125 for(int k=0; k<_nVar;k++)
1126 _GravityImplicitationMatrix[k*_nVar]=-_gravite[k];
1129 double pression=Vi[0];
1130 getDensityDerivatives( pression, T, norm_u*norm_u);
1131 for(int k=0; k<_nVar;k++)
1133 _GravityImplicitationMatrix[k*_nVar+0] =-_gravite[k]*_drho_sur_dp;
1134 _GravityImplicitationMatrix[k*_nVar+_nVar-1]=-_gravite[k]*_drho_sur_dT;
1138 if(_verbose && _nbTimeStep%_freqSave ==0)
1140 cout<<"SinglePhase::sourceVector"<<endl;
1142 for(int k=0;k<_nVar;k++)
1146 for(int k=0;k<_nVar;k++)
1150 for(int k=0;k<_nVar;k++)
1153 if(_timeScheme==Implicit)
1154 displayMatrix(_GravityImplicitationMatrix, _nVar, "Gravity implicitation matrix");
1158 void SinglePhase::pressureLossVector(PetscScalar * pressureLoss, double K, PetscScalar * Ui, PetscScalar * Vi, PetscScalar * Uj, PetscScalar * Vj)
1160 double norm_u=0, u_n=0, rho;
1161 for(int i=0;i<_Ndim;i++)
1162 u_n += _Uroe[1+i]*_vec_normal[i];
1166 for(int i=0;i<_Ndim;i++)
1167 norm_u += Vi[1+i]*Vi[1+i];
1168 norm_u=sqrt(norm_u);
1170 for(int i=0;i<_Ndim;i++)
1171 pressureLoss[1+i]=-1/2*K*rho*norm_u*Vi[1+i];
1174 for(int i=0;i<_Ndim;i++)
1175 norm_u += Vj[1+i]*Vj[1+i];
1176 norm_u=sqrt(norm_u);
1178 for(int i=0;i<_Ndim;i++)
1179 pressureLoss[1+i]=-1/2*K*rho*norm_u*Vj[1+i];
1181 pressureLoss[_nVar-1]=-1/2*K*rho*norm_u*norm_u*norm_u;
1183 if(_verbose && _nbTimeStep%_freqSave ==0)
1185 cout<<"SinglePhase::pressureLossVector K= "<<K<<endl;
1187 for(int k=0;k<_nVar;k++)
1191 for(int k=0;k<_nVar;k++)
1195 for(int k=0;k<_nVar;k++)
1199 for(int k=0;k<_nVar;k++)
1202 cout<<"pressureLoss="<<endl;
1203 for(int k=0;k<_nVar;k++)
1204 cout<<pressureLoss[k]<<", ";
1209 void SinglePhase::porosityGradientSourceVector()
1211 double u_ni=0, u_nj=0, rhoi,rhoj, pi=_Vi[0], pj=_Vj[0],pij;
1212 for(int i=0;i<_Ndim;i++) {
1213 u_ni += _Vi[1+i]*_vec_normal[i];
1214 u_nj += _Vj[1+i]*_vec_normal[i];
1216 _porosityGradientSourceVector[0]=0;
1217 rhoj=_Uj[0]/_porosityj;
1218 rhoi=_Ui[0]/_porosityi;
1219 pij=(pi+pj)/2+rhoi*rhoj/2/(rhoi+rhoj)*(u_ni-u_nj)*(u_ni-u_nj);
1220 for(int i=0;i<_Ndim;i++)
1221 _porosityGradientSourceVector[1+i]=pij*(_porosityi-_porosityj)*2/(1/_inv_dxi+1/_inv_dxj);
1222 _porosityGradientSourceVector[_nVar-1]=0;
1225 void SinglePhase::jacobian(const int &j, string nameOfGroup,double * normale)
1227 if(_verbose && _nbTimeStep%_freqSave ==0)
1228 cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
1231 for(k=0; k<_nVar*_nVar;k++)
1232 _Jcb[k] = 0;//No implicitation at this stage
1235 for(k=1; k<_nVar; k++)
1236 _idm[k] = _idm[k-1] + 1;
1237 VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);
1238 double q_n=0;//quantité de mouvement normale à la paroi
1239 for(k=0; k<_Ndim; k++)
1240 q_n+=_externalStates[(k+1)]*normale[k];
1242 // loop of boundary types
1243 if (_limitField[nameOfGroup].bcType==Wall)
1245 for(k=0; k<_nVar;k++)
1246 _Jcb[k*_nVar + k] = 1;
1247 for(k=1; k<_nVar-1;k++)
1248 for(int l=1; l<_nVar-1;l++)
1249 _Jcb[k*_nVar + l] -= 2*normale[k-1]*normale[l-1];
1251 else if (_limitField[nameOfGroup].bcType==Inlet)
1255 double v[_Ndim], ve[_Ndim], v2, ve2;
1258 for(k=1; k<_nVar; k++)
1259 _idm[k] = _idm[k-1] + 1;
1260 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1261 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1263 ve[0] = _limitField[nameOfGroup].v_x[0];
1268 ve[1] = _limitField[nameOfGroup].v_y[0];
1274 ve[2] = _limitField[nameOfGroup].v_z[0];
1279 double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,_Uj[0]);
1280 double total_energy=internal_energy+ve2/2;
1283 _Jcb[0]=v2/(2*internal_energy);
1284 for(k=0; k<_Ndim;k++)
1285 _Jcb[1+k]=-v[k]/internal_energy;
1286 _Jcb[_nVar-1]=1/internal_energy;
1288 for(int l =1;l<1+_Ndim;l++){
1289 _Jcb[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
1290 for(k=0; k<_Ndim;k++)
1291 _Jcb[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
1292 _Jcb[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
1295 _Jcb[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
1296 for(k=0; k<_Ndim;k++)
1297 _Jcb[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
1298 _Jcb[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
1301 for(k=0;k<_nVar;k++)
1306 _Jcb[_nVar]=_limitField[nameOfGroup].v_x[0];//Kieu
1307 v2 +=(_limitField[nameOfGroup].v_x[0])*(_limitField[nameOfGroup].v_x[0]);
1310 _Jcb[2*_nVar]= _limitField[nameOfGroup].v_y[0];
1311 v2 +=_limitField[nameOfGroup].v_y[0]*_limitField[nameOfGroup].v_y[0];
1313 _Jcb[3*_nVar]=_limitField[nameOfGroup].v_z[0];
1314 v2 +=_limitField[nameOfGroup].v_z[0]*_limitField[nameOfGroup].v_z[0];
1317 _Jcb[(_nVar-1)*_nVar]=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho) + v2/2;
1320 else if (_limitField[nameOfGroup].bcType==InletPressure && q_n<0){
1322 double v[_Ndim], v2=0;
1324 for(k=1; k<_nVar; k++)
1325 _idm[k] = _idm[k-1] + 1;
1326 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1328 for(k=0; k<_Ndim;k++){
1333 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _limitField[nameOfGroup].T);
1334 double rho_int = _externalStates[0];
1335 double density_ratio=rho_ext/rho_int;
1337 for(int l =1;l<1+_Ndim;l++){
1338 _Jcb[l*_nVar]=-density_ratio*v[l-1];
1339 _Jcb[l*_nVar+l]=density_ratio;
1342 _Jcb[(_nVar-1)*_nVar]=-v2*density_ratio;
1343 for(k=0; k<_Ndim;k++)
1344 _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k];
1346 // not wall, not inlet, not inletPressure
1347 else if(_limitField[nameOfGroup].bcType==Outlet || (_limitField[nameOfGroup].bcType==InletPressure && q_n>=0))
1350 double v[_Ndim], v2=0;
1352 for(k=1; k<_nVar; k++)
1353 _idm[k] = _idm[k-1] + 1;
1354 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1356 for(k=0; k<_Ndim;k++){
1361 double rho_ext=_fluides[0]->getDensity(_limitField[nameOfGroup].p, _externalStates[_nVar-1]);
1362 double rho_int = _externalStates[0];
1363 double density_ratio=rho_ext/rho_int;
1364 double internal_energy=_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho_int);
1365 double total_energy=internal_energy+v2/2;
1368 _Jcb[0]=density_ratio*(1-v2/(2*internal_energy));
1369 for(k=0; k<_Ndim;k++)
1370 _Jcb[1+k]=density_ratio*v[k]/internal_energy;
1371 _Jcb[_nVar-1]=-density_ratio/internal_energy;
1373 for(int l =1;l<1+_Ndim;l++){
1374 _Jcb[l*_nVar]=density_ratio*v2*v[l-1]/(2*internal_energy);
1375 for(k=0; k<_Ndim;k++)
1376 _Jcb[l*_nVar+1+k]=density_ratio*v[k]*v[l-1]/internal_energy;
1377 _Jcb[l*_nVar+1+k]-=density_ratio;
1378 _Jcb[l*_nVar+_nVar-1]=-density_ratio*v[l-1]/internal_energy;
1381 _Jcb[(_nVar-1)*_nVar]=density_ratio*v2*total_energy/(2*internal_energy);
1382 for(k=0; k<_Ndim;k++)
1383 _Jcb[(_nVar-1)*_nVar+1+k]=density_ratio*v[k]*total_energy/internal_energy;
1384 _Jcb[(_nVar-1)*_nVar+_nVar-1]=density_ratio*(1-total_energy/internal_energy);
1388 double cd = 1,cn=0,p0, gamma;
1389 _idm[0] = j*_nVar;// Kieu
1390 for(k=1; k<_nVar;k++)
1391 _idm[k] = _idm[k-1] + 1;
1392 VecGetValues(_conservativeVars, _nVar, _idm, _phi);
1393 VecGetValues(_primitiveVars, _nVar, _idm, _externalStates);
1395 // compute the common numerator and common denominator
1396 p0=_fluides[0]->constante("p0");
1397 gamma =_fluides[0]->constante("gamma");
1398 cn =_limitField[nameOfGroup].p +p0;
1399 cd = _phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0;
1403 for(k=1; k<_nVar-1;k++)
1404 v2+=_externalStates[k]*_externalStates[k];
1406 _JcbDiff[0] = cn*(_phi[_nVar-1] -v2 -p0)/cd;
1407 for(k=1; k<_nVar-1;k++)
1408 _JcbDiff[k]=cn*_phi[k]/cd;
1409 _JcbDiff[_nVar-1]= -cn*_phi[0]/cd;
1411 for(idim=0; idim<_Ndim;idim++)
1414 _JcbDiff[(1+idim)*_nVar]=-(v2*cn*_phi[idim+1])/(2*cd);
1415 //colonnes intermediaire
1416 for(jdim=0; jdim<_Ndim;jdim++)
1418 _JcbDiff[(1+idim)*_nVar + jdim + 1] =_externalStates[idim+1]*_phi[jdim+1];
1419 _JcbDiff[(1+idim)*_nVar + jdim + 1]*=cn/cd;
1421 //matrice identite*cn*(rhoe- p0)
1422 _JcbDiff[(1+idim)*_nVar + idim + 1] +=( cn*(_phi[0]*_fluides[0]->getInternalEnergy(_externalStates[_nVar-1],rho)-p0))/cd;
1425 _JcbDiff[(1+idim)*_nVar + _nVar-1]=-_phi[idim+1]*cn/cd;
1428 _JcbDiff[_nVar*(_nVar-1)] = -(v2*_phi[_nVar -1]*cn)/(2*cd);
1429 for(int idim=0; idim<_Ndim;idim++)
1430 _JcbDiff[_nVar*(_nVar-1)+idim+1]=_externalStates[idim+1]*_phi[_nVar -1]*cn/cd;
1431 _JcbDiff[_nVar*_nVar -1] = -(v2/2+p0)*cn/cd;
1434 else if (_limitField[nameOfGroup].bcType!=Neumann)// not wall, not inlet, not outlet
1436 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1437 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1438 _runLogFile->close();
1439 throw CdmathException("SinglePhase::jacobian: This boundary condition is not treated");
1443 //calcule la jacobienne pour une CL de diffusion
1444 void SinglePhase::jacobianDiff(const int &j, string nameOfGroup)
1446 if(_verbose && _nbTimeStep%_freqSave ==0)
1447 cout<<"Jacobienne condition limite diffusion bord "<< nameOfGroup<<endl;
1450 for(k=0; k<_nVar*_nVar;k++)
1453 if (_limitField[nameOfGroup].bcType==Wall){
1454 double v[_Ndim], ve[_Ndim], v2, ve2;
1457 for(k=1; k<_nVar; k++)
1458 _idm[k] = _idm[k-1] + 1;
1459 VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
1460 VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
1462 ve[0] = _limitField[nameOfGroup].v_x[0];
1467 ve[1] = _limitField[nameOfGroup].v_y[0];
1473 ve[2] = _limitField[nameOfGroup].v_z[0];
1479 double internal_energy=_fluides[0]->getInternalEnergy(_limitField[nameOfGroup].T,rho);
1480 double total_energy=internal_energy+ve2/2;
1483 _JcbDiff[0]=v2/(2*internal_energy);
1484 for(k=0; k<_Ndim;k++)
1485 _JcbDiff[1+k]=-v[k]/internal_energy;
1486 _JcbDiff[_nVar-1]=1/internal_energy;
1488 for(int l =1;l<1+_Ndim;l++){
1489 _JcbDiff[l*_nVar]=v2*ve[l-1]/(2*internal_energy);
1490 for(k=0; k<_Ndim;k++)
1491 _JcbDiff[l*_nVar+1+k]=-v[k]*ve[l-1]/internal_energy;
1492 _JcbDiff[l*_nVar+_nVar-1]=ve[l-1]/internal_energy;
1495 _JcbDiff[(_nVar-1)*_nVar]=v2*total_energy/(2*internal_energy);
1496 for(k=0; k<_Ndim;k++)
1497 _JcbDiff[(_nVar-1)*_nVar+1+k]=-v[k]*total_energy/internal_energy;
1498 _JcbDiff[(_nVar-1)*_nVar+_nVar-1]=total_energy/internal_energy;
1500 else if (_limitField[nameOfGroup].bcType==Outlet || _limitField[nameOfGroup].bcType==Neumann
1501 ||_limitField[nameOfGroup].bcType==Inlet || _limitField[nameOfGroup].bcType==InletPressure)
1503 for(k=0;k<_nVar;k++)
1504 _JcbDiff[k*_nVar+k]=1;
1507 cout << "group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1508 *_runLogFile<<"group named "<<nameOfGroup << " : unknown boundary condition" << endl;
1509 _runLogFile->close();
1510 throw CdmathException("SinglePhase::jacobianDiff: This boundary condition is not recognised");
1514 void SinglePhase::primToCons(const double *P, const int &i, double *W, const int &j){
1515 //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;
1517 double rho=_fluides[0]->getDensity(P[i*(_Ndim+2)], P[i*(_Ndim+2)+(_Ndim+1)]);
1518 W[j*(_Ndim+2)] = _porosityField(j)*rho;//phi*rho
1519 for(int k=0; k<_Ndim; k++)
1520 W[j*(_Ndim+2)+(k+1)] = W[j*(_Ndim+2)]*P[i*(_Ndim+2)+(k+1)];//phi*rho*u
1522 W[j*(_Ndim+2)+(_Ndim+1)] = W[j*(_Ndim+2)]*_fluides[0]->getInternalEnergy(P[i*(_Ndim+2)+ (_Ndim+1)],rho);//rho*e
1523 for(int k=0; k<_Ndim; k++)
1524 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
1527 void SinglePhase::primToConsJacobianMatrix(double *V)
1529 double pression=V[0];
1530 double temperature=V[_nVar-1];
1531 double vitesse[_Ndim];
1532 for(int idim=0;idim<_Ndim;idim++)
1533 vitesse[idim]=V[1+idim];
1535 for(int idim=0;idim<_Ndim;idim++)
1536 v2+=vitesse[idim]*vitesse[idim];
1538 double rho=_fluides[0]->getDensity(pression,temperature);
1539 double gamma=_fluides[0]->constante("gamma");
1540 double Pinf=_fluides[0]->constante("p0");
1541 double q=_fluides[0]->constante("q");
1542 double cv=_fluides[0]->constante("cv");
1544 for(int k=0;k<_nVar*_nVar; k++)
1545 _primToConsJacoMat[k]=0;
1547 if( !_useDellacherieEOS)
1549 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
1550 double e=fluide0->getInternalEnergy(temperature);
1553 _primToConsJacoMat[0]=1/((gamma-1)*(e-q));
1554 _primToConsJacoMat[_nVar-1]=-rho*cv/(e-q);
1556 for(int idim=0;idim<_Ndim;idim++)
1558 _primToConsJacoMat[_nVar+idim*_nVar]=vitesse[idim]/((gamma-1)*(e-q));
1559 _primToConsJacoMat[_nVar+idim*_nVar+1+idim]=rho;
1560 _primToConsJacoMat[_nVar+idim*_nVar+_nVar-1]=-rho*vitesse[idim]*cv/(e-q);
1562 _primToConsJacoMat[(_nVar-1)*_nVar]=E/((gamma-1)*(e-q));
1563 for(int idim=0;idim<_Ndim;idim++)
1564 _primToConsJacoMat[(_nVar-1)*_nVar+1+idim]=rho*vitesse[idim];
1565 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*cv*(1-E/(e-q));
1567 else if( _useDellacherieEOS)
1569 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
1570 double h=fluide0->getEnthalpy(temperature);
1572 double cp=_fluides[0]->constante("cp");
1574 _primToConsJacoMat[0]=gamma/((gamma-1)*(h-q));
1575 _primToConsJacoMat[_nVar-1]=-rho*cp/(h-q);
1577 for(int idim=0;idim<_Ndim;idim++)
1579 _primToConsJacoMat[_nVar+idim*_nVar]=gamma*vitesse[idim]/((gamma-1)*(h-q));
1580 _primToConsJacoMat[_nVar+idim*_nVar+1+idim]=rho;
1581 _primToConsJacoMat[_nVar+idim*_nVar+_nVar-1]=-rho*vitesse[idim]*cp/(h-q);
1583 _primToConsJacoMat[(_nVar-1)*_nVar]=gamma*H/((gamma-1)*(h-q))-1;
1584 for(int idim=0;idim<_Ndim;idim++)
1585 _primToConsJacoMat[(_nVar-1)*_nVar+1+idim]=rho*vitesse[idim];
1586 _primToConsJacoMat[(_nVar-1)*_nVar+_nVar-1]=rho*cp*(1-H/(h-q));
1590 *_runLogFile<<"SinglePhase::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie"<<endl;
1591 _runLogFile->close();
1592 throw CdmathException("SinglePhase::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
1595 if(_verbose && _nbTimeStep%_freqSave ==0)
1597 cout<<" SinglePhase::primToConsJacobianMatrix" << endl;
1598 displayVector(_Vi,_nVar," _Vi " );
1599 cout<<" Jacobienne primToCons: " << endl;
1600 displayMatrix(_primToConsJacoMat,_nVar," Jacobienne primToCons: ");
1604 void SinglePhase::consToPrim(const double *Wcons, double* Wprim,double porosity)
1607 for(int k=1;k<=_Ndim;k++)
1608 q_2 += Wcons[k]*Wcons[k];
1609 q_2 /= Wcons[0]; //phi rho u²
1610 double rhoe=(Wcons[(_Ndim+2)-1]-q_2/2)/porosity;
1611 double rho=Wcons[0]/porosity;
1612 Wprim[0] = _fluides[0]->getPressure(rhoe,rho);//pressure p
1614 cout << "pressure = "<< Wprim[0] << " < 0 " << endl;
1615 *_runLogFile<< "pressure = "<< Wprim[0] << " < 0 " << endl;
1616 _runLogFile->close();
1617 throw CdmathException("SinglePhase::consToPrim: negative pressure");
1619 for(int k=1;k<=_Ndim;k++)
1620 Wprim[k] = Wcons[k]/Wcons[0];//velocity u
1621 Wprim[(_Ndim+2)-1] = _fluides[0]->getTemperatureFromPressure(Wprim[0],Wcons[0]/porosity);
1623 if(_verbose && _nbTimeStep%_freqSave ==0)
1625 cout<<"ConsToPrim Vecteur conservatif"<<endl;
1626 for(int k=0;k<_nVar;k++)
1627 cout<<Wcons[k]<<endl;
1628 cout<<"ConsToPrim Vecteur primitif"<<endl;
1629 for(int k=0;k<_nVar;k++)
1630 cout<<Wprim[k]<<endl;
1634 void SinglePhase::entropicShift(double* n)//TO do: make sure _Vi and _Vj are well set
1636 /*Left and right values */
1637 double ul_n = 0, ul_2=0, ur_n=0, ur_2 = 0; //valeurs de vitesse normale et |u|² a droite et a gauche
1638 for(int i=0;i<_Ndim;i++)
1640 ul_n += _Vi[1+i]*n[i];
1641 ul_2 += _Vi[1+i]*_Vi[1+i];
1642 ur_n += _Vj[1+i]*n[i];
1643 ur_2 += _Vj[1+i]*_Vj[1+i];
1646 double cl = _fluides[0]->vitesseSonEnthalpie(_Vi[_Ndim+1]-ul_2/2);//vitesse du son a l'interface
1647 double cr = _fluides[0]->vitesseSonEnthalpie(_Vj[_Ndim+1]-ur_2/2);//vitesse du son a l'interface
1649 _entropicShift[0]=fabs(ul_n-cl - (ur_n-cr));
1650 _entropicShift[1]=fabs(ul_n - ur_n);
1651 _entropicShift[2]=fabs(ul_n+cl - (ur_n+cr));
1653 if(_verbose && _nbTimeStep%_freqSave ==0)
1655 cout<<"Entropic shift left eigenvalues: "<<endl;
1656 cout<<"("<< ul_n-cl<< ", " <<ul_n<<", "<<ul_n+cl << ")";
1658 cout<<"Entropic shift right eigenvalues: "<<endl;
1659 cout<<"("<< ur_n-cr<< ", " <<ur_n<<", "<<ur_n+cr << ")";
1661 cout<<"eigenvalue jumps "<<endl;
1662 cout<< _entropicShift[0] << ", " << _entropicShift[1] << ", "<< _entropicShift[2] <<endl;
1666 Vector SinglePhase::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
1667 if(_verbose && _nbTimeStep%_freqSave ==0)
1669 cout<<"SinglePhase::convectionFlux start"<<endl;
1670 cout<<"Ucons"<<endl;
1672 cout<<"Vprim"<<endl;
1676 double phirho=U(0);//phi rho
1677 Vector phiq(_Ndim);//phi rho u
1678 for(int i=0;i<_Ndim;i++)
1681 double pression=V(0);
1682 Vector vitesse(_Ndim);
1683 for(int i=0;i<_Ndim;i++)
1685 double Temperature= V(1+_Ndim);
1687 double vitessen=vitesse*normale;
1688 double rho=phirho/porosity;
1689 double e_int=_fluides[0]->getInternalEnergy(Temperature,rho);
1692 F(0)=phirho*vitessen;
1693 for(int i=0;i<_Ndim;i++)
1694 F(1+i)=phirho*vitessen*vitesse(i)+pression*normale(i)*porosity;
1695 F(1+_Ndim)=phirho*(e_int+0.5*vitesse*vitesse+pression/rho)*vitessen;
1697 if(_verbose && _nbTimeStep%_freqSave ==0)
1699 cout<<"SinglePhase::convectionFlux end"<<endl;
1700 cout<<"Flux F(U,V)"<<endl;
1707 void SinglePhase::RoeMatrixConservativeVariables(double u_n, double H,Vector velocity, double k, double K)
1709 /******** Construction de la matrice de Roe *********/
1710 //premiere ligne (masse)
1712 for(int idim=0; idim<_Ndim;idim++)
1713 _Aroe[1+idim]=_vec_normal[idim];
1716 //lignes intermadiaires(qdm)
1717 for(int idim=0; idim<_Ndim;idim++)
1720 _Aroe[(1+idim)*_nVar]=K*_vec_normal[idim] - u_n*_Uroe[1+idim];
1721 //colonnes intermediaires
1722 for(int jdim=0; jdim<_Ndim;jdim++)
1723 _Aroe[(1+idim)*_nVar + jdim + 1] = _Uroe[1+idim]*_vec_normal[jdim]-k*_vec_normal[idim]*_Uroe[1+jdim];
1725 _Aroe[(1+idim)*_nVar + idim + 1] += u_n;
1727 _Aroe[(1+idim)*_nVar + _nVar-1]=k*_vec_normal[idim];
1730 //derniere ligne (energie)
1731 _Aroe[_nVar*(_nVar-1)] = (K - H)*u_n;
1732 for(int idim=0; idim<_Ndim;idim++)
1733 _Aroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] - k*u_n*_Uroe[idim+1];
1734 _Aroe[_nVar*_nVar -1] = (1 + k)*u_n;
1736 void SinglePhase::convectionMatrixPrimitiveVariables( double rho, double u_n, double H,Vector vitesse)
1738 //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1739 //On remplit la matrice de Roe en variables primitives : F(V_L)-F(V_R)=Aroe (V_L-V_R)
1740 //EOS is more involved with primitive variables
1741 // call to getDensityDerivatives(double concentration, double pression, double temperature,double v2) needed
1742 _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1743 for(int i=0;i<_Ndim;i++)
1744 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1745 _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1746 for(int i=0;i<_Ndim;i++)
1748 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]+_vec_normal[i];
1749 for(int j=0;j<_Ndim;j++)
1750 _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1751 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1752 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1754 _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1755 for(int i=0;i<_Ndim;i++)
1756 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1757 _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1759 void SinglePhase::staggeredRoeUpwindingMatrixConservativeVariables( double u_n, double H,Vector velocity, double k, double K)
1761 //Calcul de décentrement de type décalé pour formulation de Roe
1762 if(fabs(u_n)>_precision)
1764 //premiere ligne (masse)
1766 for(int idim=0; idim<_Ndim;idim++)
1767 _absAroe[1+idim]=_vec_normal[idim];
1768 _absAroe[_nVar-1]=0;
1770 //lignes intermadiaires(qdm)
1771 for(int idim=0; idim<_Ndim;idim++)
1774 _absAroe[(1+idim)*_nVar]=-K*_vec_normal[idim] - u_n*_Uroe[1+idim];
1775 //colonnes intermediaires
1776 for(int jdim=0; jdim<_Ndim;jdim++)
1777 _absAroe[(1+idim)*_nVar + jdim + 1] = _Uroe[1+idim]*_vec_normal[jdim]+k*_vec_normal[idim]*_Uroe[1+jdim];
1779 _absAroe[(1+idim)*_nVar + idim + 1] += u_n;
1781 _absAroe[(1+idim)*_nVar + _nVar-1]=-k*_vec_normal[idim];
1784 //derniere ligne (energie)
1785 _absAroe[_nVar*(_nVar-1)] = (-K - H)*u_n;
1786 for(int idim=0; idim<_Ndim;idim++)
1787 _absAroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] + k*u_n*_Uroe[idim+1];
1788 _absAroe[_nVar*_nVar -1] = (1 - k)*u_n;
1796 for(int i=0; i<_nVar*_nVar;i++)
1797 _absAroe[i] *= signu;
1799 else//umn=0 ->centered scheme
1801 for(int i=0; i<_nVar*_nVar;i++)
1805 void SinglePhase::staggeredRoeUpwindingMatrixPrimitiveVariables(double rho, double u_n,double H, Vector vitesse)
1807 //Not used. Suppress or use in alternative implicitation in primitive variable of the staggered-roe scheme
1808 //Calcul de décentrement de type décalé pour formulation Roe
1809 _AroeImplicit[0*_nVar+0]=_drho_sur_dp*u_n;
1810 for(int i=0;i<_Ndim;i++)
1811 _AroeImplicit[0*_nVar+1+i]=rho*_vec_normal[i];
1812 _AroeImplicit[0*_nVar+1+_Ndim]=_drho_sur_dT*u_n;
1813 for(int i=0;i<_Ndim;i++)
1815 _AroeImplicit[(1+i)*_nVar+0]=_drho_sur_dp *u_n*vitesse[i]-_vec_normal[i];
1816 for(int j=0;j<_Ndim;j++)
1817 _AroeImplicit[(1+i)*_nVar+1+j]=rho*vitesse[i]*_vec_normal[j];
1818 _AroeImplicit[(1+i)*_nVar+1+i]+=rho*u_n;
1819 _AroeImplicit[(1+i)*_nVar+1+_Ndim]=_drho_sur_dT*u_n*vitesse[i];
1821 _AroeImplicit[(1+_Ndim)*_nVar+0]=(_drhoE_sur_dp+1)*u_n;
1822 for(int i=0;i<_Ndim;i++)
1823 _AroeImplicit[(1+_Ndim)*_nVar+1+i]=rho*(H*_vec_normal[i]+u_n*vitesse[i]);
1824 _AroeImplicit[(1+_Ndim)*_nVar+1+_Ndim]=_drhoE_sur_dT*u_n;
1827 Vector SinglePhase::staggeredVFFCFlux()
1829 if(_verbose && _nbTimeStep%_freqSave ==0)
1830 cout<<"SinglePhase::staggeredVFFCFlux() start"<<endl;
1832 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1834 *_runLogFile<< "SinglePhase::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding, pressure = "<< endl;
1835 _runLogFile->close();
1836 throw CdmathException("SinglePhase::staggeredVFFCFlux: staggeredVFFCFlux method should be called only for VFFC formulation and staggered upwinding");
1838 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
1842 double uijn=0, phiqn=0, uin=0, ujn=0;
1843 for(int idim=0;idim<_Ndim;idim++)
1845 uijn+=_vec_normal[idim]*_Uroe[1+idim];//URoe = rho, u, H
1846 uin +=_vec_normal[idim]*_Ui[1+idim];
1847 ujn +=_vec_normal[idim]*_Uj[1+idim];
1850 if( (uin>0 && ujn >0) || (uin>=0 && ujn <=0 && uijn>0) ) // formerly (uijn>_precision)
1852 for(int idim=0;idim<_Ndim;idim++)
1853 phiqn+=_vec_normal[idim]*_Ui[1+idim];//phi rho u n
1855 for(int idim=0;idim<_Ndim;idim++)
1856 Fij(1+idim)=phiqn*_Vi[1+idim]+_Vj[0]*_vec_normal[idim]*_porosityj;
1857 Fij(_nVar-1)=phiqn/_Ui[0]*(_Ui[_nVar-1]+_Vj[0]*sqrt(_porosityj/_porosityi));
1859 else if( (uin<0 && ujn <0) || (uin>=0 && ujn <=0 && uijn<0) ) // formerly (uijn<-_precision)
1861 for(int idim=0;idim<_Ndim;idim++)
1862 phiqn+=_vec_normal[idim]*_Uj[1+idim];//phi rho u n
1864 for(int idim=0;idim<_Ndim;idim++)
1865 Fij(1+idim)=phiqn*_Vj[1+idim]+_Vi[0]*_vec_normal[idim]*_porosityi;
1866 Fij(_nVar-1)=phiqn/_Uj[0]*(_Uj[_nVar-1]+_Vi[0]*sqrt(_porosityi/_porosityj));
1868 else//case (uin<=0 && ujn >=0) or (uin>=0 && ujn <=0 && uijn==0), apply centered scheme
1870 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar);
1871 Vector normale(_Ndim);
1872 for(int i1=0;i1<_Ndim;i1++)
1873 normale(i1)=_vec_normal[i1];
1874 for(int i1=0;i1<_nVar;i1++)
1881 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
1882 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
1883 Fij=(Fi+Fj)/2;//+_maxvploc*(Ui-Uj)/2;
1885 if(_verbose && _nbTimeStep%_freqSave ==0)
1887 cout<<"SinglePhase::staggeredVFFCFlux() endf uijn="<<uijn<<endl;
1894 void SinglePhase::staggeredVFFCMatricesConservativeVariables(double un)//vitesse normale de Roe en entrée
1896 if(_verbose && _nbTimeStep%_freqSave ==0)
1897 cout<<"SinglePhase::staggeredVFFCMatrices()"<<endl;
1899 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
1901 *_runLogFile<< "SinglePhase::staggeredVFFCMatrices: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding"<< endl;
1902 _runLogFile->close();
1903 throw CdmathException("SinglePhase::staggeredVFFCMatrices: staggeredVFFCMatrices method should be called only for VFFC formulation and staggered upwinding");
1905 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
1907 double ui_n=0, ui_2=0, uj_n=0, uj_2=0;//vitesse normale et carré du module
1908 double H;//enthalpie totale (expression particulière)
1909 consToPrim(_Ui,_Vi,_porosityi);
1910 consToPrim(_Uj,_Vj,_porosityj);
1911 for(int i=0;i<_Ndim;i++)
1913 ui_2 += _Vi[1+i]*_Vi[1+i];
1914 ui_n += _Vi[1+i]*_vec_normal[i];
1915 uj_2 += _Vj[1+i]*_Vj[1+i];
1916 uj_n += _Vj[1+i]*_vec_normal[i];
1919 double rhoi,pj, Ei, rhoj;
1920 double cj, Kj, kj;//dérivées de la pression
1921 rhoi=_Ui[0]/_porosityi;
1922 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
1924 rhoj=_Uj[0]/_porosityj;
1926 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
1927 kj = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
1928 Kj = uj_2*kj/2; //g-1/2 *|u|²
1931 double ci, Ki, ki;//dérivées de la pression
1932 Ej= _Uj[_Ndim+1]/rhoj;
1935 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
1936 ki = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
1937 Ki = ui_2*ki/2; //g-1/2 *|u|²
1941 /***********Calcul des valeurs propres ********/
1942 vector<complex<double> > vp_dist(3,0);
1943 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
1945 _maxvploc=fabs(ui_n)+cj;
1946 if(_maxvploc>_maxvp)
1949 if(_verbose && _nbTimeStep%_freqSave ==0)
1950 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
1952 /******** Construction de la matrice A^+ *********/
1953 //premiere ligne (masse)
1955 for(int idim=0; idim<_Ndim;idim++)
1956 _AroePlus[1+idim]=_vec_normal[idim];
1957 _AroePlus[_nVar-1]=0;
1959 //lignes intermadiaires(qdm)
1960 for(int idim=0; idim<_Ndim;idim++)
1963 _AroePlus[(1+idim)*_nVar]=- ui_n*_Vi[1+idim];
1964 //colonnes intermediaires
1965 for(int jdim=0; jdim<_Ndim;jdim++)
1966 _AroePlus[(1+idim)*_nVar + jdim + 1] = _Vi[1+idim]*_vec_normal[jdim];
1968 _AroePlus[(1+idim)*_nVar + idim + 1] += ui_n;
1970 _AroePlus[(1+idim)*_nVar + _nVar-1]=0;
1973 //derniere ligne (energie)
1974 _AroePlus[_nVar*(_nVar-1)] = - H*ui_n;
1975 for(int idim=0; idim<_Ndim;idim++)
1976 _AroePlus[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] ;
1977 _AroePlus[_nVar*_nVar -1] = ui_n;
1979 /******** Construction de la matrice A^- *********/
1980 //premiere ligne (masse)
1982 for(int idim=0; idim<_Ndim;idim++)
1983 _AroeMinus[1+idim]=0;
1984 _AroeMinus[_nVar-1]=0;
1986 //lignes intermadiaires(qdm)
1987 for(int idim=0; idim<_Ndim;idim++)
1990 _AroeMinus[(1+idim)*_nVar]=Kj*_vec_normal[idim] ;
1991 //colonnes intermediaires
1992 for(int jdim=0; jdim<_Ndim;jdim++)
1993 _AroeMinus[(1+idim)*_nVar + jdim + 1] = -kj*_vec_normal[idim]*_Vj[1+jdim];
1995 _AroeMinus[(1+idim)*_nVar + idim + 1] += 0;
1997 _AroeMinus[(1+idim)*_nVar + _nVar-1]=kj*_vec_normal[idim];
2000 //derniere ligne (energie)
2001 _AroeMinus[_nVar*(_nVar-1)] = Kj *ui_n;
2002 for(int idim=0; idim<_Ndim;idim++)
2003 _AroeMinus[_nVar*(_nVar-1)+idim+1]= - kj*uj_n*_Vi[idim+1];
2004 _AroeMinus[_nVar*_nVar -1] = kj*ui_n;
2006 else if(un<-_precision)
2008 /***********Calcul des valeurs propres ********/
2009 vector<complex<double> > vp_dist(3,0);
2010 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2012 _maxvploc=fabs(uj_n)+ci;
2013 if(_maxvploc>_maxvp)
2016 if(_verbose && _nbTimeStep%_freqSave ==0)
2017 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2019 /******** Construction de la matrice A^+ *********/
2020 //premiere ligne (masse)
2022 for(int idim=0; idim<_Ndim;idim++)
2023 _AroePlus[1+idim]=0;
2024 _AroePlus[_nVar-1]=0;
2026 //lignes intermadiaires(qdm)
2027 for(int idim=0; idim<_Ndim;idim++)
2030 _AroePlus[(1+idim)*_nVar]=Ki*_vec_normal[idim] ;
2031 //colonnes intermediaires
2032 for(int jdim=0; jdim<_Ndim;jdim++)
2033 _AroePlus[(1+idim)*_nVar + jdim + 1] = -ki*_vec_normal[idim]*_Vi[1+jdim];
2035 _AroePlus[(1+idim)*_nVar + idim + 1] += 0;
2037 _AroePlus[(1+idim)*_nVar + _nVar-1]=ki*_vec_normal[idim];
2040 //derniere ligne (energie)
2041 _AroePlus[_nVar*(_nVar-1)] = Ki *uj_n;
2042 for(int idim=0; idim<_Ndim;idim++)
2043 _AroePlus[_nVar*(_nVar-1)+idim+1]= - ki*ui_n*_Vj[idim+1];
2044 _AroePlus[_nVar*_nVar -1] = ki*uj_n;
2046 /******** Construction de la matrice A^- *********/
2047 //premiere ligne (masse)
2049 for(int idim=0; idim<_Ndim;idim++)
2050 _AroeMinus[1+idim]=_vec_normal[idim];
2051 _AroeMinus[_nVar-1]=0;
2053 //lignes intermadiaires(qdm)
2054 for(int idim=0; idim<_Ndim;idim++)
2057 _AroeMinus[(1+idim)*_nVar]= - uj_n*_Vj[1+idim];
2058 //colonnes intermediaires
2059 for(int jdim=0; jdim<_Ndim;jdim++)
2060 _AroeMinus[(1+idim)*_nVar + jdim + 1] = _Vj[1+idim]*_vec_normal[jdim];
2062 _AroeMinus[(1+idim)*_nVar + idim + 1] += uj_n;
2064 _AroeMinus[(1+idim)*_nVar + _nVar-1]=0;
2067 //derniere ligne (energie)
2068 _AroeMinus[_nVar*(_nVar-1)] = - H*uj_n;
2069 for(int idim=0; idim<_Ndim;idim++)
2070 _AroeMinus[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] ;
2071 _AroeMinus[_nVar*_nVar -1] = uj_n;
2073 else//case nil velocity on the interface, apply centered scheme
2075 double u_n=0, u_2=0;//vitesse normale et carré du module
2076 for(int i=0;i<_Ndim;i++)
2078 u_2 += _Uroe[1+i]*_Uroe[1+i];
2079 u_n += _Uroe[1+i]*_vec_normal[i];
2081 Vector vitesse(_Ndim);
2082 for(int idim=0;idim<_Ndim;idim++)
2083 vitesse[idim]=_Uroe[1+idim];
2086 /***********Calcul des valeurs propres ********/
2088 c = _fluides[0]->vitesseSonEnthalpie(H-u_2/2);//vitesse du son a l'interface
2089 k = _fluides[0]->constante("gamma") - 1;//A generaliser pour porosite et stephane gas law
2090 K = u_2*k/2; //g-1/2 *|u|²
2092 _maxvploc=fabs(u_n)+c;
2093 if(_maxvploc>_maxvp)
2096 /******** Construction de la matrice A^+ *********/
2097 //premiere ligne (masse)
2099 for(int idim=0; idim<_Ndim;idim++)
2100 _AroePlus[1+idim]=0;
2101 _AroePlus[_nVar-1]=0;
2103 //lignes intermadiaires(qdm)
2104 for(int idim=0; idim<_Ndim;idim++)
2107 _AroePlus[(1+idim)*_nVar]=- ui_n*_Vi[1+idim];
2108 //colonnes intermediaires
2109 for(int jdim=0; jdim<_Ndim;jdim++)
2110 _AroePlus[(1+idim)*_nVar + jdim + 1] = _Vi[1+idim]*_vec_normal[jdim]-0.5*ki*_vec_normal[idim]*_Vi[1+jdim];
2112 _AroePlus[(1+idim)*_nVar + idim + 1] += 0.5*ui_n;
2114 _AroePlus[(1+idim)*_nVar + _nVar-1]=0.5*ki*_vec_normal[idim];
2117 //derniere ligne (energie)
2118 _AroePlus[_nVar*(_nVar-1)] = 0;
2119 for(int idim=0; idim<_Ndim;idim++)
2120 _AroePlus[_nVar*(_nVar-1)+idim+1]=0 ;
2121 _AroePlus[_nVar*_nVar -1] = 0;
2123 /******** Construction de la matrice A^- *********/
2124 //premiere ligne (masse)
2126 for(int idim=0; idim<_Ndim;idim++)
2127 _AroeMinus[1+idim]=0;
2128 _AroeMinus[_nVar-1]=0;
2130 //lignes intermadiaires(qdm)
2131 for(int idim=0; idim<_Ndim;idim++)
2134 _AroeMinus[(1+idim)*_nVar]=Kj*_vec_normal[idim] - uj_n*_Vj[1+idim];
2135 //colonnes intermediaires
2136 for(int jdim=0; jdim<_Ndim;jdim++)
2137 _AroeMinus[(1+idim)*_nVar + jdim + 1] = -0.5*kj*_vec_normal[idim]*_Vj[1+jdim];
2139 _AroeMinus[(1+idim)*_nVar + idim + 1] += 0.5*uj_n;
2141 _AroeMinus[(1+idim)*_nVar + _nVar-1]=0.5*kj*_vec_normal[idim];
2144 //derniere ligne (energie)
2145 _AroeMinus[_nVar*(_nVar-1)] = 0;
2146 for(int idim=0; idim<_Ndim;idim++)
2147 _AroeMinus[_nVar*(_nVar-1)+idim+1]= 0;
2148 _AroeMinus[_nVar*_nVar -1] = 0;
2151 if(_timeScheme==Implicit)
2152 for(int i=0; i<_nVar*_nVar;i++)
2154 _AroeMinusImplicit[i] = _AroeMinus[i];
2155 _AroePlusImplicit[i] = _AroePlus[i];
2158 /******** Construction de la matrices Aroe *********/
2160 //premiere ligne (masse)
2162 for(int idim=0; idim<_Ndim;idim++)
2163 _Aroe[1+idim]=_vec_normal[idim];
2166 //lignes intermadiaires(qdm)
2167 for(int idim=0; idim<_Ndim;idim++)
2170 _Aroe[(1+idim)*_nVar]=Ki*_vec_normal[idim] - uj_n*_Uj[1+idim];
2171 //colonnes intermediaires
2172 for(int jdim=0; jdim<_Ndim;jdim++)
2173 _Aroe[(1+idim)*_nVar + jdim + 1] = _Uj[1+idim]*_vec_normal[jdim]-ki*_vec_normal[idim]*_Ui[1+jdim];
2175 _Aroe[(1+idim)*_nVar + idim + 1] += uj_n;
2177 _Aroe[(1+idim)*_nVar + _nVar-1]=ki*_vec_normal[idim];
2180 //derniere ligne (energie)
2181 _Aroe[_nVar*(_nVar-1)] = (Ki - H)*uj_n;
2182 for(int idim=0; idim<_Ndim;idim++)
2183 _Aroe[_nVar*(_nVar-1)+idim+1]=H*_vec_normal[idim] - ki*ui_n*_Uj[idim+1];
2184 _Aroe[_nVar*_nVar -1] = (1 + ki)*uj_n;
2188 void SinglePhase::staggeredVFFCMatricesPrimitiveVariables(double un)//vitesse normale de Roe en entrée
2190 if(_verbose && _nbTimeStep%_freqSave ==0)
2191 cout<<"SinglePhase::staggeredVFFCMatricesPrimitiveVariables()"<<endl;
2193 if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
2195 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding" << endl;
2196 _runLogFile->close();
2197 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: staggeredVFFCMatricesPrimitiveVariables method should be called only for VFFC formulation and staggered upwinding");
2199 else//_spaceScheme==staggered && _nonLinearFormulation==VFFC
2201 double ui_n=0., ui_2=0., uj_n=0., uj_2=0.;//vitesse normale et carré du module
2202 double H;//enthalpie totale (expression particulière)
2203 consToPrim(_Ui,_Vi,_porosityi);
2204 consToPrim(_Uj,_Vj,_porosityj);
2206 for(int i=0;i<_Ndim;i++)
2208 ui_2 += _Vi[1+i]*_Vi[1+i];
2209 ui_n += _Vi[1+i]*_vec_normal[i];
2210 uj_2 += _Vj[1+i]*_Vj[1+i];
2211 uj_n += _Vj[1+i]*_vec_normal[i];
2214 if(_verbose && _nbTimeStep%_freqSave ==0){
2215 cout <<"SinglePhase::staggeredVFFCMatricesPrimitiveVariables " << endl;
2216 cout<<"Vecteur primitif _Vi" << endl;
2217 for(int i=0;i<_nVar;i++)
2220 cout<<"Vecteur primitif _Vj" << endl;
2221 for(int i=0;i<_nVar;i++)
2226 double gamma=_fluides[0]->constante("gamma");
2227 double q=_fluides[0]->constante("q");
2229 if(fabs(un)>_precision)//non zero velocity on the interface
2231 if( !_useDellacherieEOS)
2233 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2234 double cv=fluide0->constante("cv");
2238 double rhoi,rhoj,pj, Ei, ei;
2239 double cj;//vitesse du son pour calcul valeurs propres
2240 rhoi=_Ui[0]/_porosityi;
2241 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
2244 rhoj=_Uj[0]/_porosityj;
2245 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
2247 /***********Calcul des valeurs propres ********/
2248 vector<complex<double> > vp_dist(3,0);
2249 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
2251 _maxvploc=fabs(ui_n)+cj;
2252 if(_maxvploc>_maxvp)
2255 if(_verbose && _nbTimeStep%_freqSave ==0)
2256 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
2258 /******** Construction de la matrice A^+ *********/
2259 //premiere ligne (masse)
2260 _AroePlusImplicit[0]=ui_n/((gamma-1)*(ei-q));
2261 for(int idim=0; idim<_Ndim;idim++)
2262 _AroePlusImplicit[1+idim]=rhoi*_vec_normal[idim];
2263 _AroePlusImplicit[_nVar-1]=-rhoi*ui_n*cv/(ei-q);
2265 //lignes intermadiaires(qdm)
2266 for(int idim=0; idim<_Ndim;idim++)
2269 _AroePlusImplicit[(1+idim)*_nVar]=ui_n/((gamma-1)*(ei-q))*_Vi[1+idim];
2270 //colonnes intermediaires
2271 for(int jdim=0; jdim<_Ndim;jdim++)
2272 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = rhoi*_Vi[1+idim]*_vec_normal[jdim];
2274 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += rhoi*ui_n;
2276 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoi*ui_n*cv/(ei-q)*_Vi[1+idim];
2279 //derniere ligne (energie)
2280 _AroePlusImplicit[_nVar*(_nVar-1)] = Ei*ui_n/((gamma-1)*(ei-q));
2281 for(int idim=0; idim<_Ndim;idim++)
2282 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoi*Ei+pj)*_vec_normal[idim]+rhoi*ui_n*_Vi[1+idim];
2283 _AroePlusImplicit[_nVar*_nVar -1] = rhoi*ui_n*(1-Ei/(ei-q))*cv;
2285 /******** Construction de la matrice A^- *********/
2286 //premiere ligne (masse)
2287 _AroeMinusImplicit[0]=0;
2288 for(int idim=0; idim<_Ndim;idim++)
2289 _AroeMinusImplicit[1+idim]=0;
2290 _AroeMinusImplicit[_nVar-1]=0;
2292 //lignes intermadiaires(qdm)
2293 for(int idim=0; idim<_Ndim;idim++)
2296 _AroeMinusImplicit[(1+idim)*_nVar]=_vec_normal[idim] ;
2297 //colonnes intermediaires
2298 for(int jdim=0; jdim<_Ndim;jdim++)
2299 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2301 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2303 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2306 //derniere ligne (energie)
2307 _AroeMinusImplicit[_nVar*(_nVar-1)] = ui_n;
2308 for(int idim=0; idim<_Ndim;idim++)
2309 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2310 _AroeMinusImplicit[_nVar*_nVar -1] = 0;
2312 else if(un<-_precision)
2314 double pi, Ej, rhoi, rhoj, ej;
2315 double ci;//vitesse du son pour calcul valeurs propres
2316 rhoj=_Uj[0]/_porosityj;
2317 Ej= _Uj[_Ndim+1]/rhoj;
2320 rhoi=_Ui[0]/_porosityi;
2321 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
2323 /***********Calcul des valeurs propres ********/
2324 vector<complex<double> > vp_dist(3,0);
2325 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2327 _maxvploc=fabs(uj_n)+ci;
2328 if(_maxvploc>_maxvp)
2331 if(_verbose && _nbTimeStep%_freqSave ==0)
2332 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2334 /******** Construction de la matrice A^+ *********/
2335 //premiere ligne (masse)
2336 _AroePlusImplicit[0]=0;
2337 for(int idim=0; idim<_Ndim;idim++)
2338 _AroePlusImplicit[1+idim]=0;
2339 _AroePlusImplicit[_nVar-1]=0;
2341 //lignes intermadiaires(qdm)
2342 for(int idim=0; idim<_Ndim;idim++)
2345 _AroePlusImplicit[(1+idim)*_nVar]=0;
2346 //colonnes intermediaires
2347 for(int jdim=0; jdim<_Ndim;jdim++)
2348 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2350 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2352 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2355 //derniere ligne (energie)
2356 _AroePlusImplicit[_nVar*(_nVar-1)] = uj_n;
2357 for(int idim=0; idim<_Ndim;idim++)
2358 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2359 _AroePlusImplicit[_nVar*_nVar -1] = 0;
2361 /******** Construction de la matrice A^- *********/
2362 //premiere ligne (masse)
2363 _AroeMinusImplicit[0]=uj_n/((gamma-1)*(ej-q));
2364 for(int idim=0; idim<_Ndim;idim++)
2365 _AroeMinusImplicit[1+idim]=rhoj*_vec_normal[idim];
2366 _AroeMinusImplicit[_nVar-1]=-rhoj*uj_n*cv/(ej-q);
2368 //lignes intermadiaires(qdm)
2369 for(int idim=0; idim<_Ndim;idim++)
2372 _AroeMinusImplicit[(1+idim)*_nVar]= uj_n/((gamma-1)*(ej-q))*_Vj[1+idim];
2373 //colonnes intermediaires
2374 for(int jdim=0; jdim<_Ndim;jdim++)
2375 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = rhoj*_Vj[1+idim]*_vec_normal[jdim];
2377 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += rhoj*uj_n;
2379 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoj*uj_n*cv/(ej-q)*_Vj[1+idim];
2382 //derniere ligne (energie)
2383 _AroeMinusImplicit[_nVar*(_nVar-1)] = Ej*uj_n/((gamma-1)*(ej-q));
2384 for(int idim=0; idim<_Ndim;idim++)
2385 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoj*Ej+pi)*_vec_normal[idim]+rhoj*uj_n*_Vj[1+idim];
2386 _AroeMinusImplicit[_nVar*_nVar -1] = rhoj*uj_n*(1-Ej/(ej-q))*cv;
2390 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero" << endl;
2391 _runLogFile->close();
2392 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero");
2395 else if(_useDellacherieEOS )
2397 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2398 double cp=fluide0->constante("cp");
2402 double rhoi,rhoj,pj, Ei, hi, Hi;
2403 double cj;//vitesse du son pour calcul valeurs propres
2404 rhoi=_Ui[0]/_porosityi;
2405 Ei= _Ui[_Ndim+1]/(rhoi*_porosityi);
2409 rhoj=_Uj[0]/_porosityj;
2410 cj = _fluides[0]->vitesseSonTemperature(_Vj[_Ndim+1],rhoj);
2412 /***********Calcul des valeurs propres ********/
2413 vector<complex<double> > vp_dist(3,0);
2414 vp_dist[0]=ui_n-cj;vp_dist[1]=ui_n;vp_dist[2]=ui_n+cj;
2416 _maxvploc=fabs(ui_n)+cj;
2417 if(_maxvploc>_maxvp)
2420 if(_verbose && _nbTimeStep%_freqSave ==0)
2421 cout<<"Valeurs propres "<<ui_n-cj<<" , "<<ui_n<<" , "<<ui_n+cj<<endl;
2423 /******** Construction de la matrice A^+ *********/
2424 //premiere ligne (masse)
2425 _AroePlusImplicit[0]=ui_n*gamma/((gamma-1)*(hi-q));
2426 for(int idim=0; idim<_Ndim;idim++)
2427 _AroePlusImplicit[1+idim]=rhoi*_vec_normal[idim];
2428 _AroePlusImplicit[_nVar-1]=-rhoi*ui_n*cp/(hi-q);
2430 //lignes intermadiaires(qdm)
2431 for(int idim=0; idim<_Ndim;idim++)
2434 _AroePlusImplicit[(1+idim)*_nVar]=ui_n*gamma/((gamma-1)*(hi-q))*_Vi[1+idim];
2435 //colonnes intermediaires
2436 for(int jdim=0; jdim<_Ndim;jdim++)
2437 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = rhoi*_Vi[1+idim]*_vec_normal[jdim];
2439 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += rhoi*ui_n;
2441 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoi*ui_n*cp/(hi-q)*_Vi[1+idim];
2444 //derniere ligne (energie)
2445 _AroePlusImplicit[_nVar*(_nVar-1)] = ui_n*(Hi*gamma/((gamma-1)*(hi-q))-1);
2446 for(int idim=0; idim<_Ndim;idim++)
2447 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoi*Ei+pj)*_vec_normal[idim]+rhoi*ui_n*_Vi[1+idim];
2448 _AroePlusImplicit[_nVar*_nVar -1] = rhoi*ui_n*(1-Hi/(hi-q))*cp;
2450 /******** Construction de la matrice A^- *********/
2451 //premiere ligne (masse)
2452 _AroeMinusImplicit[0]=0;
2453 for(int idim=0; idim<_Ndim;idim++)
2454 _AroeMinusImplicit[1+idim]=0;
2455 _AroeMinusImplicit[_nVar-1]=0;
2457 //lignes intermadiaires(qdm)
2458 for(int idim=0; idim<_Ndim;idim++)
2461 _AroeMinusImplicit[(1+idim)*_nVar]=_vec_normal[idim] ;
2462 //colonnes intermediaires
2463 for(int jdim=0; jdim<_Ndim;jdim++)
2464 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2466 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2468 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2471 //derniere ligne (energie)
2472 _AroeMinusImplicit[_nVar*(_nVar-1)] = ui_n;
2473 for(int idim=0; idim<_Ndim;idim++)
2474 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2475 _AroeMinusImplicit[_nVar*_nVar -1] = 0;
2477 else if(un<-_precision)
2479 double pi, Ej, rhoi,rhoj, Hj, hj;
2480 double ci;//vitesse du son pour calcul valeurs propres
2481 rhoj=_Uj[0]/_porosityj;
2482 Ej= _Uj[_Ndim+1]/rhoj;
2486 rhoi=_Ui[0]/_porosityi;
2487 ci = _fluides[0]->vitesseSonTemperature(_Vi[_Ndim+1],rhoi);
2489 /***********Calcul des valeurs propres ********/
2490 vector<complex<double> > vp_dist(3,0);
2491 vp_dist[0]=uj_n-ci;vp_dist[1]=uj_n;vp_dist[2]=uj_n+ci;
2493 _maxvploc=fabs(uj_n)+ci;
2494 if(_maxvploc>_maxvp)
2497 if(_verbose && _nbTimeStep%_freqSave ==0)
2498 cout<<"Valeurs propres "<<uj_n-ci<<" , "<<uj_n<<" , "<<uj_n+ci<<endl;
2500 /******** Construction de la matrice A^+ *********/
2501 //premiere ligne (masse)
2502 _AroePlusImplicit[0]=0;
2503 for(int idim=0; idim<_Ndim;idim++)
2504 _AroePlusImplicit[1+idim]=0;
2505 _AroePlusImplicit[_nVar-1]=0;
2507 //lignes intermadiaires(qdm)
2508 for(int idim=0; idim<_Ndim;idim++)
2511 _AroePlusImplicit[(1+idim)*_nVar]=0;
2512 //colonnes intermediaires
2513 for(int jdim=0; jdim<_Ndim;jdim++)
2514 _AroePlusImplicit[(1+idim)*_nVar + jdim + 1] = 0;
2516 _AroePlusImplicit[(1+idim)*_nVar + idim + 1] += 0;
2518 _AroePlusImplicit[(1+idim)*_nVar + _nVar-1]=0;
2521 //derniere ligne (energie)
2522 _AroePlusImplicit[_nVar*(_nVar-1)] = uj_n;
2523 for(int idim=0; idim<_Ndim;idim++)
2524 _AroePlusImplicit[_nVar*(_nVar-1)+idim+1]= 0;
2525 _AroePlusImplicit[_nVar*_nVar -1] = 0;
2527 /******** Construction de la matrice A^- *********/
2528 //premiere ligne (masse)
2529 _AroeMinusImplicit[0]=uj_n*gamma/((gamma-1)*(hj-q));
2530 for(int idim=0; idim<_Ndim;idim++)
2531 _AroeMinusImplicit[1+idim]=rhoj*_vec_normal[idim];
2532 _AroeMinusImplicit[_nVar-1]=-rhoj*uj_n*cp/(hj-q);
2534 //lignes intermadiaires(qdm)
2535 for(int idim=0; idim<_Ndim;idim++)
2538 _AroeMinusImplicit[(1+idim)*_nVar]= uj_n*gamma/((gamma-1)*(hj-q))*_Vj[1+idim];
2539 //colonnes intermediaires
2540 for(int jdim=0; jdim<_Ndim;jdim++)
2541 _AroeMinusImplicit[(1+idim)*_nVar + jdim + 1] = rhoj*_Vj[1+idim]*_vec_normal[jdim];
2543 _AroeMinusImplicit[(1+idim)*_nVar + idim + 1] += rhoj*uj_n;
2545 _AroeMinusImplicit[(1+idim)*_nVar + _nVar-1]=-rhoj*uj_n*cp/(hj-q)*_Vj[1+idim];
2548 //derniere ligne (energie)
2549 _AroeMinusImplicit[_nVar*(_nVar-1)] = uj_n*(Hj*gamma/((gamma-1)*(hj-q))-1);
2550 for(int idim=0; idim<_Ndim;idim++)
2551 _AroeMinusImplicit[_nVar*(_nVar-1)+idim+1]=(rhoj*Ej+pi)*_vec_normal[idim]+rhoj*uj_n*_Vj[1+idim];
2552 _AroeMinusImplicit[_nVar*_nVar -1] = rhoj*uj_n*(1-Hj/(hj-q))*cp;
2556 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero" << endl;
2557 _runLogFile->close();
2558 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: velocity un should be non zero");
2563 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
2564 _runLogFile->close();
2565 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
2568 else//case nil velocity on the interface, apply centered scheme
2571 primToConsJacobianMatrix(_Vj);
2572 Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
2573 primToConsJacobianMatrix(_Vi);
2574 Poly.matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
2578 void SinglePhase::applyVFRoeLowMachCorrections(bool isBord, string groupname)
2580 if(_nonLinearFormulation!=VFRoe)
2582 *_runLogFile<< "SinglePhase::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation" << endl;
2583 _runLogFile->close();
2584 throw CdmathException("SinglePhase::applyVFRoeLowMachCorrections: applyVFRoeLowMachCorrections method should be called only for VFRoe formulation");
2586 else//_nonLinearFormulation==VFRoe
2588 if(_spaceScheme==lowMach){
2590 for(int i=0;i<_Ndim;i++)
2591 u_2 += _Uroe[1+i]*_Uroe[1+i];
2592 double c = _maxvploc;//vitesse du son a l'interface
2593 double M=sqrt(u_2)/c;//Mach number
2594 _Vij[0]=M*_Vij[0]+(1-M)*(_Vi[0]+_Vj[0])/2;
2596 else if(_spaceScheme==pressureCorrection)
2597 {//order 1 : no correction, oarder 2 : correction everywhere, order 3 : correction only inside, orders 4 and 5 : special correction at boundaries
2598 if(_pressureCorrectionOrder==2 || (!isBord && _pressureCorrectionOrder==3) || (!isBord && _pressureCorrectionOrder==4) || (!isBord && _pressureCorrectionOrder==5) )
2600 double norm_uij=0, uij_n=0, ui_n=0, uj_n=0;
2601 for(int i=0;i<_Ndim;i++)
2603 norm_uij += _Uroe[1+i]*_Uroe[1+i];
2604 uij_n += _Uroe[1+i]*_vec_normal[i];
2605 ui_n += _Vi[1+i]*_vec_normal[i];
2606 uj_n += _Vj[1+i]*_vec_normal[i];
2608 norm_uij=sqrt(norm_uij);
2609 if(norm_uij>_precision)//avoid division by zero
2610 _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;
2612 _Vij[0]=(_Vi[0]+_Vj[0])/2 - _Uroe[0]*norm_uij*(uj_n-ui_n)/4;
2614 else if(_pressureCorrectionOrder==4 && isBord)
2616 else if(_pressureCorrectionOrder==5 && isBord)
2618 double g_n=0;//scalar product of gravity and normal vector
2619 for(int i=0;i<_Ndim;i++)
2620 g_n += _GravityField3d[i]*_vec_normal[i];
2621 _Vij[0]=_Vi[0]- _Ui[0]*g_n/_inv_dxi/2;
2624 else if(_spaceScheme==staggered)
2627 for(int i=0;i<_Ndim;i++)
2628 uij_n += _Uroe[1+i]*_vec_normal[i];
2630 if(uij_n>_precision){
2632 for(int i=0;i<_Ndim;i++)
2634 _Vij[_nVar-1]=_Vi[_nVar-1];
2636 else if(uij_n<-_precision){
2638 for(int i=0;i<_Ndim;i++)
2640 _Vij[_nVar-1]=_Vj[_nVar-1];
2643 _Vij[0]=(_Vi[0]+_Vi[0])/2;
2644 for(int i=0;i<_Ndim;i++)
2645 _Vij[1+i]=(_Vj[1+i]+_Vj[1+i])/2;
2646 _Vij[_nVar-1]=(_Vj[_nVar-1]+_Vj[_nVar-1])/2;
2649 primToCons(_Vij,0,_Uij,0);
2653 void SinglePhase::testConservation()
2655 double SUM, DELTA, x;
2657 for(int i=0; i<_nVar; i++)
2661 cout << "Masse totale (kg): ";
2665 cout << "Energie totale (J): ";
2667 cout << "Quantite de mouvement totale (kg.m.s^-1): ";
2673 for(int j=0; j<_Nmailles; j++)
2675 if(!_usePrimitiveVarsInNewton)
2676 VecGetValues(_conservativeVars, 1, &I, &x);//on recupere la valeur du champ
2678 VecGetValues(_primitiveVars, 1, &I, &x);//on recupere la valeur du champ
2679 SUM += x*_mesh.getCell(j).getMeasure();
2680 VecGetValues(_newtonVariation, 1, &I, &x);//on recupere la variation du champ
2681 DELTA += x*_mesh.getCell(j).getMeasure();
2684 if(fabs(SUM)>_precision)
2685 cout << SUM << ", variation relative: " << fabs(DELTA /SUM) << endl;
2687 cout << " a une somme quasi nulle, variation absolue: " << fabs(DELTA) << endl;
2691 void SinglePhase::getDensityDerivatives( double pressure, double temperature, double v2)
2693 double rho=_fluides[0]->getDensity(pressure,temperature);
2694 double gamma=_fluides[0]->constante("gamma");
2695 double q=_fluides[0]->constante("q");
2697 if( !_useDellacherieEOS)
2699 StiffenedGas* fluide0=dynamic_cast<StiffenedGas*>(_fluides[0]);
2700 double e = fluide0->getInternalEnergy(temperature);
2701 double cv=fluide0->constante("cv");
2704 _drho_sur_dp=1/((gamma-1)*(e-q));
2705 _drho_sur_dT=-rho*cv/(e-q);
2706 _drhoE_sur_dp=E/((gamma-1)*(e-q));
2707 _drhoE_sur_dT=rho*cv*(1-E/(e-q));
2709 else if(_useDellacherieEOS )
2711 StiffenedGasDellacherie* fluide0=dynamic_cast<StiffenedGasDellacherie*>(_fluides[0]);
2712 double h=fluide0->getEnthalpy(temperature);
2714 double cp=fluide0->constante("cp");
2716 _drho_sur_dp=gamma/((gamma-1)*(h-q));
2717 _drho_sur_dT=-rho*cp/(h-q);
2718 _drhoE_sur_dp=gamma*H/((gamma-1)*(h-q))-1;
2719 _drhoE_sur_dT=rho*cp*(1-H/(h-q));
2723 *_runLogFile<< "SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie" << endl;
2724 _runLogFile->close();
2725 throw CdmathException("SinglePhase::staggeredVFFCMatricesPrimitiveVariables: eos should be StiffenedGas or StiffenedGasDellacherie");
2728 if(_verbose && _nbTimeStep%_freqSave ==0)
2730 cout<<"_drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
2731 cout<<"_drhoE_sur_dp= "<<_drhoE_sur_dp<<", _drhoE_sur_dT= "<<_drhoE_sur_dT<<endl;
2734 void SinglePhase::save(){
2735 string prim(_path+"/SinglePhasePrim_");///Results
2736 string cons(_path+"/SinglePhaseCons_");
2737 string allFields(_path+"/");
2740 allFields+=_fileName;
2743 for (long i = 0; i < _Nmailles; i++){
2744 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
2745 for (int j = 0; j < _nVar; j++){
2747 VecGetValues(_primitiveVars,1,&Ii,&_VV(i,j));
2750 if(_saveConservativeField){
2751 for (long i = 0; i < _Nmailles; i++){
2752 // j = 0 : density; j = _nVar - 1 : energy j = 1,..,_nVar-2: momentum
2753 for (int j = 0; j < _nVar; j++){
2755 VecGetValues(_conservativeVars,1,&Ii,&_UU(i,j));
2758 _UU.setTime(_time,_nbTimeStep);
2760 _VV.setTime(_time,_nbTimeStep);
2762 // create mesh and component info
2763 if (_nbTimeStep ==0 || _restartWithNewFileName){
2764 string prim_suppress ="rm -rf "+prim+"_*";
2765 string cons_suppress ="rm -rf "+cons+"_*";
2767 system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
2768 system(cons_suppress.c_str());//Nettoyage des précédents calculs identiques
2770 if(_saveConservativeField){
2771 _UU.setInfoOnComponent(0,"Density_(kg/m^3)");
2772 _UU.setInfoOnComponent(1,"Momentum_x");// (kg/m^2/s)
2774 _UU.setInfoOnComponent(2,"Momentum_y");// (kg/m^2/s)
2776 _UU.setInfoOnComponent(3,"Momentum_z");// (kg/m^2/s)
2778 _UU.setInfoOnComponent(_nVar-1,"Energy_(J/m^3)");
2793 _VV.setInfoOnComponent(0,"Pressure_(Pa)");
2794 _VV.setInfoOnComponent(1,"Velocity_x_(m/s)");
2796 _VV.setInfoOnComponent(2,"Velocity_y_(m/s)");
2798 _VV.setInfoOnComponent(3,"Velocity_z_(m/s)");
2799 _VV.setInfoOnComponent(_nVar-1,"Temperature_(K)");
2815 // do not create mesh
2820 _VV.writeVTK(prim,false);
2823 _VV.writeMED(prim,false);
2829 if(_saveConservativeField){
2833 _UU.writeVTK(cons,false);
2836 _UU.writeMED(cons,false);
2844 if(_saveVelocity || _saveAllFields){
2845 for (long i = 0; i < _Nmailles; i++){
2846 // j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
2847 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
2849 int Ii = i*_nVar +1+j;
2850 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
2852 for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
2855 _Vitesse.setTime(_time,_nbTimeStep);
2856 if (_nbTimeStep ==0 || _restartWithNewFileName){
2857 _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
2858 _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
2859 _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
2864 _Vitesse.writeVTK(prim+"_Velocity");
2867 _Vitesse.writeMED(prim+"_Velocity");
2870 _Vitesse.writeCSV(prim+"_Velocity");
2878 _Vitesse.writeVTK(prim+"_Velocity",false);
2881 _Vitesse.writeMED(prim+"_Velocity",false);
2884 _Vitesse.writeCSV(prim+"_Velocity");
2892 double p,T,rho, h, vx,vy,vz;
2894 for (long i = 0; i < _Nmailles; i++){
2896 VecGetValues(_conservativeVars,1,&Ii,&rho);
2898 VecGetValues(_primitiveVars,1,&Ii,&p);
2899 Ii = i*_nVar +_nVar-1;
2900 VecGetValues(_primitiveVars,1,&Ii,&T);
2902 VecGetValues(_primitiveVars,1,&Ii,&vx);
2906 VecGetValues(_primitiveVars,1,&Ii,&vy);
2909 VecGetValues(_primitiveVars,1,&Ii,&vz);
2913 h = _fluides[0]->getEnthalpy(T,rho);
2927 _Enthalpy.setTime(_time,_nbTimeStep);
2928 _Density.setTime(_time,_nbTimeStep);
2929 _Pressure.setTime(_time,_nbTimeStep);
2930 _Temperature.setTime(_time,_nbTimeStep);
2931 _VitesseX.setTime(_time,_nbTimeStep);
2934 _VitesseY.setTime(_time,_nbTimeStep);
2936 _VitesseZ.setTime(_time,_nbTimeStep);
2938 if (_nbTimeStep ==0 || _restartWithNewFileName){
2942 _Enthalpy.writeVTK(allFields+"_Enthalpy");
2943 _Density.writeVTK(allFields+"_Density");
2944 _Pressure.writeVTK(allFields+"_Pressure");
2945 _Temperature.writeVTK(allFields+"_Temperature");
2946 _VitesseX.writeVTK(allFields+"_VelocityX");
2949 _VitesseY.writeVTK(allFields+"_VelocityY");
2951 _VitesseZ.writeVTK(allFields+"_VelocityZ");
2955 _Enthalpy.writeMED(allFields+"_Enthalpy");
2956 _Density.writeMED(allFields+"_Density");
2957 _Pressure.writeMED(allFields+"_Pressure");
2958 _Temperature.writeMED(allFields+"_Temperature");
2959 _VitesseX.writeMED(allFields+"_VelocityX");
2962 _VitesseY.writeMED(allFields+"_VelocityY");
2964 _VitesseZ.writeMED(allFields+"_VelocityZ");
2968 _Enthalpy.writeCSV(allFields+"_Enthalpy");
2969 _Density.writeCSV(allFields+"_Density");
2970 _Pressure.writeCSV(allFields+"_Pressure");
2971 _Temperature.writeCSV(allFields+"_Temperature");
2972 _VitesseX.writeCSV(allFields+"_VelocityX");
2975 _VitesseY.writeCSV(allFields+"_VelocityY");
2977 _VitesseZ.writeCSV(allFields+"_VelocityZ");
2986 _Enthalpy.writeVTK(allFields+"_Enthalpy",false);
2987 _Density.writeVTK(allFields+"_Density",false);
2988 _Pressure.writeVTK(allFields+"_Pressure",false);
2989 _Temperature.writeVTK(allFields+"_Temperature",false);
2990 _VitesseX.writeVTK(allFields+"_VelocityX",false);
2993 _VitesseY.writeVTK(allFields+"_VelocityY",false);
2995 _VitesseZ.writeVTK(allFields+"_VelocityZ",false);
2999 _Enthalpy.writeMED(allFields+"_Enthalpy",false);
3000 _Density.writeMED(allFields+"_Density",false);
3001 _Pressure.writeMED(allFields+"_Pressure",false);
3002 _Temperature.writeMED(allFields+"_Temperature",false);
3003 _VitesseX.writeMED(allFields+"_VelocityX",false);
3006 _VitesseY.writeMED(allFields+"_VelocityY",false);
3008 _VitesseZ.writeMED(allFields+"_VelocityZ",false);
3012 _Enthalpy.writeCSV(allFields+"_Enthalpy");
3013 _Density.writeCSV(allFields+"_Density");
3014 _Pressure.writeCSV(allFields+"_Pressure");
3015 _Temperature.writeCSV(allFields+"_Temperature");
3016 _VitesseX.writeCSV(allFields+"_VelocityX");
3019 _VitesseY.writeCSV(allFields+"_VelocityY");
3021 _VitesseZ.writeCSV(allFields+"_VelocityZ");
3046 if(_saveConservativeField){
3061 if(_saveVelocity || _saveAllFields){
3065 _Vitesse.writeVTK(prim+"_Velocity");
3068 _Vitesse.writeMED(prim+"_Velocity");
3071 _Vitesse.writeCSV(prim+"_Velocity");
3077 if (_restartWithNewFileName)
3078 _restartWithNewFileName=false;
3081 Field& SinglePhase::getPressureField()
3085 _Pressure=Field("Pressure",CELLS,_mesh,1);
3087 for (long i = 0; i < _Nmailles; i++){
3089 VecGetValues(_primitiveVars,1,&Ii,&_Pressure(i));
3091 _Pressure.setTime(_time,_nbTimeStep);
3096 Field& SinglePhase::getTemperatureField()
3100 _Temperature=Field("Temperature",CELLS,_mesh,1);
3102 for (long i = 0; i < _Nmailles; i++){
3103 Ii = i*_nVar +_nVar-1;
3104 VecGetValues(_primitiveVars,1,&Ii,&_Temperature(i));
3106 _Temperature.setTime(_time,_nbTimeStep);
3108 return _Temperature;
3111 Field& SinglePhase::getVelocityField()
3113 if(!_saveAllFields )
3115 _Vitesse=Field("Vitesse",CELLS,_mesh,3);
3117 for (long i = 0; i < _Nmailles; i++)
3119 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
3121 int Ii = i*_nVar +1+j;
3122 VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
3124 for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
3127 _Vitesse.setTime(_time,_nbTimeStep);
3128 _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
3129 _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
3130 _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
3136 Field& SinglePhase::getVelocityXField()
3138 if(!_saveAllFields )
3140 _VitesseX=Field("Velocity X",CELLS,_mesh,1);
3142 for (long i = 0; i < _Nmailles; i++)
3144 int Ii = i*_nVar +1;
3145 VecGetValues(_primitiveVars,1,&Ii,&_VitesseX(i));
3147 _VitesseX.setTime(_time,_nbTimeStep);
3148 _VitesseX.setInfoOnComponent(0,"Velocity_x_(m/s)");
3154 Field& SinglePhase::getDensityField()
3156 if(!_saveAllFields )
3158 _Density=Field("Density",CELLS,_mesh,1);
3160 for (long i = 0; i < _Nmailles; i++){
3162 VecGetValues(_conservativeVars,1,&Ii,&_Density(i));
3164 _Density.setTime(_time,_nbTimeStep);
3169 Field& SinglePhase::getMomentumField()//not yet managed by parameter _saveAllFields
3171 _Momentum=Field("Momentum",CELLS,_mesh,_Ndim);
3173 for (long i = 0; i < _Nmailles; i++)
3174 for (int j = 0; j < _Ndim; j++)//On récupère les composantes de qdm
3176 int Ii = i*_nVar +1+j;
3177 VecGetValues(_conservativeVars,1,&Ii,&_Momentum(i,j));
3179 _Momentum.setTime(_time,_nbTimeStep);
3184 Field& SinglePhase::getTotalEnergyField()//not yet managed by parameter _saveAllFields
3186 _TotalEnergy=Field("TotalEnergy",CELLS,_mesh,1);
3188 for (long i = 0; i < _Nmailles; i++){
3189 Ii = i*_nVar +_nVar-1;
3190 VecGetValues(_conservativeVars,1,&Ii,&_TotalEnergy(i));
3192 _TotalEnergy.setTime(_time,_nbTimeStep);
3194 return _TotalEnergy;
3197 Field& SinglePhase::getEnthalpyField()
3199 if(!_saveAllFields )
3201 _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
3204 for (long i = 0; i < _Nmailles; i++){
3206 VecGetValues(_primitiveVars,1,&Ii,&p);
3207 Ii = i*_nVar +_nVar-1;
3208 VecGetValues(_primitiveVars,1,&Ii,&T);
3210 rho=_fluides[0]->getDensity(p,T);
3211 _Enthalpy(i)=_fluides[0]->getEnthalpy(T,rho);
3213 _Enthalpy.setTime(_time,_nbTimeStep);
3219 vector<string> SinglePhase::getOutputFieldsNames()
3221 vector<string> result(8);
3223 result[0]="Pressure";
3224 result[1]="Velocity";
3225 result[2]="Temperature";
3226 result[3]="Density";
3227 result[4]="Momentum";
3228 result[5]="TotalEnergy";
3229 result[6]="Enthalpy";
3230 result[7]="VelocityX";
3235 Field& SinglePhase::getOutputField(const string& nameField )
3237 if(nameField=="pressure" || nameField=="Pressure" || nameField=="PRESSURE" || nameField=="PRESSION" || nameField=="Pression" || nameField=="pression" )
3238 return getPressureField();
3239 else if(nameField=="velocity" || nameField=="Velocity" || nameField=="VELOCITY" || nameField=="Vitesse" || nameField=="VITESSE" || nameField=="vitesse" )
3240 return getVelocityField();
3241 else if(nameField=="velocityX" || nameField=="VelocityX" || nameField=="VELOCITYX" || nameField=="VitesseX" || nameField=="VITESSEX" || nameField=="vitesseX" )
3242 return getVelocityXField();
3243 else if(nameField=="temperature" || nameField=="Temperature" || nameField=="TEMPERATURE" || nameField=="temperature" )
3244 return getTemperatureField();
3245 else if(nameField=="density" || nameField=="Density" || nameField=="DENSITY" || nameField=="Densite" || nameField=="DENSITE" || nameField=="densite" )
3246 return getDensityField();
3247 else if(nameField=="momentum" || nameField=="Momentum" || nameField=="MOMENTUM" || nameField=="Qdm" || nameField=="QDM" || nameField=="qdm" )
3248 return getMomentumField();
3249 else if(nameField=="enthalpy" || nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPIE" || nameField=="enthalpie" )
3250 return getEnthalpyField();
3251 else if(nameField=="totalenergy" || nameField=="TotalEnergy" || nameField=="TOTALENERGY" || nameField=="ENERGIETOTALE" || nameField=="EnergieTotale" || nameField=="energietotale" )
3252 return getTotalEnergyField();
3255 cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
3256 throw CdmathException("SinglePhase::getOutputField error : Unknown Field name");