CDMATH/tests/examples/transport1d/main.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 import math
   5
   6 import cdmath
   7
   8
   9 def main():
  10     a = -5.0
  11     b = 5.0
  12     nx = 1000
  13     ntmax = 1000
  14     dx = (b - a) / nx
  15     pi = 3.1415927
  16     # Transport velocity
  17     cfl = 0.5
  18     u = 3.
  19     dt = cfl * dx / u
  20
  21     my_mesh = cdmath.Mesh(a, b, nx)
  22     conc = cdmath.Field("Concentration", cdmath.CELLS, my_mesh, 1)
  23
  24     # Initial conditions
  25     sigma = math.sqrt(0.2)
  26     for i in range(my_mesh.getNumberOfCells()):
  27         x = my_mesh.getCell(i).x()
  28         conc[i] = 0.5 / (sigma * math.sqrt(2 * pi)) * math.exp(-0.5 * math.pow((x / sigma), 2))
  29         pass
  30
  31     time = 0.
  32     tmax = 3.0
  33     it = 0
  34
  35     print("MED post-treatment of the solution at T=" + str(time) + "…")
  36     output_filename = "EqTr1D"
  37     conc.setTime(time, it)
  38     conc.writeMED(output_filename)
  39     conc.writeVTK(output_filename)
  40     conc.writeCSV(output_filename)
  41     output_freq = 10
  42
  43     # Time loop
  44     while (it < ntmax and time <= tmax):
  45         print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
  46         conc[0] = conc[0] - u * dt / dx * (conc[0] - conc[my_mesh.getNumberOfCells() - 1])
  47         for j in range(1, my_mesh.getNumberOfCells()):
  48             conc[j] = conc[j] - u * dt / dx * (conc[j] - conc[j - 1])
  49             pass
  50         time += dt
  51         it += 1
  52         if (it % output_freq == 0):
  53             conc.setTime(time, it)
  54             conc.writeMED(output_filename, False)
  55             conc.writeVTK(output_filename, False)
  56             conc.writeCSV(output_filename)
  57             pass
  58         pass
  59     print("CDMATH calculation done.")
  60     return
  61
  62
  63 if __name__ == """__main__""":
  64     main()