CDMATH/tests/examples/HeatEquation1DImplicit/HeatEquation1DImplicit.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 #===============================================================================================================================
   5 # Name        : Résolution VF implicite 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 import cdmath
  22
  23 def implicitSchemeMatrix(nx,cfl):
  24     implMat=cdmath.SparseMatrixPetsc(nx,nx,3)
  25     for i in range(nx):
  26         implMat.setValue(i,(i+1)%nx,-cfl)
  27         implMat.setValue(i,i,1.+2*cfl)
  28         implMat.setValue(i,(i-1)%nx,-cfl)
  29
  30     return implMat
  31
  32 def HeatEquation1DImplicit(nx,cfl):
  33     print( "Simulation of 1D heat equation with an implicit scheme")
  34
  35     ##################### Simulation parameters
  36     a = 0.0 # space domain :  a <= x <= b
  37     b = 1.0
  38     dx = (b - a) / nx #space step
  39
  40     d = 1. # thermal diffusivity
  41     tmax = (b-a)/d # runs the simulation for 0 <= t <= tMax
  42     dt = cfl * dx *dx/ (2*d)
  43     ntmax = 100
  44
  45     x=[a+0.5*dx + i*dx for i in range(nx)]   # array of cell center (1D mesh)
  46
  47     ########################## Initial data
  48
  49     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
  50     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
  51
  52     max_initial=max(u_initial)
  53     min_initial=min(u_initial)
  54     total_var_initial = np.sum([abs(u_initial[i] - u_initial[(i-1)%nx]) for i in range(nx)])
  55
  56     time = 0.
  57     it = 0
  58     output_freq = 10
  59
  60     #Initialisation of the linear system
  61     systemMat=implicitSchemeMatrix(nx,cfl)
  62     iterGMRESMax=50
  63     precision=1.e-5
  64     Un =cdmath.Vector(nx)
  65     for i in range(nx):
  66         Un[i]=u[i]
  67     LS=cdmath.LinearSolver(systemMat,Un,iterGMRESMax, precision, "GMRES","ILU")
  68     LS.setComputeConditionNumber()
  69
  70     # Video settings
  71     FFMpegWriter = manimation.writers['ffmpeg']
  72     metadata = dict(title="Implicit scheme for heat equation, "+"CFL="+str(cfl), artist = "CEA Saclay", comment="Stable for any CFL number")
  73     writer=FFMpegWriter(fps=output_freq, metadata=metadata, codec='h264')
  74     with writer.saving(plt.figure(), "HeatEquation1D_Implicit_"+str(nx)+"Cells_CFL"+str(cfl)+".mp4", ntmax):
  75         ########################### Postprocessing initialisation
  76         # Picture frame
  77         plt.legend()
  78         plt.xlabel('x')
  79         plt.ylabel('u')
  80         plt.xlim(a,b)
  81         plt.ylim( min_initial - 0.1*(max_initial-min_initial), max_initial +  0.1*(max_initial-min_initial) )
  82         plt.title("Implicit scheme for heat equation, "+"CFL="+str(cfl))
  83         line1, = plt.plot(x, u, label='u') #new picture for video # Returns a tuple of line objects, thus the comma
  84
  85         print("Starting time loop")
  86         print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
  87         np.savetxt( "HeatEquation1D_Implicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.txt", u, delimiter="\n")
  88         writer.grab_frame()
  89         plt.savefig("HeatEquation1D_Implicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_0.png")
  90
  91         ############################# Time loop
  92         while (it < ntmax and time <= tmax):
  93             # Solve linear system
  94             for i in range(nx):
  95                 Un[i]=u[i]
  96             LS.setSndMember(Un)
  97             Un=LS.solve()
  98             if(not LS.getStatus()):
  99                 print( "Linear system did not converge ", iterGMRES, " GMRES iterations" )
 100                 raise ValueError("Pas de convergence du système linéaire");
 101             for i in range(nx):
 102                 u[i]=Un[i]
 103
 104             assert max(u) <= max_initial
 105             assert min(u) >= min_initial
 106             assert np.sum([abs(u[i] - u[(i-1)%nx]) for i in range(nx)]) <= total_var_initial
 107
 108             time += dt
 109             it += 1
 110
 111             # Postprocessing
 112             line1.set_ydata(u)
 113             writer.grab_frame()
 114             if (it % output_freq == 0):
 115                 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 116                 np.savetxt( "HeatEquation1D_Implicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".txt", u, delimiter="\n")
 117                 plt.savefig("HeatEquation1D_Implicit_"+str(nx)+"Cells_CFL"+str(cfl)+"_ResultField_"+str(it)+".png")
 118                 #plt.show()
 119
 120     print( "Exact solution minimum   : ", min(u_initial), "Numerical solution minimum   : ",  min(u) )
 121     print( "Exact solution maximum   : ", max(u_initial), "Numerical solution maximum   : ",  max(u) )
 122     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)]) )
 123     print( "l1 numerical error       : ", dx*np.sum([abs(u[i] - u_initial[i]) for i in range(nx)]) )
 124
 125     print("Simulation of heat equation with implicit scheme done.")
 126
 127     #return min, max, total variation and l1 error
 128     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)])
 129
 130
 131 if __name__ == """__main__""":
 132     if len(sys.argv) >2 :
 133         nx = int(sys.argv[1])
 134         cfl = float(sys.argv[2])
 135         HeatEquation1DImplicit(nx,cfl)
 136     else :
 137         nx = 50 # number of cells
 138         cfl = 1 # c*dt/(2*dx) <= CFL
 139         HeatEquation1DImplicit(nx,cfl)
 140