CoreFlows/examples/C/WaveSystem_FV_SphericalExplosion_MPI.cxx

   1 //============================================================================
   2 // Author      : Michael NDJINGA
   3 // Date        : November 2020
   4 // Description : multiD linear wave system run in parallel with calls to MPI
   5 //               Test used in the MPI version of the Salome platform
   6 //============================================================================
   7
   8 #include <iostream>
   9 #include <string>
  10 #include <cmath>
  11
  12 #include "Mesh.hxx"
  13 #include "Cell.hxx"
  14 #include "Face.hxx"
  15 #include "Field.hxx"
  16 #include "CdmathException.hxx"
  17
  18 #include <petscksp.h>
  19
  20 using namespace std;
  21
  22 double p0  =155e5;   //reference pressure in a pressurised nuclear vessel
  23 double c0  =700.;    //reference sound speed for water at 155 bars, 600K
  24 double rho0=p0/c0*c0;//reference density
  25 double precision=1e-5;
  26
  27 void initial_conditions_shock(Mesh my_mesh,Field& pressure_field,Field& velocity_field)
  28 {
  29     double rayon=0.15;
  30     double xcentre=0.;
  31     double ycentre=0;
  32     double zcentre=0;
  33
  34         double x, y, z;
  35         double val, valX, valY, valZ;
  36
  37     int dim    =my_mesh.getMeshDimension();
  38     int nbCells=my_mesh.getNumberOfCells();
  39
  40     for (int j=0 ; j<nbCells ; j++)
  41     {
  42         x = my_mesh.getCell(j).x() ;
  43                 if(dim>1)
  44                 {
  45                         y = my_mesh.getCell(j).y() ;
  46                         if(dim==3)
  47                                 z = my_mesh.getCell(j).z() ;
  48                 }
  49
  50         valX=(x-xcentre)*(x-xcentre);
  51                 if(dim==1)
  52                         val=sqrt(valX);
  53                 else if(dim==2)
  54                 {
  55                         valY=(y-ycentre)*(y-ycentre);
  56                         val=sqrt(valX+valY);
  57                 }
  58                 else if(dim==3)
  59                 {
  60                         valY=(y-ycentre)*(y-ycentre);
  61                         valZ=(z-zcentre)*(z-zcentre);
  62                         val=sqrt(valX+valY+valZ);
  63                 }
  64
  65                 for(int idim=0; idim<dim; idim++)
  66                         velocity_field[j,idim]=0;
  67
  68         if (val<rayon)
  69             pressure_field(j) = 155e5;
  70         else
  71             pressure_field(j) = 70e5;
  72     }
  73 }
  74
  75 void addValue( int i, int j, Matrix M, Mat * mat )
  76 {
  77     int I,J;
  78     for (int k=0; k<M.getNumberOfRows(); k++)
  79         for (int l=0; l<M.getNumberOfColumns(); l++)
  80         {
  81             I=i+k;
  82             J=j+l;
  83             MatSetValues( *mat,1, &I, 1, &J, &M(k,l), ADD_VALUES);
  84         }
  85 }
  86
  87 Matrix jacobianMatrices(Vector normal, double coeff)
  88 {
  89     int dim=normal.size();
  90     Matrix    A(dim+1,dim+1);
  91     Matrix absA(dim+1,dim+1);
  92
  93     absA(0,0)=c0*coeff;
  94     for(int i=0 ; i<dim; i++)
  95     {
  96         A(i+1,0)=      normal[i]*coeff;
  97         A(0,i+1)=c0*c0*normal[i]*coeff;
  98         for( int j=0 ; j<dim; j++)
  99             absA(i+1,j+1)=c0*normal[i]*normal[j]*coeff;
 100     }
 101     return (A - absA)*(1./2);
 102 }
 103
 104 void computeDivergenceMatrix(Mesh my_mesh, Mat * implMat, double dt)
 105 {
 106     int nbCells = my_mesh.getNumberOfCells();
 107     int dim=my_mesh.getMeshDimension();
 108     int nbComp=dim+1;
 109     Vector normal(dim);
 110
 111     Matrix idMoinsJacCL(nbComp);
 112
 113     for(int j=0; j<nbCells; j++)//On parcourt les cellules
 114     {
 115         Cell Cj = my_mesh.getCell(j);
 116         int nbFaces = Cj.getNumberOfFaces();
 117
 118         for(int k=0; k<nbFaces; k++)
 119         {
 120             int indexFace = Cj.getFacesId()[k];
 121             Face Fk = my_mesh.getFace(indexFace);
 122             for(int i =0; i<dim ; i++)
 123                 normal[i] = Cj.getNormalVector(k, i);//normale sortante
 124
 125             Matrix Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
 126
 127             int cellAutre =-1;
 128             if ( not Fk.isBorder())
 129             {
 130                 /* hypothese: La cellule d'index indexC1 est la cellule courante index j */
 131                 if (Fk.getCellsId()[0] == j)
 132                     // hypothese verifiée
 133                     cellAutre = Fk.getCellsId()[1];
 134                 else if(Fk.getCellsId()[1] == j)
 135                     // hypothese non verifiée
 136                     cellAutre = Fk.getCellsId()[0];
 137                 else
 138                     throw CdmathException("computeDivergenceMatrix: problem with mesh, unknown cell number");
 139
 140                 addValue(j*nbComp,cellAutre*nbComp,Am      ,implMat);
 141                 addValue(j*nbComp,        j*nbComp,Am*(-1.),implMat);
 142             }
 143             else
 144             {
 145                                 if( Fk.getGroupName() != "Periodic" && Fk.getGroupName() != "Neumann")//Wall boundary condition unless Periodic/Neumann specified explicitly
 146                 {
 147                                         Vector v(dim+1);
 148                     for(int i=0; i<dim; i++)
 149                         v[i+1]=normal[i];
 150                     Matrix idMoinsJacCL=v.tensProduct(v)*2;
 151
 152                     addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL,implMat);
 153                                 }
 154                 else if( Fk.getGroupName() == "Periodic")//Periodic boundary condition
 155                 {
 156                     int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
 157                     Face Fp = my_mesh.getFace(indexFP);
 158                     cellAutre = Fp.getCellsId()[0];
 159
 160                     addValue(j*nbComp,cellAutre*nbComp,Am      ,implMat);
 161                     addValue(j*nbComp,        j*nbComp,Am*(-1.),implMat);
 162                                 }
 163                 else if(Fk.getGroupName() != "Neumann")//Nothing to do for Neumann boundary condition
 164                 {
 165                     cout<< Fk.getGroupName() <<endl;
 166                     throw CdmathException("computeDivergenceMatrix: Unknown boundary condition name");
 167                                 }
 168             }
 169         }
 170     }
 171 }
 172
 173 void WaveSystem(double tmax, int ntmax, double cfl, int output_freq, const Mesh& my_mesh, const string file, int rank, int size, string resultDirectory)
 174 {
 175         /* Time iteration variables */
 176     int it=0;
 177     bool isStationary=false;
 178     double time=0.;
 179     double dt;
 180     double norm;
 181
 182         /* PETSc variables */
 183         int globalNbUnknowns;
 184         int  localNbUnknowns;
 185         int d_nnz, o_nnz;
 186         Vec Un, dUn, Un_seq;
 187     Mat divMat;
 188         VecScatter scat;
 189
 190
 191         /* Mesh parameters managed only by proc 0 */
 192     int nbCells;
 193     int dim;
 194     int nbComp;
 195         int idx;//Index where to add the vector values
 196         double value;//value to add in the vector
 197     Field pressure_field, velocity_field;
 198     std::string meshName;
 199
 200         if(rank == 0)
 201             globalNbUnknowns=my_mesh.getNumberOfCells()*(my_mesh.getMeshDimension()+1);//nbCells*nbComp
 202
 203         MPI_Bcast(&globalNbUnknowns, 1, MPI_INT, 0, MPI_COMM_WORLD);
 204
 205     /* iteration vectors */
 206         VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE    ,globalNbUnknowns,&Un);
 207         VecDuplicate (Un,&dUn);
 208         if(rank == 0)
 209                 VecCreateSeq(PETSC_COMM_SELF,globalNbUnknowns,&Un_seq);//For saving results on proc 0
 210
 211         VecScatterCreateToZero(Un,&scat,&Un_seq);
 212
 213         /* System matrix */
 214         VecGetLocalSize(Un, &localNbUnknowns);
 215
 216         if(rank == 0)
 217         {
 218                 int nbVoisinsMax = my_mesh.getMaxNbNeighbours(CELLS);
 219             d_nnz=(nbVoisinsMax+1)*(my_mesh.getMeshDimension()+1);//(nbVoisinsMax+1)*nbComp
 220             o_nnz= nbVoisinsMax   *(my_mesh.getMeshDimension()+1);//                 nbComp
 221         }
 222         MPI_Bcast(&d_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
 223         MPI_Bcast(&o_nnz, 1, MPI_INT, 0, MPI_COMM_WORLD);
 224
 225         MatCreateAIJ(PETSC_COMM_WORLD,localNbUnknowns,localNbUnknowns,globalNbUnknowns,globalNbUnknowns,d_nnz,NULL,o_nnz,NULL,&divMat);
 226
 227         if(rank == 0)
 228                 {
 229                 /* Retrieve mesh data */
 230             nbCells = my_mesh.getNumberOfCells();
 231             dim=my_mesh.getMeshDimension();
 232             nbComp=dim+1;
 233             meshName=my_mesh.getName();
 234             double dx_min=my_mesh.minRatioVolSurf();
 235                 dt = cfl * dx_min / c0;
 236
 237             /* Initial conditions */
 238             cout<<"Building the initial condition on processor 0" << endl;
 239
 240             pressure_field=Field("Pressure",CELLS,my_mesh,1) ;
 241             velocity_field=Field("Velocity",CELLS,my_mesh,dim) ;
 242             initial_conditions_shock(my_mesh,pressure_field, velocity_field);
 243
 244             cout << "Saving the solution at T=" << time <<"  on processor 0"<<endl;
 245             pressure_field.setTime(time,it);
 246             pressure_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_pressure");
 247             velocity_field.setTime(time,it);
 248             velocity_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_velocity");
 249             /* --------------------------------------------- */
 250
 251
 252             for(int k =0; k<nbCells; k++)
 253             {
 254                         idx = k*nbComp;
 255                         value=pressure_field[k];//vale to add in the vector
 256                         VecSetValues(Un,1,&idx,&value,INSERT_VALUES);
 257                         for(int idim =0; idim<dim; idim++)
 258                         {
 259                                 idx = k*nbComp+1+idim;
 260                                 value =rho0*velocity_field[k,idim];
 261                                 VecSetValues(Un,1,&idx,&value,INSERT_VALUES);
 262                         }
 263                 }
 264             computeDivergenceMatrix(my_mesh,&divMat,dt);
 265         }
 266
 267         VecAssemblyBegin(Un);
 268         VecAssemblyEnd(Un);
 269
 270         //VecView(Un, PETSC_VIEWER_STDOUT_WORLD );
 271
 272         MatAssemblyBegin(divMat, MAT_FINAL_ASSEMBLY);
 273         MatAssemblyEnd(  divMat, MAT_FINAL_ASSEMBLY);
 274
 275         //MatView(divMat,       PETSC_VIEWER_STDOUT_WORLD );
 276
 277         MPI_Bcast(&dt, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
 278
 279     /* Time loop */
 280         PetscPrintf(PETSC_COMM_WORLD,"Starting computation of the linear wave system on all processors : \n\n");
 281     while (it<ntmax && time <= tmax && !isStationary)
 282     {
 283                 MatMult(divMat,Un,dUn);
 284         VecAXPY(Un,-1,dUn);
 285
 286         time=time+dt;
 287         it=it+1;
 288
 289                 VecNorm(dUn,NORM_2,&norm);
 290         isStationary = norm<precision;
 291         /* Sauvegardes */
 292         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax )
 293         {
 294                         PetscPrintf(PETSC_COMM_WORLD,"-- Iteration: %d, Time: %f, dt: %f, saving results on processor 0 \n", it, time, dt);
 295                         VecScatterBegin(scat,Un,Un_seq,INSERT_VALUES,SCATTER_FORWARD);
 296                         VecScatterEnd(  scat,Un,Un_seq,INSERT_VALUES,SCATTER_FORWARD);
 297
 298                         if(rank == 0)
 299                         {
 300                     for(int k=0; k<nbCells; k++)
 301                     {
 302                                         idx = k*(dim+1)+0;
 303                                         VecGetValues(Un_seq,1,&idx,&value);
 304                         pressure_field[k]  =value;
 305                                         for(int idim =0; idim<dim; idim++)
 306                                         {
 307                                                 idx = k*nbComp+1+idim;
 308                                                 VecGetValues(Un_seq,1,&idx,&value);
 309                                                 velocity_field[k,idim] = value/rho0;
 310                                         }
 311                                 }
 312                     pressure_field.setTime(time,it);
 313                     pressure_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_pressure",false);
 314                     velocity_field.setTime(time,it);
 315                     velocity_field.writeMED(resultDirectory+"/WaveSystem"+to_string(dim)+"DUpwind_"+to_string(size)+"Procs_"+meshName+"_velocity",false);
 316                         }
 317                 }
 318     }
 319
 320         PetscPrintf(PETSC_COMM_WORLD,"\n End of calculations at iteration: %d, and time: %f\n", it, time);
 321     if(it>=ntmax)
 322         PetscPrintf(PETSC_COMM_WORLD, "Nombre de pas de temps maximum ntmax= %d atteint\n", ntmax);
 323     else if(isStationary)
 324         PetscPrintf(PETSC_COMM_WORLD, "Régime stationnaire atteint au pas de temps %d, t= %f\n", it, time);
 325     else
 326         PetscPrintf(PETSC_COMM_WORLD, "Temps maximum tmax= %f atteint\n", tmax);
 327
 328         VecDestroy(&Un);
 329         VecDestroy(&Un_seq);
 330         VecDestroy(&dUn);
 331         MatDestroy(&divMat);
 332 }
 333
 334 int main(int argc, char *argv[])
 335 {
 336         /* PETSc initialisation */
 337         PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
 338         PetscMPIInt    size;        /* size of communicator */
 339         PetscMPIInt    rank;        /* processor rank */
 340         MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 341         MPI_Comm_size(PETSC_COMM_WORLD,&size);
 342
 343     // Problem data
 344     double tmax=1.;
 345     int ntmax=2000;
 346     int freqSortie=10;
 347     string fileOutPut="SphericalWave";
 348     Mesh myMesh;
 349         string resultDirectory="./";
 350
 351         if(size>1)
 352                 PetscPrintf(PETSC_COMM_WORLD,"---- More than one processor detected : running a parallel simulation ----\n");
 353                 PetscPrintf(PETSC_COMM_WORLD,"---- Limited parallelism : input and output remain sequential ----\n");
 354                 PetscPrintf(PETSC_COMM_WORLD,"---- Only the matrix-vector products are done in parallel ----\n");
 355                 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");
 356
 357         if(rank == 0)
 358         {
 359             cout << "-- Starting the RESOLUTION OF THE 2D WAVE SYSTEM on "<< size <<" processors"<<endl;
 360             cout << "- Numerical scheme : Upwind explicit scheme" << endl;
 361             cout << "- Boundary conditions : WALL" << endl;
 362
 363                 /* Read or create mesh */
 364                 if(argc<2)
 365                 {
 366                     cout << "- DOMAIN : SQUARE" << endl;
 367                     cout << "- MESH : CARTESIAN, GENERATED INTERNALLY WITH CDMATH" << endl<< endl;
 368                     cout << "Construction of a cartesian mesh on processor 0" << endl;
 369                     double xinf=-0.5;
 370                     double xsup= 0.5;
 371                     double yinf=-0.5;
 372                     double ysup= 0.5;
 373                     int nx=50;
 374                     int ny=50;
 375                     myMesh=Mesh(xinf,xsup,nx,yinf,ysup,ny);
 376
 377                     double eps=precision;
 378                     myMesh.setGroupAtPlan(xsup,0,eps,"RightEdge");
 379                     myMesh.setGroupAtPlan(xinf,0,eps,"LeftEdge");
 380                     myMesh.setGroupAtPlan(yinf,1,eps,"BottomEdge");
 381                     myMesh.setGroupAtPlan(ysup,1,eps,"TopEdge");
 382                 }
 383                 else
 384                 {
 385                     cout << "- MESH:  GENERATED EXTERNALLY WITH SALOME" << endl<< endl;
 386                     cout << "Loading of a mesh named "<<argv[1]<<" on processor 0" << endl;
 387                     string filename = argv[1];
 388                     myMesh=Mesh(filename);
 389                 }
 390
 391                 /*Detect directory where to same results*/
 392                 if(argc>2)
 393                     resultDirectory = argv[2];
 394         }
 395
 396     double cfl=1./myMesh.getSpaceDimension();
 397         WaveSystem(tmax,ntmax,cfl,freqSortie,myMesh,fileOutPut, rank, size, resultDirectory);
 398
 399         if(rank == 0)
 400                 cout << "Simulation complete." << endl;
 401
 402         PetscFinalize();
 403     return 0;
 404 }