1 #include "ProblemFluid.hxx"
6 ProblemFluid::ProblemFluid(MPI_Comm comm):ProblemCoreFlows(comm)
12 _porosityFieldSet=false;
13 _pressureLossFieldSet=false;
14 _sectionFieldSet=false;
15 _GravityField3d=vector<double>(3,0);
16 _gravityReferencePoint=vector<double>(3,0);
17 _Uroe=NULL;_Udiff=NULL;_temp=NULL;_l=NULL;_vec_normal=NULL;;
20 _saveConservativeField=false;
21 _usePrimitiveVarsInNewton=false;
22 _saveInterfacialField=false;
23 _err_press_max=0; _part_imag_max=0; _nbMaillesNeg=0; _nbVpCplx=0;_minm1=1e30;_minm2=1e30;//Pour affichage paramètres diphasiques dans IterateTimeStep
25 _entropicCorrection=false;
26 _pressureCorrectionOrder=2;
27 _nonLinearFormulation=Roe;
32 SpaceScheme ProblemFluid::getSpaceScheme() const
36 void ProblemFluid::setNumericalScheme(SpaceScheme spaceScheme, TimeScheme timeScheme)
38 if( _nbTimeStep>0 && timeScheme!=_timeScheme)//This is a change of time scheme during a simulation
39 _restartWithNewTimeScheme=true;
40 _timeScheme = timeScheme;
41 _spaceScheme = spaceScheme;
44 void ProblemFluid::initialize()
48 *_runLogFile<<"!!!!!!!!ProblemFluid::initialize() set initial data first"<<endl;
50 throw CdmathException("!!!!!!!!ProblemFluid::initialize() set initial data first");
52 else if (_VV.getTypeOfField() != CELLS)
54 *_runLogFile<<"!!!!!!!!Initial data should be a field on CELLS, not NODES, neither FACES"<<endl;
56 throw CdmathException("!!!!!!!!ProblemFluid::initialize() Initial data should be a field on CELLS, not NODES, neither FACES");
58 cout << "Number of Phases = " << _nbPhases << " mesh dimension = "<<_Ndim<<" number of variables = "<<_nVar<<endl;
59 *_runLogFile << "Number of Phases = " << _nbPhases << " spaceDim= "<<_Ndim<<" number of variables= "<<_nVar<<endl;
61 /********* local arrays ****************/
62 _AroePlus = new PetscScalar[_nVar*_nVar];
63 _Diffusion = new PetscScalar[_nVar*_nVar];
64 _AroeMinus = new PetscScalar[_nVar*_nVar];
65 _Aroe = new PetscScalar[_nVar*_nVar];
66 _absAroe = new PetscScalar[_nVar*_nVar];
67 _signAroe = new PetscScalar[_nVar*_nVar];
68 _invAroe = new PetscScalar[_nVar*_nVar];
69 _AroeImplicit = new PetscScalar[_nVar*_nVar];
70 _AroeMinusImplicit = new PetscScalar[_nVar*_nVar];
71 _AroePlusImplicit = new PetscScalar[_nVar*_nVar];
72 _absAroeImplicit = new PetscScalar[_nVar*_nVar];
73 _primToConsJacoMat = new PetscScalar[_nVar*_nVar];
74 _phi = new PetscScalar[_nVar];
75 _Jcb = new PetscScalar[_nVar*_nVar];
76 _JcbDiff = new PetscScalar[_nVar*_nVar];
77 _a = new PetscScalar[_nVar*_nVar];
78 _externalStates = new PetscScalar[_nVar];
79 _Ui = new PetscScalar[_nVar];
80 _Uj = new PetscScalar[_nVar];
81 _Vi = new PetscScalar[_nVar];
82 _Vj = new PetscScalar[_nVar];
83 _Uij = new PetscScalar[_nVar];//used for VFRoe scheme
84 _Vij = new PetscScalar[_nVar];//used for VFRoe scheme
85 _Si = new PetscScalar[_nVar];
86 _Sj = new PetscScalar[_nVar];
87 _pressureLossVector= new PetscScalar[_nVar];
88 _porosityGradientSourceVector= new PetscScalar[_nVar];
90 _l = new double[_nVar];
91 _r = new double[_nVar];
92 _Udiff = new double[_nVar];
93 _temp = new double[_nVar*_nVar];
95 _idm = new int[_nVar];
96 _idn = new int[_nVar];
97 _vec_normal = new double[_Ndim];
99 for(int k=0;k<_nVar*_nVar;k++)
104 for(int k=0; k<_nVar; k++){
109 /**************** Field creation *********************/
110 if(!_heatPowerFieldSet){
111 _heatPowerField=Field("Heat Power",CELLS,_mesh,1);
112 for(int i =0; i<_Nmailles; i++)
113 _heatPowerField(i) = _heatSource;
116 _dp_over_dt=Field("dP/dt",CELLS,_mesh,1);
117 if(!_porosityFieldSet){
118 _porosityField=Field("Porosity field",CELLS,_mesh,1);
119 for(int i =0; i<_Nmailles; i++)
120 _porosityField(i) = 1;
121 _porosityFieldSet=true;
124 //conservative field used only for saving results
125 _UU=Field ("Conservative vec", CELLS, _mesh, _nVar);
126 if(_saveInterfacialField && _nonLinearFormulation==VFRoe)
128 _UUstar=Field ("Interfacial U", CELLS, _mesh, _nVar);
129 _VVstar=Field ("Interfacial V", CELLS, _mesh, _nVar);
132 //Construction des champs primitifs et conservatifs initiaux comme avant dans ParaFlow
133 double * initialFieldPrim = new double[_nVar*_Nmailles];
134 for(int i =0; i<_Nmailles;i++)
135 for(int j =0; j<_nVar;j++)
136 initialFieldPrim[i*_nVar+j]=_VV(i,j);
137 double *initialFieldCons = new double[_nVar*_Nmailles];
138 for(int i=0; i<_Nmailles; i++)
139 primToCons(initialFieldPrim, i, initialFieldCons, i);
140 for(int i =0; i<_Nmailles;i++)
141 for(int j =0; j<_nVar;j++)
142 _UU(i,j)=initialFieldCons[i*_nVar+j];
144 /**********Petsc structures: ****************/
146 //creation de la matrice
147 if(_timeScheme == Implicit)
148 MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
150 //creation des vecteurs
151 VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uext);
152 VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Vext);
153 VecCreateSeq(PETSC_COMM_SELF, _nVar, &_Uextdiff);
154 // VecCreateSeq(PETSC_COMM_SELF, _nVar*_Nmailles, &_conservativeVars);
155 VecCreate(PETSC_COMM_SELF, &_conservativeVars);//Current conservative variables at Newton iteration k between time steps n and n+1
156 VecSetSizes(_conservativeVars,PETSC_DECIDE,_nVar*_Nmailles);
157 VecSetBlockSize(_conservativeVars,_nVar);
158 VecSetFromOptions(_conservativeVars);
159 VecDuplicate(_conservativeVars, &_old);//Old conservative variables at time step n
160 VecDuplicate(_conservativeVars, &_newtonVariation);//Newton variation Uk+1-Uk to be computed between time steps n and n+1
161 VecDuplicate(_conservativeVars, &_b);//Right hand side of Newton method at iteration k between time steps n and n+1
162 VecDuplicate(_conservativeVars, &_primitiveVars);//Current primitive variables at Newton iteration k between time steps n and n+1
166 VecDuplicate(_conservativeVars, &_vecScaling);
167 VecDuplicate(_conservativeVars, &_invVecScaling);
168 VecDuplicate(_conservativeVars, &_bScaling);
170 _blockDiag = new PetscScalar[_nVar];
171 _invBlockDiag = new PetscScalar[_nVar];
175 int *indices = new int[_Nmailles];
176 for(int i=0; i<_Nmailles; i++)
178 VecSetValuesBlocked(_conservativeVars, _Nmailles, indices, initialFieldCons, INSERT_VALUES);
179 VecAssemblyBegin(_conservativeVars);
180 VecAssemblyEnd(_conservativeVars);
181 VecCopy(_conservativeVars, _old);
182 VecAssemblyBegin(_old);
183 VecAssemblyEnd(_old);
184 VecSetValuesBlocked(_primitiveVars, _Nmailles, indices, initialFieldPrim, INSERT_VALUES);
185 VecAssemblyBegin(_primitiveVars);
186 VecAssemblyEnd(_primitiveVars);
189 cout << "Variables primitives initiales : " << endl;
190 VecView(_primitiveVars, PETSC_VIEWER_STDOUT_WORLD);
192 cout<<"Variables conservatives initiales : "<<endl;
193 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_SELF);
196 delete[] initialFieldPrim;
197 delete[] initialFieldCons;
201 KSPCreate(PETSC_COMM_SELF, &_ksp);
202 KSPSetType(_ksp, _ksptype);
203 // if(_ksptype == KSPGMRES) KSPGMRESSetRestart(_ksp,10000);
204 KSPSetTolerances(_ksp,_precision,_precision,PETSC_DEFAULT,_maxPetscIts);
205 KSPGetPC(_ksp, &_pc);
206 PCSetType(_pc, _pctype);
208 // Creation du solveur de Newton de PETSc
209 if( _timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
213 // set nonlinear solver
214 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)
215 snestype = (char*)&SNESNEWTONLS;
216 else if (_nonLinearSolver == Newton_PETSC_TRUSTREGION)
217 snestype = (char*)&SNESNEWTONTR;
218 else if (_nonLinearSolver == Newton_PETSC_NGMRES)
219 snestype = (char*)&SNESNGMRES;
220 else if (_nonLinearSolver ==Newton_PETSC_ASPIN)
221 snestype = (char*)&SNESASPIN;
222 else if(_nonLinearSolver != Newton_SOLVERLAB)
224 cout << "!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" << endl;
225 *_runLogFile << "!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" << endl;
226 _runLogFile->close();
227 throw CdmathException("!!! Error : only 'Newton_PETSC_LINESEARCH', 'Newton_PETSC_TRUSTREGION', 'Newton_PETSC_NGMRES', 'Newton_PETSC_ASPIN' or 'Newton_SOLVERLAB' nonlinear solvers are acceptable !!!" );
230 SNESCreate(PETSC_COMM_WORLD, &_snes);
231 SNESSetType( _snes, snestype);
232 SNESGetLineSearch( _snes, &_linesearch);
233 if(_nonLinearSolver == Newton_PETSC_LINESEARCH_BASIC)
234 SNESLineSearchSetType( _linesearch, SNESLINESEARCHBASIC );
235 else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_BT)
236 SNESLineSearchSetType( _linesearch, SNESLINESEARCHBT );
237 else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_SECANT)
238 SNESLineSearchSetType( _linesearch, SNESLINESEARCHL2 );
239 else if(_nonLinearSolver == Newton_PETSC_LINESEARCH_NLEQERR)
240 SNESLineSearchSetType( _linesearch, SNESLINESEARCHNLEQERR );
242 PetscViewerCreate(PETSC_COMM_WORLD,&_monitorLineSearch);
243 PetscViewerSetType(_monitorLineSearch, PETSCVIEWERASCII);
245 SNESSetTolerances(_snes,_precision_Newton,_precision_Newton,_precision_Newton,_maxNewtonIts,-1);
247 SNESSetFunction(_snes,_newtonVariation,computeSnesRHS,this);
248 SNESSetJacobian(_snes,_A,_A,computeSnesJacobian,this);
251 _initializedMemory=true;
252 save();//save initial data
255 bool ProblemFluid::initTimeStep(double dt){
257 return _dt>0;//No need to call MatShift as the linear system matrix is filled at each Newton iteration (unlike linear problem)
260 bool ProblemFluid::solveTimeStep(){
261 if(_timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
262 return solveNewtonPETSc();
264 return ProblemCoreFlows::solveTimeStep();
267 bool ProblemFluid::solveNewtonPETSc()
269 if( (_nbTimeStep-1)%_freqSave ==0)
270 SNESLineSearchSetDefaultMonitor(_linesearch, _monitorLineSearch);
272 SNESLineSearchSetDefaultMonitor(_linesearch, NULL);
274 SNESSolve(_snes,NULL,_conservativeVars);
276 SNESConvergedReason reason;
277 SNESGetConvergedReason(_snes,&reason);
279 if( (_nbTimeStep-1)%_freqSave ==0)
281 if(reason == SNES_CONVERGED_FNORM_ABS )
282 cout<<"Converged with absolute norm (absolute tolerance) less than "<<_precision_Newton<<", (||F|| < atol)"<<endl;
283 else if(reason == SNES_CONVERGED_FNORM_RELATIVE )
284 cout<<"Converged because residual has decreased by a factor less than "<<_precision_Newton<<", (||F|| < rtol*||F_initial||)"<<endl;
285 else if(reason == SNES_CONVERGED_SNORM_RELATIVE )
286 cout<<"Converged with relative norm (relative tolerance) less than "<<_precision_Newton<<", (|| delta x || < stol || x ||)"<<endl;
287 else if(reason == SNES_CONVERGED_ITS )
288 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;
289 else if(reason == SNES_DIVERGED_LINEAR_SOLVE )
290 cout<<"Solving linear system failed"<<endl;
291 else if(reason == SNES_DIVERGED_LINE_SEARCH )
292 cout<<"Line search failed"<<endl;
293 else if(reason == SNES_DIVERGED_MAX_IT )
294 cout<<"Reached the maximum number of iterations"<<endl;
296 SNESGetIterationNumber(_snes, &_NEWTON_its);
297 PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n\n", _NEWTON_its);
299 *_runLogFile << "Nombre d'iterations de Newton "<< _NEWTON_its <<endl;
305 bool ProblemFluid::iterateTimeStep(bool &converged)
309 if(_NEWTON_its>0){//Pas besoin de computeTimeStep à la première iteration de Newton
311 computeTimeStep(stop);//This compute timestep is just to update the linear system. The time step was imposed before starting the Newton iterations
313 if(stop){//Le compute time step ne s'est pas bien passé
314 cout<<"ComputeTimeStep failed"<<endl;
319 computeNewtonVariation();
321 //converged=convergence des iterations
322 if(_timeScheme == Explicit)
324 else{//Implicit scheme
326 KSPGetIterationNumber(_ksp, &_PetscIts);
327 if( _MaxIterLinearSolver < _PetscIts)//save the maximum number of iterations needed during the newton scheme
328 _MaxIterLinearSolver = _PetscIts;
330 KSPConvergedReason reason;
331 KSPGetConvergedReason(_ksp,&reason);
333 if(reason<0)//solving the linear system failed
335 cout<<"Systeme lineaire : pas de convergence de Petsc. Raison PETSC numéro "<<reason<<endl;
336 cout<<"Nombre d'itérations effectuées "<< _PetscIts<<" nombre maximal Itérations autorisées "<<_maxPetscIts<<endl;
337 *_runLogFile<<"Systeme lineaire : pas de convergence de Petsc. Raison PETSC numéro "<<reason<<endl;
338 *_runLogFile<<"Nombre d'itérations effectuées "<< _PetscIts<<" nombre maximal Itérations autorisées "<<_maxPetscIts<<endl;
342 else{//solving the linear system succeeded
343 //Calcul de la variation relative Uk+1-Uk
347 for(int j=1; j<=_Nmailles; j++)
349 for(int k=0; k<_nVar; k++)
352 VecGetValues(_newtonVariation, 1, &I, &dx);
353 VecGetValues(_conservativeVars, 1, &I, &x);
354 if (fabs(x)*fabs(x)< _precision)
356 if(_erreur_rel < fabs(dx))
357 _erreur_rel = fabs(dx);
359 else if(_erreur_rel < fabs(dx/x))
360 _erreur_rel = fabs(dx/x);
364 converged = _erreur_rel <= _precision_Newton;
367 double relaxation=1;//Uk+1=Uk+relaxation*deltaU
369 VecAXPY(_conservativeVars, relaxation, _newtonVariation);
371 //mise a jour du champ primitif
374 if(_nbPhases==2 && fabs(_err_press_max) > _precision)//la pression n'a pu être calculée en diphasique à partir des variables conservatives
376 cout<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
377 *_runLogFile<<"Warning consToPrim: nbiter max atteint, erreur relative pression= "<<_err_press_max<<" precision= " <<_precision<<endl;
383 cout<<"Vecteur Ukp1-Uk "<<endl;
384 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
385 cout << "Nouvel etat courant Uk de l'iteration Newton: " << endl;
386 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_SELF);
389 if(_nbPhases==2 && (_nbTimeStep-1)%_freqSave ==0){
390 if(_minm1<-_precision || _minm2<-_precision)
392 cout<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
393 *_runLogFile<<"!!!!!!!!! WARNING masse partielle negative sur " << _nbMaillesNeg << " faces, min m1= "<< _minm1 << " , minm2= "<< _minm2<< " precision "<<_precision<<endl;
397 cout << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
398 *_runLogFile << "!!!!!!!!!!!!!!!!!!!!!!!! Complex eigenvalues on " << _nbVpCplx << " cells, max imag= " << _part_imag_max << endl;
405 double ProblemFluid::computeTimeStep(bool & stop){//dt is not known and will not contribute to the Newton scheme
407 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
409 cout << "ProblemFluid::computeTimeStep : Début calcul matrice implicite et second membre"<<endl;
412 if(_restartWithNewTimeScheme)//This is a change of time scheme during a simulation
414 if(_timeScheme == Implicit)
415 MatCreateSeqBAIJ(PETSC_COMM_SELF, _nVar, _nVar*_Nmailles, _nVar*_Nmailles, (1+_neibMaxNbCells), PETSC_NULL, &_A);
418 _restartWithNewTimeScheme=false;
420 if(_timeScheme == Implicit)
423 VecAssemblyBegin(_b);
426 std::vector< int > idCells(2);
427 PetscInt idm, idn, size = 1;
429 long nbFaces = _mesh.getNumberOfFaces();
434 for (int j=0; j<nbFaces;j++){
435 Fj = _mesh.getFace(j);
436 _isBoundary=Fj.isBorder();
437 idCells = Fj.getCellsId();
439 // If Fj is on the boundary
442 for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
444 // compute the normal vector corresponding to face j : from Ctemp1 outward
445 Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
447 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
448 {//we look for l the index of the face Fj for the cell Ctemp1
449 if (j == Ctemp1.getFacesId()[l])
451 for (int idim = 0; idim < _Ndim; ++idim)
452 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
456 }else{ // _Ndim = 1, build normal vector (bug cdmath)
457 if(!_sectionFieldSet)
459 if (Fj.x()<Ctemp1.x())
472 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
474 cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
475 for(int p=0; p<_Ndim; p++)
476 cout << _vec_normal[p] << ",";
479 nameOfGroup = Fj.getGroupName();
480 _porosityi=_porosityField(idCells[k]);
481 _porosityj=_porosityi;
482 setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
483 convectionState(idCells[k],0,true);
484 convectionMatrices();
485 diffusionStateAndMatrices(idCells[k], 0, true);
488 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
490 _inv_dxi = 1/Ctemp1.getMeasure();
492 addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
493 addDiffusionToSecondMember(idCells[k],-1,true);
494 addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
496 if(_timeScheme == Implicit){
497 for(int l=0; l<_nVar*_nVar;l++){
498 _AroeMinusImplicit[l] *= _inv_dxi;
499 _Diffusion[l] *=_inv_dxi*_inv_dxi;
502 jacobian(idCells[k],nameOfGroup,_vec_normal);
503 jacobianDiff(idCells[k],nameOfGroup);
504 if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
505 cout << "Matrice Jacobienne CL convection:" << endl;
506 for(int p=0; p<_nVar; p++){
507 for(int q=0; q<_nVar; q++)
508 cout << _Jcb[p*_nVar+q] << "\t";
512 cout << "Matrice Jacobienne CL diffusion:" << endl;
513 for(int p=0; p<_nVar; p++){
514 for(int q=0; q<_nVar; q++)
515 cout << _JcbDiff[p*_nVar+q] << "\t";
521 //calcul et insertion de A^-*Jcb
522 Polynoms::matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
523 MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
526 displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
529 for(int k=0; k<_nVar*_nVar;k++){
530 _AroeMinusImplicit[k] *= -1;
532 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
534 displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
536 //calcul et insertion de D*JcbDiff
537 Polynoms::matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
538 MatSetValuesBlocked(_A, size, &idm, size, &idm, _a, ADD_VALUES);
539 for(int k=0; k<_nVar*_nVar;k++)
541 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
544 } else if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
545 // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
546 Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
547 Ctemp2 = _mesh.getCell(idCells[1]);
549 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
550 if (j == Ctemp1.getFacesId()[l]){
551 for (int idim = 0; idim < _Ndim; ++idim)
552 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
556 }else{ // _Ndim = 1, build normal vector (bug cdmath)
557 if(!_sectionFieldSet)
559 if (Fj.x()<Ctemp1.x())
566 if(idCells[0]>idCells[1])
572 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
574 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
575 cout<<" Normal vector= ";
576 for (int idim = 0; idim < _Ndim; ++idim)
577 cout<<_vec_normal[idim]<<", ";
580 _porosityi=_porosityField(idCells[0]);
581 _porosityj=_porosityField(idCells[1]);
582 convectionState(idCells[0],idCells[1],false);
583 convectionMatrices();
584 diffusionStateAndMatrices(idCells[0], idCells[1], false);
586 // compute 1/dxi and 1/dxj
589 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
590 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
594 _inv_dxi = 1/Ctemp1.getMeasure();
595 _inv_dxj = 1/Ctemp2.getMeasure();
598 addConvectionToSecondMember(idCells[0],idCells[1], false);
599 addDiffusionToSecondMember( idCells[0],idCells[1], false);
600 addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
602 if(_timeScheme == Implicit){
603 for(int k=0; k<_nVar*_nVar;k++)
605 _AroeMinusImplicit[k] *= _inv_dxi;
606 _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);
610 //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNbCells<<endl;
611 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
612 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
615 displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
616 displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
618 for(int k=0;k<_nVar*_nVar;k++){
619 _AroeMinusImplicit[k] *= -1;
622 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
623 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
625 displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
626 displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
628 for(int k=0; k<_nVar*_nVar;k++)
630 _AroePlusImplicit[k] *= _inv_dxj;
631 _Diffusion[k] *=_inv_dxj/_inv_dxi;
633 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
634 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
636 displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
638 for(int k=0;k<_nVar*_nVar;k++){
639 _AroePlusImplicit[k] *= -1;
642 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
643 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
646 displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
649 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
650 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
651 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
652 *_runLogFile<<"Warning: treatment of a junction node"<<endl;
654 if(!_sectionFieldSet)
656 _runLogFile->close();
657 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
659 int largestSectionCellIndex=0;
660 for(int i=1;i<Fj.getNumberOfCells();i++){
661 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
662 largestSectionCellIndex=i;
664 idm = idCells[largestSectionCellIndex];
665 Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
666 _porosityi=_porosityField(idm);
668 if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
672 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
674 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
675 cout << " ; vecteur normal=(";
676 for(int p=0; p<_Ndim; p++)
677 cout << _vec_normal[p] << ",";
680 for(int i=0;i<Fj.getNumberOfCells();i++){
681 if(i != largestSectionCellIndex){
683 Ctemp2 = _mesh.getCell(idn);
684 _porosityj=_porosityField(idn);
685 convectionState(idm,idn,false);
686 convectionMatrices();
687 diffusionStateAndMatrices(idm, idn,false);
689 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
690 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
692 _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
693 _inv_dxj = 1/Ctemp2.getMeasure();
695 addConvectionToSecondMember(idm,idn, false);
696 _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
697 addDiffusionToSecondMember(idm,idn, false);
698 _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
699 addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
701 if(_timeScheme == Implicit){
702 for(int k=0; k<_nVar*_nVar;k++)
704 _AroeMinusImplicit[k] *= _inv_dxi;
705 _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
707 MatSetValuesBlocked(_A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
708 MatSetValuesBlocked(_A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
711 displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
712 displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
714 for(int k=0;k<_nVar*_nVar;k++){
715 _AroeMinusImplicit[k] *= -1;
718 MatSetValuesBlocked(_A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
719 MatSetValuesBlocked(_A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
721 displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
722 displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
724 for(int k=0; k<_nVar*_nVar;k++)
726 _AroePlusImplicit[k] *= _inv_dxj;
727 _Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
729 MatSetValuesBlocked(_A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
730 MatSetValuesBlocked(_A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
732 displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
734 for(int k=0;k<_nVar*_nVar;k++){
735 _AroePlusImplicit[k] *= -1;
738 MatSetValuesBlocked(_A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
739 MatSetValuesBlocked(_A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
742 displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
749 cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
750 _runLogFile->close();
751 throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
757 if(_timeScheme == Implicit){
758 for(int imaille = 0; imaille<_Nmailles; imaille++)
759 MatSetValuesBlocked(_A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
761 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
762 displayMatrix(_GravityImplicitationMatrix,_nVar,"Gravity matrix:");
764 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
765 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
766 if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
767 cout << "ProblemFluid::computeTimeStep : Fin calcul matrice implicite et second membre"<<endl;
768 cout << "ProblemFluid::computeTimeStep : Matrice implicite :"<<endl;
769 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
770 cout << "ProblemFluid::computeTimeStep : Second membre :"<<endl;
771 VecView(_b, PETSC_VIEWER_STDOUT_WORLD);
779 if(_nbTimeStep+1<_cfl)
780 return (_nbTimeStep+1)*_minl/_maxvp;
783 return _cfl*_minl/_maxvp;
786 void ProblemFluid::computeNewtonVariation()
790 cout<<"Vecteur courant Uk "<<endl;
791 VecView(_conservativeVars,PETSC_VIEWER_STDOUT_SELF);
794 if(_timeScheme == Explicit)
796 VecCopy(_b,_newtonVariation);
797 VecScale(_newtonVariation, _dt);
798 if(_system && (_nbTimeStep-1)%_freqSave ==0)
800 cout<<"Vecteur _newtonVariation =_b*dt"<<endl;
801 VecView(_newtonVariation,PETSC_VIEWER_STDOUT_SELF);
809 cout << "Matrice du système linéaire avant contribution delta t" << endl;
810 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
812 cout << "Second membre du système linéaire avant contribution delta t" << endl;
813 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
816 MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
817 MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
819 VecAXPY(_b, 1/_dt, _old);
820 VecAXPY(_b, -1/_dt, _conservativeVars);
823 #if PETSC_VERSION_GREATER_3_5
824 KSPSetOperators(_ksp, _A, _A);
826 KSPSetOperators(_ksp, _A, _A,SAME_NONZERO_PATTERN);
831 cout << "Matrice du système linéaire après contribution delta t" << endl;
832 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
834 cout << "Second membre du système linéaire après contribution delta t" << endl;
835 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
840 KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
843 KSPSolve(_ksp, _b, _newtonVariation);
847 computeScaling(_maxvp);
849 VecAssemblyBegin(_vecScaling);
850 VecAssemblyBegin(_invVecScaling);
851 for(int imaille = 0; imaille<_Nmailles; imaille++)
854 VecSetValuesBlocked(_vecScaling,1 , &indice, _blockDiag, INSERT_VALUES);
855 VecSetValuesBlocked(_invVecScaling,1,&indice,_invBlockDiag, INSERT_VALUES);
857 VecAssemblyEnd(_vecScaling);
858 VecAssemblyEnd(_invVecScaling);
862 cout << "Matrice avant le preconditionneur des vecteurs propres " << endl;
863 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
865 cout << "Second membre avant le preconditionneur des vecteurs propres " << endl;
866 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
869 MatDiagonalScale(_A,_vecScaling,_invVecScaling);
872 cout << "Matrice apres le preconditionneur des vecteurs propres " << endl;
873 MatView(_A,PETSC_VIEWER_STDOUT_SELF);
876 VecCopy(_b,_bScaling);
877 VecPointwiseMult(_b,_vecScaling,_bScaling);
880 cout << "Produit du second membre par le preconditionneur bloc diagonal a gauche" << endl;
881 VecView(_b, PETSC_VIEWER_STDOUT_SELF);
885 KSPSolve(_ksp,_b, _bScaling);
886 VecPointwiseMult(_newtonVariation,_invVecScaling,_bScaling);
890 cout << "solution du systeme lineaire local:" << endl;
891 VecView(_newtonVariation, PETSC_VIEWER_STDOUT_SELF);
897 void ProblemFluid::computeNewtonRHS( Vec X, Vec F_X){//dt is known and will contribute to the right hand side of the Newton scheme
899 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
901 cout << "ProblemFluid::computeNewtonRHS : Début calcul second membre pour PETSc, _dt="<<_dt<<endl;
903 cout<<"Vecteur courant Uk "<<endl;
904 VecView(X,PETSC_VIEWER_STDOUT_SELF);
908 VecCopy(X,_conservativeVars);
911 VecAssemblyBegin(_b);
914 std::vector< int > idCells(2);
915 PetscInt idm, idn, size = 1;
917 long nbFaces = _mesh.getNumberOfFaces();
922 for (int j=0; j<nbFaces;j++){
923 Fj = _mesh.getFace(j);
924 _isBoundary=Fj.isBorder();
925 idCells = Fj.getCellsId();
927 // If Fj is on the boundary
930 for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
932 // compute the normal vector corresponding to face j : from Ctemp1 outward
933 Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
935 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
936 {//we look for l the index of the face Fj for the cell Ctemp1
937 if (j == Ctemp1.getFacesId()[l])
939 for (int idim = 0; idim < _Ndim; ++idim)
940 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
944 }else{ // _Ndim = 1, build normal vector (bug cdmath)
945 if(!_sectionFieldSet)
947 if (Fj.x()<Ctemp1.x())
960 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
962 cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
963 for(int p=0; p<_Ndim; p++)
964 cout << _vec_normal[p] << ",";
967 nameOfGroup = Fj.getGroupName();
968 _porosityi=_porosityField(idCells[k]);
969 _porosityj=_porosityi;
970 setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
971 convectionState(idCells[k],0,true);
972 convectionMatrices();
973 diffusionStateAndMatrices(idCells[k], 0, true);
976 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
978 _inv_dxi = 1/Ctemp1.getMeasure();
980 addConvectionToSecondMember(idCells[k],-1,true,nameOfGroup);
981 addDiffusionToSecondMember(idCells[k],-1,true);
982 addSourceTermToSecondMember(idCells[k],(_mesh.getCell(idCells[k])).getNumberOfFaces(),-1, -1,true,j,_inv_dxi*Ctemp1.getMeasure());
984 } else if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
985 // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
986 Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
987 Ctemp2 = _mesh.getCell(idCells[1]);
989 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
990 if (j == Ctemp1.getFacesId()[l]){
991 for (int idim = 0; idim < _Ndim; ++idim)
992 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
996 }else{ // _Ndim = 1, build normal vector (bug cdmath)
997 if(!_sectionFieldSet)
999 if (Fj.x()<Ctemp1.x())
1000 _vec_normal[0] = -1;
1006 if(idCells[0]>idCells[1])
1007 _vec_normal[0] = -1;
1012 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1014 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
1015 cout<<" Normal vector= ";
1016 for (int idim = 0; idim < _Ndim; ++idim)
1017 cout<<_vec_normal[idim]<<", ";
1020 _porosityi=_porosityField(idCells[0]);
1021 _porosityj=_porosityField(idCells[1]);
1022 convectionState(idCells[0],idCells[1],false);
1023 convectionMatrices();
1024 diffusionStateAndMatrices(idCells[0], idCells[1], false);
1026 // compute 1/dxi and 1/dxj
1029 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
1030 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
1034 _inv_dxi = 1/Ctemp1.getMeasure();
1035 _inv_dxj = 1/Ctemp2.getMeasure();
1038 addConvectionToSecondMember(idCells[0],idCells[1], false);
1039 addDiffusionToSecondMember( idCells[0],idCells[1], false);
1040 addSourceTermToSecondMember(idCells[0], Ctemp1.getNumberOfFaces(),idCells[1], _mesh.getCell(idCells[1]).getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
1043 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
1044 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1045 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
1046 *_runLogFile<<"Warning: treatment of a junction node"<<endl;
1048 if(!_sectionFieldSet)
1050 _runLogFile->close();
1051 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
1053 int largestSectionCellIndex=0;
1054 for(int i=1;i<Fj.getNumberOfCells();i++){
1055 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
1056 largestSectionCellIndex=i;
1058 idm = idCells[largestSectionCellIndex];
1059 Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
1060 _porosityi=_porosityField(idm);
1062 if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
1065 _vec_normal[0] = -1;
1066 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1068 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
1069 cout << " ; vecteur normal=(";
1070 for(int p=0; p<_Ndim; p++)
1071 cout << _vec_normal[p] << ",";
1072 cout << "). "<<endl;
1074 for(int i=0;i<Fj.getNumberOfCells();i++){
1075 if(i != largestSectionCellIndex){
1077 Ctemp2 = _mesh.getCell(idn);
1078 _porosityj=_porosityField(idn);
1079 convectionState(idm,idn,false);
1080 convectionMatrices();
1081 diffusionStateAndMatrices(idm, idn,false);
1083 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1084 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
1086 _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
1087 _inv_dxj = 1/Ctemp2.getMeasure();
1089 addConvectionToSecondMember(idm,idn, false);
1090 _inv_dxi = sqrt(_sectionField(idn)/_sectionField(idm))/Ctemp1.getMeasure();
1091 addDiffusionToSecondMember(idm,idn, false);
1092 _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
1093 addSourceTermToSecondMember(idm, Ctemp1.getNumberOfFaces()*(Fj.getNumberOfCells()-1),idn, Ctemp2.getNumberOfFaces(),false,j,_inv_dxi*Ctemp1.getMeasure());
1099 cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
1100 _runLogFile->close();
1101 throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
1106 //Contribution from delta t
1107 VecAXPY(_b, 1/_dt, _old);
1108 VecAXPY(_b, -1/_dt, _conservativeVars);
1114 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1116 cout << "ProblemFluid::computeNewtonRHS : Fin calcul second membre pour PETSc"<<endl;
1117 VecView(F_X, PETSC_VIEWER_STDOUT_WORLD);
1122 int ProblemFluid::computeSnesRHS(SNES snes, Vec X, Vec F_X, void *ctx)//Prototype imposé par PETSc
1124 ProblemFluid * myProblem = (ProblemFluid *) ctx;
1125 myProblem->computeNewtonRHS( X, F_X);
1130 void ProblemFluid::computeNewtonJacobian( Vec X, Mat A){//dt is known and will contribute to the jacobian matrix of the Newton scheme
1132 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1134 cout << "ProblemFluid::computeNewtonJacobian : Début calcul Jacobienne schéma Newton pour PETSc, _dt="<<_dt<<endl;
1138 if(_timeScheme == Explicit){
1139 MatCreateConstantDiagonal(PETSC_COMM_SELF, _nVar, _nVar, _nVar*_Nmailles, _nVar*_Nmailles,1./_dt, &A);
1144 VecCopy(X,_conservativeVars);
1146 std::vector< int > idCells(2);
1147 PetscInt idm, idn, size = 1;
1149 long nbFaces = _mesh.getNumberOfFaces();
1154 for (int j=0; j<nbFaces;j++){
1155 Fj = _mesh.getFace(j);
1156 _isBoundary=Fj.isBorder();
1157 idCells = Fj.getCellsId();
1159 // If Fj is on the boundary
1162 for(int k=0;k<Fj.getNumberOfCells();k++)//there will be at most two neighours in the case of an inner wall
1164 // compute the normal vector corresponding to face j : from Ctemp1 outward
1165 Ctemp1 = _mesh.getCell(idCells[k]);//origin of the normal vector
1167 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++)
1168 {//we look for l the index of the face Fj for the cell Ctemp1
1169 if (j == Ctemp1.getFacesId()[l])
1171 for (int idim = 0; idim < _Ndim; ++idim)
1172 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
1176 }else{ // _Ndim = 1, build normal vector (bug cdmath)
1177 if(!_sectionFieldSet)
1179 if (Fj.x()<Ctemp1.x())
1180 _vec_normal[0] = -1;
1187 _vec_normal[0] = -1;
1188 else//idCells[0]==31
1192 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1194 cout << "face numero " << j << " cellule frontiere " << idCells[k] << " ; vecteur normal=(";
1195 for(int p=0; p<_Ndim; p++)
1196 cout << _vec_normal[p] << ",";
1197 cout << "). "<<endl;
1199 nameOfGroup = Fj.getGroupName();
1200 _porosityi=_porosityField(idCells[k]);
1201 _porosityj=_porosityi;
1202 setBoundaryState(nameOfGroup,idCells[k],_vec_normal);
1203 convectionState(idCells[k],0,true);
1204 convectionMatrices();
1205 diffusionStateAndMatrices(idCells[k], 0, true);
1208 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
1210 _inv_dxi = 1/Ctemp1.getMeasure();
1212 for(int l=0; l<_nVar*_nVar;l++){
1213 _AroeMinusImplicit[l] *= _inv_dxi;
1214 _Diffusion[l] *=_inv_dxi*_inv_dxi;
1217 jacobian(idCells[k],nameOfGroup,_vec_normal);
1218 jacobianDiff(idCells[k],nameOfGroup);
1219 if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
1220 cout << "Matrice Jacobienne CL convection:" << endl;
1221 for(int p=0; p<_nVar; p++){
1222 for(int q=0; q<_nVar; q++)
1223 cout << _Jcb[p*_nVar+q] << "\t";
1227 cout << "Matrice Jacobienne CL diffusion:" << endl;
1228 for(int p=0; p<_nVar; p++){
1229 for(int q=0; q<_nVar; q++)
1230 cout << _JcbDiff[p*_nVar+q] << "\t";
1236 //calcul et insertion de A^-*Jcb
1237 Polynoms::matrixProduct(_AroeMinusImplicit, _nVar, _nVar, _Jcb, _nVar, _nVar, _a);
1238 MatSetValuesBlocked(A, size, &idm, size, &idm, _a, ADD_VALUES);
1241 displayMatrix(_a, _nVar, "produit A^-*Jcb pour CL");
1244 for(int k=0; k<_nVar*_nVar;k++){
1245 _AroeMinusImplicit[k] *= -1;
1247 MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
1249 displayMatrix(_AroeMinusImplicit, _nVar,"-_AroeMinusImplicit: ");
1251 //calcul et insertion de D*JcbDiff
1252 Polynoms::matrixProduct(_Diffusion, _nVar, _nVar, _JcbDiff, _nVar, _nVar, _a);
1253 MatSetValuesBlocked(A, size, &idm, size, &idm, _a, ADD_VALUES);
1254 for(int k=0; k<_nVar*_nVar;k++)
1255 _Diffusion[k] *= -1;
1256 MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
1258 } else if (Fj.getNumberOfCells()==2 ){ // Fj is inside the domain and has two neighours (no junction)
1259 // compute the normal vector corresponding to face j : from Ctemp1 to Ctemp2
1260 Ctemp1 = _mesh.getCell(idCells[0]);//origin of the normal vector
1261 Ctemp2 = _mesh.getCell(idCells[1]);
1263 for(int l=0; l<Ctemp1.getNumberOfFaces(); l++){//we look for l the index of the face Fj for the cell Ctemp1
1264 if (j == Ctemp1.getFacesId()[l]){
1265 for (int idim = 0; idim < _Ndim; ++idim)
1266 _vec_normal[idim] = Ctemp1.getNormalVector(l,idim);
1270 }else{ // _Ndim = 1, build normal vector (bug cdmath)
1271 if(!_sectionFieldSet)
1273 if (Fj.x()<Ctemp1.x())
1274 _vec_normal[0] = -1;
1280 if(idCells[0]>idCells[1])
1281 _vec_normal[0] = -1;
1286 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1288 cout << "face numero " << j << " cellule gauche " << idCells[0] << " cellule droite " << idCells[1];
1289 cout<<" Normal vector= ";
1290 for (int idim = 0; idim < _Ndim; ++idim)
1291 cout<<_vec_normal[idim]<<", ";
1294 _porosityi=_porosityField(idCells[0]);
1295 _porosityj=_porosityField(idCells[1]);
1296 convectionState(idCells[0],idCells[1],false);
1297 convectionMatrices();
1298 diffusionStateAndMatrices(idCells[0], idCells[1], false);
1300 // compute 1/dxi and 1/dxj
1303 _inv_dxi = Fj.getMeasure()/Ctemp1.getMeasure();
1304 _inv_dxj = Fj.getMeasure()/Ctemp2.getMeasure();
1308 _inv_dxi = 1/Ctemp1.getMeasure();
1309 _inv_dxj = 1/Ctemp2.getMeasure();
1312 for(int k=0; k<_nVar*_nVar;k++)
1314 _AroeMinusImplicit[k] *= _inv_dxi;
1315 _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);
1319 //cout<<"idm= "<<idm<<"idn= "<<idn<<"nbvoismax= "<<_neibMaxNbCells<<endl;
1320 MatSetValuesBlocked(A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
1321 MatSetValuesBlocked(A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
1324 displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
1325 displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
1327 for(int k=0;k<_nVar*_nVar;k++){
1328 _AroeMinusImplicit[k] *= -1;
1329 _Diffusion[k] *= -1;
1331 MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
1332 MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
1334 displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
1335 displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
1337 for(int k=0; k<_nVar*_nVar;k++)
1339 _AroePlusImplicit[k] *= _inv_dxj;
1340 _Diffusion[k] *=_inv_dxj/_inv_dxi;
1342 MatSetValuesBlocked(A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
1343 MatSetValuesBlocked(A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
1345 displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
1347 for(int k=0;k<_nVar*_nVar;k++){
1348 _AroePlusImplicit[k] *= -1;
1349 _Diffusion[k] *= -1;
1351 MatSetValuesBlocked(A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
1352 MatSetValuesBlocked(A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
1355 displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
1357 else if( Fj.getNumberOfCells()>2 && _Ndim==1 ){//inner face with more than two neighbours
1358 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1359 cout<<"lattice mesh junction at face "<<j<<" nbvoismax= "<<_neibMaxNbCells<<endl;
1360 *_runLogFile<<"Warning: treatment of a junction node"<<endl;
1362 if(!_sectionFieldSet)
1364 _runLogFile->close();
1365 throw CdmathException("ProblemFluid::ComputeTimeStep(): pipe network requires section field");
1367 int largestSectionCellIndex=0;
1368 for(int i=1;i<Fj.getNumberOfCells();i++){
1369 if(_sectionField(idCells[i])>_sectionField(idCells[largestSectionCellIndex]))
1370 largestSectionCellIndex=i;
1372 idm = idCells[largestSectionCellIndex];
1373 Ctemp1 = _mesh.getCell(idm);//origin of the normal vector
1374 _porosityi=_porosityField(idm);
1376 if (j==15)// bug cdmath (Fj.x() > _mesh.getCell(idm).x())
1379 _vec_normal[0] = -1;
1380 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1382 cout<<"Cell with largest section has index "<< largestSectionCellIndex <<" and number "<<idm<<endl;
1383 cout << " ; vecteur normal=(";
1384 for(int p=0; p<_Ndim; p++)
1385 cout << _vec_normal[p] << ",";
1386 cout << "). "<<endl;
1388 for(int i=0;i<Fj.getNumberOfCells();i++){
1389 if(i != largestSectionCellIndex){
1391 Ctemp2 = _mesh.getCell(idn);
1392 _porosityj=_porosityField(idn);
1393 convectionState(idm,idn,false);
1394 convectionMatrices();
1395 diffusionStateAndMatrices(idm, idn,false);
1397 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1398 cout<<"Neighbour index "<<i<<" cell number "<< idn<<endl;
1400 _inv_dxi = _sectionField(idn)/_sectionField(idm)/Ctemp1.getMeasure();
1401 _inv_dxj = 1/Ctemp2.getMeasure();
1403 for(int k=0; k<_nVar*_nVar;k++)
1405 _AroeMinusImplicit[k] *= _inv_dxi;
1406 _Diffusion[k] *=_inv_dxi*2/(1/_inv_dxi+1/_inv_dxj);//use sqrt as above
1408 MatSetValuesBlocked(A, size, &idm, size, &idn, _AroeMinusImplicit, ADD_VALUES);
1409 MatSetValuesBlocked(A, size, &idm, size, &idn, _Diffusion, ADD_VALUES);
1412 displayMatrix(_AroeMinusImplicit, _nVar, "+_AroeMinusImplicit: ");
1413 displayMatrix(_Diffusion, _nVar, "+_Diffusion: ");
1415 for(int k=0;k<_nVar*_nVar;k++){
1416 _AroeMinusImplicit[k] *= -1;
1417 _Diffusion[k] *= -1;
1419 MatSetValuesBlocked(A, size, &idm, size, &idm, _AroeMinusImplicit, ADD_VALUES);
1420 MatSetValuesBlocked(A, size, &idm, size, &idm, _Diffusion, ADD_VALUES);
1422 displayMatrix(_AroeMinusImplicit, _nVar, "-_AroeMinusImplicit: ");
1423 displayMatrix(_Diffusion, _nVar, "-_Diffusion: ");
1425 for(int k=0; k<_nVar*_nVar;k++)
1427 _AroePlusImplicit[k] *= _inv_dxj;
1428 _Diffusion[k] *=_inv_dxj/_inv_dxi;//use sqrt as above
1430 MatSetValuesBlocked(A, size, &idn, size, &idn, _AroePlusImplicit, ADD_VALUES);
1431 MatSetValuesBlocked(A, size, &idn, size, &idn, _Diffusion, ADD_VALUES);
1433 displayMatrix(_AroePlusImplicit, _nVar, "+_AroePlusImplicit: ");
1435 for(int k=0;k<_nVar*_nVar;k++){
1436 _AroePlusImplicit[k] *= -1;
1437 _Diffusion[k] *= -1;
1439 MatSetValuesBlocked(A, size, &idn, size, &idm, _AroePlusImplicit, ADD_VALUES);
1440 MatSetValuesBlocked(A, size, &idn, size, &idm, _Diffusion, ADD_VALUES);
1443 displayMatrix(_AroePlusImplicit, _nVar, "-_AroePlusImplicit: ");
1449 cout<< "Face j="<<j<< " is not a boundary face and has "<<Fj.getNumberOfCells()<< " neighbour cells"<<endl;
1450 _runLogFile->close();
1451 throw CdmathException("ProblemFluid::ComputeTimeStep(): incompatible number of cells around a face");
1456 for(int imaille = 0; imaille<_Nmailles; imaille++)
1457 MatSetValuesBlocked(A, size, &imaille, size, &imaille, _GravityImplicitationMatrix, ADD_VALUES);
1459 MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1460 MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1464 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1466 cout << "ProblemFluid::computeNewtonJacobian : Fin calcul Jacobienne schéma Newton pour PETSc"<<endl;
1467 MatView(A,PETSC_VIEWER_STDOUT_SELF);
1472 int ProblemFluid::computeSnesJacobian(SNES snes, Vec X, Mat A, Mat Aapprox, void *ctx)//Propotype imposé par PETSc
1474 ProblemFluid * myProblem = (ProblemFluid *) ctx;
1475 myProblem->computeNewtonJacobian( X, A);
1482 void ProblemFluid::validateTimeStep()
1486 cout<<" Vecteur Un"<<endl;
1487 VecView(_old, PETSC_VIEWER_STDOUT_WORLD);
1488 cout<<" Vecteur Un+1"<<endl;
1489 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_WORLD);
1491 VecAXPY(_old, -1, _conservativeVars);//old contient old-courant
1493 //Calcul de la variation Un+1-Un
1498 for(int j=1; j<=_Nmailles; j++)
1500 for(int k=0; k<_nVar; k++)
1502 I = (j-1)*_nVar + k;
1503 VecGetValues(_old, 1, &I, &dx);
1504 VecGetValues(_conservativeVars, 1, &I, &x);
1505 if (fabs(x)< _precision)
1507 if(_erreur_rel < fabs(dx))
1508 _erreur_rel = fabs(dx);
1510 else if(_erreur_rel < fabs(dx/x))
1511 _erreur_rel = fabs(dx/x);
1515 _isStationary =_erreur_rel <_precision;
1517 VecCopy(_conservativeVars, _old);
1519 if(_verbose && (_nbTimeStep-1)%_freqSave ==0){
1520 if(!_usePrimitiveVarsInNewton)
1522 cout <<"Valeur propre locale max: " << _maxvp << endl;
1525 if(_nbPhases==2 && (_nbTimeStep-1)%_freqSave ==0){
1526 //Find minimum and maximum void fractions
1527 double alphamin=1e30;
1528 double alphamax=-1e30;
1529 double T, Tmax=-1e30;
1532 for(int j=0; j<_Nmailles; j++)
1534 VecGetValues(_primitiveVars, 1, &I, &x);//extract void fraction
1540 VecGetValues(_primitiveVars, 1, &J, &T);//extract void fraction
1545 cout<<"Alpha min = " << alphamin << ", Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
1547 *_runLogFile<<"Alpha min = " << alphamin << " Alpha max = " << alphamax<<", Tmax= "<<Tmax<<" K"<<endl;
1552 if ((_nbTimeStep-1)%_freqSave ==0 || _isStationary || _time>=_timeMax || _nbTimeStep>=_maxNbOfTimeStep)
1556 void ProblemFluid::abortTimeStep(){
1560 void ProblemFluid::addConvectionToSecondMember
1562 const int &j, bool isBord, string groupname
1565 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1566 cout<<"ProblemFluid::addConvectionToSecondMember start"<<endl;
1568 //extraction des valeurs
1569 for(int k=0; k<_nVar; k++)
1570 _idm[k] = _nVar*i + k;
1571 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1572 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1575 for(int k=0; k<_nVar; k++)
1576 _idn[k] = _nVar*j + k;
1577 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1578 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1581 for(int k=0; k<_nVar; k++)
1583 VecGetValues(_Uext, _nVar, _idn, _Uj);
1584 consToPrim(_Uj, _Vj,_porosityj);
1589 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1591 cout << "addConvectionToSecondMember : état i= " << i << ", _Ui=" << endl;
1592 for(int q=0; q<_nVar; q++)
1593 cout << _Ui[q] << endl;
1595 cout << "addConvectionToSecondMember : état j= " << j << ", _Uj=" << endl;
1596 for(int q=0; q<_nVar; q++)
1597 cout << _Uj[q] << endl;
1600 if(_nonLinearFormulation!=reducedRoe){//VFRoe, Roe or VFFC
1601 Vector Ui(_nVar), Uj(_nVar), Vi(_nVar), Vj(_nVar), Fi(_nVar), Fj(_nVar),Fij(_nVar);
1602 for(int i1=0;i1<_nVar;i1++)
1609 Vector normale(_Ndim);
1610 for(int i1=0;i1<_Ndim;i1++)
1611 normale(i1)=_vec_normal[i1];
1613 Matrix signAroe(_nVar);
1614 for(int i1=0;i1<_nVar;i1++)
1615 for(int i2=0;i2<_nVar;i2++)
1616 signAroe(i1,i2)=_signAroe[i1*_nVar+i2];
1618 Matrix absAroe(_nVar);
1619 for(int i1=0;i1<_nVar;i1++)
1620 for(int i2=0;i2<_nVar;i2++)
1621 absAroe(i1,i2)=_absAroe[i1*_nVar+i2];
1623 if(_nonLinearFormulation==VFRoe)//VFRoe
1625 Vector Uij(_nVar),Vij(_nVar);
1626 double porosityij=sqrt(_porosityi*_porosityj);
1628 Uij=(Ui+Uj)/2+signAroe*(Ui-Uj)/2;
1630 for(int i1=0;i1<_nVar;i1++)
1633 consToPrim(_Uij, _Vij,porosityij);
1635 applyVFRoeLowMachCorrections(isBord, groupname);
1637 for(int i1=0;i1<_nVar;i1++)
1643 Fij=convectionFlux(Uij,Vij,normale,porosityij);
1645 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1647 cout<<"Etat interfacial conservatif "<<i<<", "<<j<< endl;
1649 cout<<"Etat interfacial primitif "<<i<<", "<<j<< endl;
1655 if(_spaceScheme==staggered && _nonLinearFormulation==VFFC)//special case
1656 Fij=staggeredVFFCFlux();
1659 Fi=convectionFlux(Ui,Vi,normale,_porosityi);
1660 Fj=convectionFlux(Uj,Vj,normale,_porosityj);
1661 if(_nonLinearFormulation==VFFC)//VFFC
1662 Fij=(Fi+Fj)/2+signAroe*(Fi-Fj)/2;
1663 else if(_nonLinearFormulation==Roe)//Roe
1664 Fij=(Fi+Fj)/2+absAroe*(Ui-Uj)/2;
1666 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1668 cout<<"Flux cellule "<<i<<" = "<< endl;
1670 cout<<"Flux cellule "<<j<<" = "<< endl;
1675 for(int i1=0;i1<_nVar;i1++)
1676 _phi[i1]=-Fij(i1)*_inv_dxi;
1677 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1678 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1680 cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1681 cout<<"Flux interfacial "<<i<<", "<<j<< endl;
1683 cout << "Contribution convection à " << i << ", -Fij(i1)*_inv_dxi= "<<endl;
1684 for(int q=0; q<_nVar; q++)
1685 cout << _phi[q] << endl;
1689 for(int i1=0;i1<_nVar;i1++)
1690 _phi[i1]*=-_inv_dxj/_inv_dxi;
1691 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1692 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1694 cout << "Contribution convection à " << j << ", Fij(i1)*_inv_dxj= "<<endl;
1695 for(int q=0; q<_nVar; q++)
1696 cout << _phi[q] << endl;
1700 }else //_nonLinearFormulation==reducedRoe)
1702 for(int k=0; k<_nVar; k++)
1703 _temp[k]=(_Ui[k] - _Uj[k])*_inv_dxi;//(Ui-Uj)*_inv_dxi
1704 Polynoms::matrixProdVec(_AroeMinus, _nVar, _nVar, _temp, _phi);//phi=A^-(U_i-U_j)/dx
1705 VecSetValuesBlocked(_b, 1, _idm, _phi, ADD_VALUES);
1707 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1709 cout << "Ajout convection au 2nd membre pour les etats " << i << "," << j << endl;
1710 cout << "(Ui - Uj)*_inv_dxi= "<<endl;;
1711 for(int q=0; q<_nVar; q++)
1712 cout << _temp[q] << endl;
1714 cout << "Contribution convection à " << i << ", A^-*(Ui - Uj)*_inv_dxi= "<<endl;
1715 for(int q=0; q<_nVar; q++)
1716 cout << _phi[q] << endl;
1722 for(int k=0; k<_nVar; k++)
1723 _temp[k]*=_inv_dxj/_inv_dxi;//(Ui-Uj)*_inv_dxj
1724 Polynoms::matrixProdVec(_AroePlus, _nVar, _nVar, _temp, _phi);//phi=A^+(U_i-U_j)/dx
1725 VecSetValuesBlocked(_b, 1, _idn, _phi, ADD_VALUES);
1727 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1729 cout << "Contribution convection à " << j << ", A^+*(Ui - Uj)*_inv_dxi= "<<endl;
1730 for(int q=0; q<_nVar; q++)
1731 cout << _phi[q] << endl;
1736 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1738 cout<<"ProblemFluid::addConvectionToSecondMember end : matrices de décentrement cellules i= " << i << ", et j= " << j<< "):"<<endl;
1739 displayMatrix(_absAroe, _nVar,"Valeur absolue matrice de Roe");
1740 displayMatrix(_AroeMinus, _nVar,"Matrice _AroeMinus");
1741 displayMatrix(_AroePlus, _nVar,"Matrice _AroePlus");
1742 displayMatrix(_signAroe, _nVar,"Signe de la matrice de Roe");
1746 void ProblemFluid::addSourceTermToSecondMember
1747 ( const int i, int nbVoisinsi,
1748 const int j, int nbVoisinsj,
1749 bool isBord, int ij, double mesureFace)//To do : generalise to unstructured meshes
1751 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1752 cout<<"ProblemFluid::addSourceTerm cell i= "<<i<< " cell j= "<< j<< " isbord "<<isBord<<endl;
1755 for(int k=1; k<_nVar;k++)
1756 _idm[k] = _idm[k-1] + 1;
1757 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1758 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1759 sourceVector(_Si,_Ui,_Vi,i);
1761 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1763 cout << "Terme source S(Ui), i= " << i<<endl;
1764 for(int q=0; q<_nVar; q++)
1765 cout << _Si[q] << endl;
1769 for(int k=0; k<_nVar; k++)
1770 _idn[k] = _nVar*j + k;
1771 VecGetValues(_conservativeVars, _nVar, _idn, _Uj);
1772 VecGetValues(_primitiveVars, _nVar, _idn, _Vj);
1773 sourceVector(_Sj,_Uj,_Vj,j);
1776 for(int k=0; k<_nVar; k++)
1778 VecGetValues(_Uext, _nVar, _idn, _Uj);
1779 consToPrim(_Uj, _Vj,_porosityj);
1780 sourceVector(_Sj,_Uj,_Vj,i);
1782 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1785 cout << "Terme source S(Uj), j= " << j<<endl;
1787 cout << "Terme source S(Uj), cellule fantôme "<<endl;
1788 for(int q=0; q<_nVar; q++)
1789 cout << _Sj[q] << endl;
1793 //Compute pressure loss vector
1795 if(_pressureLossFieldSet){
1796 K=_pressureLossField(ij);
1797 pressureLossVector(_pressureLossVector, K,_Ui, _Vi, _Uj, _Vj);
1801 for(int k=0; k<_nVar;k++)
1802 _pressureLossVector[k]=0;
1804 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1806 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", K="<<K<<endl;
1807 for(int k=0; k<_nVar;k++)
1808 cout<< _pressureLossVector[k]<<", ";
1811 //Contribution of the porosityField gradient:
1813 porosityGradientSourceVector();
1815 for(int k=0; k<_nVar;k++)
1816 _porosityGradientSourceVector[k]=0;
1819 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1822 cout<<"interface i= "<<i<<", j= "<<j<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<", dxj= "<<1/_inv_dxj<<endl;
1824 cout<<"interface frontière i= "<<i<<", ij="<<ij<<", dxi= "<<1/_inv_dxi<<endl;
1825 cout<<"Gradient de porosite à l'interface"<<endl;
1826 for(int k=0; k<_nVar;k++)
1827 cout<< _porosityGradientSourceVector[k]<<", ";
1832 if(_wellBalancedCorrection){
1833 for(int k=0; k<_nVar;k++)
1834 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1835 Polynoms::matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1836 for(int k=0; k<_nVar;k++){
1837 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1838 _Sj[k]=(_phi[k]+_l[k])*mesureFace/_perimeters(j);///nbVoisinsj;
1840 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1842 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant (après décentrement) de la face (i,j), j="<<j<<endl;
1843 for(int q=0; q<_nVar; q++)
1844 cout << _Si[q] << endl;
1845 cout << "Contribution au terme source Sj de la cellule j= " << j<<" venant (après décentrement) de la face (i,j), i="<<i<<endl;
1846 for(int q=0; q<_nVar; q++)
1847 cout << _Sj[q] << endl;
1849 cout<<"ratio surface sur volume i = "<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1850 cout<<"ratio surface sur volume j = "<<mesureFace/_perimeters(j)<<" perimeter = "<< _perimeters(j)<<endl;
1855 for(int k=0; k<_nVar;k++){
1856 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i)
1857 _Sj[k]=_Sj[k]/nbVoisinsj+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(j)
1859 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1861 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,j), j="<<j<<endl;
1862 for(int q=0; q<_nVar; q++)
1863 cout << _Si[q] << endl;
1864 cout << "Contribution au terme source Sj de la cellule j = " << j<<" venant de la face (i,j), i="<<i <<endl;
1865 for(int q=0; q<_nVar; q++)
1866 cout << _Sj[q] << endl;
1871 VecSetValuesBlocked(_b, 1, _idn, _Sj, ADD_VALUES);
1873 if(_wellBalancedCorrection){
1874 for(int k=0; k<_nVar;k++)
1875 _phi[k]=(_Si[k]+_Sj[k])/2+_pressureLossVector[k]+_porosityGradientSourceVector[k];
1876 Polynoms::matrixProdVec(_signAroe, _nVar, _nVar, _phi, _l);
1877 for(int k=0; k<_nVar;k++)
1878 _Si[k]=(_phi[k]-_l[k])*mesureFace/_perimeters(i);///nbVoisinsi;
1879 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1881 cout << "Contribution au terme source Si de la cellule i= " << i<<" venant (après décentrement) de la face (i,bord)"<<endl;
1882 for(int q=0; q<_nVar; q++)
1883 cout << _Si[q] << endl;
1884 cout<<"ratio surface sur volume i ="<<mesureFace/_perimeters(i)<<" perimeter = "<< _perimeters(i)<<endl;
1890 for(int k=0; k<_nVar;k++)
1891 _Si[k]=_Si[k]/nbVoisinsi+_pressureLossVector[k]/2+_porosityGradientSourceVector[k]/2;//mesureFace/_perimeters(i);//
1892 if (_verbose && (_nbTimeStep-1)%_freqSave ==0)
1894 cout << "Contribution au terme source Si de la cellule i = " << i<<" venant de la face (i,bord) "<<endl;
1895 for(int q=0; q<_nVar; q++)
1896 cout << _Si[q] << endl;
1902 VecSetValuesBlocked(_b, 1, _idm, _Si, ADD_VALUES);
1904 if(_verbose && (_nbTimeStep-1)%_freqSave ==0 && _wellBalancedCorrection)
1905 displayMatrix( _signAroe,_nVar,"Signe matrice de Roe");
1908 void ProblemFluid::updatePrimitives()
1910 VecAssemblyBegin(_primitiveVars);
1911 for(int i=1; i<=_Nmailles; i++)
1913 _idm[0] = _nVar*( (i-1));
1914 for(int k=1; k<_nVar; k++)
1915 _idm[k] = _idm[k-1] + 1;
1917 VecGetValues(_conservativeVars, _nVar, _idm, _Ui);
1918 consToPrim(_Ui,_Vi,_porosityField(i-1));
1919 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1921 cout << "ProblemFluid::updatePrimitives() cell " << i-1 << endl;
1923 for(int q=0; q<_nVar; q++)
1924 cout << _Ui[q] << "\t";
1927 for(int q=0; q<_nVar; q++)
1928 cout << _Vi[q] << "\t";
1932 if(_nbPhases==2 && _Psat>-1e30){//Cas simulation flashing
1934 VecGetValues(_primitiveVars, 1, _idm+1, &pressure);
1935 _dp_over_dt(i-1)=(_Vi[1]-pressure)/_dt;//pn+1-pn
1938 VecSetValuesBlocked(_primitiveVars, 1, _idm, _Vi, INSERT_VALUES);
1940 VecAssemblyEnd(_primitiveVars);
1944 cout << "Nouvelles variables primitives : " << endl;
1945 VecView(_primitiveVars, PETSC_VIEWER_STDOUT_WORLD);
1950 void ProblemFluid::updateConservatives()
1952 VecAssemblyBegin(_conservativeVars);
1953 for(int i=1; i<=_Nmailles; i++)
1955 _idm[0] = _nVar*( (i-1));
1956 for(int k=1; k<_nVar; k++)
1957 _idm[k] = _idm[k-1] + 1;
1959 VecGetValues(_primitiveVars, _nVar, _idm, _Vi);
1961 primToCons(_Vi,0,_Ui,0);
1963 VecSetValuesBlocked(_conservativeVars, 1, _idm, _Ui, INSERT_VALUES);
1965 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
1967 cout << "ProblemFluid::updateConservatives() cell " << i-1 << endl;
1969 for(int q=0; q<_nVar; q++)
1970 cout << _Vi[q] << "\t";
1973 for(int q=0; q<_nVar; q++)
1974 cout << _Ui[q] << "\t";
1978 VecAssemblyEnd(_conservativeVars);
1982 cout << "Nouvelles variables conservatives : " << endl;
1983 VecView(_conservativeVars, PETSC_VIEWER_STDOUT_WORLD);
1988 vector< complex<double> > ProblemFluid::getRacines(vector< double > pol_car){
1989 int degre_polynom = pol_car.size()-1;
1991 vector< complex< double > > valeurs_propres (degre_polynom);
1992 vector< double > zeror(degre_polynom);
1993 vector< double > zeroi(degre_polynom);
1994 for(int j=0; j<(degre_polynom+1)/2; j++){//coefficients in order of decreasing powers for rpoly
1996 pol_car[j] =pol_car[degre_polynom-j];
1997 pol_car[degre_polynom-j]=tmp;
2000 //Calcul des racines du polynome
2001 roots_polynoms roots;
2002 int size_vp = roots.rpoly(&pol_car[0],degre_polynom,&zeror[0],&zeroi[0]);
2004 //On ordonne les valeurs propres
2005 if(zeror[1]<zeror[0])
2015 if(size_vp == degre_polynom) {
2016 for(int ct =0; ct<degre_polynom; ct++)
2017 valeurs_propres[ct] = complex<double> (zeror[ct],zeroi[ct]); //vp non triviales
2020 cout << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
2021 *_runLogFile << " Pb, found only " << size_vp << " eigenvalues over " << degre_polynom << endl;
2022 cout<<"Coefficients polynome caracteristique: "<<endl;
2023 for(int ct =0; ct<degre_polynom+1; ct++)
2024 cout<<pol_car[ct]<< " , " <<endl;
2026 *_runLogFile<<"getRacines computation failed"<<endl;
2027 _runLogFile->close();
2028 throw CdmathException("getRacines computation failed");
2031 return valeurs_propres;
2034 void ProblemFluid::AbsMatriceRoe(vector< complex<double> > valeurs_propres_dist)
2036 int nbVp_dist=valeurs_propres_dist.size();
2037 vector< complex< double > > y (nbVp_dist,0);
2038 for( int i=0 ; i<nbVp_dist ; i++)
2039 y[i] = Polynoms::abs_generalise(valeurs_propres_dist[i]);
2040 Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _absAroe,y);
2041 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2043 cout<< endl<<"ProblemFluid::AbsMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
2044 for(int ct =0; ct<nbVp_dist; ct++)
2045 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
2047 cout<<"ProblemFluid::AbsMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
2048 for(int ct =0; ct<nbVp_dist; ct++)
2049 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
2051 displayMatrix(_absAroe,_nVar,"Valeur absolue de la matrice de Roe");
2055 void ProblemFluid::SigneMatriceRoe(vector< complex<double> > valeurs_propres_dist)
2057 int nbVp_dist=valeurs_propres_dist.size();
2058 vector< complex< double > > y (nbVp_dist,0);
2059 for( int i=0 ; i<nbVp_dist ; i++)
2061 if(valeurs_propres_dist[i].real()>0)
2063 else if(valeurs_propres_dist[i].real()<0)
2069 Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _signAroe,y);
2070 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2072 cout<< endl<<"ProblemFluid::SigneMatriceRoe: Valeurs propres :" << nbVp_dist<<endl;
2073 for(int ct =0; ct<nbVp_dist; ct++)
2074 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
2076 cout<<"ProblemFluid::SigneMatriceRoe: Valeurs à atteindre pour l'interpolation"<<endl;
2077 for(int ct =0; ct<nbVp_dist; ct++)
2078 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
2080 displayMatrix(_signAroe,_nVar,"Signe matrice de Roe");
2083 void ProblemFluid::InvMatriceRoe(vector< complex<double> > valeurs_propres_dist)
2085 int nbVp_dist=valeurs_propres_dist.size();
2086 vector< complex< double > > y (nbVp_dist,0);
2088 for( int i=0 ; i<nbVp_dist ; i++)
2090 if(Polynoms::module(valeurs_propres_dist[i])>_precision)
2091 y[i] = 1./valeurs_propres_dist[i];
2093 y[i] = 1./_precision;
2095 Polynoms::abs_par_interp_directe(nbVp_dist,valeurs_propres_dist, _Aroe, _nVar,_precision, _invAroe,y);
2096 if(_verbose && (_nbTimeStep-1)%_freqSave ==0)
2098 cout<< endl<<"ProblemFluid::InvMatriceRoe : Valeurs propres :" << nbVp_dist<<endl;
2099 for(int ct =0; ct<nbVp_dist; ct++)
2100 cout<< "("<<valeurs_propres_dist[ct].real()<< ", " <<valeurs_propres_dist[ct].imag() <<") ";
2102 cout<<"ProblemFluid::InvMatriceRoe : Valeurs à atteindre pour l'interpolation"<<endl;
2103 for(int ct =0; ct<nbVp_dist; ct++)
2104 cout<< "("<<y[ct].real()<< ", " <<y[ct].imag() <<") ";
2106 displayMatrix(_invAroe,_nVar,"Inverse matrice de Roe");
2110 Field ProblemFluid::getConservativeField() const
2112 if(!_initializedMemory)
2114 _runLogFile->close();
2115 throw CdmathException("ProblemFluid::getConservativeField call initialize() first");
2120 Field ProblemFluid::getPrimitiveField() const
2122 if(!_initializedMemory)
2124 _runLogFile->close();
2125 throw CdmathException("ProblemFluid::getPrimitiveField call initialize() first");
2130 void ProblemFluid::terminate(){
2133 delete[] _Diffusion;
2134 delete[] _GravityImplicitationMatrix;
2135 delete[] _AroeMinus;
2140 delete[] _AroeImplicit;
2141 delete[] _AroeMinusImplicit;
2142 delete[] _AroePlusImplicit;
2143 delete[] _absAroeImplicit;
2148 delete[] _primToConsJacoMat;
2155 delete[] _externalStates;
2158 delete[] _vec_normal;
2163 if(_nonLinearFormulation==VFRoe){
2169 delete[] _pressureLossVector;
2170 delete[] _porosityGradientSourceVector;
2173 delete[] _blockDiag;
2174 delete[] _invBlockDiag;
2176 VecDestroy(&_vecScaling);
2177 VecDestroy(&_invVecScaling);
2178 VecDestroy(&_bScaling);
2181 VecDestroy(&_conservativeVars);
2182 VecDestroy(&_newtonVariation);
2184 VecDestroy(&_primitiveVars);
2187 VecDestroy(&_Uextdiff);
2191 for(int i=0;i<_nbPhases;i++)
2194 // Destruction du solveur de Newton de PETSc
2195 if(_timeScheme == Implicit && _nonLinearSolver != Newton_SOLVERLAB)
2196 SNESDestroy(&_snes);
2200 ProblemFluid::getInputFieldsNames()
2202 vector<string> result(1);
2204 result[0]="HeatPower";
2205 result[1]="Porosity";
2206 result[2]="PressureLoss";
2207 result[3]="Section";
2213 ProblemFluid::setInputField(const string& nameField, Field& inputField )
2215 if(nameField=="Porosity" || nameField=="POROSITY" || nameField=="Porosité" || nameField=="POROSITE")
2216 return setPorosityField( inputField) ;
2217 else if(nameField=="HeatPower" || nameField=="HEATPOWER" || nameField=="PuissanceThermique" || nameField=="PUISSANCETHERMIQUE" )
2218 return setHeatPowerField( inputField );
2219 else if(nameField=="PressureLoss" || nameField=="PRESSURELOSS" || nameField=="PerteDeCharge" || nameField=="PPERTEDECHARGE")
2220 return setPressureLossField( inputField) ;
2221 else if(nameField=="Section" || nameField=="SECTION" )
2222 return setSectionField( inputField );
2225 cout<<"Error : Field name "<< nameField << " is not an input field name, call getInputFieldsNames first" << endl;
2226 throw CdmathException("SinglePhase::setInputField error : Unknown Field name");
2231 ProblemFluid::setPorosityField(string fileName, string fieldName){
2232 if(!_initialDataSet)
2233 throw CdmathException("!!!!!!!! ProblemFluid::setPorosityField set initial field first");
2235 _porosityField=Field(fileName, CELLS,fieldName);
2236 _porosityField.getMesh().checkFastEquivalWith(_mesh);
2237 _porosityFieldSet=true;
2241 ProblemFluid::setPressureLossField(Field PressureLoss){
2242 if(!_initialDataSet)
2243 throw CdmathException("!!!!!!!! ProblemFluid::setPressureLossField set initial field first");
2245 PressureLoss.getMesh().checkFastEquivalWith(_mesh);
2246 _pressureLossField=PressureLoss;
2247 _pressureLossFieldSet=true;
2251 ProblemFluid::setPressureLossField(string fileName, string fieldName){
2252 if(!_initialDataSet)
2253 throw CdmathException("!!!!!!!! ProblemFluid::setPressureLossField set initial field first");
2255 _pressureLossField=Field(fileName, FACES,fieldName);
2256 _pressureLossField.getMesh().checkFastEquivalWith(_mesh);
2257 _pressureLossFieldSet=true;
2261 ProblemFluid::setSectionField(Field sectionField){
2262 if(!_initialDataSet)
2263 throw CdmathException("!!!!!!!! ProblemFluid::setSectionField set initial field first");
2265 sectionField.getMesh().checkFastEquivalWith(_mesh);
2266 _sectionField=sectionField;
2267 _sectionFieldSet=true;
2271 ProblemFluid::setSectionField(string fileName, string fieldName){
2272 if(!_initialDataSet)
2273 throw CdmathException("!!!!!!!! ProblemFluid::setSectionField set initial field first");
2275 _sectionField=Field(fileName, CELLS,fieldName);
2276 _sectionField.getMesh().checkFastEquivalWith(_mesh);
2277 _sectionFieldSet=true;
2281 ProblemFluid::setPorosityField(Field Porosity){
2282 if(!_initialDataSet)
2283 throw CdmathException("!!!!!!!! ProblemFluid::setPorosityField set initial field first");
2285 Porosity.getMesh().checkFastEquivalWith(_mesh);
2286 _porosityField=Porosity;
2287 _porosityFieldSet=true;