#include "ProblemFluid.hxx"
-#include "Fluide.h"
#include "math.h"
+#include <numeric>
using namespace std;
-ProblemFluid::ProblemFluid(void){
+ProblemFluid::ProblemFluid(MPI_Comm comm):ProblemCoreFlows(comm)
+{
_latentHeat=1e30;
_Tsat=1e30;
_Psat=-1e30;
void ProblemFluid::initialize()
{
- if(!_initialDataSet){
- *_runLogFile<<"ProblemFluid::initialize() set initial data first"<<endl;
+ if(!_initialDataSet)
+ {
+ *_runLogFile<<"!!!!!!!!ProblemFluid::initialize() set initial data first"<<endl;
+ _runLogFile->close();
+ throw CdmathException("!!!!!!!!ProblemFluid::initialize() set initial data first");
+ }
+ else if (_VV.getTypeOfField() != CELLS)
+ {
+ *_runLogFile<<"!!!!!!!!Initial data should be a field on CELLS, not NODES, neither FACES"<<endl;
_runLogFile->close();
- throw CdmathException("ProblemFluid::initialize() set initial data first");
+ throw CdmathException("!!!!!!!!ProblemFluid::initialize() Initial data should be a field on CELLS, not NODES, neither FACES");
}
cout << "Number of Phases = " << _nbPhases << " mesh dimension = "<<_Ndim<<" number of variables = "<<_nVar<<endl;
*_runLogFile << "Number of Phases = " << _nbPhases << " spaceDim= "<<_Ndim<<" number of variables= "<<_nVar<<endl;
//creation de la matrice
if(_timeScheme == Implicit)
- MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNb), PETSC_NULL, &_A);
+ MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
//creation des vecteurs
VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uext);
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)
{
int *indices = new int[_Nmailles];
- for(int i=0; i<_Nmailles; i++)
- indices[i] = i;
+ std::iota(indices, indices +_Nmailles, 0);
VecSetValuesBlocked(_conservativeVars, _Nmailles, indices, initialFieldCons, INSERT_VALUES);
VecAssemblyBegin(_conservativeVars);
VecAssemblyEnd(_conservativeVars);
delete[] initialFieldCons;
delete[] indices;
- //Linear solver
- KSPCreate(PETSC_COMM_SELF, &_ksp);
- KSPSetType(_ksp, _ksptype);
- // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
- KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
- KSPGetPC(_ksp, &_pc);
- PCSetType(_pc, _pctype);
+ createKSP();
+
+ // 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;
VecAXPY(_conservativeVars, relaxation, _newtonVariation);
- //mise a jour du champ primitif
- updatePrimitives();
+ //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+_neibMaxNbCells), 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= "<<_neibMaxNbCells<<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= "<<_neibMaxNbCells<<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= "<<_neibMaxNbCells<<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;
- 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++)
cout << endl;
}
idm = idCells[k];
- Polynoms Poly;
//calcul et insertion de A^-*Jcb
- Poly.matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
+ 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");
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
- Poly.matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
- MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
+ 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);
+ 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);
+ //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNbCells<<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: ");
_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)
- cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNb<<endl;
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
+ cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
*_runLogFile<<"Warning: treatment of a junction node"<<endl;
if(!_sectionFieldSet)
_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++)
{
for(int k=0; k<_nVar; k++)
_temp[k]=(_Ui[k] - _Uj[k])*_inv_dxi;//(Ui-Uj)*_inv_dxi
- Polynoms Poly;
- Poly.matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
+ 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;;
{
for(int k=0; k<_nVar; k++)
_temp[k]*=_inv_dxj/_inv_dxi;//(Ui-Uj)*_inv_dxj
- Poly.matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
+ 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;
cout<< _porosityGradientSourceVector[k]<<", ";
cout<<endl;
}
- Polynoms Poly;
+
if(!isBord){
if(_wellBalancedCorrection){
for(int k=0; k<_nVar;k++)
_phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
- Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
+ Polynoms::matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
for(int k=0; k<_nVar;k++){
_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++)
if(_wellBalancedCorrection){
for(int k=0; k<_nVar;k++)
_phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
- Poly.matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
+ 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 = ";
void ProblemFluid::AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist)
{
- Polynoms Poly;
int nbVp_dist=valeurs_propres_dist.size();
vector< complex< double > > y (nbVp_dist,0);
for( int i=0 ; i<nbVp_dist ; i++)
- y[i] = Poly.abs_generalise(valeurs_propres_dist[i]);
- Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ y[i] = Polynoms::abs_generalise(valeurs_propres_dist[i]);
+ Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<< endl<<"ProblemFluid::AbsMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
for(int ct =0; ct<nbVp_dist; ct++)
else
y[i] = 0;
}
- Polynoms Poly;
- Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+
+ Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<< endl<<"ProblemFluid::SigneMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
for(int ct =0; ct<nbVp_dist; ct++)
{
int nbVp_dist=valeurs_propres_dist.size();
vector< complex< double > > y (nbVp_dist,0);
- Polynoms Poly;
+
for( int i=0 ; i<nbVp_dist ; i++)
{
- if(Poly.module(valeurs_propres_dist[i])>_precision)
+ if(Polynoms::module(valeurs_propres_dist[i])>_precision)
y[i] = 1./valeurs_propres_dist[i];
else
y[i] = 1./_precision;
}
- Poly.abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
- if(_verbose && _nbTimeStep%_freqSave ==0)
+ Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
+ if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
{
cout<< endl<<"ProblemFluid::InvMatriceRoe : Valeurs propres :" << nbVp_dist<<endl;
for(int ct =0; ct<nbVp_dist; ct++)
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);
+}
+
+vector<string>
+ProblemFluid::getInputFieldsNames()
+{
+ vector<string> result(1);
+
+ result[0]="HeatPower";
+ result[1]="Porosity";
+ result[2]="PressureLoss";
+ result[3]="Section";
+
+ return result;
+}
+
+void
+ProblemFluid::setInputField(const string& nameField, Field& inputField )
+{
+ if(nameField=="Porosity" || nameField=="POROSITY" || nameField=="Porosité" || nameField=="POROSITE")
+ return setPorosityField( inputField) ;
+ else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
+ return setHeatPowerField( inputField );
+ else if(nameField=="PressureLoss" || nameField=="PRESSURELOSS" || nameField=="PerteDeCharge" || nameField=="PPERTEDECHARGE")
+ return setPressureLossField( inputField) ;
+ else if(nameField=="Section" || nameField=="SECTION" )
+ return setSectionField( inputField );
+ else
+ {
+ cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
+ throw CdmathException("SinglePhase::setInputField error : Unknown Field name");
+ }
+}
+
+void
+ProblemFluid::setPorosityField(string fileName, string fieldName){
+ if(!_initialDataSet)
+ throw CdmathException("!!!!!!!! ProblemFluid::setPorosityField set initial field first");
+
+ _porosityField=Field(fileName, CELLS,fieldName);
+ _porosityField.getMesh().checkFastEquivalWith(_mesh);
+ _porosityFieldSet=true;
+}
+
+void
+ProblemFluid::setPressureLossField(Field PressureLoss){
+ if(!_initialDataSet)
+ throw CdmathException("!!!!!!!! ProblemFluid::setPressureLossField set initial field first");
+
+ PressureLoss.getMesh().checkFastEquivalWith(_mesh);
+ _pressureLossField=PressureLoss;
+ _pressureLossFieldSet=true;
+}
+
+void
+ProblemFluid::setPressureLossField(string fileName, string fieldName){
+ if(!_initialDataSet)
+ throw CdmathException("!!!!!!!! ProblemFluid::setPressureLossField set initial field first");
+
+ _pressureLossField=Field(fileName, FACES,fieldName);
+ _pressureLossField.getMesh().checkFastEquivalWith(_mesh);
+ _pressureLossFieldSet=true;
+}
+
+void
+ProblemFluid::setSectionField(Field sectionField){
+ if(!_initialDataSet)
+ throw CdmathException("!!!!!!!! ProblemFluid::setSectionField set initial field first");
+
+ sectionField.getMesh().checkFastEquivalWith(_mesh);
+ _sectionField=sectionField;
+ _sectionFieldSet=true;
+}
+
+void
+ProblemFluid::setSectionField(string fileName, string fieldName){
+ if(!_initialDataSet)
+ throw CdmathException("!!!!!!!! ProblemFluid::setSectionField set initial field first");
+
+ _sectionField=Field(fileName, CELLS,fieldName);
+ _sectionField.getMesh().checkFastEquivalWith(_mesh);
+ _sectionFieldSet=true;
+}
+
+void
+ProblemFluid::setPorosityField(Field Porosity){
+ if(!_initialDataSet)
+ throw CdmathException("!!!!!!!! ProblemFluid::setPorosityField set initial field first");
+
+ Porosity.getMesh().checkFastEquivalWith(_mesh);
+ _porosityField=Porosity;
+ _porosityFieldSet=true;
}
