From b71fc203edda97a03592f60b301f8d21bac4a0fa Mon Sep 17 00:00:00 2001
From: michael <michael@localhost.localdomain>
Date: Mon, 14 Dec 2020 17:12:43 +0100
Subject: [PATCH] First MPI tests

---
 CoreFlows/examples/C/CMakeLists.txt           |  15 +-
 ...eSystem_2DFV_SphericalExplosion_CDMATH.cxx | 277 +++++++++++++
 ...WaveSystem_2DFV_SphericalExplosion_MPI.cxx | 374 ++++++++++++++++++
 3 files changed, 662 insertions(+), 4 deletions(-)
 create mode 100755 CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx
 create mode 100755 CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_MPI.cxx

diff --git a/CoreFlows/examples/C/CMakeLists.txt b/CoreFlows/examples/C/CMakeLists.txt
index d9a5aa4..0abbd28 100755
--- a/CoreFlows/examples/C/CMakeLists.txt
+++ b/CoreFlows/examples/C/CMakeLists.txt
@@ -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
index 0000000..24ad784
--- /dev/null
+++ b/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_CDMATH.cxx
@@ -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
index 0000000..0d1cc35
--- /dev/null
+++ b/CoreFlows/examples/C/WaveSystem_2DFV_SphericalExplosion_MPI.cxx
@@ -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;
+}
-- 
2.39.2