CDMATH/tests/examples/TransportEquation1DUpwindImplicit/1DTransportEquationUpwindImplicit.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   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 #================================================================================================================================
  12
  13
  14 import numpy as np
  15 import matplotlib
  16 matplotlib.use("Agg")
  17 import matplotlib.pyplot as plt
  18 import matplotlib.animation as manimation
  19 import sys
  20 from math import sin, pi, ceil
  21 import cdmath
  22
  23 precision=1.e-5
  24
  25 def upwindSchemeMatrix(nx,cfl):
  26     upwindMat=cdmath.SparseMatrixPetsc(nx,nx,2)
  27     for i in range(nx):
  28         upwindMat.setValue(i,i,1+cfl)
  29         upwindMat.setValue(i,(i-1)%nx,-cfl)
  30
  31     return upwindMat
  32
  33 def Transport1DUpwindImplicit(nx,cfl):
  34     print "Simulation of 1D transport equation with implicit upwind scheme"
  35
  36     ##################### Simulation parameters
  37     a = 0.0 # space domain :  a <= x <= b
  38     b = 1.0
  39     dx = (b - a) / nx #space step
  40
  41     c = 0.25 # advection velocity
  42     tmax = (b-a)/c # runs the simulation for 0 <= t <= tMax
  43     dt = cfl * dx / c
  44     ntmax = ceil(tmax/dt)
  45
  46     if(cfl>nx):
  47         raise("Impossible to run this simulation with cfl>nx. Choose another value for nx or cfl.")
  48
  49     x=[a+0.5*dx + i*dx for i in range(nx)]   # array of cell center (1D mesh)
  50
  51     ########################## Initial data
  52
  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
  55
  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)])
  59
  60     time = 0.
  61     it = 0
  62     output_freq = 10
  63
  64     #Linear system initialisation
  65     systemMat=upwindSchemeMatrix(nx,cfl)
  66     iterGMRESMax=50
  67     Un =cdmath.Vector(nx)
  68     for i in range(nx):
  69         Un[i]=u[i]
  70     LS=cdmath.LinearSolver(systemMat,Un,iterGMRESMax, precision, "GMRES","ILU")
  71
  72     # Video settings
  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
  78         # Picture frame
  79         plt.legend()
  80         plt.xlabel('x')
  81         plt.ylabel('u')
  82         plt.xlim(a,b)
  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
  86
  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")
  90         writer.grab_frame()
  91         plt.savefig("TransportEquation_UpwindImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.png")
  92
  93         ############################# Time loop
  94         while (it < ntmax and time <= tmax):
  95             # Solve linear system
  96             for i in range(nx):
  97                 Un[i]=u[i]
  98             LS.setSndMember(Un)
  99             Un=LS.solve()
 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");
 103             for i in range(nx):
 104                 u[i]=Un[i]
 105
 106             time += dt
 107             it += 1
 108
 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
 112
 113             # Postprocessing
 114             line1.set_ydata(u)
 115             writer.grab_frame()
 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")
 120                 #plt.show()
 121
 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)])
 126
 127     print("Simulation of transport equation with implicit upwind scheme done.")
 128
 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)])
 131
 132
 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)
 138     else :
 139         nx = 50 # number of cells
 140         cfl = 0.99 # c*dt/dx <= CFL
 141         Transport1DUpwindImplicit(nx,cfl)
 142