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 implicite 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 #================================================================================================================================
17 import matplotlib.pyplot as plt
18 import matplotlib.animation as manimation
20 from math import sin, pi, ceil
25 def upwindSchemeMatrix(nx,cfl):
26 upwindMat=cdmath.SparseMatrixPetsc(nx,nx,2)
28 upwindMat.setValue(i,i,1+cfl)
29 upwindMat.setValue(i,(i-1)%nx,-cfl)
33 def Transport1DUpwindImplicit(nx,cfl):
34 print( "Simulation of 1D transport equation with implicit upwind scheme" )
36 ##################### Simulation parameters
37 a = 0.0 # space domain : a <= x <= b
39 dx = (b - a) / nx #space step
41 c = 0.25 # advection velocity
42 tmax = (b-a)/c # runs the simulation for 0 <= t <= tMax
47 raise("Impossible to run this simulation with cfl>nx. Choose another value for nx or cfl.")
49 x=[a+0.5*dx + i*dx for i in range(nx)] # array of cell center (1D mesh)
51 ########################## Initial data
53 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
54 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
56 max_initial=max(u_initial)
57 min_initial=min(u_initial)
58 total_var_initial = np.sum([abs(u_initial[i] - u_initial[(i-1)%nx]) for i in range(nx)])
64 #Linear system initialisation
65 systemMat=upwindSchemeMatrix(nx,cfl)
70 LS=cdmath.LinearSolver(systemMat,Un,iterGMRESMax, precision, "GMRES","ILU")
73 FFMpegWriter = manimation.writers['ffmpeg']
74 metadata = dict(title="Upwind implicit scheme for transport equation, "+"CFL="+str(cfl), artist = "CEA Saclay", comment="Stable for any CFL>0")
75 writer=FFMpegWriter(fps=output_freq, metadata=metadata, codec='h264')
76 with writer.saving(plt.figure(), "1DTransportEquation_UpwindImplicit_"+str(nx)+"CELLS_CFL"+str(cfl)+".mp4", ntmax):
77 ########################### Postprocessing initialisation
83 plt.ylim( min_initial - 0.1*(max_initial-min_initial), max_initial + 0.1*(max_initial-min_initial) )
84 plt.title("Upwind implicit scheme for transport equation, "+"CFL="+str(cfl))
85 line1, = plt.plot(x, u, label='u') #new picture for video # Returns a tuple of line objects, thus the comma
87 print("Starting time loop")
88 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
89 np.savetxt( "TransportEquation_UpwindImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.txt", u, delimiter="\n")
91 plt.savefig("TransportEquation_UpwindImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.png")
93 ############################# Time loop
94 while (it < ntmax and time <= tmax):
100 if(not LS.getStatus()):
101 print( "Linear system did not converge ", iterGMRES, " GMRES iterations")
102 raise ValueError("Pas de convergence du système linéaire");
109 assert max(u) <= max_initial
110 assert min(u) >= min_initial
111 assert np.sum([abs(u[i] - u[(i-1)%nx]) for i in range(nx)]) <= total_var_initial
116 if (it % output_freq == 0):
117 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
118 np.savetxt( "TransportEquation_UpwindImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".txt", u, delimiter="\n")
119 plt.savefig("TransportEquation_UpwindImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".png")
122 print( "Exact solution minimum : ", min(u_initial), "Numerical solution minimum : ", min(u) )
123 print( "Exact solution maximum : ", max(u_initial), "Numerical solution maximum : ", max(u) )
124 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)]) )
125 print( "l1 numerical error : ", dx*np.sum([abs(u[i] - u_initial[i]) for i in range(nx)]) )
127 print("Simulation of transport equation with implicit upwind scheme done.")
129 #return min, max, total variation and l1 error
130 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)])
133 if __name__ == """__main__""":
134 if len(sys.argv) >2 :
135 nx = int(sys.argv[1])
136 cfl = float(sys.argv[2])
137 Transport1DUpwindImplicit(nx,cfl)
139 nx = 50 # number of cells
140 cfl = 0.99 # c*dt/dx <= CFL
141 Transport1DUpwindImplicit(nx,cfl)