CDMATH/tests/examples/BurgersEquation1D/1DBurgersEquation_shock_formation.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 du schéma de Godunov explicite avec des schémas alternatifs sur un maillage 1D régulier
   9 #               Premier schéma alternatif = Godunov sur l'équation conservative \partial_t u^3 + 1/3 \partial_x u^3 = 0
  10 #               Deuxième schéma alternatif = Upwind sur l'équation non conservative \partial_t u + u \partial_x u = 0
  11 #               Donnée initiale continue générant une onde de choc
  12 #               Conditions aux limites de Neumann
  13 #               Conclusion : 1) le choc ne peut être capturé qu'avec un schéma conservatif 2) la vitesse du choc dépend de la formulation conservative employée
  14 #                       Création et sauvegarde du champ résultant et des figures
  15 #               Génération d'une video sauvegardée dans un fichier .mp4
  16 #================================================================================================================================
  17
  18
  19 from math import sin, cos, pi, sqrt
  20 import numpy as np
  21 from copy import deepcopy
  22 import matplotlib
  23 matplotlib.use("Agg")
  24 import matplotlib.pyplot as plt
  25 import matplotlib.animation as manimation
  26
  27 def Flux_Godunov(u_l, u_r):
  28     if (u_l==u_r):
  29         flux = 0.5*u_l*u_l
  30     elif (u_l<0 and 0<u_r):
  31         flux = 0.;
  32     elif (u_l<u_r):
  33         flux = min(0.5*u_l*u_l,0.5*u_r*u_r);
  34     elif (u_l>u_r):
  35         flux = max(0.5*u_l*u_l,0.5*u_r*u_r);
  36
  37     return flux
  38
  39 def Flux_Godunov2(u_l, u_r):
  40     if (u_l==u_r):
  41         flux = 1./3.*u_l*u_l*u_l
  42     elif (u_l*u_l<u_r*u_r):
  43         flux = min(1./3.*u_l*u_l*u_l,1./3.*u_r*u_r*u_r)
  44     elif (u_l*u_l>u_r*u_r):
  45         flux = max(1./3.*u_l*u_l*u_l,1./3.*u_r*u_r*u_r)
  46     else:
  47         print "u_l ",u_l, " u_r ", u_r
  48     return flux
  49
  50 def Du_ncsv(u_l, u_i, u_r):
  51     if (u_i<0):
  52         Du= u_r-u_i
  53     elif (0<u_i):
  54         Du= u_i-u_l
  55     else:
  56         Du= u_r-u_l/2
  57
  58     return Du
  59
  60 def Burgers1D():
  61     print( "Simulation of 1D Burgers' equation with various explicit schemes" )
  62
  63     ##################### Simulation parameters
  64     a = 0.0 # space domain :  a <= x <= b
  65     b = 1.
  66     nx=100
  67     dx = (b - a) / nx #space step
  68
  69     tmax = 1.5 # runs the simulation for 0 <= t <= tMax
  70     ntmax=100
  71     cfl=0.95
  72
  73     x=[a+0.5*dx + i*dx for i in range(nx)]   # array of cell center (1D mesh)
  74
  75     ########################## Initial data
  76     xa=0.1
  77     xb=0.5
  78     u_initial = [ (xi<xa)+(xa<=xi)*(xi<=xb)*(np.cos(np.pi*(xa-xi)/(xa-xb))+1.0)*0.5  for xi in x];# to be used with a=0, b=1
  79     u_godunov = [ (xi<xa)+(xa<=xi)*(xi<=xb)*(np.cos(np.pi*(xa-xi)/(xa-xb))+1.0)*0.5  for xi in x];# to be used with a=0, b=1
  80     u_ncsv    = [ (xi<xa)+(xa<=xi)*(xi<=xb)*(np.cos(np.pi*(xa-xi)/(xa-xb))+1.0)*0.5  for xi in x];# to be used with a=0, b=1
  81     u_csv2    = deepcopy(u_initial)# to be used with a=0, b=1
  82     Unp1      = [0]*nx
  83     Unp1_ncsv = [0]*nx
  84     Unp1_csv2 = [0]*nx
  85
  86     max_initial=max(u_initial)
  87     min_initial=min(u_initial)
  88
  89     time = 0.
  90     it = 0
  91     output_freq = 10
  92
  93     # Video settings
  94     FFMpegWriter = manimation.writers['ffmpeg']
  95     metadata = dict(title="Finite volumes schemes for Burgers equation", artist = "CEA Saclay", comment="Shock formation")
  96     writer=FFMpegWriter(fps=output_freq, metadata=metadata, codec='h264')
  97     with writer.saving(plt.figure(), "1DBurgersEquation_FV"+".mp4", ntmax):
  98         ########################### Postprocessing initialisation
  99         # Picture frame
 100         plt.xlabel('x')
 101         plt.ylabel('u')
 102         plt.xlim(a,b)
 103         plt.ylim( min_initial - 0.1*(max_initial-min_initial), max_initial +  0.3*(max_initial-min_initial) )
 104         plt.title("Finite volume schemes for Burgers' equation on " + str(nx) + " cells")
 105         line1, = plt.plot(x, u_godunov, label='Conservative (Godunov) scheme') #new picture for video # Returns a tuple of line objects, thus the comma
 106         line2, = plt.plot(x, u_ncsv,    label='Non conservative scheme') #new picture for video # Returns a tuple of line objects, thus the comma
 107         line3, = plt.plot(x, u_csv2,    label='Alternative conservative scheme') #new picture for video # Returns a tuple of line objects, thus the comma
 108         plt.legend()
 109
 110         print("Starting time loop")
 111         print("-- Iter: " + str(it) + ", Time: " + str(time) )
 112         np.savetxt("BurgersEquation_FVGodunov_ResultField_0"+".txt", u_godunov, delimiter="\n")
 113         np.savetxt("BurgersEquation_FVNonCons_ResultField_0"+".txt", u_ncsv, delimiter="\n")
 114         np.savetxt("BurgersEquation_FVAltCons_ResultField_0"+".txt", u_csv2, delimiter="\n")
 115         writer.grab_frame()
 116         plt.savefig("BurgersEquation_FV_Shock_ResultField_0"+".png")
 117
 118         ############################# Time loop
 119         while (it < ntmax and time <= tmax):
 120             upw = max(np.abs(u_godunov))
 121             dt = cfl * dx / upw
 122             # Loop on all cells
 123             for i in range(0,nx):
 124                 if (i==0):
 125                     flux_iminus = 0.5*u_godunov[0]*u_godunov[0]#Flux at the left Neumann boundary
 126                 if (i==nx-1):
 127                     flux_iplus  = 0.5*u_godunov[nx-1]*u_godunov[nx-1]#Flux at the right Neumann boundary
 128                 else:
 129                     flux_iplus = Flux_Godunov(u_godunov[i],u_godunov[i+1])
 130                 pass
 131                 Unp1[i] = u_godunov[i] - dt/dx*(flux_iplus-flux_iminus);
 132                 flux_iminus = flux_iplus;
 133             u_godunov = Unp1
 134
 135             for i in range(0,nx):
 136                 if (i==0):
 137                     Du = Du_ncsv(u_ncsv[0],    u_ncsv[0],    u_ncsv[1])# Neumann boundary condition
 138                 elif (i==nx-1):
 139                     Du = Du_ncsv(u_ncsv[nx-2], u_ncsv[nx-1], u_ncsv[nx-1])# Neumann boundary condition
 140                 else:
 141                     Du = Du_ncsv(u_ncsv[i-1],    u_ncsv[i],    u_ncsv[i+1])
 142                 pass
 143                 Unp1_ncsv[i] = u_ncsv[i] - dt/dx*u_ncsv[i]*Du
 144             u_ncsv = Unp1_ncsv
 145
 146             for i in range(0,nx):
 147                 if (i==0):
 148                     flux_iminus = 1./3.*u_csv2[0]*u_csv2[0]*u_csv2[0]#Flux at the left Neumann boundary
 149                 if (i==nx-1):
 150                     flux_iplus  = 1./3.*u_csv2[nx-1]*u_csv2[nx-1]*u_csv2[nx-1]#Flux at the right Neumann boundary
 151                 else:
 152                     flux_iplus = Flux_Godunov2(u_csv2[i],u_csv2[i+1])
 153                 pass
 154                 Unp1_csv2[i] = sqrt(u_csv2[i]*u_csv2[i] - dt/dx*(flux_iplus-flux_iminus))
 155                 flux_iminus = flux_iplus;
 156             u_csv2 = Unp1_csv2
 157
 158             time += dt
 159             it += 1
 160
 161             # Postprocessing
 162             line1.set_ydata(u_godunov)
 163             line2.set_ydata(u_ncsv)
 164             line3.set_ydata(u_csv2)
 165             writer.grab_frame()
 166             if (it % output_freq == 0):
 167                 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 168                 np.savetxt( "BurgersEquation_FVGodunov_ResultField_"+str(it)+".txt", u_godunov, delimiter="\n")
 169                 np.savetxt( "BurgersEquation_FVNonCons_ResultField_"+str(it)+".txt", u_ncsv,    delimiter="\n")
 170                 np.savetxt( "BurgersEquation_FVAltCons_ResultField_"+str(it)+".txt", u_csv2,    delimiter="\n")
 171                 plt.savefig("BurgersEquation_FV_Shock_ResultField_"+str(it)+".png")
 172                 #plt.show()
 173
 174     print("Simulation of shock formation on Burgers' equation with finite volume schemes done.")
 175
 176     return
 177
 178 if __name__ == """__main__""":
 179     Burgers1D()