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