]> SALOME platform Git repositories - tools/solverlab.git/commitdiff
Salome HOME
First MPI tests
authormichael <michael@localhost.localdomain>
Mon, 14 Dec 2020 16:12:43 +0000 (17:12 +0100)
committermichael <michael@localhost.localdomain>
Mon, 14 Dec 2020 16:12:43 +0000 (17:12 +0100)
CoreFlows/examples/C/CMakeLists.txt
CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx [new file with mode: 0755]
CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_MPI.cxx [new file with mode: 0755]

index d9a5aa4cf06594631a4db9f4be9283cac3ea2960..0abbd28adbec0b893d6371c27c5c4b1094045658 100755 (executable)
@@ -7,9 +7,6 @@ INCLUDE_DIRECTORIES(
 )
 
 
-SET(_extra_lib_CoreFlows CoreFlowsLibs )
-
-
 if(CMAKE_COMPILER_IS_GNUCXX)
     if (CMAKE_BUILD_TYPE STREQUAL Debug)
     include(CodeCoverage)
@@ -44,7 +41,7 @@ function(CreateTestExecAndInstall SourceTestFile libList)
 endfunction(CreateTestExecAndInstall)
 
 
-set( libs_for_tests ${_extra_lib_CoreFlows} )
+set( libs_for_tests CoreFlowsLibs )
 
 # copy tests resources (med files etc.) into the build directory
 file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/../resources DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
