VecView(_primitiveVars, PETSC_VIEWER_STDOUT_SELF);
}
- if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
+ if(_nbPhases==2 && (_nbTimeStep-1)%_freqSave ==0){
if(_minm1<-_precision || _minm2<-_precision)
{
cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
{
VecCopy(_b,_newtonVariation);
VecScale(_newtonVariation, _dt);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
VecGetValues(_Uext, _nVar, _idm, _Uj);
else
VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"DriftModel::convectionState Left state cell " << i<< ": "<<endl;
for(int k =0; k<_nVar; k++)
rj = sqrt(_Uj[0]);
_Uroe[0] = ri*rj/sqrt(_porosityi*_porosityj); //moyenne geometrique des densites
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout << "Densité moyenne Roe gauche " << i << ": " << ri*ri << ", droite " << j << ": " << rj*rj << "->" << _Uroe[0] << endl;
xi = _Ui[1]/_Ui[0];//cm gauche
xj = _Uj[1]/_Uj[0];//cm droite
_Uroe[1] = (xi*ri + xj*rj)/(ri + rj);//moyenne de Roe des concentrations
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout << "Concentration de Roe gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[1] << endl;
for(int k=0;k<_Ndim;k++){
xi = _Ui[k+2];//phi rho u gauche
xj = _Uj[k+2];//phi rho u droite
_Uroe[2+k] = (xi/ri + xj/rj)/(ri + rj);
//"moyenne" des vitesses
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout << "Vitesse de Roe composante "<< k<<" gauche " << i << ": " << xi/(ri*ri) << ", droite " << j << ": " << xj/(rj*rj) << "->" << _Uroe[k+2] << endl;
}
// H = (rho E + p)/rho
xi = (xi + pi)/_Ui[0];
xj = (xj + pj)/_Uj[0];
_Uroe[_nVar-1] = (ri*xi + rj*xj)/(ri + rj);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout << "Enthalpie totale de Roe H gauche " << i << ": " << xi << ", droite " << j << ": " << xj << "->" << _Uroe[_nVar-1] << endl;
// Moyenne de Roe de Tg et Td
Ii = i*_nVar+_nVar-1;
Ii = j*_nVar+_nVar-1;
VecGetValues(_primitiveVars, 1, &Ii, &xj);
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"Convection interfacial state"<<endl;
for(int k=0;k<_nVar;k++)
else
VecGetValues(_conservativeVars, _nVar, _idm, _Uj);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "DriftModel::diffusionStateAndMatrices cellule gauche" << i << endl;
cout << "Ui = ";
for(int k=0; k<_nVar; k++)
_Udiff[k] = (_Ui[k]/_porosityi+_Uj[k]/_porosityj)/2;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "DriftModel::diffusionStateAndMatrices conservative diffusion state" << endl;
cout << "_Udiff = ";
}
_Diffusion[_nVar*_nVar-1]=-lambda/(m_v*C_v+m_l*C_l);
/*Affichages */
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Matrice de diffusion D, pour le couple (" << i << "," << j<< "):" << endl;
displayMatrix( _Diffusion,_nVar," Matrice de diffusion ");
double porosityj=_porosityField(j);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "setBoundaryState for group "<< nameOfGroup<< ", inner cell j= "<<j << " face unit normal vector "<<endl;
for(k=0; k<_Ndim; k++){
}
_externalStates[_nVar-1] = _externalStates[1]*_fluides[0]->getInternalEnergy(T,rho_v)+(_externalStates[0]-_externalStates[1])*_fluides[1]->getInternalEnergy(T,rho_l) + _externalStates[0]*v2/2;
}
- else if(_nbTimeStep%_freqSave ==0)
+ else if((_nbTimeStep-1)%_freqSave ==0)
cout<< "Warning : fluid going out through inlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition"<<endl;
_idm[0] = 0;
Tm=_limitField[nameOfGroup].T;
}
else{
- if(_nbTimeStep%_freqSave ==0)
+ if((_nbTimeStep-1)%_freqSave ==0)
cout<< "Warning : fluid going out through inletPressure boundary "<<nameOfGroup<<". Applying Neumann boundary condition for concentration, velocity and temperature"<<endl;
concentration=_Vj[0];
Tm=_Vj[_nVar-1];
VecAssemblyEnd(_Uextdiff);
}
else if (_limitField[nameOfGroup].bcType==Outlet){
- if(q_n<=0 && _nbTimeStep%_freqSave ==0)
+ if(q_n<=0 && (_nbTimeStep-1)%_freqSave ==0)
cout<< "Warning : fluid going in through outlet boundary "<<nameOfGroup<<". Applying Neumann boundary condition for concentration, velocity and temperature"<<endl;
//Computation of the hydrostatic contribution : scalar product between gravity vector and position vector
//entree: URoe = rho, cm, u, H
//sortie: matrices Roe+ et Roe- +Roe si scheme centre
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"DriftModel::convectionMatrices()"<<endl;
double umn=0, u_2=0; //valeur de u.normale et |u|²
double pression= getMixturePressure(cm, rhom, Tm);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"Roe state rhom="<<rhom<<" cm= "<<cm<<" Hm= "<<Hm<<" Tm= "<<Tm<<" pression "<<pression<<endl;
getMixturePressureDerivatives( m_v, m_l, pression, Tm);//EOS is involved to express pressure jump and sound speed
throw CdmathException("Calcul matrice de Roe: vitesse du son complexe");
}
double am=sqrt(_kappa*hm+_khi+cm*_ksi);//vitesse du son du melange
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"DriftModel::convectionMatrices, sound speed am= "<<am<<endl;
//On remplit la matrice de Roe
_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])+1*_entropicShift[i];
- Poly.abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
+ y[i] = Polynoms::abs_generalise(vp_dist[i])+1*_entropicShift[i];
+ Polynoms::abs_par_interp_directe(3,vp_dist, _Aroe, _nVar,_precision, _absAroe,y);
if( _spaceScheme ==pressureCorrection)
{
for(int idim=0;idim<_Ndim; idim++)
_Vij[2+idim]=_Uroe[2+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++)
_AroePlusImplicit[i] = _AroePlus[i];
}
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
displayMatrix(_Aroe, _nVar,"Matrice de Roe");
displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
}
}
- if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0 && _timeScheme==Implicit)
{
displayMatrix(_AroeMinusImplicit, _nVar,"Matrice _AroeMinusImplicit");
displayMatrix(_AroePlusImplicit, _nVar,"Matrice _AroePlusImplicit");
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 ==lowMach) || (_spaceScheme ==pressureCorrection))//upwind sans entropic
SigneMatriceRoe( vp_dist);
throw CdmathException("DriftModel::convectionMatrices: well balanced option not treated");
}
- if(_verbose && _nbTimeStep%_freqSave ==0 && _timeScheme==Implicit)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0 && _timeScheme==Implicit)
displayMatrix(_signAroe, _nVar,"Signe de la matrice de Roe");
}
_idm[k] = _nVar*i + k;
VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution diffusion: variables primitives maille " << i<<endl;
for(int q=0; q<_nVar; q++)
VecGetValues(_Uextdiff, _nVar, _idn, _Uj);
consToPrim(_Uj,_Vj,1);
}
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution diffusion: variables primitives maille " <<j <<endl;
for(int q=0; q<_nVar; q++)
_idm[0] = i;
VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution diffusion au 2nd membre pour la maille " << i << ": "<<endl;
for(int q=0; q<_nVar; q++)
VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution diffusion au 2nd membre pour la maille " << j << ": "<<endl;
for(int q=0; q<_nVar; q++)
}
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"DriftModel::sourceVector"<<endl;
cout<<"Ui="<<endl;
pressureLoss[_nVar-1]=-K*phirho*norm_u*norm_u*norm_u;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"DriftModel::pressureLossVector K= "<<K<<endl;
cout<<"pressure loss vector phirho= "<< phirho<<" norm_u= "<<norm_u<<endl;
void DriftModel::jacobian(const int &j, string nameOfGroup,double * normale)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"Jacobienne condition limite convection bord "<< nameOfGroup<<endl;
int k;
//calcule la jacobienne pour une CL de diffusion
void DriftModel::jacobianDiff(const int &j, string nameOfGroup)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"Jacobienne condition limite diffusion bord "<< nameOfGroup<<endl;
for(int k=0; k<_Ndim; k++)
W[j*_nVar+_nVar-1] += W[j*_nVar]*P[i*_nVar+(k+2)]*P[i*_nVar+(k+2)]*0.5;//phi rhom e_m + phi rho u*u
- if(_verbose && _nbTimeStep%_freqSave ==0){
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
cout<<"DriftModel::primToCons: T="<<temperature<<", pression= "<<pression<<endl;
cout<<"rhov= "<<rho_v<<", rhol= "<<rho_l<<", rhom= "<<W[j*(_nVar)]<<" e_v="<<e_v<<" e_l="<<e_l<<endl;
}
throw CdmathException("DriftModel::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<" DriftModel::primToConsJacobianMatrix" << endl;
displayVector(_Vi,_nVar," _Vi " );
for(int k=0;k<_Ndim;k++)
Wprim[k+2] = Wcons[k+2]/Wcons[0];
Wprim[_nVar-1] = Temperature;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"ConsToPrim Vecteur conservatif"<<endl;
for(int k=0;k<_nVar;k++)
double delta= b*b-4*a*c;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"DriftModel::getMixturePressure: a= "<<a<<", b= "<<b<<", c= "<<c<<", delta= "<<delta<<endl;
if(delta<0){
}
- if(_verbose && _nbTimeStep%_freqSave ==0){
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
cout<<"DriftModel::getMixturePressureAndTemperature: a= "<<a<<", b= "<<b<<", c= "<<c<<", delta= "<<delta<<endl;
cout<<"pressure= "<<pression<<", temperature= "<<temperature<<endl;
}
_kappa=temp2/denom;
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"rho_l= "<<rho_l<<", temp1= "<<temp1<<", temp2= "<<temp2<<", denom= "<<denom<<endl;
cout<<"khi= "<<_khi<<", ksi= "<<_ksi<<", kappa= "<<_kappa<<endl;
double rhom=_Ui[0];
double cm=_Vi[0];
double Hm=(_Ui[_nVar-1]+_Vi[1])/rhom;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"Entropic shift left state rhom="<<rhom<<" cm= "<<cm<<"Hm= "<<Hm<<endl;
double Tm=_Vi[_nVar-1];
double hm=Hm-0.5*ul_2;
rhom=_Uj[0];
cm=_Vj[0];
Hm=(_Uj[_nVar-1]+_Vj[1])/rhom;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"Entropic shift right state rhom="<<rhom<<" cm= "<<cm<<"Hm= "<<Hm<<endl;
Tm=_Vj[_nVar-1];
hm=Hm-0.5*ul_2;
_entropicShift[1]=abs(umnl - umnr);
_entropicShift[2]=abs(umnl+aml - (umnr+amr));
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"Entropic shift left eigenvalues: "<<endl;
cout<<"("<< umnl-aml<< ", " <<umnl<<", "<<umnl+aml << ")";
}
Vector DriftModel::convectionFlux(Vector U,Vector V, Vector normale, double porosity){
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"DriftModel::convectionFlux start"<<endl;
cout<<"Ucons"<<endl;
F(2+i)=m_v*vitesse_vn*vitesse_v(i)+m_l*vitesse_ln*vitesse_l(i)+pression*normale(i)*porosity;
F(2+_Ndim)=m_v*(e_v+0.5*vitesse_v*vitesse_v+pression/rho_v)*vitesse_vn+m_l*(e_l+0.5*vitesse_l*vitesse_l+pression/rho_l)*vitesse_ln;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"DriftModel::convectionFlux end"<<endl;
cout<<"Flux F(U,V)"<<endl;
Vector DriftModel::staggeredVFFCFlux()
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"DriftModel::staggeredVFFCFlux()start"<<endl;
if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
Fij=(1+theta)/2*F1+(1-theta)/2*F2;
*/
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"DriftModel::staggeredVFFCFlux() end uijn="<<uijn<<endl;
cout<<Fij<<endl;
void DriftModel::staggeredVFFCMatricesConservativeVariables(double u_mn)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"DriftModel::staggeredVFFCMatricesConservativeVariables()"<<endl;
if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
if(u_mn>_precision)
{
Hm=Emi+pj/rhomi;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Hm= "<<Hm<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
/***********Calcul des valeurs propres ********/
}
double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
//On remplit les valeurs propres
else if(u_mn<-_precision)
{
Hm=Emj+pi/rhomj;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Hm= "<<Hm<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
/***********Calcul des valeurs propres ********/
throw CdmathException("staggeredVFFCMatricesConservativeVariables: vitesse du son complexe");
}
double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
//On remplit les valeurs propres
void DriftModel::staggeredVFFCMatricesPrimitiveVariables(double u_mn)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"DriftModel::staggeredVFFCMatricesPrimitiveVariables()"<<endl;
if(_spaceScheme!=staggered || _nonLinearFormulation!=VFFC)
if(u_mn>_precision)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Emi= "<<Emi<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
/***********Calcul des valeurs propres ********/
}
double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
//On remplit les valeurs propres
}
else if(u_mn<-_precision)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Emj= "<<Emj<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
/***********Calcul des valeurs propres ********/
throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
}
double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
//On remplit les valeurs propres
if(u_mn>_precision)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"VFFC Staggered state rhomi="<<rhomi<<" cmi= "<<cmi<<" Hmi= "<<Hmi<<" Tmi= "<<Tmi<<" pj= "<<pj<<endl;
/***********Calcul des valeurs propres ********/
}
double amj=sqrt(_kappa*hmj+_khi+cmj*_ksi);//vitesse du son du melange
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", amj= "<<amj<<endl;
//On remplit les valeurs propres
}
else if(u_mn<-_precision)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"VFFC Staggered state rhomj="<<rhomj<<" cmj= "<<cmj<<" Hmj= "<<Hmj<<" Tmj= "<<Tmj<<" pi= "<<pi<<endl;
/***********Calcul des valeurs propres ********/
throw CdmathException("staggeredVFFCMatricesPrimitiveVariables: vitesse du son complexe");
}
double ami=sqrt(_kappa*hmi+_khi+cmi*_ksi);//vitesse du son du melange
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"_khi= "<<_khi<<", _kappa= "<< _kappa << ", _ksi= "<<_ksi <<", ami= "<<ami<<endl;
//On remplit les valeurs propres
}
else//case nil velocity on the interface, multiply by jacobian matrix
{
- 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);
}
}
}
else
throw CdmathException("DriftModel::primToConsJacobianMatrix: eos should be StiffenedGas or StiffenedGasDellacherie");
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"_drho_sur_dcv= "<<_drho_sur_dcv<<", _drho_sur_dp= "<<_drho_sur_dp<<", _drho_sur_dT= "<<_drho_sur_dT<<endl;
cout<<"_drhocv_sur_dcv= "<<_drhocv_sur_dcv<<", _drhocv_sur_dp= "<<_drhocv_sur_dp<<", _drhocv_sur_dT= "<<_drhocv_sur_dT<<endl;
#include "ProblemFluid.hxx"
-#include "Fluide.h"
#include "math.h"
using namespace std;
VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Vext);
VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uextdiff);
// VecCreateSeq(PETSC_COMM_SELF, _nVar*_Nmailles, &_conservativeVars);
- VecCreate(PETSC_COMM_SELF, &_conservativeVars);
+ VecCreate(PETSC_COMM_SELF, &_conservativeVars);//Current conservative variables at Newton iteration k between time steps n and n+1
VecSetSizes(_conservativeVars,PETSC_DECIDE,_nVar*_Nmailles);
VecSetBlockSize(_conservativeVars,_nVar);
VecSetFromOptions(_conservativeVars);
- VecDuplicate(_conservativeVars, &_old);
- VecDuplicate(_conservativeVars, &_newtonVariation);
- VecDuplicate(_conservativeVars, &_b);
- VecDuplicate(_conservativeVars, &_primitiveVars);
+ VecDuplicate(_conservativeVars, &_old);//Old conservative variables at time step n
+ VecDuplicate(_conservativeVars, &_newtonVariation);//Newton variation Uk+1-Uk to be computed between time steps n and n+1
+ VecDuplicate(_conservativeVars, &_b);//Right hand side of Newton method at iteration k between time steps n and n+1
+ VecDuplicate(_conservativeVars, &_primitiveVars);//Current primitive variables at Newton iteration k between time steps n and n+1
if(_isScaling)
{
KSPGetPC(_ksp, &_pc);
PCSetType(_pc, _pctype);
+ // Creation du solveur de Newton de PETSc
+ if( _timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
+ {
+ SNESType snestype;
+
+ // set nonlinear solver
+ if (_nonLinearSolver == Newton_PETSC_LINESEARCH || _nonLinearSolver == Newton_PETSC_LINESEARCH_BASIC || _nonLinearSolver == Newton_PETSC_LINESEARCH_BT || _nonLinearSolver == Newton_PETSC_LINESEARCH_SECANT || _nonLinearSolver == Newton_PETSC_LINESEARCH_NLEQERR)
+ snestype = (char*)&SNESNEWTONLS;
+ else if (_nonLinearSolver == Newton_PETSC_TRUSTREGION)
+ snestype = (char*)&SNESNEWTONTR;
+ else if (_nonLinearSolver == Newton_PETSC_NGMRES)
+ snestype = (char*)&SNESNGMRES;
+ else if (_nonLinearSolver ==Newton_PETSC_ASPIN)
+ snestype = (char*)&SNESASPIN;
+ else if(_nonLinearSolver != Newton_SOLVERLAB)
+ {
+ cout << "!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" << endl;
+ *_runLogFile << "!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" << endl;
+ _runLogFile->close();
+ throw CdmathException("!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" );
+ }
+
+ SNESCreate(PETSC_COMM_WORLD, &_snes);
+ SNESSetType( _snes, snestype);
+ SNESGetLineSearch( _snes, &_linesearch);
+ if(_nonLinearSolver == Newton_PETSC_LINESEARCH_BASIC)
+ SNESLineSearchSetType( _linesearch, SNESLINESEARCHBASIC );
+ else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_BT)
+ SNESLineSearchSetType( _linesearch, SNESLINESEARCHBT );
+ else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_SECANT)
+ SNESLineSearchSetType( _linesearch, SNESLINESEARCHL2 );
+ else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_NLEQERR)
+ SNESLineSearchSetType( _linesearch, SNESLINESEARCHNLEQERR );
+
+ PetscViewerCreate(PETSC_COMM_WORLD,&_monitorLineSearch);
+ PetscViewerSetType(_monitorLineSearch, PETSCVIEWERASCII);
+
+ SNESSetTolerances(_snes,_precision_Newton,_precision_Newton,_precision_Newton,_maxNewtonIts,-1);
+
+ SNESSetFunction(_snes,_newtonVariation,computeSnesRHS,this);
+ SNESSetJacobian(_snes,_A,_A,computeSnesJacobian,this);
+ }
+
_initializedMemory=true;
save();//save initial data
}
return _dt>0;//No need to call MatShift as the linear system matrix is filled at each Newton iteration (unlike linear problem)
}
+bool ProblemFluid::solveTimeStep(){
+ if(_timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
+ return solveNewtonPETSc();
+ else
+ return ProblemCoreFlows::solveTimeStep();
+}
+
+bool ProblemFluid::solveNewtonPETSc()
+{
+ if( (_nbTimeStep-1)%_freqSave ==0)
+ SNESLineSearchSetDefaultMonitor(_linesearch, _monitorLineSearch);
+ else
+ SNESLineSearchSetDefaultMonitor(_linesearch, NULL);
+
+ SNESSolve(_snes,NULL,_conservativeVars);
+
+ SNESConvergedReason reason;
+ SNESGetConvergedReason(_snes,&reason);
+
+ if( (_nbTimeStep-1)%_freqSave ==0)
+ {
+ if(reason == SNES_CONVERGED_FNORM_ABS )
+ cout<<"Converged with absolute norm (absolute tolerance) less than "<<_precision_Newton<<", (||F|| < atol)"<<endl;
+ else if(reason == SNES_CONVERGED_FNORM_RELATIVE )
+ cout<<"Converged because residual has decreased by a factor less than "<<_precision_Newton<<", (||F|| < rtol*||F_initial||)"<<endl;
+ else if(reason == SNES_CONVERGED_SNORM_RELATIVE )
+ cout<<"Converged with relative norm (relative tolerance) less than "<<_precision_Newton<<", (|| delta x || < stol || x ||)"<<endl;
+ else if(reason == SNES_CONVERGED_ITS )
+ cout<<"SNESConvergedSkip() was chosen as the convergence test; thus the usual convergence criteria have not been checked and may or may not be satisfied"<<endl;
+ else if(reason == SNES_DIVERGED_LINEAR_SOLVE )
+ cout<<"Solving linear system failed"<<endl;
+ else if(reason == SNES_DIVERGED_LINE_SEARCH )
+ cout<<"Line search failed"<<endl;
+ else if(reason == SNES_DIVERGED_MAX_IT )
+ cout<<"Reached the maximum number of iterations"<<endl;
+
+ SNESGetIterationNumber(_snes, &_NEWTON_its);
+ PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n\n", _NEWTON_its);
+ *_runLogFile <<endl;
+ *_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its <<endl;
+ }
+
+ return reason>0;
+}
+
bool ProblemFluid::iterateTimeStep(bool &converged)
{
bool stop=false;
if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
_maxvp=0.;
- computeTimeStep(stop);//This compute timestep is just to update the linear system. The time step was imposed befor starting the Newton iterations
+ computeTimeStep(stop);//This compute timestep is just to update the linear system. The time step was imposed before starting the Newton iterations
}
if(stop){//Le compute time step ne s'est pas bien passé
cout<<"ComputeTimeStep failed"<<endl;
double relaxation=1;//Uk+1=Uk+relaxation*deltaU
- VecAXPY(_conservativeVars, relaxation, _newtonVariation);
+ VecAXPY(_conservativeVars, relaxation, _newtonVariation);
+
+ //mise a jour du champ primitif
+ updatePrimitives();
+
+ if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
+ {
+ cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
+ *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
+ converged=false;
+ return false;
+ }
+ if(_system)
+ {
+ cout<<"Vecteur Ukp1-Uk "<<endl;
+ VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
+ cout << "Nouvel etat courant Uk de l'iteration Newton: " << endl;
+ VecView(_conservativeVars, PETSC_VIEWER_STDOUT_SELF);
+ }
+
+ if(_nbPhases==2 && (_nbTimeStep-1)%_freqSave ==0){
+ if(_minm1<-_precision || _minm2<-_precision)
+ {
+ cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
+ *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
+ }
+
+ if (_nbVpCplx>0){
+ cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
+ *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
+ }
+ }
+
+ return true;
+}
+
+double ProblemFluid::computeTimeStep(bool & stop){//dt is not known and will not contribute to the Newton scheme
+
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "ProblemFluid::computeTimeStep : Début calcul matrice implicite et second membre"<<endl;
+ cout << endl;
+ }
+ if(_restartWithNewTimeScheme)//This is a change of time scheme during a simulation
+ {
+ if(_timeScheme == Implicit)
+ MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
+ else
+ MatDestroy(&_A);
+ _restartWithNewTimeScheme=false;
+ }
+ if(_timeScheme == Implicit)
+ MatZeroEntries(_A);
+
+ VecAssemblyBegin(_b);
+ VecZeroEntries(_b);
+
+ std::vector< int > idCells(2);
+ PetscInt idm, idn, size = 1;
+
+ long nbFaces = _mesh.getNumberOfFaces();
+ Face Fj;
+ Cell Ctemp1,Ctemp2;
+ string nameOfGroup;
+
+ for (int j=0; j<nbFaces;j++){
+ Fj = _mesh.getFace(j);
+ _isBoundary=Fj.isBorder();
+ idCells = Fj.getCellsId();
+
+ // If Fj is on the boundary
+ if (_isBoundary)
+ {
+ for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
+ {
+ // compute the normal vector corresponding to face j : from Ctemp1 outward
+ Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
+ if (_Ndim >1){
+ for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
+ {//we look for l the index of the face Fj for the cell Ctemp1
+ if (j == Ctemp1.getFacesId()[l])
+ {
+ for (int idim = 0; idim < _Ndim; ++idim)
+ _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
+ break;
+ }
+ }
+ }else{ // _Ndim = 1, build normal vector (bug cdmath)
+ if(!_sectionFieldSet)
+ {
+ if (Fj.x()<Ctemp1.x())
+ _vec_normal[0] = -1;
+ else
+ _vec_normal[0] = 1;
+ }
+ else
+ {
+ if(idCells[0]==0)
+ _vec_normal[0] = -1;
+ else//idCells[0]==31
+ _vec_normal[0] = 1;
+ }
+ }
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << _vec_normal[p] << ",";
+ cout << "). "<<endl;
+ }
+ nameOfGroup = Fj.getGroupName();
+ _porosityi=_porosityField(idCells[k]);
+ _porosityj=_porosityi;
+ setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
+ convectionState(idCells[k],0,true);
+ convectionMatrices();
+ diffusionStateAndMatrices(idCells[k], 0, true);
+ // compute 1/dxi
+ if (_Ndim > 1)
+ _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
+ else
+ _inv_dxi = 1/Ctemp1.getMeasure();
+
+ addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
+ addDiffusionToSecondMember(idCells[k],-1,true);
+ addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
+
+ if(_timeScheme == Implicit){
+ for(int l=0; l<_nVar*_nVar;l++){
+ _AroeMinusImplicit[l] *= _inv_dxi;
+ _Diffusion[l] *=_inv_dxi*_inv_dxi;
+ }
+
+ jacobian(idCells[k],nameOfGroup,_vec_normal);
+ jacobianDiff(idCells[k],nameOfGroup);
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
+ cout << "Matrice Jacobienne CL convection:" << endl;
+ for(int p=0; p<_nVar; p++){
+ for(int q=0; q<_nVar; q++)
+ cout << _Jcb[p*_nVar+q] << "\t";
+ cout << endl;
+ }
+ cout << endl;
+ cout << "Matrice Jacobienne CL diffusion:" << endl;
+ for(int p=0; p<_nVar; p++){
+ for(int q=0; q<_nVar; q++)
+ cout << _JcbDiff[p*_nVar+q] << "\t";
+ cout << endl;
+ }
+ cout << endl;
+ }
+ idm = idCells[k];
+ //calcul et insertion de A^-*Jcb
+ Polynoms::matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
+ MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
+
+ if(_verbose)
+ displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
+
+ //insertion de -A^-
+ for(int k=0; k<_nVar*_nVar;k++){
+ _AroeMinusImplicit[k] *= -1;
+ }
+ MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
+ if(_verbose)
+ displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
+
+ //calcul et insertion de D*JcbDiff
+ Polynoms::matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
+ MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
+ for(int k=0; k<_nVar*_nVar;k++)
+ _Diffusion[k] *= -1;
+ MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
+ }
+ }
+ } else if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
+ // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
+ Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
+ Ctemp2 = _mesh.getCell(idCells[1]);
+ if (_Ndim >1){
+ for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
+ if (j == Ctemp1.getFacesId()[l]){
+ for (int idim = 0; idim < _Ndim; ++idim)
+ _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
+ break;
+ }
+ }
+ }else{ // _Ndim = 1, build normal vector (bug cdmath)
+ if(!_sectionFieldSet)
+ {
+ if (Fj.x()<Ctemp1.x())
+ _vec_normal[0] = -1;
+ else
+ _vec_normal[0] = 1;
+ }
+ else
+ {
+ if(idCells[0]>idCells[1])
+ _vec_normal[0] = -1;
+ else
+ _vec_normal[0] = 1;
+ }
+ }
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
+ cout<<" Normal vector= ";
+ for (int idim = 0; idim < _Ndim; ++idim)
+ cout<<_vec_normal[idim]<<", ";
+ cout<<endl;
+ }
+ _porosityi=_porosityField(idCells[0]);
+ _porosityj=_porosityField(idCells[1]);
+ convectionState(idCells[0],idCells[1],false);
+ convectionMatrices();
+ diffusionStateAndMatrices(idCells[0], idCells[1], false);
+
+ // compute 1/dxi and 1/dxj
+ if (_Ndim > 1)
+ {
+ _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
+ _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
+ }
+ else
+ {
+ _inv_dxi = 1/Ctemp1.getMeasure();
+ _inv_dxj = 1/Ctemp2.getMeasure();
+ }
+
+ addConvectionToSecondMember(idCells[0],idCells[1], false);
+ addDiffusionToSecondMember( idCells[0],idCells[1], false);
+ addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
+
+ if(_timeScheme == Implicit){
+ for(int k=0; k<_nVar*_nVar;k++)
+ {
+ _AroeMinusImplicit[k] *= _inv_dxi;
+ _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);
+ }
+ idm = idCells[0];
+ idn = idCells[1];
+ //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNb<<endl;
+ MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
+
+ if(_verbose){
+ displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
+ displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
+ }
+ for(int k=0;k<_nVar*_nVar;k++){
+ _AroeMinusImplicit[k] *= -1;
+ _Diffusion[k] *= -1;
+ }
+ MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
+ if(_verbose){
+ displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
+ displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
+ }
+ for(int k=0; k<_nVar*_nVar;k++)
+ {
+ _AroePlusImplicit[k] *= _inv_dxj;
+ _Diffusion[k] *=_inv_dxj/_inv_dxi;
+ }
+ MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
+ if(_verbose)
+ displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
+
+ for(int k=0;k<_nVar*_nVar;k++){
+ _AroePlusImplicit[k] *= -1;
+ _Diffusion[k] *= -1;
+ }
+ MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
+
+ if(_verbose)
+ displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
+ }
+ }
+ else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
+ *_runLogFile<<"Warning: treatment of a junction node"<<endl;
+
+ if(!_sectionFieldSet)
+ {
+ _runLogFile->close();
+ throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
+ }
+ int largestSectionCellIndex=0;
+ for(int i=1;i<Fj.getNumberOfCells();i++){
+ if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
+ largestSectionCellIndex=i;
+ }
+ idm = idCells[largestSectionCellIndex];
+ Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
+ _porosityi=_porosityField(idm);
+
+ if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
+ _vec_normal[0] = 1;
+ else//j==16
+ _vec_normal[0] = -1;
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
+ cout << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << _vec_normal[p] << ",";
+ cout << "). "<<endl;
+ }
+ for(int i=0;i<Fj.getNumberOfCells();i++){
+ if(i != largestSectionCellIndex){
+ idn = idCells[i];
+ Ctemp2 = _mesh.getCell(idn);
+ _porosityj=_porosityField(idn);
+ convectionState(idm,idn,false);
+ convectionMatrices();
+ diffusionStateAndMatrices(idm, idn,false);
+
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
+
+ _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
+ _inv_dxj = 1/Ctemp2.getMeasure();
+
+ addConvectionToSecondMember(idm,idn, false);
+ _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
+ addDiffusionToSecondMember(idm,idn, false);
+ _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
+ addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
+
+ if(_timeScheme == Implicit){
+ for(int k=0; k<_nVar*_nVar;k++)
+ {
+ _AroeMinusImplicit[k] *= _inv_dxi;
+ _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
+ }
+ MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
+
+ if(_verbose){
+ displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
+ displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
+ }
+ for(int k=0;k<_nVar*_nVar;k++){
+ _AroeMinusImplicit[k] *= -1;
+ _Diffusion[k] *= -1;
+ }
+ MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
+ if(_verbose){
+ displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
+ displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
+ }
+ for(int k=0; k<_nVar*_nVar;k++)
+ {
+ _AroePlusImplicit[k] *= _inv_dxj;
+ _Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
+ }
+ MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
+ if(_verbose)
+ displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
+
+ for(int k=0;k<_nVar*_nVar;k++){
+ _AroePlusImplicit[k] *= -1;
+ _Diffusion[k] *= -1;
+ }
+ MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
+
+ if(_verbose)
+ displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
+ }
+ }
+ }
+ }
+ else
+ {
+ cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
+ _runLogFile->close();
+ throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
+ }
+
+ }
+ VecAssemblyEnd(_b);
+
+ if(_timeScheme == Implicit){
+ for(int imaille = 0; imaille<_Nmailles; imaille++)
+ MatSetValuesBlocked(_A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
+
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ displayMatrix(_GravityImplicitationMatrix,_nVar,"Gravity matrix:");
+
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
+ cout << "ProblemFluid::computeTimeStep : Fin calcul matrice implicite et second membre"<<endl;
+ cout << "ProblemFluid::computeTimeStep : Matrice implicite :"<<endl;
+ MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ cout << "ProblemFluid::computeTimeStep : Second membre :"<<endl;
+ VecView(_b, PETSC_VIEWER_STDOUT_WORLD);
+ cout << endl;
+ }
+ }
+
+ stop=false;
+
+ /*
+ if(_nbTimeStep+1<_cfl)
+ return (_nbTimeStep+1)*_minl/_maxvp;
+ else
+ */
+ return _cfl*_minl/_maxvp;
+}
+
+void ProblemFluid::computeNewtonVariation()
+{
+ if(_system)
+ {
+ cout<<"Vecteur courant Uk "<<endl;
+ VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+ if(_timeScheme == Explicit)
+ {
+ VecCopy(_b,_newtonVariation);
+ VecScale(_newtonVariation, _dt);
+ if(_system && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
+ VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+ }
+ else
+ {
+ if(_system)
+ {
+ cout << "Matrice du système linéaire avant contribution delta t" << endl;
+ MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ cout << "Second membre du système linéaire avant contribution delta t" << endl;
+ VecView(_b, PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+ MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
+
+ VecAXPY(_b, 1/_dt, _old);
+ VecAXPY(_b, -1/_dt, _conservativeVars);
+ MatShift(_A, 1/_dt);
+
+#if PETSC_VERSION_GREATER_3_5
+ KSPSetOperators(_ksp, _A, _A);
+#else
+ KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
+#endif
+
+ if(_system)
+ {
+ cout << "Matrice du système linéaire après contribution delta t" << endl;
+ MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ cout << "Second membre du système linéaire après contribution delta t" << endl;
+ VecView(_b, PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+
+ if(_conditionNumber)
+ KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
+ if(!_isScaling)
+ {
+ KSPSolve(_ksp, _b, _newtonVariation);
+ }
+ else
+ {
+ computeScaling(_maxvp);
+ int indice;
+ VecAssemblyBegin(_vecScaling);
+ VecAssemblyBegin(_invVecScaling);
+ for(int imaille = 0; imaille<_Nmailles; imaille++)
+ {
+ indice = imaille;
+ VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
+ VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
+ }
+ VecAssemblyEnd(_vecScaling);
+ VecAssemblyEnd(_invVecScaling);
+
+ if(_system)
+ {
+ cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
+ MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
+ VecView(_b, PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+ MatDiagonalScale(_A,_vecScaling,_invVecScaling);
+ if(_system)
+ {
+ cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
+ MatView(_A,PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+ VecCopy(_b,_bScaling);
+ VecPointwiseMult(_b,_vecScaling,_bScaling);
+ if(_system)
+ {
+ cout << "Produit du second membre par le preconditionneur bloc diagonal a gauche" << endl;
+ VecView(_b, PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+
+ KSPSolve(_ksp,_b, _bScaling);
+ VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
+ }
+ if(_system)
+ {
+ cout << "solution du systeme lineaire local:" << endl;
+ VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+ }
+}
+
+void ProblemFluid::computeNewtonRHS( Vec X, Vec F_X){//dt is known and will contribute to the right hand side of the Newton scheme
+
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "ProblemFluid::computeNewtonRHS : Début calcul second membre pour PETSc, _dt="<<_dt<<endl;
+ cout << endl;
+ cout<<"Vecteur courant Uk "<<endl;
+ VecView(X,PETSC_VIEWER_STDOUT_SELF);
+ cout << endl;
+ }
+
+ VecCopy(X,_conservativeVars);
+ updatePrimitives();
+
+ VecAssemblyBegin(_b);
+ VecZeroEntries(_b);
+
+ std::vector< int > idCells(2);
+ PetscInt idm, idn, size = 1;
+
+ long nbFaces = _mesh.getNumberOfFaces();
+ Face Fj;
+ Cell Ctemp1,Ctemp2;
+ string nameOfGroup;
+
+ for (int j=0; j<nbFaces;j++){
+ Fj = _mesh.getFace(j);
+ _isBoundary=Fj.isBorder();
+ idCells = Fj.getCellsId();
+
+ // If Fj is on the boundary
+ if (_isBoundary)
+ {
+ for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
+ {
+ // compute the normal vector corresponding to face j : from Ctemp1 outward
+ Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
+ if (_Ndim >1){
+ for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
+ {//we look for l the index of the face Fj for the cell Ctemp1
+ if (j == Ctemp1.getFacesId()[l])
+ {
+ for (int idim = 0; idim < _Ndim; ++idim)
+ _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
+ break;
+ }
+ }
+ }else{ // _Ndim = 1, build normal vector (bug cdmath)
+ if(!_sectionFieldSet)
+ {
+ if (Fj.x()<Ctemp1.x())
+ _vec_normal[0] = -1;
+ else
+ _vec_normal[0] = 1;
+ }
+ else
+ {
+ if(idCells[0]==0)
+ _vec_normal[0] = -1;
+ else//idCells[0]==31
+ _vec_normal[0] = 1;
+ }
+ }
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << _vec_normal[p] << ",";
+ cout << "). "<<endl;
+ }
+ nameOfGroup = Fj.getGroupName();
+ _porosityi=_porosityField(idCells[k]);
+ _porosityj=_porosityi;
+ setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
+ convectionState(idCells[k],0,true);
+ convectionMatrices();
+ diffusionStateAndMatrices(idCells[k], 0, true);
+ // compute 1/dxi
+ if (_Ndim > 1)
+ _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
+ else
+ _inv_dxi = 1/Ctemp1.getMeasure();
+
+ addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
+ addDiffusionToSecondMember(idCells[k],-1,true);
+ addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
+ }
+ } else if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
+ // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
+ Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
+ Ctemp2 = _mesh.getCell(idCells[1]);
+ if (_Ndim >1){
+ for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
+ if (j == Ctemp1.getFacesId()[l]){
+ for (int idim = 0; idim < _Ndim; ++idim)
+ _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
+ break;
+ }
+ }
+ }else{ // _Ndim = 1, build normal vector (bug cdmath)
+ if(!_sectionFieldSet)
+ {
+ if (Fj.x()<Ctemp1.x())
+ _vec_normal[0] = -1;
+ else
+ _vec_normal[0] = 1;
+ }
+ else
+ {
+ if(idCells[0]>idCells[1])
+ _vec_normal[0] = -1;
+ else
+ _vec_normal[0] = 1;
+ }
+ }
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
+ cout<<" Normal vector= ";
+ for (int idim = 0; idim < _Ndim; ++idim)
+ cout<<_vec_normal[idim]<<", ";
+ cout<<endl;
+ }
+ _porosityi=_porosityField(idCells[0]);
+ _porosityj=_porosityField(idCells[1]);
+ convectionState(idCells[0],idCells[1],false);
+ convectionMatrices();
+ diffusionStateAndMatrices(idCells[0], idCells[1], false);
+
+ // compute 1/dxi and 1/dxj
+ if (_Ndim > 1)
+ {
+ _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
+ _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
+ }
+ else
+ {
+ _inv_dxi = 1/Ctemp1.getMeasure();
+ _inv_dxj = 1/Ctemp2.getMeasure();
+ }
+
+ addConvectionToSecondMember(idCells[0],idCells[1], false);
+ addDiffusionToSecondMember( idCells[0],idCells[1], false);
+ addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
+
+ }
+ else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
+ *_runLogFile<<"Warning: treatment of a junction node"<<endl;
+
+ if(!_sectionFieldSet)
+ {
+ _runLogFile->close();
+ throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
+ }
+ int largestSectionCellIndex=0;
+ for(int i=1;i<Fj.getNumberOfCells();i++){
+ if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
+ largestSectionCellIndex=i;
+ }
+ idm = idCells[largestSectionCellIndex];
+ Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
+ _porosityi=_porosityField(idm);
+
+ if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
+ _vec_normal[0] = 1;
+ else//j==16
+ _vec_normal[0] = -1;
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
+ cout << " ; vecteur normal=(";
+ for(int p=0; p<_Ndim; p++)
+ cout << _vec_normal[p] << ",";
+ cout << "). "<<endl;
+ }
+ for(int i=0;i<Fj.getNumberOfCells();i++){
+ if(i != largestSectionCellIndex){
+ idn = idCells[i];
+ Ctemp2 = _mesh.getCell(idn);
+ _porosityj=_porosityField(idn);
+ convectionState(idm,idn,false);
+ convectionMatrices();
+ diffusionStateAndMatrices(idm, idn,false);
- //mise a jour du champ primitif
- updatePrimitives();
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
- if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
- {
- cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
- *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
- converged=false;
- return false;
- }
- if(_system)
- {
- cout<<"Vecteur Ukp1-Uk "<<endl;
- VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
- cout << "Nouvel etat courant Uk de l'iteration Newton: " << endl;
- VecView(_conservativeVars, PETSC_VIEWER_STDOUT_SELF);
- }
+ _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
+ _inv_dxj = 1/Ctemp2.getMeasure();
- if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
- if(_minm1<-_precision || _minm2<-_precision)
+ addConvectionToSecondMember(idm,idn, false);
+ _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
+ addDiffusionToSecondMember(idm,idn, false);
+ _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
+ addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
+ }
+ }
+ }
+ else
{
- cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
- *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
+ cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
+ _runLogFile->close();
+ throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
}
- if (_nbVpCplx>0){
- cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
- *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
- }
}
- _minm1=1e30;
- _minm2=1e30;
- _nbMaillesNeg=0;
- _nbVpCplx =0;
- _part_imag_max=0;
+
+ //Contribution from delta t
+ VecAXPY(_b, 1/_dt, _old);
+ VecAXPY(_b, -1/_dt, _conservativeVars);
- return true;
+ VecAssemblyEnd(_b);
+ VecCopy(_b,F_X);
+ VecScale(F_X,-1.);
+
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ {
+ cout << "ProblemFluid::computeNewtonRHS : Fin calcul second membre pour PETSc"<<endl;
+ VecView(F_X, PETSC_VIEWER_STDOUT_WORLD);
+ cout << endl;
+ }
}
-double ProblemFluid::computeTimeStep(bool & stop){
- VecZeroEntries(_b);
+int ProblemFluid::computeSnesRHS(SNES snes, Vec X, Vec F_X, void *ctx)//Prototype imposé par PETSc
+{
+ ProblemFluid * myProblem = (ProblemFluid *) ctx;
+ myProblem->computeNewtonRHS( X, F_X);
+
+ return 0;
+}
- if(_restartWithNewTimeScheme)//This is a change of time scheme during a simulation
+void ProblemFluid::computeNewtonJacobian( Vec X, Mat A){//dt is known and will contribute to the jacobian matrix of the Newton scheme
+
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
- if(_timeScheme == Implicit)
- MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
- else
- MatDestroy(&_A);
- _restartWithNewTimeScheme=false;
+ cout << "ProblemFluid::computeNewtonJacobian : Début calcul Jacobienne schéma Newton pour PETSc, _dt="<<_dt<<endl;
+ cout << endl;
}
- if(_timeScheme == Implicit)
- MatZeroEntries(_A);
- VecAssemblyBegin(_b);
- VecAssemblyBegin(_conservativeVars);
- std::vector< int > idCells;
+ if(_timeScheme == Explicit){
+ MatCreateConstantDiagonal(PETSC_COMM_SELF, _nVar, _nVar, _nVar*_Nmailles, _nVar*_Nmailles,1./_dt, &A);
+ return ;
+ }
+
+ MatZeroEntries(A);
+ VecCopy(X,_conservativeVars);
+
+ std::vector< int > idCells(2);
PetscInt idm, idn, size = 1;
long nbFaces = _mesh.getNumberOfFaces();
_vec_normal[0] = 1;
}
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
for(int p=0; p<_Ndim; p++)
else
_inv_dxi = 1/Ctemp1.getMeasure();
- addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
- addDiffusionToSecondMember(idCells[k],-1,true);
- addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
-
- if(_timeScheme == Implicit){
for(int l=0; l<_nVar*_nVar;l++){
_AroeMinusImplicit[l] *= _inv_dxi;
_Diffusion[l] *=_inv_dxi*_inv_dxi;
jacobian(idCells[k],nameOfGroup,_vec_normal);
jacobianDiff(idCells[k],nameOfGroup);
- if(_verbose && _nbTimeStep%_freqSave ==0){
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
cout << "Matrice Jacobienne CL convection:" << endl;
for(int p=0; p<_nVar; p++){
for(int q=0; q<_nVar; q++)
idm = idCells[k];
//calcul et insertion de A^-*Jcb
Polynoms::matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idm, _a, ADD_VALUES);
if(_verbose)
displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
for(int k=0; k<_nVar*_nVar;k++){
_AroeMinusImplicit[k] *= -1;
}
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
if(_verbose)
displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
//calcul et insertion de D*JcbDiff
Polynoms::matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idm, _a, ADD_VALUES);
for(int k=0; k<_nVar*_nVar;k++)
_Diffusion[k] *= -1;
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
}
- }
} else if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
// compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
_vec_normal[0] = 1;
}
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
cout<<" Normal vector= ";
_inv_dxj = 1/Ctemp2.getMeasure();
}
- addConvectionToSecondMember( idCells[0],idCells[1], false);
- addDiffusionToSecondMember( idCells[0],idCells[1], false);
- addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
-
- if(_timeScheme == Implicit){
for(int k=0; k<_nVar*_nVar;k++)
{
_AroeMinusImplicit[k] *= _inv_dxi;
idm = idCells[0];
idn = idCells[1];
//cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNb<<endl;
- MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
- MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
if(_verbose){
displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
_AroeMinusImplicit[k] *= -1;
_Diffusion[k] *= -1;
}
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
if(_verbose){
displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
_AroePlusImplicit[k] *= _inv_dxj;
_Diffusion[k] *=_inv_dxj/_inv_dxi;
}
- MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
- MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
if(_verbose)
displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
_AroePlusImplicit[k] *= -1;
_Diffusion[k] *= -1;
}
- MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
- MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
if(_verbose)
displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
- }
}
else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
*_runLogFile<<"Warning: treatment of a junction node"<<endl;
_vec_normal[0] = 1;
else//j==16
_vec_normal[0] = -1;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
cout << " ; vecteur normal=(";
convectionMatrices();
diffusionStateAndMatrices(idm, idn,false);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
_inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
_inv_dxj = 1/Ctemp2.getMeasure();
- addConvectionToSecondMember(idm,idn, false);
- _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
- addDiffusionToSecondMember(idm,idn, false);
- _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
- addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
-
- if(_timeScheme == Implicit){
for(int k=0; k<_nVar*_nVar;k++)
{
_AroeMinusImplicit[k] *= _inv_dxi;
_Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
}
- MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
- MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
if(_verbose){
displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
_AroeMinusImplicit[k] *= -1;
_Diffusion[k] *= -1;
}
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
if(_verbose){
displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
_AroePlusImplicit[k] *= _inv_dxj;
_Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
}
- MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
- MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
if(_verbose)
displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
_AroePlusImplicit[k] *= -1;
_Diffusion[k] *= -1;
}
- MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
- MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
+ MatSetValuesBlocked(A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
if(_verbose)
displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
}
}
- }
}
else
{
}
}
- VecAssemblyEnd(_conservativeVars);
- VecAssemblyEnd(_b);
-
- if(_timeScheme == Implicit){
- for(int imaille = 0; imaille<_Nmailles; imaille++)
- MatSetValuesBlocked(_A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
-
- if(_verbose && _nbTimeStep%_freqSave ==0)
- displayMatrix(_GravityImplicitationMatrix,_nVar,"Gravity matrix:");
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
- if(_verbose && _nbTimeStep%_freqSave ==0){
- cout << "Fin ProblemFluid.cxx : matrice implicite"<<endl;
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- }
- }
+ for(int imaille = 0; imaille<_Nmailles; imaille++)
+ MatSetValuesBlocked(A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
- stop=false;
+ MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
- /*
- if(_nbTimeStep+1<_cfl)
- return (_nbTimeStep+1)*_minl/_maxvp;
- else
- */
- return _cfl*_minl/_maxvp;
-}
+ MatShift(A, 1/_dt);
-void ProblemFluid::computeNewtonVariation()
-{
- if(_verbose)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
- cout<<"Vecteur courant Uk "<<endl;
- VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
+ cout << "ProblemFluid::computeNewtonJacobian : Fin calcul Jacobienne schéma Newton pour PETSc"<<endl;
+ MatView(A,PETSC_VIEWER_STDOUT_SELF);
cout << endl;
}
- if(_timeScheme == Explicit)
- {
- VecCopy(_b,_newtonVariation);
- VecScale(_newtonVariation, _dt);
- if(_verbose && _nbTimeStep%_freqSave ==0)
- {
- cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
- VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- }
- }
- else
- {
- if(_verbose)
- {
- cout << "Matrice du système linéaire avant contribution delta t" << endl;
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- cout << "Second membre du système linéaire avant contribution delta t" << endl;
- VecView(_b, PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- }
- MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
-
- VecAXPY(_b, 1/_dt, _old);
- VecAXPY(_b, -1/_dt, _conservativeVars);
- MatShift(_A, 1/_dt);
-
-#if PETSC_VERSION_GREATER_3_5
- KSPSetOperators(_ksp, _A, _A);
-#else
- KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
-#endif
-
- if(_verbose)
- {
- cout << "Matrice du système linéaire après contribution delta t" << endl;
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- cout << "Second membre du système linéaire après contribution delta t" << endl;
- VecView(_b, PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- }
+}
- if(_conditionNumber)
- KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
- if(!_isScaling)
- {
- KSPSolve(_ksp, _b, _newtonVariation);
- }
- else
- {
- computeScaling(_maxvp);
- int indice;
- VecAssemblyBegin(_vecScaling);
- VecAssemblyBegin(_invVecScaling);
- for(int imaille = 0; imaille<_Nmailles; imaille++)
- {
- indice = imaille;
- VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
- VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
- }
- VecAssemblyEnd(_vecScaling);
- VecAssemblyEnd(_invVecScaling);
+int ProblemFluid::computeSnesJacobian(SNES snes, Vec X, Mat A, Mat Aapprox, void *ctx)//Propotype imposé par PETSc
+{
+ ProblemFluid * myProblem = (ProblemFluid *) ctx;
+ myProblem->computeNewtonJacobian( X, A);
- if(_system)
- {
- cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
- VecView(_b, PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- }
- MatDiagonalScale(_A,_vecScaling,_invVecScaling);
- if(_system)
- {
- cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
- MatView(_A,PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- }
- VecCopy(_b,_bScaling);
- VecPointwiseMult(_b,_vecScaling,_bScaling);
- if(_system)
- {
- cout << "Produit du second membre par le preconditionneur bloc diagonal a gauche" << endl;
- VecView(_b, PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- }
+ Aapprox = A;
- KSPSolve(_ksp,_b, _bScaling);
- VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
- }
- if(_verbose)
- {
- cout << "solution du systeme lineaire local:" << endl;
- VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
- cout << endl;
- }
- }
+ return 0;
}
void ProblemFluid::validateTimeStep()
VecCopy(_conservativeVars, _old);
- if(_verbose && _nbTimeStep%_freqSave ==0){
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
if(!_usePrimitiveVarsInNewton)
testConservation();
cout <<"Valeur propre locale max: " << _maxvp << endl;
}
- if(_nbPhases==2 && _nbTimeStep%_freqSave ==0){
+ if(_nbPhases==2 && (_nbTimeStep-1)%_freqSave ==0){
//Find minimum and maximum void fractions
double alphamin=1e30;
double alphamax=-1e30;
_time+=_dt;
_nbTimeStep++;
- if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
+ if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
save();
}
const int &j, bool isBord, string groupname
)
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"ProblemFluid::addConvectionToSecondMember start"<<endl;
//extraction des valeurs
_idm[0] = i;
_idn[0] = j;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "addConvectionToSecondMember : état i= " << i << ", _Ui=" << endl;
for(int q=0; q<_nVar; q++)
Fij=convectionFlux(Uij,Vij,normale,porosityij);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"Etat interfacial conservatif "<<i<<", "<<j<< endl;
cout<<Uij<<endl;
else if(_nonLinearFormulation==Roe)//Roe
Fij=(Fi+Fj)/2+absAroe*(Ui-Uj)/2;
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"Flux cellule "<<i<<" = "<< endl;
cout<<Fi<<endl;
for(int i1=0;i1<_nVar;i1++)
_phi[i1]=-Fij(i1)*_inv_dxi;
VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
cout<<"Flux interfacial "<<i<<", "<<j<< endl;
for(int i1=0;i1<_nVar;i1++)
_phi[i1]*=-_inv_dxj/_inv_dxi;
VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution convection à " << j << ", Fij(i1)*_inv_dxj= "<<endl;
for(int q=0; q<_nVar; q++)
Polynoms::matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
cout << "(Ui - Uj)*_inv_dxi= "<<endl;;
Polynoms::matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution convection à " << j << ", A^+*(Ui - Uj)*_inv_dxi= "<<endl;
for(int q=0; q<_nVar; q++)
}
}
}
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"ProblemFluid::addConvectionToSecondMember end : matrices de décentrement cellules i= " << i << ", et j= " << j<< "):"<<endl;
displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
const int j, int nbVoisinsj,
bool isBord, int ij, double mesureFace)//To do : generalise to unstructured meshes
{
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
cout<<"ProblemFluid::addSourceTerm cell i= "<<i<< " cell j= "<< j<< " isbord "<<isBord<<endl;
_idm[0] = i*_nVar;
VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
sourceVector(_Si,_Ui,_Vi,i);
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Terme source S(Ui), i= " << i<<endl;
for(int q=0; q<_nVar; q++)
consToPrim(_Uj, _Vj,_porosityj);
sourceVector(_Sj,_Uj,_Vj,i);
}
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
if(!isBord)
cout << "Terme source S(Uj), j= " << j<<endl;
for(int k=0; k<_nVar;k++)
_pressureLossVector[k]=0;
}
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", K="<<K<<endl;
for(int k=0; k<_nVar;k++)
_porosityGradientSourceVector[k]=0;
}
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
if(!isBord)
cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<", dxj= "<<1/_inv_dxj<<endl;
_Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
_Sj[k]=(_phi[k]+_l[k])*mesureFace/_perimeters(j);///nbVoisinsj;
}
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution au terme source Si de la cellule i= " << i<<" venant (après décentrement) de la face (i,j), j="<<j<<endl;
for(int q=0; q<_nVar; q++)
_Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i)
_Sj[k]=_Sj[k]/nbVoisinsj+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(j)
}
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,j), j="<<j<<endl;
for(int q=0; q<_nVar; q++)
Polynoms::matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
for(int k=0; k<_nVar;k++)
_Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution au terme source Si de la cellule i= " << i<<" venant (après décentrement) de la face (i,bord)"<<endl;
for(int q=0; q<_nVar; q++)
{
for(int k=0; k<_nVar;k++)
_Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i);//
- if (_verbose && _nbTimeStep%_freqSave ==0)
+ if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,bord) "<<endl;
for(int q=0; q<_nVar; q++)
_idm[0] = i;
VecSetValuesBlocked(_b, 1, _idm, _Si, ADD_VALUES);
- if(_verbose && _nbTimeStep%_freqSave ==0 && _wellBalancedCorrection)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0 && _wellBalancedCorrection)
displayMatrix( _signAroe,_nVar,"Signe matrice de Roe");
}
VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
consToPrim(_Ui,_Vi,_porosityField(i-1));
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "ProblemFluid::updatePrimitives() cell " << i-1 << endl;
cout << "Ui = ";
_idm[0] = i-1;
VecSetValuesBlocked(_conservativeVars, 1, _idm, _Ui, INSERT_VALUES);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout << "ProblemFluid::updateConservatives() cell " << i-1 << endl;
cout << "Vi = ";
KSPDestroy(&_ksp);
for(int i=0;i<_nbPhases;i++)
delete _fluides[i];
+
+ // Destruction du solveur de Newton de PETSc
+ if(_timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
+ SNESDestroy(&_snes);
}
