CDMATH/tests/examples/HeatEquation1DExplicit/HeatEquation1DExplicit.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 #===============================================================================================================================
   5 # Name        : Résolution VF de l'équation de la chaleur 1D \partial_t u = d \partial_xx u avec conditions aux limites périodiques
   6 # Author      : Michaël Ndjinga
   7 # Copyright   : CEA Saclay 2019
   8 # Description : 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 from math import sin, pi, ceil
  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 copy import deepcopy
  22
  23 def HeatEquation1DExplicit(nx,cfl):
  24     print "Simulation of 1D heat equation with an explicit scheme"
  25
  26     ##################### Simulation parameters
  27     a = 0.0 # space domain :  a <= x <= b
  28     b = 1.0
  29     dx = (b - a) / nx #space step
  30
  31     d = 1. # thermal diffusivity
  32     tmax = (b-a)/d # runs the simulation for 0 <= t <= tMax
  33     dt = cfl * dx *dx/ (2*d)
  34     ntmax = 100
  35
  36     x=[a+0.5*dx + i*dx for i in range(nx)]   # array of cell center (1D mesh)
  37
  38     ########################## Initial data
  39
  40     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
  41     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
  42
  43     max_initial=max(u_initial)
  44     min_initial=min(u_initial)
  45     total_var_initial = np.sum([abs(u_initial[i] - u_initial[(i-1)%nx]) for i in range(nx)])
  46
  47     time = 0.
  48     it = 0
  49     output_freq = 10
  50
  51     # Video settings
  52     FFMpegWriter = manimation.writers['ffmpeg']
  53     metadata = dict(title="Explicit scheme for heat equation", artist = "CEA Saclay", comment="CFL="+str(cfl)+", Stable if CFL<1")
  54     writer=FFMpegWriter(fps=output_freq, metadata=metadata, codec='h264')
  55     with writer.saving(plt.figure(), "HeatEquation1D_Explicit_"+str(nx)+"Cells_CFL"+str(cfl)+".mp4", ntmax):
  56         ########################### Postprocessing initialisation
  57         # Picture frame
  58         plt.legend()
  59         plt.xlabel('x')
  60         plt.ylabel('u')
  61         plt.xlim(a,b)
  62         plt.ylim( min_initial - 0.1*(max_initial-min_initial), max_initial +  0.1*(max_initial-min_initial) )
  63         plt.title("Explicit scheme for the heat equation, CFL="+str(cfl))
  64         line1, = plt.plot(x, u, label='u') #new picture for video # Returns a tuple of line objects, thus the comma
  65
  66         print("Starting time loop")
  67         print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
  68         np.savetxt( "HeatEquation1D_Explicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.txt", u, delimiter="\n")
  69         writer.grab_frame()
  70         plt.savefig("HeatEquation1D_Explicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.png")
  71
  72         ############################# Time loop
  73         while (it < ntmax and time <= tmax):
  74             un=deepcopy(u)
  75             for i in range(nx):
  76                 u[i] = un[i] + d * dt / (dx*dx) * (un[(i+1)%nx] -2*un[i] + un[(i-1)%nx])/2
  77
  78             time += dt
  79             it += 1
  80
  81             if cfl<1 :
  82                 assert max(u) <= max_initial
  83                 assert min(u) >= min_initial
  84                 assert np.sum([abs(u[i] - u[(i-1)%nx]) for i in range(nx)]) <= total_var_initial
  85
  86             # Postprocessing
  87             line1.set_ydata(u)
  88             writer.grab_frame()
  89             if (it % output_freq == 0):
  90                 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
  91                 np.savetxt( "HeatEquation1D_Explicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".txt", u, delimiter="\n")
  92                 plt.savefig("HeatEquation1D_Explicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".png")
  93                 #plt.show()
  94
  95     print "Exact solution minimum   : ", min(u_initial), "Numerical solution minimum   : ",  min(u)
  96     print "Exact solution maximum   : ", max(u_initial), "Numerical solution maximum   : ",  max(u)
  97     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)])
  98     print "l1 numerical error       : ", dx*np.sum([abs(u[i] - u_initial[i]) for i in range(nx)])
  99
 100     print("Simulation of heat equation with explicit scheme done.")
 101
 102     #return min, max, total variation and l1 error
 103     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)])
 104
 105
 106 if __name__ == """__main__""":
 107     if len(sys.argv) >2 :
 108         nx = int(sys.argv[1])
 109         cfl = float(sys.argv[2])
 110         HeatEquation1DExplicit(nx,cfl)
 111     else :
 112         nx = 50 # number of cells
 113         cfl = 1 # c*dt/(2*dx) <= CFL
 114         HeatEquation1DExplicit(nx,cfl)
 115