CDMATH/tests/examples/EulerSystem_Shock/EulerSystemUpwind/EulerSystemUpwind.py

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