CDMATH/tests/examples/WaveSystem1DUpwind_RiemannProblem/WaveSystem1DUpwind_RiemannProblem.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 import cdmath
   5 import numpy as np
   6 import matplotlib
   7 matplotlib.use("Agg")
   8 import matplotlib.pyplot as plt
   9 import matplotlib.animation as manimation
  10 import sys
  11
  12 p0=155.e5#reference pressure in a pressurised nuclear vessel
  13 c0=700.#reference sound speed for water at 155 bars
  14 rho0=p0/c0*c0#reference density
  15 precision=1e-5
  16
  17 def initial_conditions_Riemann_problem(a,b,nx):
  18     print( "Initial data Riemann problem" )
  19
  20     dx = (b - a) / nx #space step
  21     x=[a+0.5*dx + i*dx for i in range(nx)]   # array of cell center (1D mesh)
  22
  23     u_initial = [ 0 ]*nx
  24     p_initial = [ (xi<(a+b)/2)*p0 + (xi>=(a+b)/2)*p0/2  for xi in x]
  25
  26     return p_initial, u_initial
  27
  28 def jacobianMatrices(coeff,scaling):
  29     dim=1
  30     A=cdmath.Matrix(dim+1,dim+1)
  31     absA=cdmath.Matrix(dim+1,dim+1)
  32
  33     absA[0,0]=c0*coeff
  34     for i in range(dim):
  35         for j in range(dim):
  36             absA[i+1,j+1]=c0*coeff
  37         if( scaling==0):
  38             A[0,i+1]=c0*c0*coeff
  39             A[i+1,0]=      coeff
  40         elif( scaling==1):
  41             A[0,i+1]=      coeff
  42             A[i+1,0]=      coeff
  43         else:
  44             A[0,i+1]=   c0*coeff
  45             A[i+1,0]=   c0*coeff
  46
  47     return A,absA
  48
  49
  50 def computeUpwindDivergenceMatrix(a,b,nx,nbVoisinsMax,dt,scaling):
  51     nbCells = nx
  52     dx=(b-a)/nx
  53     dim=1
  54     nbComp=dim+1
  55     normal=cdmath.Vector(dim)
  56
  57     implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
  58
  59     A,absA= jacobianMatrices(dt/dx,scaling)
  60     for j in range(nbCells):#On parcourt les cellules
  61         if ( j==0) :
  62             implMat.addValue(j*nbComp,(j+1)*nbComp,(A-absA)/2)
  63             implMat.addValue(j*nbComp,    j*nbComp,(A-absA)/2*(-1.))
  64         elif ( j==nbCells-1) :
  65             implMat.addValue(j*nbComp,    j*nbComp,(A+absA)/2)
  66             implMat.addValue(j*nbComp,(j-1)*nbComp,(A+absA)/2*(-1.))
  67         else :
  68             implMat.addValue(j*nbComp,(j+1)*nbComp,(A-absA)/2)
  69             implMat.addValue(j*nbComp,    j*nbComp,(A-absA)/2*(-1.))
  70
  71             implMat.addValue(j*nbComp,    j*nbComp,(A+absA)/2)
  72             implMat.addValue(j*nbComp,(j-1)*nbComp,(A+absA)/2*(-1.))
  73
  74     return implMat
  75
  76 def WaveSystemVF(ntmax, tmax, cfl, a,b,nx, output_freq, meshName,scaling):
  77     dim=1
  78     nbCells = nx
  79
  80     dt = 0.
  81     time = 0.
  82     it=0;
  83     isStationary=False
  84
  85     dx=(b-a)/nx
  86     dt = cfl * dx / c0
  87
  88     nbVoisinsMax=2
  89
  90     #iteration vectors
  91     Un =cdmath.Vector(nbCells*(dim+1))
  92     dUn=cdmath.Vector(nbCells*(dim+1))
  93
  94     Un_implicite =cdmath.Vector(nbCells*(dim+1))
  95     dUn_implicite=cdmath.Vector(nbCells*(dim+1))
  96
  97     # Initial conditions #
  98     print("Construction of the initial condition …")
  99     pressure_field, velocity_field = initial_conditions_Riemann_problem(a,b,nx)
 100     pressure_field_implicite, velocity_field_implicite = initial_conditions_Riemann_problem(a,b,nx)
 101     max_initial_p=max(pressure_field)
 102     min_initial_p=min(pressure_field)
 103     max_initial_v=max(velocity_field)
 104     min_initial_v=min(velocity_field)
 105
 106     for k in range(nbCells):
 107         Un[k*(dim+1)+0] =     pressure_field[k]
 108         Un[k*(dim+1)+1] =rho0*velocity_field[k]
 109         Un_implicite[k*(dim+1)+0] =     pressure_field_implicite[k]
 110         Un_implicite[k*(dim+1)+1] =rho0*velocity_field_implicite[k]
 111
 112     # Video settings
 113     FFMpegWriter = manimation.writers['ffmpeg']
 114     metadata = dict(title="Finite volumes schemes for the 2D Wave System", artist = "CEA Saclay", comment="Shock propagation")
 115     writer=FFMpegWriter(fps=10, metadata=metadata, codec='h264')
 116     with writer.saving(plt.figure(), "2DWaveSystem_Upwind"+".mp4", ntmax):
 117         #sauvegarde de la donnée initiale
 118         plt.xlabel('x (m)')
 119         plt.ylabel('Pressure -Pa)')
 120         plt.xlim(a,b)
 121         plt.ylim( min_initial_p - 0.1*(max_initial_p-min_initial_p), max_initial_p +  0.1*(max_initial_p-min_initial_p) )
 122         plt.title("Riemann problem for Wave system on " + str(nx) + " cells")
 123         line1, = plt.plot([a+0.5*dx + i*dx for i in range(nx)], pressure_field, label='Explicit upwind scheme') #new picture for video # Returns a tuple of line objects, thus the comma
 124         line3, = plt.plot([a+0.5*dx + i*dx for i in range(nx)], pressure_field_implicite, label='Implicit upwind scheme') #new picture for video # Returns a tuple of line objects, thus the comma
 125         plt.legend()
 126         writer.grab_frame()
 127         np.savetxt( "WaveSystem"+str(dim)+"DUpwindExplicit"+meshName+"_pressure"+"_0"+".txt", pressure_field, delimiter="\n")
 128         np.savetxt( "WaveSystem"+str(dim)+"DUpwindExplicit"+meshName+"_velocity"+"_0"+".txt", velocity_field, delimiter="\n")
 129         np.savetxt( "WaveSystem"+str(dim)+"DUpwindImplicit"+meshName+"_pressure"+"_0"+".txt", pressure_field_implicite, delimiter="\n")
 130         np.savetxt( "WaveSystem"+str(dim)+"DUpwindImplicit"+meshName+"_velocity"+"_0"+".txt", velocity_field_implicite, delimiter="\n")
 131         plt.savefig("WaveSystem"+str(dim)+"DUpwind"        +meshName+"_pressure"+"_0"+".png")
 132
 133         divMat=computeUpwindDivergenceMatrix(a,b,nx,nbVoisinsMax,dt,scaling)
 134         divMat_implicit=computeUpwindDivergenceMatrix(a,b,nx,nbVoisinsMax,dt,scaling)
 135
 136         iterGMRESMax=50
 137
 138         divMat_implicit.diagonalShift(1)#only after  filling all coefficients
 139         if( scaling==0):
 140             LS=cdmath.LinearSolver(divMat_implicit,Un,iterGMRESMax, precision, "GMRES","ILU")
 141         else:
 142             LS=cdmath.LinearSolver(divMat_implicit,Vn,iterGMRESMax, precision, "GMRES","ILU")
 143
 144         print("Starting computation of the linear wave system with an upwind scheme …")
 145
 146         # Starting time loop
 147         while (it<ntmax and time <= tmax and not isStationary):
 148             dUn=divMat*Un
 149             Un-=dUn
 150
 151             dUn_implicite=Un_implicite.deepCopy()
 152             LS.setSndMember(Un_implicite)
 153             Un_implicite=LS.solve();
 154             if(not LS.getStatus()):
 155                 print( "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations" )
 156                 raise ValueError("Pas de convergence du système linéaire");
 157             dUn_implicite-=Un_implicite
 158
 159             for k in range(nbCells):
 160                 pressure_field[k] = Un[k*(dim+1)+0]
 161                 velocity_field[k] = Un[k*(dim+1)+1] / rho0
 162                 pressure_field_implicite[k] = Un_implicite[k*(dim+1)+0]
 163                 velocity_field_implicite[k] = Un_implicite[k*(dim+1)+1] / rho0
 164
 165             line1.set_ydata(pressure_field)
 166             line3.set_ydata(pressure_field_implicite)
 167             writer.grab_frame()
 168
 169             time=time+dt;
 170             it=it+1;
 171
 172             #Sauvegardes
 173             if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 174                 print( "-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
 175                 print( "Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations" )
 176
 177                 np.savetxt("WaveSystem" +str(dim)+"DUpwindExplicit"+meshName+"_pressure"+str(it)+".txt", pressure_field          , delimiter="\n")
 178                 np.savetxt("WaveSystem" +str(dim)+"DUpwindExplicit"+meshName+"_velocity"+str(it)+".txt", velocity_field          , delimiter="\n")
 179                 np.savetxt("WaveSystem" +str(dim)+"DUpwindImplicit"+meshName+"_pressure"+str(it)+".txt", pressure_field          , delimiter="\n")
 180                 np.savetxt("WaveSystem" +str(dim)+"DUpwindImplicit"+meshName+"_velocity"+str(it)+".txt", velocity_field          , delimiter="\n")
 181                 plt.savefig("WaveSystem"+str(dim)+"DUpwind"        +meshName+"_pressure"+str(it)+".png")
 182
 183                 print()
 184     print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
 185
 186     if(it>=ntmax):
 187         print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint" )
 188         return
 189     elif(isStationary):
 190         print( "Régime stationnaire atteint au pas de temps ", it, ", t= ", time )
 191         print( "------------------------------------------------------------------------------------")
 192
 193         np.savetxt( "WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat.txt", pressure_field, delimiter="\n")
 194         np.savetxt( "WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat.txt", velocity_field, delimiter="\n")
 195         plt.savefig("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat.png")
 196
 197         return
 198     else:
 199         print( "Temps maximum Tmax= ", tmax, " atteint")
 200         return
 201
 202
 203 def solve( a,b,nx, meshName, scaling, meshType, cfl):
 204
 205     print( "Resolution of the Wave system in dimension 1 on "+str(nx)+ " cells, upwind scheme")
 206     print( "Initial data : ", "Riemann problem")
 207     print( "Boundary conditions : ", "Neumann")
 208     print( "Mesh name : ",meshName , ", ", nx, " cells")
 209
 210     # Problem data
 211     tmax = 10000.
 212     ntmax = 50
 213     output_freq = 1
 214
 215     WaveSystemVF(ntmax, tmax, cfl, a,b,nx, output_freq, meshName, scaling)
 216
 217     return
 218
 219
 220 if __name__ == """__main__""":
 221     a=0.
 222     b=1.
 223     nx=100
 224     cfl=0.99
 225     scaling=0
 226     solve( a,b,nx,"SquareRegularSquares",scaling,"RegularSquares",cfl)