@@ -99,4 +96,14 @@ CreateTestExecAndInstall(StationaryDiffusionEquation_2DFV_StructuredSquares.cxx
 CreateTestExecAndInstall(StationaryDiffusionEquation_3DEF_StructuredTetrahedra.cxx  "${libs_for_tests}" )
 CreateTestExecAndInstall(StationaryDiffusionEquation_3DFV_StructuredTetrahedra.cxx  "${libs_for_tests}" )
 CreateTestExecAndInstall(testEOS.cxx  "${libs_for_tests}" )
+CreateTestExecAndInstall(WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx  "${libs_for_tests}" )
+
+if( SOLVERLAB_WITH_MPI )
+  add_executable(WaveSystem_2DFV_SphericalExplosion_MPI.exe WaveSystem_2DFV_SphericalExplosion_MPI.cxx)                    # compilation of the testxxx.exe 
+  target_link_libraries(WaveSystem_2DFV_SphericalExplosion_MPI.exe CoreFlowsLibs  ${MPI_LIBRARY})              # provide required lib for testxxx.exe 
+  install(TARGETS WaveSystem_2DFV_SphericalExplosion_MPI.exe DESTINATION share/examples)
+  add_test(NAME WaveSystem_2DFV_SphericalExplosion_MPI_SEQ COMMAND "${MPIEXEC}" "-n" "1" "./WaveSystem_2DFV_SphericalExplosion_MPI.exe")     # adding a ctest Test
+  add_test(NAME WaveSystem_2DFV_SphericalExplosion_MPI_PAR COMMAND "${MPIEXEC}" "-n" "2" "./WaveSystem_2DFV_SphericalExplosion_MPI.exe")     # adding a ctest Test
+endif( SOLVERLAB_WITH_MPI )
 
diff --git a/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx b/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx
new file mode 100755 (executable)
index 0000000..24ad784
--- /dev/null
@@ -0,0 +1,277 @@
+//============================================================================
+// Author      : Michael NDJINGA
+// Date        : November 2020
+// Description : 2D linear wave system
+//============================================================================
+
+#include <iostream>
+#include <string>
+#include <cmath>
+
+#include "Mesh.hxx"
+#include "Cell.hxx"
+#include "Face.hxx"
+#include "Field.hxx"
+#include "SparseMatrixPetsc.hxx"
+#include "CdmathException.hxx"
+
+using namespace std;
+
+double p0  =155e5;   //reference pressure in a pressurised nuclear vessel
+double c0  =700.;    //reference sound speed for water at 155 bars, 600K
+double rho0=p0/c0*c0;//reference density
+double precision=1e-5;
+
+void initial_conditions_shock(Mesh my_mesh,Field& pressure_field,Field& velocity_field)
+{
+    double rayon=0.35;
+    double xcentre=0.;
+    double ycentre=0;
+
+    int dim    =my_mesh.getMeshDimension();
+    int nbCells=my_mesh.getNumberOfCells();
+
+    for (int j=0 ; j<nbCells ; j++)
+    {
+        double x = my_mesh.getCell(j).x() ;
+        double y = my_mesh.getCell(j).y() ;
+        double valX=(x-xcentre)*(x-xcentre);
+        double valY=(y-ycentre)*(y-ycentre);
+        double val=sqrt(valX+valY);
+
+               for(int idim=0; idim<dim; idim++)
+                       velocity_field[j,idim]=0;
+                       
+        if (val<rayon)
+            pressure_field(j) = 155e5;
+        else
+            pressure_field(j) = 70e5;
+    }
+}
+
+Matrix jacobianMatrices(Vector normal, double coeff)
+{
+    int dim=normal.size();
+    Matrix    A(dim+1,dim+1);
+    Matrix absA(dim+1,dim+1);
+
+    absA(0,0)=c0*coeff;
+    for(int i=0 ; i<dim; i++)
+    {
+        A(i+1,0)=      normal[i]*coeff;
+        A(0,i+1)=c0*c0*normal[i]*coeff;
+        for( int j=0 ; j<dim; j++)
+            absA(i+1,j+1)=c0*normal[i]*normal[j]*coeff;
+    }
+    return (A - absA)*(1./2);
+}
+    
+SparseMatrixPetsc computeDivergenceMatrix(Mesh my_mesh, int nbVoisinsMax, double dt)
+{
+    int nbCells = my_mesh.getNumberOfCells();
+    int dim=my_mesh.getMeshDimension();
+    int nbComp=dim+1;
+    Vector normal(dim);
+
+    SparseMatrixPetsc implMat(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp);
+
+    Matrix idMoinsJacCL(nbComp);
+    
+    for(int j=0; j<nbCells; j++)//On parcourt les cellules
+    {
+        Cell Cj = my_mesh.getCell(j);
+        int nbFaces = Cj.getNumberOfFaces();
+
+        for(int k=0; k<nbFaces; k++)
+        {
+            int indexFace = Cj.getFacesId()[k];
+            Face Fk = my_mesh.getFace(indexFace);
+            for(int i =0; i<dim ; i++)
+                normal[i] = Cj.getNormalVector(k, i);//normale sortante
+
+            Matrix Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
+
+            int cellAutre =-1;
+            if ( not Fk.isBorder())
+            {
+                /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
+                if (Fk.getCellsId()[0] == j) 
+                    // hypothese verifiée 
+                    cellAutre = Fk.getCellsId()[1];
+                else if(Fk.getCellsId()[1] == j) 
+                    // hypothese non verifiée 
+                    cellAutre = Fk.getCellsId()[0];
+                else
+                    throw CdmathException("computeDivergenceMatrix: problem with mesh, unknown cell number");
+                    
+                implMat.addValue(j*nbComp,cellAutre*nbComp,Am);
+                implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.));
+            }
+            else 
+            {    
+                               if( Fk.getGroupName() != "Periodic" && Fk.getGroupName() != "Neumann")//Wall boundary condition unless Periodic/Neumann specified explicitly
+                {
+                                       Vector v(dim+1);
+                    for(int i=0; i<dim; i++)
+                        v[i+1]=normal[i];
+                    Matrix idMoinsJacCL=v.tensProduct(v)*2;
+                    
+                    implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL);
+                               }
+                else if( Fk.getGroupName() == "Periodic")//Periodic boundary condition
+                {
+                    int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
+                    Face Fp = my_mesh.getFace(indexFP);
+                    cellAutre = Fp.getCellsId()[0];
+                    
+                    implMat.addValue(j*nbComp,cellAutre*nbComp,Am);
+                    implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.));
+                               }
+                else if(Fk.getGroupName() != "Neumann")//Nothing to do for Neumann boundary condition
+                {
+                    cout<< Fk.getGroupName() <<endl;
+                    throw CdmathException("computeDivergenceMatrix: Unknown boundary condition name");
+                               }
+            }
+        }   
+    }     
+    return implMat;
+}
+
+void WaveSystem2D(double tmax, int ntmax, double cfl, int output_freq, const Mesh& my_mesh, const string file)
+{
+       /* Retrieve mesh data */
+    int dim=my_mesh.getMeshDimension();
+    int nbCells = my_mesh.getNumberOfCells();
+    std::string meshName=my_mesh.getName();
+    int nbVoisinsMax=my_mesh.getMaxNbNeighbours(CELLS);
+    double dx_min=my_mesh.minRatioVolSurf();
+
+
+    /* Initial conditions */
+    cout << "Construction of the initial condition" << endl;
+    
+    Field pressure_field("Pressure",CELLS,my_mesh,1) ;
+    Field velocity_field("Velocity",CELLS,my_mesh,1) ;
+    initial_conditions_shock(my_mesh,pressure_field, velocity_field);
+
+    /* iteration vectors */
+    Vector  Un(nbCells*(dim+1));
+    Vector dUn(nbCells*(dim+1));
+    
+    for(int k =0; k<nbCells; k++){
+        Un[k*(dim+1)+0] =     pressure_field[k];
+        Un[k*(dim+1)+1] =rho0*velocity_field[k,0];
+        if(dim>=2)
+            Un[k + 2*nbCells] = rho0*velocity_field[k,1] ;
+            if(dim==3)
+                Un[k + 3*nbCells] = rho0*velocity_field[k,2];
+                       }
+
+    /*
+     * MED output of the initial condition at t=0 and iter = 0
+     */
+    int it=0;
+    bool isStationary=false;
+    double time=0.;
+    double dt = cfl * dx_min / c0;
+    
+    cout << "Saving the solution at T=" << time << endl;
+    pressure_field.setTime(time,it);
+    pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure");
+    velocity_field.setTime(time,it);
+    velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity");
+    /* --------------------------------------------- */
+
+    SparseMatrixPetsc divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt);
+
+    /* Time loop */
+    cout<< "Starting computation of the linear wave system with an explicit UPWIND scheme" << endl;
+    while (it<ntmax && time <= tmax && ! isStationary)
+    {
+        dUn=divMat*Un;
+        Un-=dUn;
+        
+        time=time+dt;
+        it=it+1;
+         isStationary = dUn.norm()<precision;
+        /* Sauvegardes */
+        if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax)
+        {
+            cout<<"-- Iteration: " << it << ", Time: " << time << ", dt: " << dt<<endl;
+
+            for(int k=0; k<nbCells; k++)
+            {
+                pressure_field[k]  =Un[k*(dim+1)+0];
+                velocity_field[k,0]=Un[k*(dim+1)+1]/rho0;
+                if(dim>1)
+                {
+                    velocity_field[k,1]=Un[k*(dim+1)+2]/rho0;
+                    if(dim>2)
+                        velocity_field[k,2]=Un[k*(dim+1)+3]/rho0;
+                               }
+                       }
+            pressure_field.setTime(time,it);
+            pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure",false);
+            velocity_field.setTime(time,it);
+            velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity",false);
+               }
+    }
+    cout<<"End of calculation -- Iteration: " << it << ", Time: "<< time<< ", dt: " << dt<<endl;
+
+    if(it>=ntmax)
+        cout<< "Nombre de pas de temps maximum ntmax= "<< ntmax<< " atteint"<<endl;
+    else if(isStationary)
+        cout<< "Régime stationnaire atteint au pas de temps "<< it<< ", t= "<< time<<endl;       
+    else
+        cout<< "Temps maximum Tmax= "<< tmax<< " atteint"<<endl;
+}
+int main(int argc, char *argv[])
+{
+           cout << "-- Starting the RESOLUTION OF THE 2D WAVE SYSTEM"<<endl;
+           cout << "- Numerical scheme : Upwind explicit scheme" << endl;
+           cout << "- Boundary conditions : WALL" << endl;
+
+    // Problem data
+    double cfl=0.49;
+    double tmax=1.;
+    double ntmax=3;//20000;
+    int freqSortie=100;
+    string fileOutPut="SphericalWave";
+
+       if(argc<2)
+       {
+                   cout << "- DOMAIN : SQUARE" << endl;
+                   cout << "- MESH : CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << endl<< endl;
+                   cout << "Construction of a cartesian mesh" << endl;
+           double xinf=0.0;
+           double xsup=1.0;
+           double yinf=0.0;
+           double ysup=1.0;
+           int nx=50;
+           int ny=50;
+           Mesh myMesh(xinf,xsup,nx,yinf,ysup,ny);
+           double eps=1.E-10;
+           myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
+           myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
+           myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
+           myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
+
+               WaveSystem2D(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut);
+       }
+       else
+       {
+           cout << "- MESH:  GENERATED EXTERNALLY WITH SALOME" << endl;
+           cout << "Loading of a mesh named "<<argv[1] << endl;
+           string filename = argv[1];
+           Mesh myMesh(filename);
+
+               WaveSystem2D(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut);
+       }
+       
+    cout << "Simulation complete." << endl;
+
+    return 0;
+}
diff --git a/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_MPI.cxx b/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_MPI.cxx
new file mode 100755 (executable)
index 0000000..0d1cc35
--- /dev/null
@@ -0,0 +1,374 @@
+//============================================================================
+// Author      : Michael NDJINGA
+// Date        : November 2020
+// Description : 2D linear wave system
+//============================================================================
+
+#include <iostream>
+#include <string>
+#include <cmath>
+
+#include "Mesh.hxx"
+#include "Cell.hxx"
+#include "Face.hxx"
+#include "Field.hxx"
+#include "CdmathException.hxx"
+
+#include <petscksp.h>
+
+using namespace std;
+
+double p0  =155e5;   //reference pressure in a pressurised nuclear vessel
+double c0  =700.;    //reference sound speed for water at 155 bars, 600K
+double rho0=p0/c0*c0;//reference density
+double precision=1e-5;
+
+void initial_conditions_shock(Mesh my_mesh,Field& pressure_field,Field& velocity_field)
+{
+    double rayon=0.35;
+    double xcentre=0.;
+    double ycentre=0;
+
+    int dim    =my_mesh.getMeshDimension();
+    int nbCells=my_mesh.getNumberOfCells();
+
+    for (int j=0 ; j<nbCells ; j++)
+    {
+        double x = my_mesh.getCell(j).x() ;
+        double y = my_mesh.getCell(j).y() ;
+        double valX=(x-xcentre)*(x-xcentre);
+        double valY=(y-ycentre)*(y-ycentre);
+        double val=sqrt(valX+valY);
+
+               for(int idim=0; idim<dim; idim++)
+                       velocity_field[j,idim]=0;
+                       
+        if (val<rayon)
+            pressure_field(j) = 155e5;
+        else
+            pressure_field(j) = 70e5;
+    }
+}
+
+void addValue( int i, int j, Matrix M, Mat * mat )
+{
+    int I,J;
+    for (int k=0; k<M.getNumberOfRows(); k++)
+        for (int l=0; l<M.getNumberOfColumns(); l++)
+        {
+            I=i+k;
+            J=j+l;
+            MatSetValues( *mat,1, &I, 1, &J, &M(k,l), ADD_VALUES);
+        }
+}
+
+Matrix jacobianMatrices(Vector normal, double coeff)
+{
+    int dim=normal.size();
+    Matrix    A(dim+1,dim+1);
+    Matrix absA(dim+1,dim+1);
+
+    absA(0,0)=c0*coeff;
+    for(int i=0 ; i<dim; i++)
+    {
+        A(i+1,0)=      normal[i]*coeff;
+        A(0,i+1)=c0*c0*normal[i]*coeff;
+        for( int j=0 ; j<dim; j++)
+            absA(i+1,j+1)=c0*normal[i]*normal[j]*coeff;
+    }
+    return (A - absA)*(1./2);
+}
+    
+void computeDivergenceMatrix(Mesh my_mesh, Mat * implMat, double dt)
+{
+    int nbCells = my_mesh.getNumberOfCells();
+    int dim=my_mesh.getMeshDimension();
+    int nbComp=dim+1;
+    Vector normal(dim);
+
+    Matrix idMoinsJacCL(nbComp);
+    
+    for(int j=0; j<nbCells; j++)//On parcourt les cellules
+    {
+        Cell Cj = my_mesh.getCell(j);
+        int nbFaces = Cj.getNumberOfFaces();
+
+        for(int k=0; k<nbFaces; k++)
+        {
+            int indexFace = Cj.getFacesId()[k];
+            Face Fk = my_mesh.getFace(indexFace);
+            for(int i =0; i<dim ; i++)
+                normal[i] = Cj.getNormalVector(k, i);//normale sortante
+
+            Matrix Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
+
+            int cellAutre =-1;
+            if ( not Fk.isBorder())
+            {
+                /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
+                if (Fk.getCellsId()[0] == j) 
+                    // hypothese verifiée 
+                    cellAutre = Fk.getCellsId()[1];
+                else if(Fk.getCellsId()[1] == j) 
+                    // hypothese non verifiée 
+                    cellAutre = Fk.getCellsId()[0];
+                else
+                    throw CdmathException("computeDivergenceMatrix: problem with mesh, unknown cell number");
+                    
+                addValue(j*nbComp,cellAutre*nbComp,Am      ,implMat);
+                addValue(j*nbComp,        j*nbComp,Am*(-1.),implMat);
+            }
+            else 
+            {    
+                               if( Fk.getGroupName() != "Periodic" && Fk.getGroupName() != "Neumann")//Wall boundary condition unless Periodic/Neumann specified explicitly
+                {
+                                       Vector v(dim+1);
+                    for(int i=0; i<dim; i++)
+                        v[i+1]=normal[i];
+                    Matrix idMoinsJacCL=v.tensProduct(v)*2;
+                    
+                    addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL,implMat);
+                               }
+                else if( Fk.getGroupName() == "Periodic")//Periodic boundary condition
+                {
+                    int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
+                    Face Fp = my_mesh.getFace(indexFP);
+                    cellAutre = Fp.getCellsId()[0];
+                    
+                    addValue(j*nbComp,cellAutre*nbComp,Am      ,implMat);
+                    addValue(j*nbComp,        j*nbComp,Am*(-1.),implMat);
+                               }
+                else if(Fk.getGroupName() != "Neumann")//Nothing to do for Neumann boundary condition
+                {
+                    cout<< Fk.getGroupName() <<endl;
+                    throw CdmathException("computeDivergenceMatrix: Unknown boundary condition name");
+                               }
+            }
+        }   
+    }     
+}
+
+void WaveSystem2D(double tmax, int ntmax, double cfl, int output_freq, const Mesh& my_mesh, const string file, int rank, int size)
+{
+       /* Time iteration variables */
+    int it=0;
+    bool isStationary=false;
+    double time=0.;
+    double dt;
+    double norm;
+           
+       /* PETSc variables */
+       int globalNbUnknowns;
+       int  localNbUnknowns;
+       int d_nnz, o_nnz;
+       Vec Un, dUn, Un_seq;
+    Mat divMat;
+       VecScatter scat;
+
+       
+       /* Mesh parameters managed only by proc 0 */
+    int nbCells;
+    int dim;
+    int nbComp;            
+       int idx;//Index where to add the vector values
+       double value;//value to add in the vector       
+    Field pressure_field, velocity_field;
+    std::string meshName;
+
+       if(rank == 0)
+           globalNbUnknowns=my_mesh.getNumberOfCells()*(my_mesh.getMeshDimension()+1);//nbCells*nbComp
+       
+       MPI_Bcast(&globalNbUnknowns, 1, MPI_INT, 0, MPI_COMM_WORLD);
+    /* iteration vectors */
+       VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE    ,globalNbUnknowns,&Un);
+       VecDuplicate (Un,&dUn);
+       if(rank == 0)
+               VecCreateSeq(PETSC_COMM_SELF,globalNbUnknowns,&Un_seq);//For saving results on proc 0
+
+       VecScatterCreateToZero(Un,&scat,&Un_seq);
+
+       /* System matrix */
+       VecGetLocalSize(Un, &localNbUnknowns);
+
+       if(rank == 0)
+       {
+           d_nnz=(my_mesh.getMaxNbNeighbours(CELLS)+1)*(my_mesh.getMeshDimension()+1);//(nbVoisinsMax+1)*nbComp
+           o_nnz=                                       my_mesh.getMeshDimension()+1 ;//                 nbComp
+       }
+       MPI_Bcast(&d_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
+       MPI_Bcast(&o_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
+
+       MatCreateAIJ(PETSC_COMM_WORLD,localNbUnknowns,localNbUnknowns,globalNbUnknowns,globalNbUnknowns,d_nnz,NULL,o_nnz,NULL,&divMat);
+
+       if(rank == 0)
+               {
+               /* Retrieve mesh data */
+           nbCells = my_mesh.getNumberOfCells();
+           dim=my_mesh.getMeshDimension();
+           nbComp=dim+1;           
+           meshName=my_mesh.getName();
+           double dx_min=my_mesh.minRatioVolSurf();
+               dt = cfl * dx_min / c0;
+               
+           /* Initial conditions */
+           cout<<"Building the initial condition on processor 0" << endl;
+           
+           pressure_field=Field("Pressure",CELLS,my_mesh,1) ;
+           velocity_field=Field("Velocity",CELLS,my_mesh,1) ;
+           initial_conditions_shock(my_mesh,pressure_field, velocity_field);
+       
+           cout << "Saving the solution at T=" << time <<"  on processor 0"<<endl;
+           pressure_field.setTime(time,it);
+           pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure");
+           velocity_field.setTime(time,it);
+           velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity");
+           /* --------------------------------------------- */
+       
+
+           for(int k =0; k<nbCells; k++)
+           {
+                       idx = k*nbComp;
+                       value=pressure_field[k];//vale to add in the vector
+                       VecSetValues(Un,1,&idx,&value,INSERT_VALUES);
+                       for(int idim =0; idim<dim; idim++)
+                       {
+                               idx = k*nbComp+1+idim;
+                               value =rho0*velocity_field[k,idim];
+                               VecSetValues(Un,1,&idx,&value,INSERT_VALUES);
+                       }
+               }
+           computeDivergenceMatrix(my_mesh,&divMat,dt);
+       }               
+
+       VecAssemblyBegin(Un);
+       VecAssemblyEnd(Un);
+       
+       //VecView(Un, PETSC_VIEWER_STDOUT_WORLD );
+
+       MatAssemblyBegin(divMat, MAT_FINAL_ASSEMBLY);
+       MatAssemblyEnd(  divMat, MAT_FINAL_ASSEMBLY);
+
+       //MatView(divMat,       PETSC_VIEWER_STDOUT_WORLD );
+
+       MPI_Bcast(&dt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+
+    /* Time loop */
+       PetscPrintf(PETSC_COMM_WORLD,"Starting computation of the linear wave system on all processors : \n\n");
+    while (it<ntmax && time <= tmax && !isStationary)
+    {
+               MatMult(divMat,Un,dUn);
+        VecAXPY(Un,-1,dUn);
+        
+        time=time+dt;
+        it=it+1;
+               VecNorm(dUn,NORM_2,&norm);
+        isStationary = norm<precision;
+        /* Sauvegardes */
+        if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax )
+        {
+                       PetscPrintf(PETSC_COMM_WORLD,"-- Iteration: %d, Time: %f, dt: %f, saving results on processor 0 \n", it, time, dt);
+                       VecScatterBegin(scat,Un,Un_seq,INSERT_VALUES,SCATTER_FORWARD);
+                       VecScatterEnd(  scat,Un,Un_seq,INSERT_VALUES,SCATTER_FORWARD);
+
+                       if(rank == 0)
+                       {
+                   for(int k=0; k<nbCells; k++)
+                   {
+                                       idx = k*(dim+1)+0;
+                                       VecGetValues(Un_seq,1,&idx,&value);
+                       pressure_field[k]  =value;
+                                       for(int idim =0; idim<dim; idim++)
+                                       {
+                                               idx = k*nbComp+1+idim;
+                                               VecGetValues(Un_seq,1,&idx,&value);
+                                               velocity_field[k,idim] = value/rho0;
+                                       }
+                               }
+                   pressure_field.setTime(time,it);
+                   pressure_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_pressure",false);
+                   velocity_field.setTime(time,it);
+                   velocity_field.writeVTK("WaveSystem"+to_string(dim)+"DUpwind"+meshName+"_velocity",false);
+                       }
+               }
+    }
+
+       PetscPrintf(PETSC_COMM_WORLD,"\n End of calculations at iteration: %d, and time: %f\n", it, time);
+    if(it>=ntmax)
+        PetscPrintf(PETSC_COMM_WORLD, "Nombre de pas de temps maximum ntmax= %d atteint\n", ntmax);
+    else if(isStationary)
+        PetscPrintf(PETSC_COMM_WORLD, "Régime stationnaire atteint au pas de temps %d, t= %f\n", it, time);       
+    else
+        PetscPrintf(PETSC_COMM_WORLD, "Temps maximum tmax= %f atteint\n", tmax);
+
+       VecDestroy(&Un);
+       VecDestroy(&Un_seq);
+       VecDestroy(&dUn);
+       MatDestroy(&divMat);
+}
+int main(int argc, char *argv[])
+{
+       /* PETSc initialisation */
+       PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
+       PetscMPIInt    size;        /* size of communicator */
+       PetscMPIInt    rank;        /* processor rank */
+       MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
+       MPI_Comm_size(PETSC_COMM_WORLD,&size);
+       
+    // Problem data
+    double cfl=0.49;
+    double tmax=1.;
+    int ntmax=80;//20000;
+    int freqSortie=10;
+    string fileOutPut="SphericalWave";
+    Mesh myMesh;
+
+       if(size>1)
+               PetscPrintf(PETSC_COMM_WORLD,"---- More than one processor detected : running a parallel simulation ----\n");
+               PetscPrintf(PETSC_COMM_WORLD,"---- Limited parallelism : input and output remain sequential ----\n");
+               PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrix-vector products are done in parallel ----\n");
+               PetscPrintf(PETSC_COMM_WORLD,"---- Processor 0 is in charge of building the mesh, saving the results, filling and then distributing the matrix to other processors.\n\n");
+               
+       if(rank == 0)
+       {
+           cout << "-- Starting the RESOLUTION OF THE 2D WAVE SYSTEM on "<< size <<" processors"<<endl;
+           cout << "- Numerical scheme : Upwind explicit scheme" << endl;
+           cout << "- Boundary conditions : WALL" << endl;
+       
+               if(argc<2)
+               {
+                   cout << "- DOMAIN : SQUARE" << endl;
+                   cout << "- MESH : CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << endl<< endl;
+                   cout << "Construction of a cartesian mesh on processor 0" << endl;
+                   double xinf=0.0;
+                   double xsup=1.0;
+                   double yinf=0.0;
+                   double ysup=1.0;
+                   int nx=50;
+                   int ny=50;
+                   myMesh=Mesh(xinf,xsup,nx,yinf,ysup,ny);
+                   
+                   double eps=1.E-10;
+                   myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
+                   myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
+                   myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
+                   myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
+               }
+               else
+               {
+                   cout << "- MESH:  GENERATED EXTERNALLY WITH SALOME" << endl<< endl;
+                   cout << "Loading of mesh named "<<argv[1]<<" on processor 0" << endl;
+                   string filename = argv[1];
+                   myMesh=Mesh(filename);
+               }
+       }
+       WaveSystem2D(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut, rank, size);
+
+       if(rank == 0)
+               cout << "Simulation complete." << endl;
+
+       PetscFinalize();
+    return 0;
+}