Check memory is initialized in output field functions

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

diff --git a/CoreFlows/Models/src/FiveEqsTwoFluid.cxx b/CoreFlows/Models/src/FiveEqsTwoFluid.cxx
index 6bc3a36e23063090a4d72033677a968c9a2c4d54..8a165fa820fe98fdcd90348e4f37b98f46b1759c 100755 (executable)
--- a/CoreFlows/Models/src/FiveEqsTwoFluid.cxx
+++ b/CoreFlows/Models/src/FiveEqsTwoFluid.cxx
@@ -95,6 +95,7 @@ void FiveEqsTwoFluid::convectionState( const long &i, const long &j, const bool
                VecGetValues(_Uext, _nVar, _idm, _Uj);
        else
                VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
+
        if(_verbose && _nbTimeStep%_freqSave ==0)
        {
                cout<<"Convection Left state cell " << i<< ": "<<endl;
@@ -832,6 +833,7 @@ void FiveEqsTwoFluid::convectionMatrices()
                _Aroe[(_Ndim+1)*_nVar+1+idim]=0;
                _Aroe[(_Ndim+1)*_nVar+1+idim+_Ndim+1]=_vec_normal[idim];
        }
+
        // Take into account the convection term in the momentum eqts
        for(int idim=0; idim<_Ndim;idim++)
                for (int jdim=0; jdim<_Ndim;jdim++){
@@ -844,6 +846,7 @@ void FiveEqsTwoFluid::convectionMatrices()
                        _Aroe[(2+_Ndim+idim)*_nVar+_Ndim+1+jdim+1] += m_u2[idim]*_vec_normal[jdim];
                        _Aroe[(2+_Ndim+idim)*_nVar+_Ndim+1+idim+1] += m_u2[jdim]*_vec_normal[jdim];
                }
+
        // update \Delta alpha
        for (int idim=0; idim<_Ndim; idim++){
                _Aroe[ (1+idim)*_nVar] += dpi1*varrho_2*(1-m_alp1)*inv_a2_2*_vec_normal[idim];
@@ -866,6 +869,7 @@ void FiveEqsTwoFluid::convectionMatrices()
                        _Aroe[(_Ndim+1)*_nVar+ (1+idim)*_nVar+i] += (1-alpha)*varrho_2*(-m_alp1*m_rho2*b1*Delta_e1[i] -(1-m_alp1)*m_rho1*b2*Delta_e2[i])*_vec_normal[idim];
                }
        }
+
        // last row (total energy)
        for (int i=0; i<_nVar; i++){
                _Aroe[(2*_Ndim+2)*_nVar +i] += A5[i];
@@ -897,7 +901,6 @@ void FiveEqsTwoFluid::convectionMatrices()
                }
        }
        /******** Compute the eigenvalues and eigenvectors of Roe Matrix (using lapack)*********/
-       Polynoms Poly;
        dgeev_(jobvl, jobvr,  &_nVar,
                        Aroe,&LDA,egvaReal,egvaImag, egVectorL,
                        &LDVL,egVectorR,
@@ -908,7 +911,7 @@ void FiveEqsTwoFluid::convectionMatrices()
                *_runLogFile<<"FiveEqsTwoFluid::convectionMatrices: error dgeev_ : argument "<<-info<<" invalid"<<endl;
                throw CdmathException("FiveEqsTwoFluid::convectionMatrices: dgeev_ unsuccessful computation of the eigenvalues ");
        }
-       else if(info <0)
+       else if(info >0)
        {
                cout<<"Warning FiveEqsTwoFluid::convectionMatrices: dgeev_ did not compute all the eigenvalues, trying Rusanov scheme "<<endl;
                cout<<"Converged eigenvalues are ";
@@ -948,7 +951,7 @@ void FiveEqsTwoFluid::convectionMatrices()
                        valeurs_propres[j] = complex<double>(egvaReal[j],egvaImag[j]);
                }
 
-               taille_vp =Poly.new_tri_selectif(valeurs_propres,valeurs_propres.size(),_precision);
+               taille_vp =Polynoms::new_tri_selectif(valeurs_propres,valeurs_propres.size(),_precision);
 
                valeurs_propres_dist=vector< complex< double > >(taille_vp);
                for( int i=0 ; i<taille_vp ; i++)
