From 3a28c4bd2b715f2b77af249d781eac4a0d0bb3c3 Mon Sep 17 00:00:00 2001 From: michael Date: Mon, 17 May 2021 16:57:14 +0200 Subject: [PATCH] Changed message display during simulation --- CoreFlows/Models/inc/ProblemFluid.hxx | 51 +- CoreFlows/Models/src/DiffusionEquation.cxx | 6 +- CoreFlows/Models/src/DriftModel.cxx | 132 +-- CoreFlows/Models/src/FiveEqsTwoFluid.cxx | 32 +- CoreFlows/Models/src/IsothermalTwoFluid.cxx | 40 +- CoreFlows/Models/src/ProblemCoreFlows.cxx | 4 +- CoreFlows/Models/src/ProblemFluid.cxx | 1176 +++++++++++++++---- CoreFlows/Models/src/TransportEquation.cxx | 8 +- 8 files changed, 1081 insertions(+), 368 deletions(-) diff --git a/CoreFlows/Models/inc/ProblemFluid.hxx b/CoreFlows/Models/inc/ProblemFluid.hxx index 15f31a4..f0100f4 100755 --- a/CoreFlows/Models/inc/ProblemFluid.hxx +++ b/CoreFlows/Models/inc/ProblemFluid.hxx @@ -142,6 +142,14 @@ public : * */ virtual void validateTimeStep(); + /** \fn solveTimeStep + * \brief calcule les valeurs inconnues au pas de temps +1 . + * \details c'est une fonction virtuelle, qui surcharge celle de la classe ProblemCoreFlows + * @param void + * \return Renvoie false en cas de problème durant le calcul (valeurs non physiques..) + * */ + virtual bool solveTimeStep();// + /* Boundary conditions */ /** \fn setNeumannBoundaryCondition * \brief adds a new boundary condition of type Neumann @@ -244,12 +252,6 @@ public : return _nbPhases; }; - /** \fn computeNewtonVariation - * \brief Builds and solves the linear system to obtain the variation Ukp1-Uk in a Newton scheme - * @param void - * */ - virtual void computeNewtonVariation(); - /** \fn testConservation * \brief Teste et affiche la conservation de masse et de la quantité de mouvement * \Details la fonction est virtuelle pure, on la surcharge dans chacun des modèles @@ -495,6 +497,11 @@ protected : /** the formulation used to compute the non viscous fluxes **/ NonLinearFormulation _nonLinearFormulation; + /** PETSc nonlinear solver and line search **/ + SNES _snes; + SNESLineSearch _linesearch; + PetscViewer _monitorLineSearch; + map _limitField; /** boolean used to specify that an entropic correction should be used **/ @@ -542,6 +549,38 @@ protected : bool _isBoundary;// la face courante est elle une face de bord ? double _maxvploc; + bool solveNewtonPETSc();//Use PETSc Newton methods to solve time step + + /** \fn computeNewtonVariation + * \brief Builds and solves the linear system to obtain the variation Ukp1-Uk in a Newton scheme + * @param void + * */ + virtual void computeNewtonVariation(); + + /** \fn computeNewtonRHS + * \brief Builds the right hand side F_X(X) of the linear system in the Newton method to obtain the variation Ukp1-Uk + * @param void + * */ + void computeNewtonRHS( Vec X, Vec F_X); + + /** \fn computeSnesRHS + * \brief Static function calling computeNewtonRHS to match PETSc nonlinear solver (SNES) structure + * @param void + * */ + static int computeSnesRHS(SNES snes, Vec X, Vec F_X, void *ctx); + + /** \fn computeNewtonJacobian + * \brief Static function calling computeNewtonJacobian to match PETSc nonlinear solver (SNES) structure + * @param void + * */ + void computeNewtonJacobian( Vec X, Mat A); + + /** \fn computeSnesJacobian + * \brief Builds the matrix A(X) of the linear system in the Newton method to obtain the variation Ukp1-Uk + * @param void + * */ + static int computeSnesJacobian(SNES snes, Vec X, Mat A, Mat Aapprox, void *ctx); + /** \fn convectionState * \brief calcule l'etat de Roe entre deux cellules voisinnes * \Details c'ets une fonction virtuelle, on la surcharge dans chacun des modèles diff --git a/CoreFlows/Models/src/DiffusionEquation.cxx b/CoreFlows/Models/src/DiffusionEquation.cxx index 59cad62..c93fef8 100755 --- a/CoreFlows/Models/src/DiffusionEquation.cxx +++ b/CoreFlows/Models/src/DiffusionEquation.cxx @@ -569,7 +569,7 @@ bool DiffusionEquation::initTimeStep(double dt){ _dt = dt; - if(_verbose && _nbTimeStep%_freqSave ==0) + if(_verbose && (_nbTimeStep-1)%_freqSave ==0) MatView(_A,PETSC_VIEWER_STDOUT_SELF); return _dt>0; @@ -692,12 +692,12 @@ void DiffusionEquation::validateTimeStep() VecCopy(_Tk, _Tn); VecCopy(_Tk, _Tkm1); - if(_verbose && _nbTimeStep%_freqSave ==0) + if(_verbose && (_nbTimeStep-1)%_freqSave ==0) cout <<"Valeur propre locale max: " << _maxvp << endl; _time+=_dt; _nbTimeStep++; - if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep) + if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep) save(); } diff --git a/CoreFlows/Models/src/DriftModel.cxx b/CoreFlows/Models/src/DriftModel.cxx index 7df3613..baa394d 100755 --- a/CoreFlows/Models/src/DriftModel.cxx +++ b/CoreFlows/Models/src/DriftModel.cxx @@ -195,7 +195,7 @@ bool DriftModel::iterateTimeStep(bool &converged) 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<" << _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 @@ -413,7 +413,7 @@ void DriftModel::convectionState( const long &i, const long &j, const bool &IsBo 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; @@ -428,7 +428,7 @@ void DriftModel::convectionState( const long &i, const long &j, const bool &IsBo 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"<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 "<_precision) { Hm=Emi+pj/rhomi; - if(_verbose && _nbTimeStep%_freqSave ==0) + if(_verbose && (_nbTimeStep-1)%_freqSave ==0) cout<<"VFFC Staggered state rhomi="<_precision) { - if(_verbose && _nbTimeStep%_freqSave ==0) + if(_verbose && (_nbTimeStep-1)%_freqSave ==0) cout<<"VFFC Staggered state rhomi="< >(taille_vp); for( int i=0 ; i > vp_right = getRacines(pol_car); int taille_vp_right = Polynoms::new_tri_selectif(vp_right,vp_right.size(),_precision); - if(_verbose && _nbTimeStep%_freqSave ==0) + if(_verbose && (_nbTimeStep-1)%_freqSave ==0) { cout<<"Entropic shift right eigenvalues: "<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 } @@ -207,13 +249,58 @@ bool ProblemFluid::initTimeStep(double dt){ 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)"<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"< _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<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"< idCells(2); + PetscInt idm, idn, size = 1; + + long nbFaces = _mesh.getNumberOfFaces(); + Face Fj; + Cell Ctemp1,Ctemp2; + string nameOfGroup; + + for (int j=0; j1){ + for(int l=0; l 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; lidCells[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< 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= "<2 && _Ndim==1 ){//inner face with more than two neighbours + if(_verbose && (_nbTimeStep-1)%_freqSave ==0) + cout<<"lattice mesh junction at face "<close(); + throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field"); + } + int largestSectionCellIndex=0; + for(int i=1;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 "<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"< idCells(2); + PetscInt idm, idn, size = 1; + + long nbFaces = _mesh.getNumberOfFaces(); + Face Fj; + Cell Ctemp1,Ctemp2; + string nameOfGroup; + + for (int j=0; j1){ + for(int l=0; l 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; lidCells[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< 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 "<close(); + throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field"); + } + int largestSectionCellIndex=0; + for(int i=1;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 "< _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<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"<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< 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(); @@ -374,7 +1181,7 @@ double ProblemFluid::computeTimeStep(bool & stop){ _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++) @@ -394,11 +1201,6 @@ double ProblemFluid::computeTimeStep(bool & stop){ 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; @@ -406,7 +1208,7 @@ double ProblemFluid::computeTimeStep(bool & stop){ 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++) @@ -425,7 +1227,7 @@ double ProblemFluid::computeTimeStep(bool & stop){ 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"); @@ -434,18 +1236,17 @@ double ProblemFluid::computeTimeStep(bool & stop){ 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 @@ -474,7 +1275,7 @@ double ProblemFluid::computeTimeStep(bool & stop){ _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= "; @@ -500,11 +1301,6 @@ double ProblemFluid::computeTimeStep(bool & stop){ _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; @@ -513,8 +1309,8 @@ double ProblemFluid::computeTimeStep(bool & stop){ idm = idCells[0]; idn = idCells[1]; //cout<<"idm= "<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 "<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() @@ -834,13 +1508,13 @@ 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; @@ -867,7 +1541,7 @@ void ProblemFluid::validateTimeStep() _time+=_dt; _nbTimeStep++; - if (_nbTimeStep%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep) + if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep) save(); } @@ -880,7 +1554,7 @@ void ProblemFluid::addConvectionToSecondMember const int &j, bool isBord, string groupname ) { - if(_verbose && _nbTimeStep%_freqSave ==0) + if(_verbose && (_nbTimeStep-1)%_freqSave ==0) cout<<"ProblemFluid::addConvectionToSecondMember start"<