CDMATH/tests/examples/WaveSystem_Shock/WaveSystemCentered/WaveSystemCentered.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 centré (implicite) sur un maillage général
  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
  83     for i in range(dim):
  84         A[i+1,0]=normal[i]*coeff
  85         A[0,i+1]=c0*c0*normal[i]*coeff
  86
  87     return A/2
  88
  89
  90 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt):
  91     nbCells = my_mesh.getNumberOfCells()
  92     dim=my_mesh.getMeshDimension()
  93     nbComp=dim+1
  94     normal=cdmath.Vector(dim)
  95
  96     implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
  97
  98     idMoinsJacCL=cdmath.Matrix(nbComp)
  99
 100     v0=cdmath.Vector(dim)
 101     for i in range(dim) :
 102         v0[i] = 1.
 103
 104     for j in range(nbCells):#On parcourt les cellules
 105         Cj = my_mesh.getCell(j)
 106         nbFaces = Cj.getNumberOfFaces();
 107
 108         for k in range(nbFaces) :
 109             indexFace = Cj.getFacesId()[k];
 110             Fk = my_mesh.getFace(indexFace);
 111             for i in range(dim) :
 112                 normal[i] = Cj.getNormalVector(k, i);#normale sortante
 113
 114             signun=sign(normal*v0)
 115             Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure(),signun);
 116
 117             cellAutre =-1
 118             if ( not Fk.isBorder()) :
 119                 # hypothese: La cellule d'index indexC1 est la cellule courante index j */
 120                 if (Fk.getCellsId()[0] == j) :
 121                     # hypothese verifiée
 122                     cellAutre = Fk.getCellsId()[1];
 123                 elif(Fk.getCellsId()[1] == j) :
 124                     # hypothese non verifiée
 125                     cellAutre = Fk.getCellsId()[0];
 126                 else :
 127                     raise ValueError("computeFluxes: problem with mesh, unknown cel number")
 128
 129                 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
 130                 implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.))
 131             else  :
 132                 if( Fk.getGroupName() != "Periodic" and Fk.getGroupName() != "Neumann"):#Wall boundary condition unless Periodic/Neumann specified explicitly
 133                     v=cdmath.Vector(dim+1)
 134                     for i in range(dim) :
 135                         v[i+1]=normal[i]
 136                     idMoinsJacCL=v.tensProduct(v)*2
 137
 138                     implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL)
 139
 140                 elif( Fk.getGroupName() == "Periodic"):#Periodic boundary condition
 141                     indexFP=my_mesh.getIndexFacePeriodic(indexFace)
 142                     Fp = my_mesh.getFace(indexFP)
 143                     cellAutre = Fp.getCellsId()[0]
 144
 145                     implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
 146                     implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.))
 147                 elif(Fk.getGroupName() != "Neumann"):#Nothing to do for Neumann boundary condition
 148                     print Fk.getGroupName()
 149                     raise ValueError("computeFluxes: Unknown boundary condition name");
 150
 151     return implMat
 152
 153 def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution):
 154     dim=my_mesh.getMeshDimension()
 155     nbCells = my_mesh.getNumberOfCells()
 156     meshName=my_mesh.getName()
 157
 158     dt = 0.
 159     time = 0.
 160     it=0;
 161     isStationary=False;
 162
 163     nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
 164     iterGMRESMax=50
 165
 166     #iteration vectors
 167     Un=cdmath.Vector(nbCells*(dim+1))
 168     dUn=cdmath.Vector(nbCells*(dim+1))
 169
 170     # Initial conditions #
 171     print("Construction of the initial condition …")
 172     if(dim==1 or filename.find("square")>-1 or filename.find("Square")>-1 or filename.find("cube")>-1 or filename.find("Cube")>-1):
 173         pressure_field, velocity_field = initial_conditions_shock(my_mesh,False)
 174     elif(filename.find("disk")>-1 or filename.find("Disk")>-1):
 175         pressure_field, velocity_field = initial_conditions_shock(my_mesh,True)
 176     else:
 177         print "Mesh name : ", filename
 178         raise ValueError("Mesh name should contain substring square, cube or disk")
 179
 180     for k in range(nbCells):
 181         Un[k*(dim+1)+0] =      pressure_field[k]
 182         Un[k*(dim+1)+1] = rho0*velocity_field[k,0]
 183         if(dim>=2):
 184             Un[k + 2*nbCells] = rho0*velocity_field[k,1] # value on the bottom face
 185             if(dim==3):
 186                 Un[k + 3*nbCells] = rho0*velocity_field[k,2]
 187
 188     #sauvegarde de la donnée initiale
 189     pressure_field.setTime(time,it);
 190     pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure");
 191     velocity_field.setTime(time,it);
 192     velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity");
 193     #Postprocessing : save 2D picture
 194     PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_initial")
 195     PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_initial")
 196
 197     dx_min=my_mesh.minRatioVolSurf()
 198
 199     dt = cfl * dx_min / c0
 200
 201     divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
 202
 203     # Add the identity matrix on the diagonal
 204     for j in range(nbCells*(dim+1)):
 205         divMat.addValue(j,j,1.)
 206     LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","LU")
 207
 208     print("Starting computation of the linear wave system with a centered scheme …")
 209
 210     # Starting time loop
 211     while (it<ntmax and time <= tmax and not isStationary):
 212         dUn=Un.deepCopy()
 213         LS.setSndMember(Un)
 214         Un=LS.solve();
 215         cvgceLS=LS.getStatus();
 216         iterGMRES=LS.getNumberOfIter();
 217         if(not cvgceLS):
 218             print "Linear system did not converge ", iterGMRES, " GMRES iterations"
 219             raise ValueError("Pas de convergence du système linéaire");
 220         dUn-=Un
 221
 222         maxVector=dUn.maxVector(dim+1)
 223         isStationary= maxVector[0]/p0<precision and maxVector[1]/rho0<precision and maxVector[2]/rho0<precision;
 224         if(dim==3):
 225             isStationary=isStationary and maxVector[3]/rho0<precision
 226         time=time+dt;
 227         it=it+1;
 228
 229         #Sauvegardes
 230         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 231             print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
 232             print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0
 233             print "Linear system converged in ", iterGMRES, " GMRES iterations"
 234
 235             for k in range(nbCells):
 236                 pressure_field[k]=Un[k*(dim+1)+0]
 237                 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
 238                 if(dim>1):
 239                     velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
 240                     if(dim>2):
 241                         velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
 242
 243             pressure_field.setTime(time,it);
 244             pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure",False);
 245             velocity_field.setTime(time,it);
 246             velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity",False);
 247
 248     print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
 249     print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0
 250     print
 251
 252     if(it>=ntmax):
 253         print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
 254     elif(isStationary):
 255         print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
 256         print "------------------------------------------------------------------------------------"
 257
 258         pressure_field.setTime(time,0);
 259         pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_Stat");
 260         velocity_field.setTime(time,0);
 261         velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_Stat");
 262
 263         #Postprocessing : save 2D picture
 264         PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_Stat")
 265         PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_Stat")
 266
 267     else:
 268         print "Temps maximum Tmax= ", tmax, " atteint"
 269
 270
 271 def solve(my_mesh,meshName,resolution):
 272     print "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension()
 273     print "Numerical method : implicit centered"
 274     print "Initial data : spherical wave"
 275     print "Wall boundary conditions"
 276     print "Mesh name : ",meshName , my_mesh.getNumberOfCells(), " cells"
 277
 278     # Problem data
 279     tmax = 1000.
 280     ntmax = 100
 281     cfl = 1./my_mesh.getSpaceDimension()
 282     output_freq = 1
 283
 284     WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution)
 285
 286 def solve_file( filename,meshName, resolution):
 287     my_mesh = cdmath.Mesh(filename+".med")
 288
 289     return solve(my_mesh, filename+str(my_mesh.getNumberOfCells()),resolution)
 290
 291
 292 if __name__ == """__main__""":
 293     if len(sys.argv) >1 :
 294         filename=sys.argv[1]
 295         my_mesh = cdmath.Mesh(filename)
 296         solve(my_mesh,filename,100)
 297     else :
 298         raise ValueError("WaveSystemUpwind.py expects a mesh file name")