CDMATH/tests/examples/TransportEquation1DUpwindExplicit/1DTransportEquationUpwindExplicit.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 #===============================================================================================================================
   5 # Name        : Résolution VF de l'équation du transport 1D \partial_t u + c \partial_x u = 0 avec conditions aux limites périodiques
   6 # Author      : Michaël Ndjinga, Katia Ait Ameur
   7 # Copyright   : CEA Saclay 2018
   8 # Description : Utilisation du schéma upwind explicite sur un maillage 1D régulier
   9 #                       Création et sauvegarde du champ résultant et des figures
  10 #               Génération d'une video sauvegardée dans un fichier .mp4
  11 #================================================================================================================================
  12
  13
  14 from math import sin, pi, ceil
  15 import numpy as np
  16 import matplotlib
  17 matplotlib.use("Agg")
  18 import matplotlib.pyplot as plt
  19 import matplotlib.animation as manimation
  20 import sys
  21 from copy import deepcopy
  22
  23 def Transport1DUpwindExplicit(nx,cfl):
  24     print( "Simulation of 1D transport equation with explicit upwind scheme" )
  25
  26     ##################### Simulation parameters
  27     a = 0.0 # space domain :  a <= x <= b
  28     b = 1.0
  29     dx = (b - a) / nx #space step
  30
  31     c = 0.25 # advection velocity
  32     tmax = (b-a)/c # runs the simulation for 0 <= t <= tMax
  33     dt = cfl * dx / c
  34     ntmax = ceil(tmax/dt)
  35
  36     x=[a+0.5*dx + i*dx for i in range(nx)]   # array of cell center (1D mesh)
  37
  38     ########################## Initial data
  39
  40     u_initial = [ 0.5*(1+sin(4*pi*xi-pi*.5))*int(xi<0.5)*int(0<xi) + int(0.6<xi)*int(xi<0.85)  for xi in x];# to be used with a=0, b=1
  41     u         = [ 0.5*(1+sin(4*pi*xi-pi*.5))*int(xi<0.5)*int(0<xi) + int(0.6<xi)*int(xi<0.85)  for xi in x];# to be used with a=0, b=1
  42
  43     max_initial=max(u_initial)
  44     min_initial=min(u_initial)
  45     total_var_initial = np.sum([abs(u_initial[i] - u_initial[(i-1)%nx]) for i in range(nx)])
  46
  47     time = 0.
  48     it = 0
  49     output_freq = 10
  50
  51     # Video settings
  52     FFMpegWriter = manimation.writers['ffmpeg']
  53     metadata = dict(title="Upwind explicit scheme for transport equation, CFL="+str(cfl), artist = "CEA Saclay", comment="Stable if CFL<1")
  54     writer=FFMpegWriter(fps=output_freq, metadata=metadata, codec='h264')
  55     with writer.saving(plt.figure(), "1DTransportEquation_UpwindExplicit_"+str(nx)+"Cells_CFL"+str(cfl)+".mp4", ntmax):
  56         ########################### Postprocessing initialisation
  57         # Picture frame
  58         plt.legend()
  59         plt.xlabel('x')
  60         plt.ylabel('u')
  61         plt.xlim(a,b)
  62         plt.ylim( min_initial - 0.1*(max_initial-min_initial), max_initial +  0.1*(max_initial-min_initial) )
  63         plt.title("Upwind explicit scheme for transport equation, CFL="+str(cfl))
  64         line1, = plt.plot(x, u, label='u') #new picture for video # Returns a tuple of line objects, thus the comma
  65
  66         print("Starting time loop")
  67         print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
  68         np.savetxt( "TransportEquation_UpwindExplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.txt", u, delimiter="\n")
  69         writer.grab_frame()
  70         plt.savefig("TransportEquation_UpwindExplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.png")
  71
  72         ############################# Time loop
  73         while (it < ntmax and time <= tmax):
  74             un=deepcopy(u)
  75             for i in range(nx):
  76                 u[i] = un[i] - c * dt / dx * (un[i] - un[(i-1)%nx])
  77
  78             time += dt
  79             it += 1
  80
  81             if cfl<1 :
  82                 assert max(u) <= max_initial
  83                 assert min(u) >= min_initial
  84                 assert np.sum([abs(u[i] - u[(i-1)%nx]) for i in range(nx)]) <= total_var_initial + 1e-15#Tiny margin required
  85
  86             # Postprocessing
  87             line1.set_ydata(u)
  88             writer.grab_frame()
  89             if (it % output_freq == 0):
  90                 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
  91                 np.savetxt( "TransportEquation_UpwindExplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".txt", u, delimiter="\n")
  92                 plt.savefig("TransportEquation_UpwindExplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".png")
  93                 #plt.show()
  94
  95     print( "Exact solution minimum   : ", min(u_initial), "Numerical solution minimum   : ",  min(u) )
  96     print( "Exact solution maximum   : ", max(u_initial), "Numerical solution maximum   : ",  max(u) )
  97     print( "Exact solution variation : ", np.sum([abs(u_initial[i] - u_initial[(i-1)%nx]) for i in range(nx)]), "Numerical solution variation : ",  np.sum([abs(u[i] - u[(i-1)%nx]) for i in range(nx)]) )
  98     print( "l1 numerical error       : ", dx*np.sum([abs(u[i] - u_initial[i]) for i in range(nx)]) )
  99
 100     print("Simulation of transport equation with explicit upwind scheme done.")
 101
 102     #return min, max, total variation and l1 error
 103     return min(u), max(u), np.sum([abs(u[i] - u[(i-1)%nx]) for i in range(nx)]), dx*np.sum([abs(u[i] - u_initial[i]) for i in range(nx)])
 104
 105
 106 if __name__ == """__main__""":
 107     if len(sys.argv) >2 :
 108         nx = int(sys.argv[1])
 109         cfl = float(sys.argv[2])
 110         Transport1DUpwindExplicit(nx,cfl)
 111     else :
 112         nx = 50 # number of cells
 113         cfl = 0.99 # c*dt/dx <= CFL
 114         Transport1DUpwindExplicit(nx,cfl)
 115