CDMATH/tests/examples/WaveSystem2DUpwind_RiemannProblem/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 droite
  11 #               Utilisation du schéma upwind explicite ou implicite sur un maillage général
  12 #               Initialisation par une discontinuité verticale
  13 #               Conditions aux limites de Neumann
  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_RiemannProblem(my_mesh):
  30     print( "Initial data : Riemann problem" )
  31     dim     = my_mesh.getMeshDimension()
  32     nbCells = my_mesh.getNumberOfCells()
  33
  34     xcentre = 0
  35
  36     pressure_field = cdmath.Field("Pressure",            cdmath.CELLS, my_mesh, 1)
  37     velocity_field = cdmath.Field("Velocity",            cdmath.CELLS, my_mesh, 3)
  38
  39     for i in range(nbCells):
  40         x = my_mesh.getCell(i).x()
  41
  42         velocity_field[i,0] = 0
  43         velocity_field[i,1] = 0
  44         velocity_field[i,2] = 0
  45
  46         if x < xcentre:
  47             pressure_field[i] = p0
  48             pass
  49         else:
  50             pressure_field[i] = p0/2
  51             pass
  52         pass
  53
  54     return pressure_field, velocity_field
  55
  56 def jacobianMatrices(normal,coeff):
  57     dim=normal.size()
  58     A=cdmath.Matrix(dim+1,dim+1)
  59     absA=cdmath.Matrix(dim+1,dim+1)
  60
  61     absA[0,0]=c0*coeff
  62     for i in range(dim):
  63         A[i+1,0]=      normal[i]*coeff
  64         A[0,i+1]=c0*c0*normal[i]*coeff
  65         for j in range(dim):
  66             absA[i+1,j+1]=c0*normal[i]*normal[j]*coeff
  67
  68     return (A - absA)/2
  69
  70 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt):
  71     nbCells = my_mesh.getNumberOfCells()
  72     dim=my_mesh.getMeshDimension()
  73     nbComp=dim+1
  74     normal=cdmath.Vector(dim)
  75
  76     implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
  77
  78     idMoinsJacCL=cdmath.Matrix(nbComp)
  79
  80     for j in range(nbCells):#On parcourt les cellules
  81         Cj = my_mesh.getCell(j)
  82         nbFaces = Cj.getNumberOfFaces();
  83
  84         for k in range(nbFaces) :
  85             indexFace = Cj.getFacesId()[k];
  86             Fk = my_mesh.getFace(indexFace);
  87             for i in range(dim) :
  88                 normal[i] = Cj.getNormalVector(k, i);#normale sortante
  89
  90             Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
  91
  92             cellAutre =-1
  93             if ( not Fk.isBorder()) :
  94                 # hypothese: La cellule d'index indexC1 est la cellule courante index j */
  95                 if (Fk.getCellsId()[0] == j) :
  96                     # hypothese verifiée
  97                     cellAutre = Fk.getCellsId()[1];
  98                 elif(Fk.getCellsId()[1] == j) :
  99                     # hypothese non verifiée
 100                     cellAutre = Fk.getCellsId()[0];
 101                 else :
 102                     raise ValueError("computeFluxes: problem with mesh, unknown cell number")
 103
 104                 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
 105                 implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.))
 106
 107     return implMat
 108
 109 def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution, isImplicit):
 110     dim=my_mesh.getMeshDimension()
 111     nbCells = my_mesh.getNumberOfCells()
 112     meshName=my_mesh.getName()
 113
 114     dt = 0.
 115     time = 0.
 116     it=0;
 117     isStationary=False;
 118     nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
 119
 120     SumFluxes = cdmath.Field("Fluxes", cdmath.CELLS, my_mesh, dim+1)
 121
 122     # Initial conditions #
 123     print("Construction of the initial condition …")
 124     pressure_field, velocity_field = initial_conditions_RiemannProblem(my_mesh)
 125
 126     #iteration vectors
 127     Un =cdmath.Vector(nbCells*(dim+1))
 128     dUn=cdmath.Vector(nbCells*(dim+1))
 129
 130     for k in range(nbCells):
 131         Un[k*(dim+1)+0] =     pressure_field[k]
 132         Un[k*(dim+1)+1] =rho0*velocity_field[k,0]
 133         Un[k*(dim+1)+2] =rho0*velocity_field[k,1]
 134         if(dim==3):
 135             Un[k*(dim+1)+3] =rho0*velocity_field[k,2]
 136
 137     #sauvegarde de la donnée initiale
 138     pressure_field.setTime(time,it);
 139     pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_pressure");
 140     velocity_field.setTime(time,it);
 141     velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_velocity");
 142
 143     dx_min=my_mesh.minRatioVolSurf()
 144
 145     dt = cfl * dx_min / c0
 146
 147     divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
 148     if( isImplicit):
 149         #Adding the identity matrix on the diagonal
 150         divMat.diagonalShift(1)#only after  filling all coefficients
 151
 152         iterGMRESMax=50
 153         LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","ILU")
 154
 155         LS.setComputeConditionNumber()
 156
 157     print("Starting computation of the linear wave system with an explicit UPWIND scheme …")
 158
 159     # Starting time loop
 160     while (it<ntmax and time <= tmax and not isStationary):
 161         if(isImplicit):
 162             dUn=Un.deepCopy()
 163             LS.setSndMember(Un)
 164             Un=LS.solve();
 165             if(not LS.getStatus()):
 166                 print( "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations" )
 167                 raise ValueError("Pas de convergence du système linéaire");
 168             dUn-=Un
 169
 170         else:
 171             dUn=divMat*Un
 172             Un-=dUn
 173
 174         time=time+dt;
 175         it=it+1;
 176
 177          #Sauvegardes
 178         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 179             print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
 180
 181             for k in range(nbCells):
 182                 pressure_field[k]  =Un[k*(dim+1)+0]
 183                 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
 184                 if(dim>1):
 185                     velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
 186                     if(dim>2):
 187                         velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
 188
 189             pressure_field.setTime(time,it);
 190             pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_pressure",False);
 191             velocity_field.setTime(time,it);
 192             velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_velocity",False);
 193
 194     print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 195     print()
 196
 197     if(it>=ntmax):
 198         print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
 199     elif(isStationary):
 200         print( "Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
 201
 202         pressure_field.setTime(time,0);
 203         pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_pressure_Stat");
 204         velocity_field.setTime(time,0);
 205         velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_velocity_Stat");
 206
 207         #Postprocessing : Extraction of the diagonal data
 208         diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,1,0],[1,0,0], resolution)
 209         diag_data_vel  =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,1,0],[1,0,0], resolution)
 210         #Postprocessing : save 2D picture
 211         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")
 212         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")
 213
 214     else:
 215         print( "Temps maximum Tmax= ", tmax, " atteint")
 216
 217
 218 def solve(my_mesh,filename,resolution, isImplicit):
 219     print( "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension() )
 220     print( "Numerical method : upwind")
 221     print( "Initial data : single straight discontinuity (Riemann problem)")
 222     print( "Neumann boundary conditions")
 223     print( "Mesh name : ",filename , my_mesh.getNumberOfCells(), " cells")
 224
 225     # Problem data
 226     tmax = 1.
 227     ntmax = 50
 228     cfl = 1./my_mesh.getSpaceDimension()
 229     output_freq = 1
 230
 231     WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution, isImplicit)
 232
 233 def solve_file( filename,resolution, isImplicit):
 234     my_mesh = cdmath.Mesh(filename+".med")
 235     solve(my_mesh, filename,resolution, isImplicit)
 236
 237 if __name__ == """__main__""":
 238     if len(sys.argv) >2 :
 239         filename=sys.argv[1]
 240         isImplicit=bool(int(sys.argv[2]))
 241         my_mesh = cdmath.Mesh(filename)
 242         solve(my_mesh,filename,100, isImplicit)
 243     else :
 244         raise ValueError("WaveSystemUpwind.py expects a mesh file name and a boolean (isImplicit)")