CDMATH/tests/examples/TransportEquation1DCenteredImplicit/1DTransportEquationCenteredImplicit.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 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 #================================================================================================================================
  13
  14
  15 import numpy as np
  16 import matplotlib
  17 matplotlib.use("Agg")
  18 import matplotlib.pyplot as plt
  19 import matplotlib.animation as manimation
  20 import sys
  21 from math import sin, pi, ceil
  22 import cdmath
  23
  24 def centeredSchemeMatrix(nx,cfl):
  25     centeredMat=cdmath.SparseMatrixPetsc(nx,nx,3)
  26     for i in range(nx):
  27         centeredMat.setValue(i,(i+1)%nx,cfl/2)
  28         centeredMat.setValue(i,i,1.)
  29         centeredMat.setValue(i,(i-1)%nx,-cfl/2.)
  30
  31     return centeredMat
  32
  33 def Transport1DCenteredImplicit(nx,cfl):
  34     print "Simulation of 1D transport equation with implicit centered 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=centeredSchemeMatrix(nx,cfl)
  66     iterGMRESMax=50
  67     precision=1.e-5
  68     Un =cdmath.Vector(nx)
  69     for i in range(nx):
  70         Un[i]=u[i]
  71     LS=cdmath.LinearSolver(systemMat,Un,iterGMRESMax, precision, "GMRES","ILU")
  72     LS.setComputeConditionNumber()
  73
  74     # Video settings
  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
  80         # Picture frame
  81         plt.legend()
  82         plt.xlabel('x')
  83         plt.ylabel('u')
  84         plt.xlim(a,b)
  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
  88
  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")
  92         writer.grab_frame()
  93         plt.savefig("TransportEquation_CenteredImplicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.png")
  94
  95         ############################# Time loop
  96         while (it < ntmax and time <= tmax):
  97             # Solve linear system
  98             for i in range(nx):
  99                 Un[i]=u[i]
 100             LS.setSndMember(Un)
 101             Un=LS.solve()
 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");
 105             for i in range(nx):
 106                 u[i]=Un[i]
 107
 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
 114
 115             time += dt
 116             it += 1
 117
 118             # Postprocessing
 119             line1.set_ydata(u)
 120             writer.grab_frame()
 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")
 125                 #plt.show()
 126
 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)])
 131
 132     print("Simulation of transport equation with implicit centered scheme done.")
 133
 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)])
 136
 137
 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)
 143     else :
 144         nx = 50 # number of cells
 145         cfl = 0.99 # c*dt/dx <= CFL
 146         Transport1DCenteredImplicit(nx,cfl)
 147