@@ -1128,7 +1131,7 @@ void FiveEqsTwoFluid::convectionMatrices()
                }
                vector< complex< double > > y (taille_vp,0);
                for( int i=0 ; i<taille_vp ; i++)
-                       y[i] = Poly.abs_generalise(valeurs_propres_dist[i]);
+                       y[i] = Polynoms::abs_generalise(valeurs_propres_dist[i]);
 
                if(_entropicCorrection)
                {
@@ -1152,7 +1155,7 @@ void FiveEqsTwoFluid::convectionMatrices()
                        for( int i=0 ; i<taille_vp ; i++)
                                cout<<y[i] <<", "<<endl;
                }
-               Poly.abs_par_interp_directe(taille_vp,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
+               Polynoms::abs_par_interp_directe(taille_vp,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
 
                if( _spaceScheme ==pressureCorrection){
                        for( int i=0 ; i<_Ndim ; i++)
@@ -1175,7 +1178,7 @@ void FiveEqsTwoFluid::convectionMatrices()
 
        if(_entropicCorrection || _spaceScheme ==pressureCorrection || _spaceScheme ==lowMach){
                InvMatriceRoe( valeurs_propres_dist);
-               Poly.matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
+               Polynoms::matrixProduct(_absAroe, _nVar, _nVar, _invAroe, _nVar, _nVar, _signAroe);
        }
        else if (_spaceScheme==upwind)//upwind sans entropic
                SigneMatriceRoe( valeurs_propres_dist);
@@ -1353,6 +1356,7 @@ void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *
                _idm[k] = _idm[k-1] + 1;
 
        VecGetValues(_conservativeVars, _nVar, _idm, _externalStates);//On initialise l'état fantôme avec l'état interne
+       VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
 
        for(k=0; k<_Ndim; k++){
                q1_n+=_externalStates[(k+1)]*normale[k];
@@ -1373,7 +1377,6 @@ void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *
        if (_limitField[nameOfGroup].bcType==Wall){
 
                //Pour la convection, inversion du sens des vitesses
-               VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
                for(k=0; k<_Ndim; k++){
                        _externalStates[(k+1)]-= 2*q1_n*normale[k];
                        _externalStates[(k+1+1+_Ndim)]-= 2*q2_n*normale[k];
@@ -1419,7 +1422,6 @@ void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *
                for(k=1; k<_nVar; k++)
                        _idm[k] = _idm[k-1] + 1;
 
-               VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
                _idm[0] = 0;
                for(k=1; k<_nVar; k++)
                        _idm[k] = _idm[k-1] + 1;
@@ -1441,7 +1443,6 @@ void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *
                for(k=1; k<_nVar; k++)
                        _idm[k] = _idm[k-1] + 1;
 
-               VecGetValues(_primitiveVars, _nVar, _idm, _Vj);
                double alpha=_limitField[nameOfGroup].alpha;//void fraction outside
                double pression=_Vj[1];//pressure inside
                double T=_limitField[nameOfGroup].T;//temperature outside
@@ -1471,7 +1472,6 @@ void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *
                }
                _externalStates[_nVar-1] = _externalStates[0]      *(_fluides[0]->getInternalEnergy(T,rho_v) + v1_2/2)
                                                                                                                                                                                                                                                                                                                                        +_externalStates[1+_Ndim]*(_fluides[1]->getInternalEnergy(T,rho_l) + v2_2/2);
-
                // _Vj external primitives
                _Vj[0] = alpha;
                _Vj[_nVar-1] = T;
@@ -1506,7 +1506,6 @@ void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *
                hydroPress*=_externalStates[0]+_externalStates[_Ndim];//multiplication by rho the total density
 
                //Building the external state
-               VecGetValues(_primitiveVars, _nVar, _idm,_Vj);
                double alpha=_limitField[nameOfGroup].alpha;
                double pression=_limitField[nameOfGroup].p+hydroPress;
                double T=_limitField[nameOfGroup].T;
@@ -1523,7 +1522,6 @@ void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *
                }
                _externalStates[_nVar-1]=    alpha *rho_v*(_fluides[0]->getInternalEnergy(T,rho_v)+v1_2/2)
                                                                                                                                                                                                                                                                                                                                                +(1-alpha)*rho_l*(_fluides[1]->getInternalEnergy(T,rho_l)+v2_2/2);
-
                // _Vj external primitives
                _Vj[0] = alpha;
                _Vj[1] = pression;
@@ -1558,7 +1556,6 @@ void FiveEqsTwoFluid::setBoundaryState(string nameOfGroup, const int &j,double *
                hydroPress*=_externalStates[0]+_externalStates[_Ndim];//multiplication by rho the total density
 
                //Building the external state
-               VecGetValues(_primitiveVars, _nVar, _idm,_Vj);
                double pression_int=_Vj[1];
                double pression_ext=_limitField[nameOfGroup].p+hydroPress;
                double T=_Vj[_nVar-1];
@@ -1923,10 +1920,9 @@ void FiveEqsTwoFluid::entropicShift(double* n, double& vpcorr0, double& vpcorr1)
                        eigValuesLeft[j] = complex<double>(egvaReal[j],egvaImag[j]);
                        eigValuesRight[j] = complex<double>(egvaReal[j],egvaImag[j]);
                }
-               Polynoms Poly;
-               int sizeLeft =  Poly.new_tri_selectif(eigValuesLeft, eigValuesLeft.size(), _precision);
-               int sizeRight =  Poly.new_tri_selectif(eigValuesRight, eigValuesRight.size(), _precision);
-               if (_verbose && _nbTimeStep%_freqSave ==0)
+               int sizeLeft =  Polynoms::new_tri_selectif(eigValuesLeft, eigValuesLeft.size(), _precision);
+               int sizeRight =  Polynoms::new_tri_selectif(eigValuesRight, eigValuesRight.size(), _precision);
+               if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
                {
                        cout<<" Eigenvalue of JacoMat Left: " << endl;
                        for(int i=0; i<sizeLeft; i++)
@@ -1998,10 +1994,9 @@ void FiveEqsTwoFluid::entropicShift(double* n)
        for(int j=0; j<_nVar; j++){
                eigValuesRight[j] = complex<double>(egvaReal[j],egvaImag[j]);
        }
-       Polynoms Poly;
-       int sizeLeft =  Poly.new_tri_selectif(eigValuesLeft, eigValuesLeft.size(), _precision);
-       int sizeRight =  Poly.new_tri_selectif(eigValuesRight, eigValuesRight.size(), _precision);
-       if (_verbose && _nbTimeStep%_freqSave ==0)
+       int sizeLeft =  Polynoms::new_tri_selectif(eigValuesLeft, eigValuesLeft.size(), _precision);
+       int sizeRight =  Polynoms::new_tri_selectif(eigValuesRight, eigValuesRight.size(), _precision);
+       if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
        {
                cout<<" Eigenvalue of JacoMat Left: " << endl;
                for(int i=0; i<sizeLeft; i++)
@@ -2262,7 +2257,7 @@ void FiveEqsTwoFluid::save(){
        }
        _VV.setTime(_time,_nbTimeStep+1);
 
-       if (_nbTimeStep ==0){
+       if (_nbTimeStep ==0 || _restartWithNewFileName){
                string prim_suppress ="rm -rf "+prim+"_*";
                string cons_suppress ="rm -rf "+cons+"_*";
                system(prim_suppress.c_str());//Nettoyage des précédents calculs identiques
@@ -2370,7 +2365,7 @@ void FiveEqsTwoFluid::save(){
                }
                _Vitesse1.setTime(_time,_nbTimeStep);
                _Vitesse2.setTime(_time,_nbTimeStep);
-               if (_nbTimeStep ==0){
+               if (_nbTimeStep ==0 || _restartWithNewFileName){                
                        _Vitesse1.setInfoOnComponent(0,"Velocity_x_(m/s)");
                        _Vitesse1.setInfoOnComponent(1,"Velocity_y_(m/s)");
                        _Vitesse1.setInfoOnComponent(2,"Velocity_z_(m/s)");
@@ -2413,4 +2408,7 @@ void FiveEqsTwoFluid::save(){
                        }
                }
        }
+
+       if (_restartWithNewFileName)
+               _restartWithNewFileName=false;
 }