CDMATH/tests/examples/WaveSystem_Shock/WaveSystemPStag/WaveSystemPStag.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 parois
  14 #                       Création et sauvegarde du champ résultant et des figures
  15 #================================================================================================================================
  16
  17
  18 from math import sqrt
  19 from numpy import sign
  20 import cdmath
  21 import PV_routines
  22 import VTK_routines
  23 import sys
  24
  25 p0=155e5#reference pressure in a pressurised nuclear vessel
  26 c0=700.#reference sound speed for water at 155 bars, 600K
  27 rho0=p0/c0*c0#reference density
  28 precision=1e-5
  29
  30 def initial_conditions_shock(my_mesh, isCircle):
  31     print "Initial data : Spherical wave"
  32     dim     = my_mesh.getMeshDimension()
  33     nbCells = my_mesh.getNumberOfCells()
  34
  35     rayon = 0.15
  36     if(not isCircle):
  37         xcentre = 0.5
  38         ycentre = 0.5
  39         zcentre = 0.5
  40     else:
  41         xcentre = 0.
  42         ycentre = 0.
  43         zcentre = 0.
  44
  45     pressure_field = cdmath.Field("Pressure",            cdmath.CELLS, my_mesh, 1)
  46     velocity_field = cdmath.Field("Velocity",            cdmath.CELLS, my_mesh, 3)
  47
  48     for i in range(nbCells):
  49         velocity_field[i,0] = 0
  50         velocity_field[i,1] = 0
  51         velocity_field[i,2] = 0
  52
  53         x = my_mesh.getCell(i).x()
  54         valX = (x - xcentre) * (x - xcentre)
  55
  56         if(dim==1):
  57             val =  sqrt(valX)
  58         if(dim==2):
  59             y = my_mesh.getCell(i).y()
  60             valY = (y - ycentre) * (y - ycentre)
  61             val =  sqrt(valX + valY)
  62         if(dim==3):
  63             y = my_mesh.getCell(i).y()
  64             z = my_mesh.getCell(i).z()
  65             valY = (y - ycentre) * (y - ycentre)
  66             valZ = (z - zcentre) * (z - zcentre)
  67             val =  sqrt(valX + valY + valZ)
  68
  69         if val < rayon:
  70             pressure_field[i] = p0
  71             pass
  72         else:
  73             pressure_field[i] = p0/2
  74             pass
  75         pass
  76
  77     return pressure_field, velocity_field
  78
  79 def jacobianMatrices(normal, coeff, signun):
  80     dim=normal.size()
  81     A=cdmath.Matrix(dim+1,dim+1)
  82     absA=cdmath.Matrix(dim+1,dim+1)
  83
  84     for i in range(dim):
  85         A[i+1,0]=normal[i]*coeff
  86         absA[i+1,0]=-signun*A[i+1,0]
  87         A[0,i+1]=c0*c0*normal[i]*coeff
  88         absA[0,i+1]=signun*A[0,i+1]
  89
  90     return (A-absA)/2
  91
  92
  93 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt):
  94     nbCells = my_mesh.getNumberOfCells()
  95     dim=my_mesh.getMeshDimension()
  96     nbComp=dim+1
  97     normal=cdmath.Vector(dim)
  98
  99     implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
 100
 101     idMoinsJacCL=cdmath.Matrix(nbComp)
 102
 103     v0=cdmath.Vector(dim)
 104     for i in range(dim) :
 105         v0[i] = 1.
 106
 107     for j in range(nbCells):#On parcourt les cellules
 108         Cj = my_mesh.getCell(j)
 109         nbFaces = Cj.getNumberOfFaces();
 110
 111         for k in range(nbFaces) :
 112             indexFace = Cj.getFacesId()[k];
 113             Fk = my_mesh.getFace(indexFace);
 114             for i in range(dim) :
 115                 normal[i] = Cj.getNormalVector(k, i);#normale sortante
 116
 117             signun=sign(normal*v0)
 118             Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure(),signun);
 119
 120             cellAutre =-1
 121             if ( not Fk.isBorder()) :
 122                 # hypothese: La cellule d'index indexC1 est la cellule courante index j */
 123                 if (Fk.getCellsId()[0] == j) :
 124                     # hypothese verifiée
 125                     cellAutre = Fk.getCellsId()[1];
 126                 elif(Fk.getCellsId()[1] == j) :
 127                     # hypothese non verifiée
 128                     cellAutre = Fk.getCellsId()[0];
 129                 else :
 130                     raise ValueError("computeFluxes: problem with mesh, unknown cel number")
 131
 132                 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
 133                 implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.))
 134             else  :
 135                 if( Fk.getGroupName() != "Periodic" and Fk.getGroupName() != "Neumann"):#Wall boundary condition unless Periodic/Neumann specified explicitly
 136                     v=cdmath.Vector(dim+1)
 137                     for i in range(dim) :
 138                         v[i+1]=normal[i]
 139                     idMoinsJacCL=v.tensProduct(v)*2
 140
 141                     implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL)
 142
 143                 elif( Fk.getGroupName() == "Periodic"):#Periodic boundary condition
 144                     indexFP=my_mesh.getIndexFacePeriodic(indexFace)
 145                     Fp = my_mesh.getFace(indexFP)
 146                     cellAutre = Fp.getCellsId()[0]
 147
 148                     implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
 149                     implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.))
 150                 elif(Fk.getGroupName() != "Neumann"):#Nothing to do for Neumann boundary condition
 151                     print Fk.getGroupName()
 152                     raise ValueError("computeFluxes: Unknown boundary condition name");
 153
 154     return implMat
 155
 156 def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution):
 157     dim=my_mesh.getMeshDimension()
 158     nbCells = my_mesh.getNumberOfCells()
 159     meshName=my_mesh.getName()
 160
 161     dt = 0.
 162     time = 0.
 163     it=0;
 164     isStationary=False;
 165
 166     nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
 167     iterGMRESMax=50
 168
 169     #iteration vectors
 170     Un=cdmath.Vector(nbCells*(dim+1))
 171     dUn=cdmath.Vector(nbCells*(dim+1))
 172
 173     # Initial conditions #
 174     print("Construction of the initial condition …")
 175     if(dim==1 or filename.find("square")>-1 or filename.find("Square")>-1 or filename.find("cube")>-1 or filename.find("Cube")>-1):
 176         pressure_field, velocity_field = initial_conditions_shock(my_mesh,False)
 177     elif(filename.find("disk")>-1 or filename.find("Disk")>-1):
 178         pressure_field, velocity_field = initial_conditions_shock(my_mesh,True)
 179     else:
 180         print "Mesh name : ", filename
 181         raise ValueError("Mesh name should contain substring square, cube or disk")
 182
 183     for k in range(nbCells):
 184         Un[k*(dim+1)+0] =      pressure_field[k]
 185         Un[k*(dim+1)+1] = rho0*velocity_field[k,0]
 186         if(dim>=2):
 187             Un[k + 2*nbCells] = rho0*velocity_field[k,1] # value on the bottom face
 188             if(dim==3):
 189                 Un[k + 3*nbCells] = rho0*velocity_field[k,2]
 190
 191     #sauvegarde de la donnée initiale
 192     pressure_field.setTime(time,it);
 193     pressure_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure");
 194     velocity_field.setTime(time,it);
 195     velocity_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity");
 196     #Postprocessing : save 2D picture
 197     PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure_initial")
 198     PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity_initial")
 199
 200     dx_min=my_mesh.minRatioVolSurf()
 201
 202     dt = cfl * dx_min / c0
 203
 204     divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
 205
 206     # Add the identity matrix on the diagonal
 207     for j in range(nbCells*(dim+1)):
 208         divMat.addValue(j,j,1.)
 209     LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","LU")
 210
 211     print("Starting computation of the linear wave system with an pseudo staggered scheme …")
 212
 213     # Starting time loop
 214     while (it<ntmax and time <= tmax and not isStationary):
 215         dUn=Un.deepCopy()
 216         LS.setSndMember(Un)
 217         Un=LS.solve();
 218         cvgceLS=LS.getStatus();
 219         iterGMRES=LS.getNumberOfIter();
 220         if(not cvgceLS):
 221             print "Linear system did not converge ", iterGMRES, " GMRES iterations"
 222             raise ValueError("Pas de convergence du système linéaire");
 223         dUn-=Un
 224
 225         maxVector=dUn.maxVector(dim+1)
 226         isStationary= maxVector[0]/p0<precision and maxVector[1]/rho0<precision and maxVector[2]/rho0<precision;
 227         if(dim==3):
 228             isStationary=isStationary and maxVector[3]/rho0<precision
 229         time=time+dt;
 230         it=it+1;
 231
 232         #Sauvegardes
 233         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 234             print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
 235             print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0
 236             print "Linear system converged in ", iterGMRES, " GMRES iterations"
 237
 238             for k in range(nbCells):
 239                 pressure_field[k]=Un[k*(dim+1)+0]
 240                 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
 241                 if(dim>1):
 242                     velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
 243                     if(dim>2):
 244                         velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
 245
 246             pressure_field.setTime(time,it);
 247             pressure_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure",False);
 248             velocity_field.setTime(time,it);
 249             velocity_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity",False);
 250
 251     print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
 252     print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0
 253     print
 254
 255     if(it>=ntmax):
 256         print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
 257     elif(isStationary):
 258         print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
 259         print "------------------------------------------------------------------------------------"
 260
 261         pressure_field.setTime(time,0);
 262         pressure_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure_Stat");
 263         velocity_field.setTime(time,0);
 264         velocity_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity_Stat");
 265
 266         #Postprocessing : save 2D picture
 267         PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure_Stat")
 268         PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity_Stat")
 269
 270     else:
 271         print "Temps maximum Tmax= ", tmax, " atteint"
 272
 273
 274 def solve(my_mesh,meshName,resolution):
 275     print "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension()
 276     print "Numerical method : implicit pseudo staggered"
 277     print "Initial data : spherical wave"
 278     print "Wall boundary conditions"
 279     print "Mesh name : ",meshName , my_mesh.getNumberOfCells(), " cells"
 280
 281     # Problem data
 282     tmax = 1000.
 283     ntmax = 100
 284     cfl = 1./my_mesh.getSpaceDimension()
 285     output_freq = 1
 286
 287     WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution)
 288
 289 def solve_file( filename,meshName, resolution):
 290     my_mesh = cdmath.Mesh(filename+".med")
 291
 292     return solve(my_mesh, filename+str(my_mesh.getNumberOfCells()),resolution)
 293
 294
 295 if __name__ == """__main__""":
 296     if len(sys.argv) >1 :
 297         filename=sys.argv[1]
 298         my_mesh = cdmath.Mesh(filename)
 299         solve(my_mesh,filename,100)
 300     else :
 301         raise ValueError("WaveSystemUpwind.py expects a mesh file name")