Check memory is initialized in output field functions

[tools/solverlab.git] / CoreFlows / Models / src / SinglePhase.cxx

diff --git a/CoreFlows/Models/src/SinglePhase.cxx b/CoreFlows/Models/src/SinglePhase.cxx
index fd4dc87544237abf9370f6e10c6f32e0bdd0dea5..acf6566f731c9bf9d459049869d27a8711dd10b5 100755 (executable)
--- a/CoreFlows/Models/src/SinglePhase.cxx
+++ b/CoreFlows/Models/src/SinglePhase.cxx
@@ -69,6 +69,7 @@ void SinglePhase::initialize(){
                _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)
                {
@@ -2885,7 +2886,7 @@ void SinglePhase::save(){
 
        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;
@@ -2913,17 +2914,24 @@ void SinglePhase::save(){
                        _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)
                {
@@ -2939,6 +2947,7 @@ void SinglePhase::save(){
                                _Density.writeVTK(allFields+"_Density");
                                _Pressure.writeVTK(allFields+"_Pressure");
                                _Temperature.writeVTK(allFields+"_Temperature");
+                               _MachNumber.writeVTK(allFields+"_MachNumber");
                                _VitesseX.writeVTK(allFields+"_VelocityX");
                                if(_Ndim>1)
                                {
@@ -2952,6 +2961,7 @@ void SinglePhase::save(){
                                _Density.writeMED(allFields+"_Density");
                                _Pressure.writeMED(allFields+"_Pressure");
                                _Temperature.writeMED(allFields+"_Temperature");
+                               _MachNumber.writeMED(allFields+"_MachNumber");
                                _VitesseX.writeMED(allFields+"_VelocityX");
                                if(_Ndim>1)
                                {
@@ -2965,6 +2975,7 @@ void SinglePhase::save(){
                                _Density.writeCSV(allFields+"_Density");
                                _Pressure.writeCSV(allFields+"_Pressure");
                                _Temperature.writeCSV(allFields+"_Temperature");
+                               _MachNumber.writeCSV(allFields+"_MachNumber");
                                _VitesseX.writeCSV(allFields+"_VelocityX");
                                if(_Ndim>1)
                                {
@@ -2983,6 +2994,7 @@ void SinglePhase::save(){
                                _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)
                                {
@@ -2996,6 +3008,7 @@ void SinglePhase::save(){
                                _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)
                                {
@@ -3009,6 +3022,7 @@ void SinglePhase::save(){
                                _Density.writeCSV(allFields+"_Density");
                                _Pressure.writeCSV(allFields+"_Pressure");
                                _Temperature.writeCSV(allFields+"_Temperature");
+                               _MachNumber.writeCSV(allFields+"_MachNumber");
                                _VitesseX.writeCSV(allFields+"_VelocityX");
                                if(_Ndim>1)
                                {
@@ -3076,6 +3090,9 @@ void SinglePhase::save(){
 
 Field& SinglePhase::getPressureField()
 {
+       if(!_initializedMemory)
+               throw CdmathException("SinglePhase::getPressureField, Call initialize first");
+
        if(!_saveAllFields)
        {
                _Pressure=Field("Pressure",CELLS,_mesh,1);
@@ -3091,6 +3108,9 @@ Field& SinglePhase::getPressureField()
 
 Field& SinglePhase::getTemperatureField()
 {
+       if(!_initializedMemory)
+               throw CdmathException("SinglePhase::getTemperatureField, Call initialize first");
+
        if(!_saveAllFields)
        {
                _Temperature=Field("Temperature",CELLS,_mesh,1);
@@ -3106,6 +3126,9 @@ Field& SinglePhase::getTemperatureField()
 
 Field& SinglePhase::getVelocityField()
 {
+       if(!_initializedMemory)
+               throw CdmathException("SinglePhase::getVelocityField, Call initialize first");
+
        if(!_saveAllFields )
        {
                _Vitesse=Field("Vitesse",CELLS,_mesh,3);
@@ -3129,8 +3152,46 @@ Field& SinglePhase::getVelocityField()
        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);
@@ -3147,8 +3208,59 @@ Field& SinglePhase::getVelocityXField()
        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);
@@ -3164,6 +3276,9 @@ Field& SinglePhase::getDensityField()
 
 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++)
@@ -3179,6 +3294,9 @@ Field& SinglePhase::getMomentumField()//not yet managed by parameter _saveAllFie
 
 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++){
@@ -3192,6 +3310,9 @@ Field& SinglePhase::getTotalEnergyField()//not yet managed by parameter _saveAll
 
 Field& SinglePhase::getEnthalpyField()
 {
+       if(!_initializedMemory)
+               throw CdmathException("SinglePhase::getEnthalpyField, Call initialize first");
+
        if(!_saveAllFields )
        {
                _Enthalpy=Field("Enthalpy",CELLS,_mesh,1);
@@ -3236,6 +3357,10 @@ Field& SinglePhase::getOutputField(const string& nameField )
                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" )
@@ -3248,7 +3373,7 @@ Field& SinglePhase::getOutputField(const string& nameField )
                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");
     }
 }