_dragCoeffs=vector<double>(1,0);
_fluides.resize(1);
_useDellacherieEOS=useDellacherieEOS;
+ _saveAllFields=false;
if(pEstimate==around1bar300K){//EOS at 1 bar and 300K
_Tref=300;
_gravite[i+1]=_GravityField3d[i];
_GravityImplicitationMatrix = new PetscScalar[_nVar*_nVar];
- if(_saveVelocity)
+ if(_saveVelocity || _saveAllFields)
_Vitesse=Field("Velocity",CELLS,_mesh,3);//Forcement en dimension 3 pour le posttraitement des lignes de courant
+ if(_saveAllFields)
+ {
+ _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
+ _Pressure=Field("Pressure",CELLS,_mesh,1);
+ _Density=Field("Density",CELLS,_mesh,1);
+ _Temperature=Field("Temperature",CELLS,_mesh,1);
+ _MachNumber=Field("MachNumber",CELLS,_mesh,1);
+ _VitesseX=Field("Velocity x",CELLS,_mesh,1);
+ if(_Ndim>1)
+ {
+ _VitesseY=Field("Velocity y",CELLS,_mesh,1);
+ if(_Ndim>2)
+ _VitesseZ=Field("Velocity z",CELLS,_mesh,1);
+ }
+ }
+
if(_entropicCorrection)
_entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
bool SinglePhase::iterateTimeStep(bool &converged)
{
if(_timeScheme == Explicit || !_usePrimitiveVarsInNewton)
- ProblemFluid::iterateTimeStep(converged);
+ return ProblemFluid::iterateTimeStep(converged);
else
{
bool stop=false;
}
}
void SinglePhase::convectionState( const long &i, const long &j, const bool &IsBord){
-
+ //First conservative state then further down we will compute interface (Roe) state and then compute primitive state
_idm[0] = _nVar*i;
for(int k=1; k<_nVar; k++)
_idm[k] = _idm[k-1] + 1;
_runLogFile->close();
throw CdmathException("densite negative, arret de calcul");
}
+
+ //Computation of Roe state
PetscScalar ri, rj, xi, xj, pi, pj;
PetscInt Ii;
ri = sqrt(_Ui[0]);//racine carre de phi_i rho_i
for(int k=0;k<_nVar;k++)
cout<< _Uroe[k]<<" , "<<endl;
}
+
+ //Extraction of primitive states
+ _idm[0] = _nVar*i; // Kieu
+ for(int k=1; k<_nVar; k++)
+ _idm[k] = _idm[k-1] + 1;
+
+ VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
+ if (_verbose && _nbTimeStep%_freqSave ==0)
+ {
+ cout << "Convection state: variables primitives maille " << i<<endl;
+ for(int q=0; q<_nVar; q++)
+ cout << _Vi[q] << endl;
+ cout << endl;
+ }
+
+ if(!IsBord ){
+ for(int k=0; k<_nVar; k++)
+ _idn[k] = _nVar*j + k;
+
+ VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
+ }
+ else
+ {
+ for(int k=0; k<_nVar; k++)
+ _idn[k] = k;
+
+ VecGetValues(_Uextdiff, _nVar, _idn, _phi);
+ consToPrim(_phi,_Vj,1);
+ }
+
+ if (_verbose && _nbTimeStep%_freqSave ==0)
+ {
+ cout << "Convection state: variables primitives maille " <<j <<endl;
+ for(int q=0; q<_nVar; q++)
+ cout << _Vj[q] << endl;
+ cout << endl;
+ }
}
void SinglePhase::diffusionStateAndMatrices(const long &i,const long &j, const bool &IsBord){
_idm[k] = _idm[k-1] + 1;
VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
- _idm[0] = _nVar*j;
- for(int k=1; k<_nVar; k++)
- _idm[k] = _idm[k-1] + 1;
+ for(int k=0; k<_nVar; k++)
+ _idn[k] = k;
if(IsBord)
- VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
+ VecGetValues(_Uextdiff, _nVar, _idn, _Uj);
else
VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
}
_Diffusion[_nVar*_nVar-1]=-lambda/(_Udiff[0]*Cv);
}
-
}
void SinglePhase::setBoundaryState(string nameOfGroup, const int &j,double *normale){
_idm[0] = _nVar*j;
_entropicShift=vector<double>(3,0);//at most 3 distinct eigenvalues
vector< complex< double > > y (3,0);
- Polynoms Poly;
for( int i=0 ; i<3 ; i++)
- y[i] = Poly.abs_generalise(vp_dist[i])+_entropicShift[i];
- Poly.abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
+ y[i] = Polynoms::abs_generalise(vp_dist[i])+_entropicShift[i];
+ Polynoms::abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
if( _spaceScheme ==pressureCorrection)
for( int i=0 ; i<_Ndim ; i++)
for(int idim=0;idim<_Ndim; idim++)
_Vij[1+idim]=_Uroe[1+idim];
primToConsJacobianMatrix(_Vij);
- Polynoms Poly;
- Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
- Poly.matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
+ Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
+ Polynoms::matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
}
else
for(int i=0; i<_nVar*_nVar;i++)
if(_entropicCorrection)
{
InvMatriceRoe( vp_dist);
- Polynoms Poly;
- Poly.matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
+ Polynoms::matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
}
else if (_spaceScheme==upwind || (_spaceScheme ==pressureCorrection ) || (_spaceScheme ==lowMach ))//upwind sans entropic
SigneMatriceRoe( vp_dist);
bool isBord)
{
- double lambda=_fluides[0]->getConductivity(_Udiff[_nVar-1]);
- double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
+ double lambda = _fluides[0]->getConductivity(_Udiff[_nVar-1]);
+ double mu = _fluides[0]->getViscosity(_Udiff[_nVar-1]);
+
+ if(isBord )
+ lambda=max(lambda,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
if(lambda==0 && mu ==0 && _heatTransfertCoeff==0)
return;
- //extraction des valeurs
- _idm[0] = _nVar*i; // Kieu
+ _idm[0] = 0;
for(int k=1; k<_nVar; k++)
_idm[k] = _idm[k-1] + 1;
-
- VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if(isBord)
{
- cout << "Calcul diffusion: variables primitives maille " << i<<endl;
- for(int q=0; q<_nVar; q++)
- cout << _Vi[q] << endl;
- cout << endl;
- }
-
- if(!isBord ){
- for(int k=0; k<_nVar; k++)
- _idn[k] = _nVar*j + k;
-
- VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
- }
- else
- {
- lambda=max(lambda,_heatTransfertCoeff);//wall nucleate boing -> larger heat transfer
- for(int k=0; k<_nVar; k++)
- _idn[k] = k;
-
- VecGetValues(_Uextdiff, _nVar, _idn, _phi);
- consToPrim(_phi,_Vj,1);
+ VecGetValues(_Uextdiff, _nVar, _idm, _Uj);
+ _inv_dxj=_inv_dxi;
}
- if (_verbose && _nbTimeStep%_freqSave ==0)
- {
- cout << "Calcul diffusion: variables primitives maille " <<j <<endl;
- for(int q=0; q<_nVar; q++)
- cout << _Vj[q] << endl;
- cout << endl;
- }
//on n'a pas de contribution sur la masse
_phi[0]=0;
//contribution visqueuse sur la quantite de mouvement
for(int k=1; k<_nVar-1; k++)
_phi[k] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*mu*(_porosityj*_Vj[k] - _porosityi*_Vi[k]);
-
//contribution visqueuse sur l'energie
_phi[_nVar-1] = _inv_dxi*2/(1/_inv_dxi+1/_inv_dxj)*lambda*(_porosityj*_Vj[_nVar-1] - _porosityi*_Vi[_nVar-1]);
-
+
_idm[0] = i;
VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
//On change de signe pour l'autre contribution
for(int k=0; k<_nVar; k++)
_phi[k] *= -_inv_dxj/_inv_dxi;
- _idn[0] = j;
+ _idn[0] = j;
VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
if(_verbose && _nbTimeStep%_freqSave ==0)
{
ur_2 += _Vj[1+i]*_Vj[1+i];
}
-
double cl = _fluides[0]->vitesseSonEnthalpie(_Vi[_Ndim+1]-ul_2/2);//vitesse du son a l'interface
double cr = _fluides[0]->vitesseSonEnthalpie(_Vj[_Ndim+1]-ur_2/2);//vitesse du son a l'interface
}
else//case nil velocity on the interface, apply centered scheme
{
- Polynoms Poly;
primToConsJacobianMatrix(_Vj);
- Poly.matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
+ Polynoms::matrixProduct(_AroeMinus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroeMinusImplicit);
primToConsJacobianMatrix(_Vi);
- Poly.matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
+ Polynoms::matrixProduct(_AroePlus, _nVar, _nVar, _primToConsJacoMat, _nVar, _nVar, _AroePlusImplicit);
}
}
}
void SinglePhase::save(){
string prim(_path+"/SinglePhasePrim_");///Results
string cons(_path+"/SinglePhaseCons_");
+ string allFields(_path+"/");
prim+=_fileName;
cons+=_fileName;
+ allFields+=_fileName;
PetscInt Ii;
for (long i = 0; i < _Nmailles; i++){
_VV.setTime(_time,_nbTimeStep);
// create mesh and component info
- if (_nbTimeStep ==0){
+ if (_nbTimeStep ==0 || _restartWithNewFileName){
string prim_suppress ="rm -rf "+prim+"_*";
string cons_suppress ="rm -rf "+cons+"_*";
_VV.writeCSV(prim);
break;
}
+
}
// do not create mesh
else{
}
}
}
- if(_saveVelocity){
+ if(_saveVelocity || _saveAllFields){
for (long i = 0; i < _Nmailles; i++){
// j = 0 : pressure; j = _nVar - 1: temperature; j = 1,..,_nVar-2: velocity
for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
_Vitesse(i,j)=0;
}
_Vitesse.setTime(_time,_nbTimeStep);
- if (_nbTimeStep ==0){
+ if (_nbTimeStep ==0 || _restartWithNewFileName){
_Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
_Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
_Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
}
}
}
+
+ if(_saveAllFields)
+ {
+ double p,T,rho, h, vx,vy,vz,v2;
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++){
+ Ii = i*_nVar;
+ VecGetValues(_conservativeVars,1,&Ii,&rho);
+ Ii = i*_nVar;
+ VecGetValues(_primitiveVars,1,&Ii,&p);
+ Ii = i*_nVar +_nVar-1;
+ VecGetValues(_primitiveVars,1,&Ii,&T);
+ Ii = i*_nVar + 1;
+ VecGetValues(_primitiveVars,1,&Ii,&vx);
+ if(_Ndim>1)
+ {
+ Ii = i*_nVar + 2;
+ VecGetValues(_primitiveVars,1,&Ii,&vy);
+ if(_Ndim>2){
+ Ii = i*_nVar + 3;
+ VecGetValues(_primitiveVars,1,&Ii,&vz);
+ }
+ }
+
+ h = _fluides[0]->getEnthalpy(T,rho);
+
+ _Enthalpy(i)=h;
+ _Density(i)=rho;
+ _Pressure(i)=p;
+ _Temperature(i)=T;
+ _VitesseX(i)=vx;
+ v2=vx*vx;
+ if(_Ndim>1)
+ {
+ _VitesseY(i)=vy;
+ v2+=vy*vy;
+ if(_Ndim>2)
+ {
+ _VitesseZ(i)=vz;
+ v2+=vz*vz;
+ }
+ }
+ _MachNumber(i)=sqrt(v2)/_fluides[0]->vitesseSonEnthalpie(h);
+ }
+ _Enthalpy.setTime(_time,_nbTimeStep);
+ _Density.setTime(_time,_nbTimeStep);
+ _Pressure.setTime(_time,_nbTimeStep);
+ _Temperature.setTime(_time,_nbTimeStep);
+ _MachNumber.setTime(_time,_nbTimeStep);
+ _VitesseX.setTime(_time,_nbTimeStep);
+ if(_Ndim>1)
+ {
+ _VitesseY.setTime(_time,_nbTimeStep);
+ if(_Ndim>2)
+ _VitesseZ.setTime(_time,_nbTimeStep);
+ }
+ if (_nbTimeStep ==0 || _restartWithNewFileName){
+ switch(_saveFormat)
+ {
+ case VTK :
+ _Enthalpy.writeVTK(allFields+"_Enthalpy");
+ _Density.writeVTK(allFields+"_Density");
+ _Pressure.writeVTK(allFields+"_Pressure");
+ _Temperature.writeVTK(allFields+"_Temperature");
+ _MachNumber.writeVTK(allFields+"_MachNumber");
+ _VitesseX.writeVTK(allFields+"_VelocityX");
+ if(_Ndim>1)
+ {
+ _VitesseY.writeVTK(allFields+"_VelocityY");
+ if(_Ndim>2)
+ _VitesseZ.writeVTK(allFields+"_VelocityZ");
+ }
+ break;
+ case MED :
+ _Enthalpy.writeMED(allFields+"_Enthalpy");
+ _Density.writeMED(allFields+"_Density");
+ _Pressure.writeMED(allFields+"_Pressure");
+ _Temperature.writeMED(allFields+"_Temperature");
+ _MachNumber.writeMED(allFields+"_MachNumber");
+ _VitesseX.writeMED(allFields+"_VelocityX");
+ if(_Ndim>1)
+ {
+ _VitesseY.writeMED(allFields+"_VelocityY");
+ if(_Ndim>2)
+ _VitesseZ.writeMED(allFields+"_VelocityZ");
+ }
+ break;
+ case CSV :
+ _Enthalpy.writeCSV(allFields+"_Enthalpy");
+ _Density.writeCSV(allFields+"_Density");
+ _Pressure.writeCSV(allFields+"_Pressure");
+ _Temperature.writeCSV(allFields+"_Temperature");
+ _MachNumber.writeCSV(allFields+"_MachNumber");
+ _VitesseX.writeCSV(allFields+"_VelocityX");
+ if(_Ndim>1)
+ {
+ _VitesseY.writeCSV(allFields+"_VelocityY");
+ if(_Ndim>2)
+ _VitesseZ.writeCSV(allFields+"_VelocityZ");
+ }
+ break;
+ }
+ }
+ else{
+ switch(_saveFormat)
+ {
+ case VTK :
+ _Enthalpy.writeVTK(allFields+"_Enthalpy",false);
+ _Density.writeVTK(allFields+"_Density",false);
+ _Pressure.writeVTK(allFields+"_Pressure",false);
+ _Temperature.writeVTK(allFields+"_Temperature",false);
+ _MachNumber.writeVTK(allFields+"_MachNumber",false);
+ _VitesseX.writeVTK(allFields+"_VelocityX",false);
+ if(_Ndim>1)
+ {
+ _VitesseY.writeVTK(allFields+"_VelocityY",false);
+ if(_Ndim>2)
+ _VitesseZ.writeVTK(allFields+"_VelocityZ",false);
+ }
+ break;
+ case MED :
+ _Enthalpy.writeMED(allFields+"_Enthalpy",false);
+ _Density.writeMED(allFields+"_Density",false);
+ _Pressure.writeMED(allFields+"_Pressure",false);
+ _Temperature.writeMED(allFields+"_Temperature",false);
+ _MachNumber.writeMED(allFields+"_MachNumber",false);
+ _VitesseX.writeMED(allFields+"_VelocityX",false);
+ if(_Ndim>1)
+ {
+ _VitesseY.writeMED(allFields+"_VelocityY",false);
+ if(_Ndim>2)
+ _VitesseZ.writeMED(allFields+"_VelocityZ",false);
+ }
+ break;
+ case CSV :
+ _Enthalpy.writeCSV(allFields+"_Enthalpy");
+ _Density.writeCSV(allFields+"_Density");
+ _Pressure.writeCSV(allFields+"_Pressure");
+ _Temperature.writeCSV(allFields+"_Temperature");
+ _MachNumber.writeCSV(allFields+"_MachNumber");
+ _VitesseX.writeCSV(allFields+"_VelocityX");
+ if(_Ndim>1)
+ {
+ _VitesseY.writeCSV(allFields+"_VelocityY");
+ if(_Ndim>2)
+ _VitesseZ.writeCSV(allFields+"_VelocityZ");
+ }
+ break;
+ }
+ }
+ }
+
if(_isStationary)
{
prim+="_Stat";
}
}
- if(_saveVelocity){
+ if(_saveVelocity || _saveAllFields){
switch(_saveFormat)
{
case VTK :
}
}
}
+
+ if (_restartWithNewFileName)
+ _restartWithNewFileName=false;
+}
+
+Field& SinglePhase::getPressureField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getPressureField, Call initialize first");
+
+ if(!_saveAllFields)
+ {
+ _Pressure=Field("Pressure",CELLS,_mesh,1);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++){
+ Ii = i*_nVar;
+ VecGetValues(_primitiveVars,1,&Ii,&_Pressure(i));
+ }
+ _Pressure.setTime(_time,_nbTimeStep);
+ }
+ return _Pressure;
+}
+
+Field& SinglePhase::getTemperatureField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getTemperatureField, Call initialize first");
+
+ if(!_saveAllFields)
+ {
+ _Temperature=Field("Temperature",CELLS,_mesh,1);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++){
+ Ii = i*_nVar +_nVar-1;
+ VecGetValues(_primitiveVars,1,&Ii,&_Temperature(i));
+ }
+ _Temperature.setTime(_time,_nbTimeStep);
+ }
+ return _Temperature;
+}
+
+Field& SinglePhase::getVelocityField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getVelocityField, Call initialize first");
+
+ if(!_saveAllFields )
+ {
+ _Vitesse=Field("Vitesse",CELLS,_mesh,3);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++)
+ {
+ for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
+ {
+ int Ii = i*_nVar +1+j;
+ VecGetValues(_primitiveVars,1,&Ii,&_Vitesse(i,j));
+ }
+ for (int j = _Ndim; j < 3; j++)//On met à zero les composantes de vitesse si la dimension est <3
+ _Vitesse(i,j)=0;
+ }
+ _Vitesse.setTime(_time,_nbTimeStep);
+ _Vitesse.setInfoOnComponent(0,"Velocity_x_(m/s)");
+ _Vitesse.setInfoOnComponent(1,"Velocity_y_(m/s)");
+ _Vitesse.setInfoOnComponent(2,"Velocity_z_(m/s)");
+ }
+
+ return _Vitesse;
+}
+
+Field& SinglePhase::getMachNumberField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getMachNumberField, Call initialize first");
+
+ if(!_saveAllFields )
+ {
+ _MachNumber=Field("Mach number",CELLS,_mesh,1);
+ int Ii;
+ double p,T,rho,h, temp, u2=0;
+ for (long i = 0; i < _Nmailles; i++){
+ Ii = i*_nVar;
+ VecGetValues(_primitiveVars,1,&Ii,&p);
+ Ii = i*_nVar +_nVar-1;
+ VecGetValues(_primitiveVars,1,&Ii,&T);
+
+ for (int j = 0; j < _Ndim; j++)//On récupère les composantes de vitesse
+ {
+ int Ii = i*_nVar +1+j;
+ VecGetValues(_primitiveVars,1,&Ii,&temp);
+ u2+=temp*temp;
+ }
+
+ rho=_fluides[0]->getDensity(p,T);
+ h =_fluides[0]->getEnthalpy(T,rho);
+ _MachNumber[i] =sqrt(u2)/_fluides[0]->vitesseSonEnthalpie(h);
+ //cout<<"u="<<sqrt(u2)<<", c= "<<_fluides[0]->vitesseSonEnthalpie(h)<<", MachNumberField[i] = "<<MachNumberField[i] <<endl;
+ }
+ _MachNumber.setTime(_time,_nbTimeStep);
+ }
+ //cout<<", MachNumberField = "<<MachNumberField <<endl;
+
+ return _MachNumber;
+}
+
+Field& SinglePhase::getVelocityXField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getVelocityXField, Call initialize first");
+
+ if(!_saveAllFields )
+ {
+ _VitesseX=Field("Velocity X",CELLS,_mesh,1);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++)
+ {
+ int Ii = i*_nVar +1;
+ VecGetValues(_primitiveVars,1,&Ii,&_VitesseX(i));
+ }
+ _VitesseX.setTime(_time,_nbTimeStep);
+ _VitesseX.setInfoOnComponent(0,"Velocity_x_(m/s)");
+ }
+
+ return _VitesseX;
+}
+
+Field& SinglePhase::getVelocityYField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getVelocityYField, Call initialize first");
+
+ if(_Ndim<2)
+ throw CdmathException("SinglePhase::getVelocityYField() error : dimension should be at least 2");
+ else
+ if(!_saveAllFields )
+ {
+ _VitesseY=Field("Velocity Y",CELLS,_mesh,1);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++)
+ {
+ int Ii = i*_nVar +2;
+ VecGetValues(_primitiveVars,1,&Ii,&_VitesseY(i));
+ }
+ _VitesseY.setTime(_time,_nbTimeStep);
+ _VitesseY.setInfoOnComponent(0,"Velocity_y_(m/s)");
+ }
+
+ return _VitesseY;
+}
+
+Field& SinglePhase::getVelocityZField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getVelocityZField, Call initialize first");
+
+ if(_Ndim<3)
+ throw CdmathException("SinglePhase::getvelocityZField() error : dimension should be 3");
+ else
+ if(!_saveAllFields )
+ {
+ _VitesseZ=Field("Velocity Z",CELLS,_mesh,1);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++)
+ {
+ int Ii = i*_nVar +3;
+ VecGetValues(_primitiveVars,1,&Ii,&_VitesseZ(i));
+ }
+ _VitesseZ.setTime(_time,_nbTimeStep);
+ _VitesseZ.setInfoOnComponent(0,"Velocity_z_(m/s)");
+ }
+
+ return _VitesseZ;
+}
+
+Field& SinglePhase::getDensityField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getDensityField, Call initialize first");
+
+ if(!_saveAllFields )
+ {
+ _Density=Field("Density",CELLS,_mesh,1);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++){
+ Ii = i*_nVar;
+ VecGetValues(_conservativeVars,1,&Ii,&_Density(i));
+ }
+ _Density.setTime(_time,_nbTimeStep);
+ }
+ return _Density;
+}
+
+Field& SinglePhase::getMomentumField()//not yet managed by parameter _saveAllFields
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getMomentumField, Call initialize first");
+
+ _Momentum=Field("Momentum",CELLS,_mesh,_Ndim);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++)
+ for (int j = 0; j < _Ndim; j++)//On récupère les composantes de qdm
+ {
+ int Ii = i*_nVar +1+j;
+ VecGetValues(_conservativeVars,1,&Ii,&_Momentum(i,j));
+ }
+ _Momentum.setTime(_time,_nbTimeStep);
+
+ return _Momentum;
+}
+
+Field& SinglePhase::getTotalEnergyField()//not yet managed by parameter _saveAllFields
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getTotalEnergyField, Call initialize first");
+
+ _TotalEnergy=Field("TotalEnergy",CELLS,_mesh,1);
+ int Ii;
+ for (long i = 0; i < _Nmailles; i++){
+ Ii = i*_nVar +_nVar-1;
+ VecGetValues(_conservativeVars,1,&Ii,&_TotalEnergy(i));
+ }
+ _TotalEnergy.setTime(_time,_nbTimeStep);
+
+ return _TotalEnergy;
+}
+
+Field& SinglePhase::getEnthalpyField()
+{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getEnthalpyField, Call initialize first");
+
+ if(!_saveAllFields )
+ {
+ _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
+ int Ii;
+ double p,T,rho;
+ for (long i = 0; i < _Nmailles; i++){
+ Ii = i*_nVar;
+ VecGetValues(_primitiveVars,1,&Ii,&p);
+ Ii = i*_nVar +_nVar-1;
+ VecGetValues(_primitiveVars,1,&Ii,&T);
+
+ rho=_fluides[0]->getDensity(p,T);
+ _Enthalpy(i)=_fluides[0]->getEnthalpy(T,rho);
+ }
+ _Enthalpy.setTime(_time,_nbTimeStep);
+ }
+
+ return _Enthalpy;
+}
+
+vector<string> SinglePhase::getOutputFieldsNames()
+{
+ vector<string> result(8);
+
+ result[0]="Pressure";
+ result[1]="Velocity";
+ result[2]="Temperature";
+ result[3]="Density";
+ result[4]="Momentum";
+ result[5]="TotalEnergy";
+ result[6]="Enthalpy";
+ result[7]="VelocityX";
+
+ return result;
+}
+
+Field& SinglePhase::getOutputField(const string& nameField )
+{
+ if(nameField=="pressure" || nameField=="Pressure" || nameField=="PRESSURE" || nameField=="PRESSION" || nameField=="Pression" || nameField=="pression" )
+ return getPressureField();
+ else if(nameField=="velocity" || nameField=="Velocity" || nameField=="VELOCITY" || nameField=="Vitesse" || nameField=="VITESSE" || nameField=="vitesse" )
+ return getVelocityField();
+ else if(nameField=="velocityX" || nameField=="VelocityX" || nameField=="VELOCITYX" || nameField=="VitesseX" || nameField=="VITESSEX" || nameField=="vitesseX" )
+ return getVelocityXField();
+ else if(nameField=="velocityY" || nameField=="VelocityY" || nameField=="VELOCITYY" || nameField=="VitesseY" || nameField=="VITESSEY" || nameField=="vitesseY" )
+ return getVelocityYField();
+ else if(nameField=="velocityZ" || nameField=="VelocityZ" || nameField=="VELOCITYZ" || nameField=="VitesseZ" || nameField=="VITESSEZ" || nameField=="vitesseZ" )
+ return getVelocityZField();
+ else if(nameField=="temperature" || nameField=="Temperature" || nameField=="TEMPERATURE" || nameField=="temperature" )
+ return getTemperatureField();
+ else if(nameField=="density" || nameField=="Density" || nameField=="DENSITY" || nameField=="Densite" || nameField=="DENSITE" || nameField=="densite" )
+ return getDensityField();
+ else if(nameField=="momentum" || nameField=="Momentum" || nameField=="MOMENTUM" || nameField=="Qdm" || nameField=="QDM" || nameField=="qdm" )
+ return getMomentumField();
+ else if(nameField=="enthalpy" || nameField=="Enthalpy" || nameField=="ENTHALPY" || nameField=="Enthalpie" || nameField=="ENTHALPIE" || nameField=="enthalpie" )
+ return getEnthalpyField();
+ else if(nameField=="totalenergy" || nameField=="TotalEnergy" || nameField=="TOTALENERGY" || nameField=="ENERGIETOTALE" || nameField=="EnergieTotale" || nameField=="energietotale" )
+ return getTotalEnergyField();
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first to check" << endl;
+ throw CdmathException("SinglePhase::getOutputField error : Unknown Field name");
+ }
}
