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 centré implicite sur un maillage 1D régulier
9 # Schéma à 3 points implicite
10 # Création et sauvegarde du champ résultant et des figures
11 # Génération d'une video sauvegardée dans un fichier .mp4
12 #================================================================================================================================
18 import matplotlib.pyplot as plt
19 import matplotlib.animation as manimation
21 from math import sin, pi, ceil
24 def centeredSchemeMatrix(nx,cfl):
25 centeredMat=cdmath.SparseMatrixPetsc(nx,nx,3)
27 centeredMat.setValue(i,(i+1)%nx,cfl/2)
28 centeredMat.setValue(i,i,1.)
29 centeredMat.setValue(i,(i-1)%nx,-cfl/2.)
33 def Transport1DCenteredImplicit(nx,cfl):
34 print( "Simulation of 1D transport equation with implicit centered 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=centeredSchemeMatrix(nx,cfl)
71 LS=cdmath.LinearSolver(systemMat,Un,iterGMRESMax, precision, "GMRES","ILU")
72 LS.setComputeConditionNumber()
75 FFMpegWriter = manimation.writers['ffmpeg']
76 metadata = dict(title="Centered implicit scheme for transport equation, "+"CFL="+str(cfl), artist = "CEA Saclay", comment="Stable for any CFL>0")
77 writer=FFMpegWriter(fps=output_freq, metadata=metadata, codec='h264')
78 with writer.saving(plt.figure(), "1DTransportEquation_CenteredImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+".mp4", ntmax):
79 ########################### Postprocessing initialisation
85 plt.ylim( min_initial - 0.1*(max_initial-min_initial), max_initial + 0.1*(max_initial-min_initial) )
86 plt.title("Centered implicit scheme for transport equation, "+"CFL="+str(cfl))
87 line1, = plt.plot(x, u, label='u') #new picture for video # Returns a tuple of line objects, thus the comma
89 print("Starting time loop")
90 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
91 np.savetxt( "TransportEquation_CenteredImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.txt", u, delimiter="\n")
93 plt.savefig("TransportEquation_CenteredImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.png")
95 ############################# Time loop
96 while (it < ntmax and time <= tmax):
102 if(not LS.getStatus()):
103 print( "Linear system did not converge ", iterGMRES, " GMRES iterations" )
104 raise ValueError("Pas de convergence du système linéaire");
108 if ( max(u) > max_initial ):
109 print( "-- Iter: " + str(it) + " max principle violated : max(t) > max(0) : max(t)= ",max(u), " max(0)= ", max_initial )
110 if ( min(u) < min_initial ):
111 print( "-- Iter: " + str(it) + " min principle violated : min(t) < min(0) : min(t)= ",min(u), " min(0)= ", min_initial )
112 if ( np.sum([abs(u[i] - u[(i-1)%nx]) for i in range(nx)]) > total_var_initial ):
113 print( "-- Iter: " + str(it) + " total variation increased : var(t) > var(0) : var(t)= ", np.sum([abs(u[i] - u[(i-1)%nx]) for i in range(nx)]), " var(0)= ", total_var_initial )
121 if (it % output_freq == 0):
122 print( "-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
123 np.savetxt( "TransportEquation_CenteredImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".txt", u, delimiter="\n")
124 plt.savefig("TransportEquation_CenteredImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".png")
127 print( "Exact solution minimum : ", min(u_initial), "Numerical solution minimum : ", min(u) )
128 print( "Exact solution maximum : ", max(u_initial), "Numerical solution maximum : ", max(u) )
129 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)]) )
130 print( "l1 numerical error : ", dx*np.sum([abs(u[i] - u_initial[i]) for i in range(nx)]) )
132 print("Simulation of transport equation with implicit centered scheme done.")
134 #return min, max, total variation and l1 error
135 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)])
138 if __name__ == """__main__""":
139 if len(sys.argv) >2 :
140 nx = int(sys.argv[1])
141 cfl = float(sys.argv[2])
142 Transport1DCenteredImplicit(nx,cfl)
144 nx = 50 # number of cells
145 cfl = 0.99 # c*dt/dx <= CFL
146 Transport1DCenteredImplicit(nx,cfl)