CDMATH/tests/examples/TransportEquation/TransportEquationUpwind.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 from math import sin, cos, pi, sqrt
   5 import time, json
   6 import cdmath
   7 import PV_routines
   8 import VTK_routines
   9
  10 precision=1e-5
  11
  12 def initial_conditions_transport_equation(my_mesh):
  13     print "Initial_data","Shock"
  14     dim     = my_mesh.getMeshDimension()
  15     nbCells = my_mesh.getNumberOfCells()
  16
  17     initial_field = cdmath.Field("unknown", cdmath.CELLS, my_mesh, 1)
  18
  19     rayon = 0.15
  20     xcentre = 0.5
  21     ycentre = 0.5
  22     zcentre = 0.5
  23
  24     for i in range(nbCells):
  25         x = my_mesh.getCell(i).x()
  26
  27         valX = (x - xcentre) * (x - xcentre)
  28
  29         if(dim==1):
  30             val=sqrt(valX)
  31         if(dim==2):
  32             y = my_mesh.getCell(i).y()
  33             valY = (y - ycentre) * (y - ycentre)
  34             val =  sqrt(valX + valY)
  35         if(dim==3):
  36             y = my_mesh.getCell(i).y()
  37             z = my_mesh.getCell(i).z()
  38             valY = (y - ycentre) * (y - ycentre)
  39             valZ = (z - zcentre) * (z - zcentre)
  40             val =  sqrt(valX + valY + valZ)
  41
  42         if val < rayon:
  43             initial_field[i] = 1
  44             pass
  45         else:
  46             initial_field[i] = 0
  47             pass
  48         pass
  49
  50     return initial_field
  51
  52 def upwinding_coeff(normal, coeff, velocity):
  53     dim=normal.size()
  54
  55     if(velocity*normal<0.):
  56         return velocity*normal*coeff
  57     else:
  58         return 0.
  59
  60
  61 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,test_bc,velocity):
  62     nbCells = my_mesh.getNumberOfCells()
  63     dim=my_mesh.getMeshDimension()
  64     nbComp=1
  65     normal=cdmath.Vector(dim)
  66
  67     implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
  68
  69     for j in range(nbCells):#On parcourt les cellules
  70         Cj = my_mesh.getCell(j)
  71         nbFaces = Cj.getNumberOfFaces();
  72
  73         for k in range(nbFaces) :
  74             indexFace = Cj.getFacesId()[k];
  75             Fk = my_mesh.getFace(indexFace);
  76             for i in range(dim) :
  77                 normal[i] = Cj.getNormalVector(k, i);#normale sortante
  78
  79             Am=upwinding_coeff( normal,dt*Fk.getMeasure()/Cj.getMeasure(),velocity);
  80
  81             cellAutre =-1
  82             if ( not Fk.isBorder()) :
  83                 # hypothese: La cellule d'index indexC1 est la cellule courante index j */
  84                 if (Fk.getCellsId()[0] == j) :
  85                     # hypothese verifiée
  86                     cellAutre = Fk.getCellsId()[1];
  87                 elif(Fk.getCellsId()[1] == j) :
  88                     # hypothese non verifiée
  89                     cellAutre = Fk.getCellsId()[0];
  90                 else :
  91                     raise ValueError("computeFluxes: problem with mesh, unknown cel number")
  92
  93                 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
  94                 implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.))
  95             else  :
  96                 if( test_bc=="Periodic" and Fk.getGroupName() != "Neumann"):#Periodic boundary condition unless Wall/Neumann specified explicitly
  97                     indexFP = my_mesh.getIndexFacePeriodic(indexFace, my_mesh.getName()== "squareWithBrickWall", my_mesh.getName()== "squareWithHexagons")
  98                     Fp = my_mesh.getFace(indexFP)
  99                     cellAutre = Fp.getCellsId()[0]
 100
 101                     implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
 102                     implMat.addValue(j*nbComp,        j*nbComp,Am*(-1.))
 103                 elif(test_bc!="Neumann" and Fk.getGroupName() != "Neumann"):#Nothing to do for Neumann boundary condition
 104                     print Fk.getGroupName()
 105                     raise ValueError("computeFluxes: Unknown boundary condition name");
 106
 107     return implMat
 108
 109 def TransportEquationVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution,test_bc,velocity,isImplicit):
 110     dim=my_mesh.getMeshDimension()
 111     nbCells = my_mesh.getNumberOfCells()
 112
 113     dt = 0.
 114     time = 0.
 115     it=0;
 116     isStationary=False
 117
 118     nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
 119
 120     #iteration vectors
 121     Un =cdmath.Vector(nbCells)
 122     dUn=cdmath.Vector(nbCells)
 123
 124     # Initial conditions #
 125     print("Construction of the initial condition …")
 126     stat_field = initial_conditions_transport_equation(my_mesh)
 127     for k in range(nbCells):
 128         Un[k] =     stat_field[k]
 129     unknown_field = initial_conditions_transport_equation(my_mesh)
 130     S = cdmath.Vector(nbCells)#source term is zero
 131
 132     total_mass_initial=unknown_field.integral()#For conservation test later
 133
 134     #sauvegarde de la donnée initiale
 135     unknown_field.writeVTK("TransportEquation"+str(dim)+"DUpwind"+meshName);
 136
 137     dx_min=my_mesh.minRatioVolSurf()
 138
 139     dt = cfl * dx_min / velocity.norm()
 140
 141     divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,test_bc,velocity)
 142     if(isImplicit):
 143         #Adding the identity matrix on the diagonal
 144         divMat.diagonalShift(1)#only after  filling all coefficients
 145         #for j in range(nbCells*(dim+1)):
 146         #    divMat.addValue(j,j,1)
 147
 148         iterGMRESMax=50
 149         LS=cdmath.LinearSolver(divMat,Un+S*dt,iterGMRESMax, precision, "GMRES","ILU")
 150         LS.setComputeConditionNumber()
 151
 152     print("Starting computation of the linear transport equation with an Upwind scheme …")
 153
 154     # Starting time loop
 155     while (it<ntmax and time <= tmax and not isStationary):
 156         if(isImplicit):
 157             dUn=Un.deepCopy()
 158             LS.setSndMember(Un+S*dt)
 159             Un=LS.solve();
 160             if(not LS.getStatus()):
 161                 print "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations"
 162                 raise ValueError("Pas de convergence du système linéaire");
 163             dUn-=Un
 164
 165         else:
 166             dUn=divMat*Un-S*dt
 167             Un-=dUn
 168
 169         maxVector=dUn.maxVector(1)
 170         isStationary= maxVector[0]<precision
 171
 172         time=time+dt;
 173         it=it+1;
 174
 175         #Sauvegardes
 176         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 177             print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
 178             print "Variation temporelle relative : ", maxVector[0]
 179             if(isImplicit):
 180                 print "Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations"
 181
 182             for k in range(nbCells):
 183                 unknown_field[k]  =Un[k]
 184
 185             unknown_field.setTime(time,it);
 186             unknown_field.writeVTK("TransportEquation"+str(dim)+"DUpwind"+meshName,False);
 187
 188     print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
 189     print "Variation temporelle relative : ", maxVector[0]
 190
 191     if(it>=ntmax):
 192         print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
 193         raise ValueError("Maximum number of time steps reached : Stationary state not found !!!!!!!")
 194     elif(isStationary):
 195         print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
 196         if(test_bc=="Periodic"):
 197             print "Mass loss: ", (total_mass_initial-unknown_field.integral()).norm(), " precision required= ", precision
 198             assert (total_mass_initial-unknown_field.integral()).norm()<precision
 199         print "------------------------------------------------------------------------------------"
 200
 201         unknown_field.setTime(time,0);
 202         unknown_field.writeVTK("TransportEquation"+str(dim)+"DUpwind"+meshName+"_Stat");
 203
 204         #Postprocessing : Extraction of the diagonal data
 205         if(dim==2):
 206             diag_data_u=VTK_routines.Extract_field_data_over_line_to_numpyArray(unknown_field,[0,1,0],[1,0,0], resolution)
 207         elif(dim==3):
 208             diag_data_u=VTK_routines.Extract_field_data_over_line_to_numpyArray(unknown_field,[0,0,0],[1,1,1], resolution)
 209         #Postprocessing : save 2D picture
 210         PV_routines.Save_PV_data_to_picture_file("TransportEquation"+str(dim)+"DUpwind"+meshName+"_Stat"+'_0.vtu',"unknown",'CELLS',"TransportEquation"+str(dim)+"DUpwind"+meshName+"_Stat")
 211
 212         return nbCells, time, it, unknown_field.getNormEuclidean().max(), diag_data_u
 213     else:
 214         print "Temps maximum Tmax= ", tmax, " atteint"
 215         raise ValueError("Maximum time reached : Stationary state not found !!!!!!!")
 216
 217
 218 def solve(my_mesh, meshName, resolution, meshType, cfl, test_bc):
 219     print "Resolution of the Transport Equation in dimension ", my_mesh.getMeshDimension()
 220     print "Numerical method : ", "Upwind"
 221     print "Initial data : ", "Spherical shock"
 222     print "Mesh name : ",meshName , ", ", my_mesh.getNumberOfCells(), " cells"
 223
 224     # Problem data
 225     tmax = 10000.
 226     ntmax = 30000
 227     output_freq = 1000
 228
 229     isImplicit=True
 230
 231     dim=my_mesh.getMeshDimension()
 232     velocity=cdmath.Vector(dim)
 233     for i in range(dim) :
 234         velocity[i] = 1
 235
 236     nbCells, t_final, ndt_final, max_unknown, diag_data_u = TransportEquationVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution, test_bc, velocity,isImplicit)
 237
 238     return nbCells, t_final, ndt_final, max_unknown, diag_data_u
 239
 240 def solve_file( filename,meshName, resolution,meshType, cfl, test_bc):
 241     my_mesh = cdmath.Mesh(filename+".med")
 242
 243     return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution, meshType, cfl, test_bc)
 244
 245
 246 if __name__ == """__main__""":
 247     M1=cdmath.Mesh(0.,1.,20,0.,1.,20,0)
 248     cfl=0.5
 249     solve(M1,"SquareRegularTriangles",100,"Regular triangles",cfl,"Periodic")
 250