CDMATH/tests/examples/BurgersEquation1D/1DBurgersEquation_rarefaction_wave.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 #===============================================================================================================================
   5 # Name        : Résolution VF de l'équation du Burgers 1D \partial_t u + 1/2 \partial_x u^2 = 0
   6 # Author      : Michaël Ndjinga
   7 # Copyright   : CEA Saclay 2019
   8 # Description : Comparaison des schémas Upwind et de Godunov explicite sur un maillage 1D régulier
   9 #               Donnée initiale discontinue générant une onde de raréfaction
  10 #               Conditions aux limites de Neumann
  11 #               Conclusion : le schéma Upwind n'est pas entropique, contrairement au schéma de Godunov
  12 #                       Création et sauvegarde du champ résultant et des figures
  13 #               Génération d'une video sauvegardée dans un fichier .mp4
  14 #================================================================================================================================
  15
  16
  17 import numpy as np
  18 from copy import deepcopy
  19 import matplotlib
  20 matplotlib.use("Agg")
  21 import matplotlib.pyplot as plt
  22 import matplotlib.animation as manimation
  23
  24 def Flux_upwind(u_l, u_r):
  25     if ((u_l+u_r)/2>0):
  26         flux = 0.5*u_l*u_l
  27     elif ((u_l+u_r)/2<0):
  28         flux = 0.5*u_r*u_r
  29     else:
  30         flux = (0.5*u_l*u_l + 0.5*u_r*u_r)/2
  31
  32     return flux
  33
  34 def Flux_Godunov(u_l, u_r):
  35     if (u_l==u_r):
  36         flux = 0.5*u_l*u_l
  37     elif (u_l<0 and 0<u_r):
  38         flux = 0.;
  39     elif (u_l<u_r):
  40         flux = min(0.5*u_l*u_l,0.5*u_r*u_r);
  41     elif (u_l>u_r):
  42         flux = max(0.5*u_l*u_l,0.5*u_r*u_r);
  43
  44     return flux
  45
  46 def Burgers1D():
  47     print( "Simulation of 1D Burgers' equation with upwind and godunov explicit schemes" )
  48
  49     ##################### Simulation parameters
  50     a = -1. # space domain :  a <= x <= b
  51     b = 1.
  52     nx=50
  53     dx = (b - a) / nx #space step
  54
  55     tmax = 4 # runs the simulation for 0 <= t <= tMax
  56     ntmax=150
  57     cfl=0.95
  58
  59     x=[a+0.5*dx + i*dx for i in range(nx)]   # array of cell center (1D mesh)
  60
  61     ########################## Initial data
  62     u_initial = [ (xi<(a+b)/2)*(-1) + (xi>(a+b)/2)*(1.)  for xi in x]
  63     u_upwind  = [ (xi<(a+b)/2)*(-1) + (xi>(a+b)/2)*(1.)  for xi in x]
  64     u_godunov = [ (xi<(a+b)/2)*(-1) + (xi>(a+b)/2)*(1.)  for xi in x]
  65     u_exact   = [ (xi<(a+b)/2)*(-1) + (xi>(a+b)/2)*(1.)  for xi in x]
  66
  67     Unp1      = [0]*nx
  68     Unp1_godunov = [0]*nx
  69
  70     max_initial=max(u_initial)
  71     min_initial=min(u_initial)
  72
  73     time = 0.
  74     it = 0
  75     output_freq = 10
  76
  77     # Video settings
  78     FFMpegWriter = manimation.writers['ffmpeg']
  79     metadata = dict(title="Finite volumes schemes for Burgers equation", artist = "CEA Saclay", comment="Rarefaction wave")
  80     writer=FFMpegWriter(fps=output_freq, metadata=metadata, codec='h264')
  81     with writer.saving(plt.figure(), "1DBurgersEquation_FV"+".mp4", ntmax):
  82         ########################### Postprocessing initialisation
  83         # Picture frame
  84         plt.xlabel('x')
  85         plt.ylabel('u')
  86         plt.xlim(a,b)
  87         plt.ylim( min_initial - 0.1*(max_initial-min_initial), max_initial +  0.5*(max_initial-min_initial) )
  88         plt.title('Finite volume schemes for Burgers equation')
  89         line1, = plt.plot(x, u_upwind,  label='Conservative (Upwind) scheme') #new picture for video # Returns a tuple of line objects, thus the comma
  90         line2, = plt.plot(x, u_godunov, label='Conservative (Godunov) scheme') #new picture for video # Returns a tuple of line objects, thus the comma
  91         line3, = plt.plot(x, u_exact,   label='Exact solution') #new picture for video # Returns a tuple of line objects, thus the comma
  92         plt.legend()
  93
  94         print("Starting time loop")
  95         print("-- Iter: " + str(it) + ", Time: " + str(time) )
  96         np.savetxt("BurgersEquation_FV_Upwind_ResultField_0" +".txt", u_upwind,  delimiter="\n")
  97         np.savetxt("BurgersEquation_FV_Godunov_ResultField_0"+".txt", u_godunov, delimiter="\n")
  98         np.savetxt("BurgersEquation_FV_Exact_ResultField_0"+".txt",   u_exact,   delimiter="\n")
  99         writer.grab_frame()
 100         plt.savefig("BurgersEquation_FV_Rarefaction_ResultField_0"+".png")
 101
 102         ############################# Time loop
 103         while (it < ntmax and time <= tmax):
 104             upw = max(np.abs(u_upwind))
 105             dt = cfl * dx / upw
 106             # Loop on all cells
 107             for i in range(0,nx):
 108                 if (i==0):
 109                     flux_iminus = 0.5*u_upwind[0]*u_upwind[0]#Flux at the left Neumann boundary
 110                 if (i==nx-1):
 111                     flux_iplus  = 0.5*u_upwind[nx-1]*u_upwind[nx-1]#Flux at the right Neumann boundary
 112                 else:
 113                     flux_iplus = Flux_upwind(u_upwind[i],u_upwind[i+1])
 114                 pass
 115                 Unp1[i] = u_upwind[i] - dt/dx*(flux_iplus-flux_iminus);
 116                 flux_iminus = flux_iplus;
 117             u_upwind = Unp1
 118
 119             for i in range(0,nx):
 120                 if (i==0):
 121                     flux_iminus = 0.5*u_godunov[0]*u_godunov[0]#Flux at the left Neumann boundary
 122                 if (i==nx-1):
 123                     flux_iplus  = 0.5*u_godunov[nx-1]*u_godunov[nx-1]#Flux at the right Neumann boundary
 124                 else:
 125                     flux_iplus = Flux_Godunov(u_godunov[i],u_godunov[i+1])
 126                 pass
 127                 Unp1_godunov[i] = u_godunov[i] - dt/dx*(flux_iplus-flux_iminus);
 128                 flux_iminus = flux_iplus;
 129             u_godunov = Unp1_godunov
 130
 131             time += dt
 132             it += 1
 133
 134             u_exact = [ (xi<-time)*(-1) + (xi>time)*(1.) + (-time<=xi and xi<=time)*(xi/time) for xi in x]
 135
 136             # Postprocessing
 137             line1.set_ydata(u_upwind)
 138             line2.set_ydata(u_godunov)
 139             line3.set_ydata(u_exact)
 140             writer.grab_frame()
 141
 142             if (it % output_freq == 0):
 143                 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 144                 np.savetxt( "BurgersEquation_Upwind_ResultField_" +str(it)+".txt", u_upwind,  delimiter="\n")
 145                 np.savetxt( "BurgersEquation_Godunov_ResultField_"+str(it)+".txt", u_godunov, delimiter="\n")
 146                 np.savetxt( "BurgersEquation_Exact_ResultField_"+str(it)+".txt", u_exact, delimiter="\n")
 147                 plt.savefig("BurgersEquation_FV_Rarefaction_ResultField_"+str(it)+".png")
 148                 #plt.show()
 149
 150     print("Simulation of rarefaction wave on Burgers' equation with finite volume schemes done.")
 151
 152     return
 153
 154 if __name__ == """__main__""":
 155     Burgers1D()