CDMATH/tests/examples/WaveSystem_Stationary/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 : Test de préservation d'un état stationnaire
  11 #               Utilisation du schéma centré (implicite) sur un maillage général
  12 #               Initialisation par un vortex stationnaire
  13 #               Conditions aux limites parois
  14 #                       Création et sauvegarde du champ résultant et des figures
  15 #================================================================================================================================
  16
  17
  18 from math import sin, cos, pi, 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_disk_vortex(my_mesh):
  31     print( "Disk vortex initial data")
  32     dim     = my_mesh.getMeshDimension()
  33     nbCells = my_mesh.getNumberOfCells()
  34
  35     if(dim!=2):
  36         raise ValueError("Wave system on disk : mesh dimension should be 2")
  37
  38     pressure_field = cdmath.Field("Pressure",            cdmath.CELLS, my_mesh, 1)
  39     velocity_field = cdmath.Field("Velocity",            cdmath.CELLS, my_mesh, 3)
  40
  41     for i in range(nbCells):
  42         x = my_mesh.getCell(i).x()
  43         y = my_mesh.getCell(i).y()
  44
  45         pressure_field[i] = p0
  46
  47         velocity_field[i,0] = -y
  48         velocity_field[i,1] =  x
  49         velocity_field[i,2] = 0
  50
  51     return pressure_field, velocity_field
  52
  53 def initial_conditions_square_vortex(my_mesh):
  54     print( "Initial data : Square vortex (Constant pressure, divergence free velocity)")
  55     dim     = my_mesh.getMeshDimension()
  56     nbCells = my_mesh.getNumberOfCells()
  57
  58     pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
  59     velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
  60
  61     for i in range(nbCells):
  62         x = my_mesh.getCell(i).x()
  63         y = my_mesh.getCell(i).y()
  64         z = my_mesh.getCell(i).z()
  65
  66         pressure_field[i] = p0
  67         if(dim==1):
  68             velocity_field[i,0] = 1
  69             velocity_field[i,1] = 0
  70             velocity_field[i,2] = 0
  71         elif(dim==2):
  72             velocity_field[i,0] =  sin(pi*x)*cos(pi*y)
  73             velocity_field[i,1] = -sin(pi*y)*cos(pi*x)
  74             velocity_field[i,2] = 0
  75         elif(dim==3):
  76             velocity_field[i,0] =    sin(pi*x)*cos(pi*y)*cos(pi*z)
  77             velocity_field[i,1] =    sin(pi*y)*cos(pi*x)*cos(pi*z)
  78             velocity_field[i,2] = -2*sin(pi*z)*cos(pi*x)*cos(pi*y)
  79
  80     return pressure_field, velocity_field
  81
  82 def jacobianMatrices(normal, coeff, signun):
  83     dim=normal.size()
  84     A=cdmath.Matrix(dim+1,dim+1)
  85
  86     for i in range(dim):
  87         A[i+1,0]=normal[i]*coeff
  88         A[0,i+1]=c0*c0*normal[i]*coeff
  89
  90     return A/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(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_square_vortex(my_mesh)
 177     elif(filename.find("disk")>-1 or filename.find("Disk")>-1):
 178         pressure_field, velocity_field = initial_conditions_disk_vortex(my_mesh)
 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         Un[k*(dim+1)+2] = rho0*velocity_field[k,1]
 187         if(dim==3):
 188             Un[k*(dim+1)+3] = rho0*velocity_field[k,2]
 189
 190     #sauvegarde de la donnée initiale
 191     pressure_field.setTime(time,it);
 192     pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure");
 193     velocity_field.setTime(time,it);
 194     velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity");
 195     #Postprocessing : save 2D picture
 196     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")
 197     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")
 198
 199     dx_min=my_mesh.minRatioVolSurf()
 200
 201     dt = cfl * dx_min / c0
 202
 203     divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
 204
 205     # Add the identity matrix on the diagonal
 206     for j in range(nbCells*(dim+1)):
 207         divMat.addValue(j,j,1.)
 208     LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","LU")
 209
 210     print("Starting computation of the linear wave system with a centered scheme …")
 211
 212     # Starting time loop
 213     while (it<ntmax and time <= tmax and not isStationary):
 214         dUn=Un.deepCopy()
 215         LS.setSndMember(Un)
 216         Un=LS.solve();
 217         cvgceLS=LS.getStatus();
 218         iterGMRES=LS.getNumberOfIter();
 219         if(not cvgceLS):
 220             print( "Linear system did not converge ", iterGMRES, " GMRES iterations")
 221             raise ValueError("Pas de convergence du système linéaire");
 222         dUn-=Un
 223
 224         maxVector=dUn.maxVector(dim+1)
 225         isStationary= maxVector[0]/p0<precision and maxVector[1]/rho0<precision and maxVector[2]/rho0<precision;
 226         if(dim==3):
 227             isStationary=isStationary and maxVector[3]/rho0<precision
 228         time=time+dt;
 229         it=it+1;
 230
 231         #Sauvegardes
 232         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 233             print( "-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
 234             print( "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0 )
 235             print( "Linear system converged in ", iterGMRES, " GMRES iterations")
 236
 237             for k in range(nbCells):
 238                 pressure_field[k]=Un[k*(dim+1)+0]
 239                 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
 240                 if(dim>1):
 241                     velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
 242                     if(dim>2):
 243                         velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
 244
 245             pressure_field.setTime(time,it);
 246             pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure",False);
 247             velocity_field.setTime(time,it);
 248             velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity",False);
 249
 250     print( "-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
 251     print( "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0 )
 252     print()
 253
 254     if(it>=ntmax):
 255         print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
 256     elif(isStationary):
 257         print( "Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
 258         print( "------------------------------------------------------------------------------------")
 259
 260         pressure_field.setTime(time,0);
 261         pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_Stat");
 262         velocity_field.setTime(time,0);
 263         velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_Stat");
 264
 265         #Postprocessing : save 2D picture
 266         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")
 267         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")
 268
 269     else:
 270         print( "Temps maximum Tmax= ", tmax, " atteint")
 271
 272
 273 def solve(my_mesh,meshName,resolution):
 274     print( "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension() )
 275     print( "Numerical method : implicit centered" )
 276     print( "Initial data : stationary solution (constant pressure, divergence free velocity)" )
 277     print( "Wall boundary conditions" )
 278     print( "Mesh name : ",meshName , my_mesh.getNumberOfCells(), " cells" )
 279
 280     # Problem data
 281     tmax = 1000.
 282     ntmax = 100
 283     cfl = 1./my_mesh.getSpaceDimension()
 284     output_freq = 100
 285
 286     WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution)
 287
 288 def solve_file( filename,meshName, resolution):
 289     my_mesh = cdmath.Mesh(filename+".med")
 290
 291     return solve(my_mesh, filename+str(my_mesh.getNumberOfCells()),resolution)
 292
 293
 294 if __name__ == """__main__""":
 295     if len(sys.argv) >1 :
 296         filename=sys.argv[1]
 297         my_mesh = cdmath.Mesh(filename)
 298         solve(my_mesh,filename,100)
 299     else :
 300         raise ValueError("WaveSystemUpwind.py expects a mesh file name")