CDMATH/tests/examples/WaveSystem_Stationary/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 : Test de préservation d'un état stationnaire
  11 #               Utilisation du schéma upwind explicite ou implicite sur un maillage général
  12 #               Initialisation par une 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 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_disk_vortex(my_mesh):
  30     print( "Disk vortex initial data")
  31     dim     = my_mesh.getMeshDimension()
  32     nbCells = my_mesh.getNumberOfCells()
  33
  34     if(dim!=2):
  35         raise ValueError("Wave system on disk : mesh dimension should be 2")
  36
  37     pressure_field = cdmath.Field("Pressure",            cdmath.CELLS, my_mesh, 1)
  38     velocity_field = cdmath.Field("Velocity",            cdmath.CELLS, my_mesh, 3)
  39
  40     for i in range(nbCells):
  41         x = my_mesh.getCell(i).x()
  42         y = my_mesh.getCell(i).y()
  43
  44         pressure_field[i] = p0
  45
  46         velocity_field[i,0] = -y
  47         velocity_field[i,1] =  x
  48         velocity_field[i,2] = 0
  49
  50     return pressure_field, velocity_field
  51
  52 def initial_conditions_square_vortex(my_mesh):
  53     print( "Initial data : Square vortex (Constant pressure, divergence free velocity)" )
  54     dim     = my_mesh.getMeshDimension()
  55     nbCells = my_mesh.getNumberOfCells()
  56
  57     pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
  58     velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
  59
  60     for i in range(nbCells):
  61         x = my_mesh.getCell(i).x()
  62         y = my_mesh.getCell(i).y()
  63         z = my_mesh.getCell(i).z()
  64
  65         pressure_field[i] = p0
  66         if(dim==1):
  67             velocity_field[i,0] = 1
  68             velocity_field[i,1] = 0
  69             velocity_field[i,2] = 0
  70         elif(dim==2):
  71             velocity_field[i,0] =  sin(pi*x)*cos(pi*y)
  72             velocity_field[i,1] = -sin(pi*y)*cos(pi*x)
  73             velocity_field[i,2] = 0
  74         elif(dim==3):
  75             velocity_field[i,0] =    sin(pi*x)*cos(pi*y)*cos(pi*z)
  76             velocity_field[i,1] =    sin(pi*y)*cos(pi*x)*cos(pi*z)
  77             velocity_field[i,2] = -2*sin(pi*z)*cos(pi*x)*cos(pi*y)
  78
  79     return pressure_field, velocity_field
  80
  81 def jacobianMatrices(normal,coeff):
  82     dim=normal.size()
  83     A=cdmath.Matrix(dim+1,dim+1)
  84     absA=cdmath.Matrix(dim+1,dim+1)
  85
  86     absA[0,0]=c0*coeff
  87     for i in range(dim):
  88         A[i+1,0]=      normal[i]*coeff
  89         A[0,i+1]=c0*c0*normal[i]*coeff
  90         for j in range(dim):
  91             absA[i+1,j+1]=c0*normal[i]*normal[j]*coeff
  92
  93     return (A - absA)/2
  94
  95 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt):
  96     nbCells = my_mesh.getNumberOfCells()
  97     dim=my_mesh.getMeshDimension()
  98     nbComp=dim+1
  99     normal=cdmath.Vector(dim)
 100
 101     implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
 102
 103     idMoinsJacCL=cdmath.Matrix(nbComp)
 104
 105     for j in range(nbCells):#On parcourt les cellules
 106         Cj = my_mesh.getCell(j)
 107         nbFaces = Cj.getNumberOfFaces();
 108
 109         for k in range(nbFaces) :
 110             indexFace = Cj.getFacesId()[k];
 111             Fk = my_mesh.getFace(indexFace);
 112             for i in range(dim) :
 113                 normal[i] = Cj.getNormalVector(k, i);#normale sortante
 114
 115             Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
 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 cell 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, isImplicit):
 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     nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
 163
 164     # Initial conditions #
 165     print("Construction of the initial condition …")
 166     if(filename.find("square")>-1 or filename.find("Square")>-1 or filename.find("cube")>-1 or filename.find("Cube")>-1):
 167         pressure_field, velocity_field = initial_conditions_square_vortex(my_mesh)
 168     elif(filename.find("disk")>-1 or filename.find("Disk")>-1):
 169         pressure_field, velocity_field = initial_conditions_disk_vortex(my_mesh)
 170     else:
 171         print( "Mesh name : ", filename)
 172         raise ValueError("Mesh name should contain substring square, cube or disk")
 173
 174     #iteration vectors
 175     Un =cdmath.Vector(nbCells*(dim+1))
 176     dUn=cdmath.Vector(nbCells*(dim+1))
 177
 178     for k in range(nbCells):
 179         Un[k*(dim+1)+0] =     pressure_field[k]
 180         Un[k*(dim+1)+1] =rho0*velocity_field[k,0]
 181         Un[k*(dim+1)+2] =rho0*velocity_field[k,1]
 182         if(dim==3):
 183             Un[k*(dim+1)+3] =rho0*velocity_field[k,2]
 184
 185     #sauvegarde de la donnée initiale
 186     pressure_field.setTime(time,it);
 187     pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure");
 188     velocity_field.setTime(time,it);
 189     velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity");
 190
 191     dx_min=my_mesh.minRatioVolSurf()
 192
 193     dt = cfl * dx_min / c0
 194
 195     divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
 196     if( isImplicit):
 197         #Adding the identity matrix on the diagonal
 198         divMat.diagonalShift(1)#only after  filling all coefficients
 199
 200         iterGMRESMax=50
 201         LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","ILU")
 202
 203         LS.setComputeConditionNumber()
 204
 205     print("Starting computation of the linear wave system with an UPWIND scheme …")
 206
 207     # Starting time loop
 208     while (it<ntmax and time <= tmax and not isStationary):
 209         if(isImplicit):
 210             dUn=Un.deepCopy()
 211             LS.setSndMember(Un)
 212             Un=LS.solve();
 213             if(not LS.getStatus()):
 214                 print( "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations")
 215                 raise ValueError("Pas de convergence du système linéaire");
 216             dUn-=Un
 217
 218         else:
 219             dUn=divMat*Un
 220             Un-=dUn
 221
 222         time=time+dt;
 223         it=it+1;
 224
 225         #Sauvegardes
 226         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 227             print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 228
 229             for k in range(nbCells):
 230                 pressure_field[k]  =Un[k*(dim+1)+0]
 231                 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
 232                 if(dim>1):
 233                     velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
 234                     if(dim>2):
 235                         velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
 236
 237             pressure_field.setTime(time,it);
 238             pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure",False);
 239             velocity_field.setTime(time,it);
 240             velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity",False);
 241
 242     print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 243     print()
 244
 245     if(it>=ntmax):
 246         print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
 247     elif(isStationary):
 248         print( "Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
 249
 250         pressure_field.setTime(time,0);
 251         pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat");
 252         velocity_field.setTime(time,0);
 253         velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat");
 254
 255         #Postprocessing : Extraction of the diagonal data
 256         diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,1,0],[1,0,0], resolution)
 257         diag_data_vel  =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,1,0],[1,0,0], resolution)
 258         #Postprocessing : save 2D picture
 259         PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat")
 260         PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat")
 261
 262     else:
 263         print( "Temps maximum Tmax= ", tmax, " atteint")
 264
 265
 266 def solve(my_mesh,filename,resolution, isImplicit):
 267     print( "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension())
 268     print( "Numerical method : upwind")
 269     print( "Initial data : stationary solution (constant pressure, divergence free velocity)")
 270     print( "Wall boundary conditions")
 271     print( "Mesh name : ",filename , my_mesh.getNumberOfCells(), " cells")
 272
 273     # Problem data
 274     tmax = 1.
 275     ntmax = 100
 276     cfl = 1./my_mesh.getSpaceDimension()
 277     output_freq = 100
 278
 279     WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution, isImplicit)
 280
 281 def solve_file( filename,resolution):
 282     my_mesh = cdmath.Mesh(filename+".med")
 283     solve(my_mesh, filename,resolution, isImplicit)
 284
 285 if __name__ == """__main__""":
 286     if len(sys.argv) >2 :
 287         filename=sys.argv[1]
 288         isImplicit=bool(int(sys.argv[2]))
 289         my_mesh = cdmath.Mesh(filename)
 290         solve(my_mesh,filename,100, isImplicit)
 291     else :
 292         raise ValueError("WaveSystemUpwind.py expects a mesh file name and a boolean (isImplicit)")