CDMATH/tests/examples/WaveSystem_Shock/WaveSystemStaggered/WaveSystemStaggered.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 #===============================================================================================================================
   5 # Name        : Résolution VF du système des ondes 2D sans terme source
   6 #                \partial_t p + c^2 \div q = 0
   7 #                \partial_t q +    \grad p = 0
   8 # Author      : Michaël Ndjinga
   9 # Copyright   : CEA Saclay 2019
  10 # Description : Propagation d'une onde de choc sphérique
  11 #               Utilisation du schéma MAC sur un maillage cartésien
  12 #               Initialisation par une surpression sphérique
  13 #               Conditions aux limites périodiques
  14 #                       Création et sauvegarde du champ résultant et des figures
  15 #================================================================================================================================
  16
  17
  18 from math import sqrt
  19 import cdmath
  20 import PV_routines
  21 import VTK_routines
  22 import sys
  23
  24 p0=155e5#reference pressure in a pressurised nuclear vessel
  25 c0=700.#reference sound speed for water at 155 bars, 600K
  26 rho0=p0/c0*c0#reference density
  27 precision=1e-5
  28
  29 def initial_conditions_shock(my_mesh, isCircle):
  30     print( "Initial data : Spherical wave")
  31     dim     = my_mesh.getMeshDimension()
  32     nbCells = my_mesh.getNumberOfCells()
  33
  34     rayon = 0.15
  35     if(not isCircle):
  36         xcentre = 0.5
  37         ycentre = 0.5
  38         zcentre = 0.5
  39     else:
  40         xcentre = 0.
  41         ycentre = 0.
  42         zcentre = 0.
  43
  44     pressure_field = cdmath.Field("Pressure",            cdmath.CELLS, my_mesh, 1)
  45     velocity_field = cdmath.Field("Velocity",            cdmath.CELLS, my_mesh, 3)
  46
  47     for i in range(nbCells):
  48         velocity_field[i,0] = 0
  49         velocity_field[i,1] = 0
  50         velocity_field[i,2] = 0
  51
  52         x = my_mesh.getCell(i).x()
  53         valX = (x - xcentre) * (x - xcentre)
  54
  55         if(dim==1):
  56             val =  sqrt(valX)
  57         if(dim==2):
  58             y = my_mesh.getCell(i).y()
  59             valY = (y - ycentre) * (y - ycentre)
  60             val =  sqrt(valX + valY)
  61         if(dim==3):
  62             y = my_mesh.getCell(i).y()
  63             z = my_mesh.getCell(i).z()
  64             valY = (y - ycentre) * (y - ycentre)
  65             valZ = (z - zcentre) * (z - zcentre)
  66             val =  sqrt(valX + valY + valZ)
  67
  68         if val < rayon:
  69             pressure_field[i] = p0
  70             pass
  71         else:
  72             pressure_field[i] = p0/2
  73             pass
  74         pass
  75
  76     return pressure_field, velocity_field
  77
  78 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling):
  79     nbCells = my_mesh.getNumberOfCells()
  80     dim=my_mesh.getMeshDimension()
  81     nbComp=dim+1
  82
  83     if(not my_mesh.isStructured()):
  84         raise ValueError("WaveSystemStaggered: the mesh should be structured");
  85
  86     NxNyNz=my_mesh.getCellGridStructure()
  87     DxDyDz=my_mesh.getDXYZ()
  88
  89     implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
  90
  91     idMoinsJacCL=cdmath.Matrix(nbComp)
  92
  93     if( dim == 1) :
  94         nx=NxNyNz[0]
  95         dx=DxDyDz[0]
  96
  97         if( scaling==0 ):
  98             for k in range(nbCells):
  99                 implMat.addValue(k,1*nbCells +  k      , -c0*c0*dt/dx)
 100                 implMat.addValue(k,1*nbCells + (k+1)%nx,  c0*c0*dt/dx)
 101
 102                 implMat.addValue(  1*nbCells +  k      ,k,  dt/dx)
 103                 implMat.addValue(  1*nbCells + (k+1)%nx,k, -dt/dx)
 104         else : # scaling >0
 105             for k in range(nbCells):
 106                 implMat.addValue(k,1*nbCells +  k      , -c0*dt/dx)
 107                 implMat.addValue(k,1*nbCells + (k+1)%nx,  c0*dt/dx)
 108
 109                 implMat.addValue(  1*nbCells +  k      ,k,  c0*dt/dx)
 110                 implMat.addValue(  1*nbCells + (k+1)%nx,k, -c0*dt/dx)
 111
 112     elif( dim == 2) :# k = j*nx+i
 113         nx=NxNyNz[0]
 114         ny=NxNyNz[1]
 115         dx=DxDyDz[0]
 116         dy=DxDyDz[1]
 117
 118         if( scaling==0 ):
 119             for k in range(nbCells):
 120                 i = k % nx
 121                 j = k //nx
 122
 123                 implMat.addValue(k,1*nbCells + j*nx +  i      ,   -c0*c0*dt/dx)
 124                 implMat.addValue(k,1*nbCells + j*nx + (i+1)%nx,    c0*c0*dt/dx)
 125
 126                 implMat.addValue(k,2*nbCells +   j       *nx + i, -c0*c0*dt/dy)
 127                 implMat.addValue(k,2*nbCells + ((j+1)%ny)*nx + i,  c0*c0*dt/dy)
 128
 129                 implMat.addValue(  1*nbCells + j*nx +  i      ,  k,  dt/dx)
 130                 implMat.addValue(  1*nbCells + j*nx + (i+1)%nx,  k, -dt/dx)
 131
 132                 implMat.addValue(  2*nbCells +   j       *nx + i,k,  dt/dy)
 133                 implMat.addValue(  2*nbCells + ((j+1)%ny)*nx + i,k, -dt/dy)
 134
 135         else :# scaling >0
 136             for k in range(nbCells):
 137                 i = k % nx
 138                 j = k //nx
 139
 140                 implMat.addValue(k,1*nbCells + j*nx +  i      ,   -c0*dt/dx)
 141                 implMat.addValue(k,1*nbCells + j*nx + (i+1)%nx,    c0*dt/dx)
 142
 143                 implMat.addValue(k,2*nbCells +   j       *nx + i, -c0*dt/dy)
 144                 implMat.addValue(k,2*nbCells + ((j+1)%ny)*nx + i,  c0*dt/dy)
 145
 146                 implMat.addValue(  1*nbCells + j*nx +  i      ,  k,  c0*dt/dx)
 147                 implMat.addValue(  1*nbCells + j*nx + (i+1)%nx,  k, -c0*dt/dx)
 148
 149                 implMat.addValue(  2*nbCells +   j       *nx + i,k,  c0*dt/dy)
 150                 implMat.addValue(  2*nbCells + ((j+1)%ny)*nx + i,k, -c0*dt/dy)
 151
 152     elif( dim == 3) :# k = l*nx*ny+j*nx+i
 153         nx=NxNyNz[0]
 154         ny=NxNyNz[1]
 155         nz=NxNyNz[2]
 156         dx=DxDyDz[0]
 157         dy=DxDyDz[1]
 158         dz=DxDyDz[2]
 159
 160         if( scaling==0 ):
 161             for k in range(nbCells):
 162                 i =  k % nx
 163                 j = (k //nx)%ny
 164                 l =  k //(nx*ny)
 165
 166                 implMat.addValue(k,1*nbCells + l*nx*ny + j*nx +  i      ,  -c0*c0*dt/dx)
 167                 implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + (i+1)%nx,   c0*c0*dt/dx)
 168
 169                 implMat.addValue(k,2*nbCells + l*nx*ny +   j       *nx + i, -c0*c0*dt/dy)
 170                 implMat.addValue(k,2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i,  c0*c0*dt/dy)
 171
 172                 implMat.addValue(k,3*nbCells +   l*nx*ny        + j*nx + i, -c0*c0*dt/dz)
 173                 implMat.addValue(k,3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i,  c0*c0*dt/dz)
 174
 175                 implMat.addValue(  1*nbCells + l*nx*ny + j*nx +  i      ,  k,  dt/dx)
 176                 implMat.addValue(  1*nbCells + l*nx*ny + j*nx + (i+1)%nx,  k, -dt/dx)
 177
 178                 implMat.addValue(  2*nbCells + l*nx*ny +   j       *nx + i,k,  dt/dy)
 179                 implMat.addValue(  2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i,k, -dt/dy)
 180
 181                 implMat.addValue(  3*nbCells +   l*nx*ny        + j*nx + i,k,  dt/dz)
 182                 implMat.addValue(  3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i,k, -dt/dz)
 183
 184         else:# scaling >0
 185             for k in range(nbCells):
 186                 i =  k % nx
 187                 j = (k //nx)%ny
 188                 l =  k //(nx*ny)
 189
 190                 implMat.addValue(k,1*nbCells + l*nx*ny + j*nx +  i      ,  -c0*dt/dx)
 191                 implMat.addValue(k,1*nbCells + l*nx*ny + j*nx + (i+1)%nx,   c0*dt/dx)
 192
 193                 implMat.addValue(k,2*nbCells + l*nx*ny +   j       *nx + i, -c0*dt/dy)
 194                 implMat.addValue(k,2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i,  c0*dt/dy)
 195
 196                 implMat.addValue(k,3*nbCells +   l*nx*ny        + j*nx + i, -c0*dt/dz)
 197                 implMat.addValue(k,3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i,  c0*dt/dz)
 198
 199                 implMat.addValue(  1*nbCells + l*nx*ny + j*nx +  i      ,  k,  c0*dt/dx)
 200                 implMat.addValue(  1*nbCells + l*nx*ny + j*nx + (i+1)%nx,  k, -c0*dt/dx)
 201
 202                 implMat.addValue(  2*nbCells + l*nx*ny +   j       *nx + i,k,  c0*dt/dy)
 203                 implMat.addValue(  2*nbCells + l*nx*ny + ((j+1)%ny)*nx + i,k, -c0*dt/dy)
 204
 205                 implMat.addValue(  3*nbCells +   l*nx*ny        + j*nx + i,k,  c0*dt/dz)
 206                 implMat.addValue(  3*nbCells + ((l+1)%nz)*nx*ny + j*nx + i,k, -c0*dt/dz)
 207
 208     return implMat
 209
 210 def WaveSystemStaggered(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution):
 211     dim=my_mesh.getMeshDimension()
 212     nbCells = my_mesh.getNumberOfCells()
 213
 214     dt = 0.
 215     time = 0.
 216     it=0;
 217     isStationary=False;
 218
 219     scaling=0
 220
 221     nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
 222     iterGMRESMax=50
 223
 224     #iteration vectors
 225     Un=cdmath.Vector(nbCells*(dim+1))
 226     dUn=cdmath.Vector(nbCells*(dim+1))
 227
 228     # Initial conditions #
 229     print("Construction of the initial condition …")
 230     if(dim==1 or meshName.find("square")>-1 or meshName.find("Square")>-1 or meshName.find("cube")>-1 or meshName.find("Cube")>-1):
 231         pressure_field, velocity_field = initial_conditions_shock(my_mesh,False)
 232     elif(meshName.find("disk")>-1 or meshName.find("Disk")>-1):
 233         pressure_field, velocity_field = initial_conditions_shock(my_mesh,True)
 234     else:
 235         print( "Mesh name : ", meshName )
 236         raise ValueError("Mesh name should contain substring square, cube or disk")
 237
 238     for k in range(nbCells):
 239         Un[k + 0*nbCells] =      pressure_field[k]
 240         Un[k + 1*nbCells] = rho0*velocity_field[k,0] # value on the left face
 241         if(dim>=2):
 242             Un[k + 2*nbCells] = rho0*velocity_field[k,1] # value on the bottom face
 243             if(dim==3):
 244                 Un[k + 3*nbCells] = rho0*velocity_field[k,2]
 245
 246     if( scaling>0):
 247         Vn = Un.deepCopy()
 248         for k in range(nbCells):
 249             Vn[k] = Vn[k]/c0
 250
 251     #sauvegarde de la donnée initiale
 252     pressure_field.setTime(time,it);
 253     pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure");
 254     velocity_field.setTime(time,it);
 255     velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity");
 256     #Postprocessing : save 2D picture
 257     PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_initial")
 258     PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_initial")
 259
 260     dx_min=my_mesh.minRatioVolSurf()
 261
 262     dt = cfl * dx_min / c0
 263
 264     divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,scaling)
 265
 266     #Add the identity matrix on the diagonal
 267     for j in range(nbCells*(dim+1)):
 268         divMat.addValue(j,j,1)
 269
 270     LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","LU")
 271
 272     print("Starting computation of the linear wave system with an pseudo staggered scheme …")
 273
 274     # Starting time loop
 275     while (it<ntmax and time <= tmax and not isStationary):
 276         dUn=Un.deepCopy()
 277         LS.setSndMember(Un)
 278         Un=LS.solve();
 279         cvgceLS=LS.getStatus();
 280         iterGMRES=LS.getNumberOfIter();
 281         if(not cvgceLS):
 282             print( "Linear system did not converge ", iterGMRES, " GMRES iterations")
 283             raise ValueError("Pas de convergence du système linéaire");
 284         dUn-=Un
 285
 286         max_dp=0 ;        max_dq=0
 287         for k in range(nbCells):
 288             max_dp = max(max_dp,abs(dUn[k]))
 289             for i in range(dim):
 290                 max_dq=max(max_dq,abs(dUn[k+(1+i)*nbCells]))
 291
 292         isStationary= max_dp/p0<precision and max_dq/rho0<precision
 293
 294         time=time+dt;
 295         it=it+1;
 296
 297         #Sauvegardes
 298         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 299             print( "-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
 300             print( "Variation temporelle relative : pressure ", max_dp/p0 ,", velocity ", max_dq/rho0 )
 301             print( "Linear system converged in ", iterGMRES, " GMRES iterations" )
 302
 303             for k in range(nbCells):
 304                 pressure_field[k]=Un[k]
 305                 velocity_field[k,0]=Un[k+1*nbCells]/rho0
 306                 if(dim>1):
 307                     velocity_field[k,1]=Un[k+2*nbCells]/rho0
 308                     if(dim>2):
 309                         velocity_field[k,2]=Un[k+3*nbCells]/rho0
 310
 311             pressure_field.setTime(time,it);
 312             pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure",False);
 313             velocity_field.setTime(time,it);
 314             velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity",False);
 315
 316     print( "-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
 317     print( "Variation temporelle relative : pressure ", max_dp/p0 ,", velocity ", max_dq/rho0 )
 318     print()
 319
 320     if(it>=ntmax):
 321         print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
 322     elif(isStationary):
 323         print( "Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
 324         print( "------------------------------------------------------------------------------------")
 325
 326         pressure_field.setTime(time,0);
 327         pressure_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat");
 328         velocity_field.setTime(time,0);
 329         velocity_field.writeVTK("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat");
 330
 331         #Postprocessing : save 2D picture
 332         PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_pressure_Stat")
 333         PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DStaggered"+meshName+"_velocity_Stat")
 334
 335     else:
 336         print( "Temps maximum Tmax= ", tmax, " atteint")
 337
 338
 339 def solve(my_mesh,meshName,resolution):
 340     print( "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension() )
 341     print( "Numerical method : staggered scheme" )
 342     print( "Initial data : Spherical wave" )
 343     print( "Periodic boundary conditions" )
 344     print( "Mesh name : ",meshName , my_mesh.getNumberOfCells(), " cells" )
 345
 346     # Problem data
 347     tmax = 1000.
 348     ntmax = 100
 349     cfl = 1./my_mesh.getSpaceDimension()
 350     output_freq = 1
 351
 352     WaveSystemStaggered(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution)
 353
 354 def solve_file( filename,meshName, resolution):
 355     my_mesh = cdmath.Mesh(filename+".med")
 356
 357     return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution)
 358
 359
 360 if __name__ == """__main__""":
 361     if len(sys.argv) >1 :
 362         my_mesh = cdmath.Mesh(sys.argv[1])
 363         solve(my_mesh,my_mesh.getName(),100)
 364     else :
 365         nx=50
 366         my_mesh = cdmath.Mesh(0,1,nx,0,1,nx)
 367         solve(my_mesh,"square",100)