_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)
{
}
}
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);
}
}
}
if(_saveAllFields)
{
- double p,T,rho, h, vx,vy,vz;
+ double p,T,rho, h, vx,vy,vz,v2;
int Ii;
for (long i = 0; i < _Nmailles; i++){
Ii = i*_nVar;
_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)
{
_Density.writeVTK(allFields+"_Density");
_Pressure.writeVTK(allFields+"_Pressure");
_Temperature.writeVTK(allFields+"_Temperature");
+ _MachNumber.writeVTK(allFields+"_MachNumber");
_VitesseX.writeVTK(allFields+"_VelocityX");
if(_Ndim>1)
{
_Density.writeMED(allFields+"_Density");
_Pressure.writeMED(allFields+"_Pressure");
_Temperature.writeMED(allFields+"_Temperature");
+ _MachNumber.writeMED(allFields+"_MachNumber");
_VitesseX.writeMED(allFields+"_VelocityX");
if(_Ndim>1)
{
_Density.writeCSV(allFields+"_Density");
_Pressure.writeCSV(allFields+"_Pressure");
_Temperature.writeCSV(allFields+"_Temperature");
+ _MachNumber.writeCSV(allFields+"_MachNumber");
_VitesseX.writeCSV(allFields+"_VelocityX");
if(_Ndim>1)
{
_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)
{
_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)
{
_Density.writeCSV(allFields+"_Density");
_Pressure.writeCSV(allFields+"_Pressure");
_Temperature.writeCSV(allFields+"_Temperature");
+ _MachNumber.writeCSV(allFields+"_MachNumber");
_VitesseX.writeCSV(allFields+"_VelocityX");
if(_Ndim>1)
{
Field& SinglePhase::getPressureField()
{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getPressureField, Call initialize first");
+
if(!_saveAllFields)
{
_Pressure=Field("Pressure",CELLS,_mesh,1);
Field& SinglePhase::getTemperatureField()
{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getTemperatureField, Call initialize first");
+
if(!_saveAllFields)
{
_Temperature=Field("Temperature",CELLS,_mesh,1);
Field& SinglePhase::getVelocityField()
{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getVelocityField, Call initialize first");
+
if(!_saveAllFields )
{
_Vitesse=Field("Vitesse",CELLS,_mesh,3);
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);
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);
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++)
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++){
Field& SinglePhase::getEnthalpyField()
{
+ if(!_initializedMemory)
+ throw CdmathException("SinglePhase::getEnthalpyField, Call initialize first");
+
if(!_saveAllFields )
{
_Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
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 getTotalEnergyField();
else
{
- cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first" << endl;
+ cout<<"Error : Field name "<< nameField << " does not exist, call getOutputFieldsNames first to check" << endl;
throw CdmathException("SinglePhase::getOutputField error : Unknown Field name");
}
}
