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 #================================================================================================================================
19 from math import sin, cos, pi, sqrt
21 from copy import deepcopy
24 import matplotlib.pyplot as plt
25 import matplotlib.animation as manimation
27 def Flux_Godunov(u_l, u_r):
30 elif (u_l<0 and 0<u_r):
33 flux = min(0.5*u_l*u_l,0.5*u_r*u_r);
35 flux = max(0.5*u_l*u_l,0.5*u_r*u_r);
39 def Flux_Godunov2(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)
47 print( "u_l ",u_l, " u_r ", u_r )
50 def Du_ncsv(u_l, u_i, u_r):
61 print( "Simulation of 1D Burgers' equation with various explicit schemes" )
63 ##################### Simulation parameters
64 a = 0.0 # space domain : a <= x <= b
67 dx = (b - a) / nx #space step
69 tmax = 1.5 # runs the simulation for 0 <= t <= tMax
73 x=[a+0.5*dx + i*dx for i in range(nx)] # array of cell center (1D mesh)
75 ########################## Initial data
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
86 max_initial=max(u_initial)
87 min_initial=min(u_initial)
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
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
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")
116 plt.savefig("BurgersEquation_FV_Shock_ResultField_0"+".png")
118 ############################# Time loop
119 while (it < ntmax and time <= tmax):
120 upw = max(np.abs(u_godunov))
123 for i in range(0,nx):
125 flux_iminus = 0.5*u_godunov[0]*u_godunov[0]#Flux at the left Neumann boundary
127 flux_iplus = 0.5*u_godunov[nx-1]*u_godunov[nx-1]#Flux at the right Neumann boundary
129 flux_iplus = Flux_Godunov(u_godunov[i],u_godunov[i+1])
131 Unp1[i] = u_godunov[i] - dt/dx*(flux_iplus-flux_iminus);
132 flux_iminus = flux_iplus;
135 for i in range(0,nx):
137 Du = Du_ncsv(u_ncsv[0], u_ncsv[0], u_ncsv[1])# Neumann boundary condition
139 Du = Du_ncsv(u_ncsv[nx-2], u_ncsv[nx-1], u_ncsv[nx-1])# Neumann boundary condition
141 Du = Du_ncsv(u_ncsv[i-1], u_ncsv[i], u_ncsv[i+1])
143 Unp1_ncsv[i] = u_ncsv[i] - dt/dx*u_ncsv[i]*Du
146 for i in range(0,nx):
148 flux_iminus = 1./3.*u_csv2[0]*u_csv2[0]*u_csv2[0]#Flux at the left Neumann boundary
150 flux_iplus = 1./3.*u_csv2[nx-1]*u_csv2[nx-1]*u_csv2[nx-1]#Flux at the right Neumann boundary
152 flux_iplus = Flux_Godunov2(u_csv2[i],u_csv2[i+1])
154 Unp1_csv2[i] = sqrt(u_csv2[i]*u_csv2[i] - dt/dx*(flux_iplus-flux_iminus))
155 flux_iminus = flux_iplus;
162 line1.set_ydata(u_godunov)
163 line2.set_ydata(u_ncsv)
164 line3.set_ydata(u_csv2)
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")
174 print("Simulation of shock formation on Burgers' equation with finite volume schemes done.")
178 if __name__ == """__main__""":