CDMATH/tests/examples/transport2d_s/main.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 import math
   5
   6 import cdmath
   7
   8
   9 def initial_conditions(my_mesh):
  10     rayon = 0.15
  11     xcentre = 0.25
  12     ycentre = 0.25
  13     y_field = cdmath.Field("Y field", cdmath.CELLS, my_mesh, 1)
  14     nbCells = my_mesh.getNumberOfCells()
  15     for j in range(nbCells):
  16         x = my_mesh.getCell(j).x()
  17         y = my_mesh.getCell(j).y()
  18         valX = (x - xcentre) * (x - xcentre)
  19         valY = (y - ycentre) * (y - ycentre)
  20         val = math.sqrt(valX + valY)
  21         if val < rayon:
  22             y_field[j] = 1.0
  23             pass
  24         else:
  25             y_field[j] = 0.0
  26             pass
  27         pass
  28     return y_field
  29
  30
  31 def sigma_flux(VitesseX, VitesseY, cfl, y_field, indexFacesPerio):
  32     # Fluxes calculation #
  33     SumFlux = cdmath.Field("Fluxes", cdmath.CELLS, y_field.getMesh(), 1)
  34     my_mesh = y_field.getMesh()
  35     nbCells = my_mesh.getNumberOfCells()
  36     normU = math.sqrt(VitesseX * VitesseX + VitesseY * VitesseY)
  37     for j in range(nbCells):
  38         Cj = my_mesh.getCell(j)
  39         nbFace = Cj.getNumberOfFaces()
  40         SumF = 0.0
  41         minlengthFk = 1.E30
  42         for k in range(nbFace):
  43             indexFace = Cj.getFacesId()[k]
  44             Fk = my_mesh.getFace(indexFace)
  45             NormalX = Cj.getNormalVector(k, 0)
  46             NormalY = Cj.getNormalVector(k, 1)
  47             LengthFk = Fk.getMeasure()
  48             UN = VitesseX * NormalX + VitesseY * NormalY
  49             minlengthFk = min(minlengthFk, LengthFk / abs(UN))
  50             minlengthFk = min(minlengthFk, LengthFk / abs(VitesseX))
  51             minlengthFk = min(minlengthFk, LengthFk / abs(VitesseY))
  52             conc = 0.0
  53             cellCourante = j
  54             cellAutre = -1
  55             if (not Fk.isBorder()):
  56                 indexC1 = Fk.getCellsId()[0]
  57                 indexC2 = Fk.getCellsId()[1]
  58                 # hypothese: La cellule d'index indexC1 est la cellule courante index j #
  59                 if (indexC1 == j):
  60                     # hypothese verifie #
  61                     cellCourante = indexC1
  62                     cellAutre = indexC2
  63                     pass
  64                 elif (indexC2 == j):
  65                     # hypothese non verifie #
  66                     cellCourante = indexC2
  67                     cellAutre = indexC1
  68                     pass
  69                 # definir la cellule gauche et droite par le prduit vitesse * normale sortante
  70                 # si u*n>0 : rien a faire sinon inverser la gauche et la droite
  71                 if (UN > 1.E-15):
  72                     conc = y_field[cellCourante]
  73                     pass
  74                 else:
  75                     conc = y_field[cellAutre]
  76                     pass
  77                 pass
  78             else:
  79                 # conditions aux limites neumann homogene #
  80                 if (Fk.getGroupName() == "LeftEdge" or Fk.getGroupName() == "RightEdge"):
  81                     if (UN > 1.E-15):
  82                         conc = y_field[cellCourante]
  83                         pass
  84                     else:
  85                         conc = 0.0
  86                         pass
  87                     pass
  88                 # conditions aux limites periodiques #
  89                 if (Fk.getGroupName() == "BottomEdge" or Fk.getGroupName() == "TopEdge"):
  90                     indexFP = indexFacesPerio[indexFace]
  91                     # une autre manière de recuperer l'index de la face periodique #
  92                     # int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
  93                     Fp = my_mesh.getFace(indexFP)
  94                     indexCp = Fp.getCellsId()[0]
  95                     if (UN > 1.E-15):
  96                         conc = y_field[cellCourante]
  97                         pass
  98                     else:
  99                         conc = y_field[indexCp]
 100                         pass
 101                     pass
 102                 pass
 103             SumF = SumF + UN * LengthFk * conc
 104             pass
 105         dt = cfl * minlengthFk / normU
 106         SumFlux[j] = dt * SumF / Cj.getMeasure()
 107         pass
 108     return dt, SumFlux
 109
 110
 111 def EquationTransport2D(tmax, VitesseX, VitesseY, cfl, output_freq, my_mesh, file):
 112
 113     # Initial conditions #
 114     print("Construction of the initial condition …")
 115     y_field = initial_conditions(my_mesh)
 116     #
 117     #  MED output of the initial conditions at t=0 and iteration = 0
 118     #
 119
 120     iteration = 0
 121     time = 0.
 122     print("Saving the solution at T=" + str(time) + "…")
 123     y_field.setTime(time, iteration)
 124     y_field.writeMED(file)
 125     y_field.writeVTK(file)
 126     y_field.writeCSV(file)
 127
 128     # Time loop #
 129     print("Resolution of the transport equation with an UPWIND scheme …")
 130     ntmax = 3
 131     indexFacesPerio = my_mesh.getIndexFacePeriodic()
 132     dt = 0.
 133     while (iteration < ntmax and time <= tmax):
 134         dt, SumFlux = sigma_flux(VitesseX, VitesseY, cfl, y_field, indexFacesPerio)
 135         print("-- Iter: " + str(iteration) + ", Time: " + str(time) + ", dt: " + str(dt))
 136
 137         # Advancing time step #
 138         y_field -= SumFlux
 139         time += dt
 140         iteration += 1
 141         # Output every output_freq iterations
 142         if (iteration % output_freq == 0):
 143             y_field.setTime(time, iteration)
 144             y_field.writeMED(file, False)
 145             y_field.writeVTK(file, False)
 146             y_field.writeCSV(file)
 147             pass
 148         pass
 149     return
 150
 151
 152 def main():
 153     print("RESOLUTION OF THE 2D TRANSPORT EQUATION:")
 154     print("- DOMAIN: SQUARE [0,1]x[0,1]")
 155     print("- MESH: CARTESIAN, INTERNAL GENERATION WITH CDMATH")
 156     print("- PERIODIC BC ON TOP AND BOTTOM")
 157     print("- HOMOGENEOUS NEUMANN BC ON LEFT AND RIGHT")
 158
 159     # Problem data
 160     cfl = 0.4
 161     VitesseX = 1.0
 162     VitesseY = 1.0
 163     tmax = 1.
 164     output_freq = 10
 165
 166     print("Construction of a cartesian mesh …")
 167     xinf = 0.0
 168     xsup = 1.0
 169     yinf = 0.0
 170     ysup = 1.0
 171     nx = 100
 172     ny = 100
 173     my_mesh = cdmath.Mesh(xinf, xsup, nx, yinf, ysup, ny)
 174     eps = 1.E-10
 175     my_mesh.setGroupAtPlan(xsup, 0, eps, "RightEdge")
 176     my_mesh.setGroupAtPlan(xinf, 0, eps, "LeftEdge")
 177     my_mesh.setGroupAtPlan(yinf, 1, eps, "BottomEdge")
 178     my_mesh.setGroupAtPlan(ysup, 1, eps, "TopEdge")
 179     fileOutPutCart = "Exercice1PyTest"
 180     EquationTransport2D(tmax, VitesseX, VitesseY, cfl, output_freq, my_mesh, fileOutPutCart)
 181     print("CDMATH calculation done.")
 182     return
 183
 184
 185 if __name__ == """__main__""":
 186     main()