SALOME platform Git repositories - tools/solverlab.git/blob

   1 #!/usr/bin/env python3
   2 # -*-coding:utf-8 -*
   3
   4 # Isothermal Euler system
   5 # d rho/d t + d q/d x =0
   6 # d q/d t + d (q^2/rho+p)/d x = 0, where q = rho*u and p = c^2*rho
   7 # UU = (rho,q) : conservative variable
   8 # Scheme : Conservative stagerred scheme (Ait Ameur et al)
   9 # Author :  Katia Ait Ameur
  10 # Date : November 2020
  11
  12 import cdmath
  13 import numpy as np
  14 import matplotlib
  15 matplotlib.use("Agg")
  16 import matplotlib.pyplot as plt
  17 import matplotlib.animation as manimation
  18 from math import sqrt
  19
  20 c0=330.#reference sound speed for water at 15 bars
  21 precision=1e-5
  22 rho0=1
  23
  24 def initial_conditions_Riemann_problem(a,b,nx):
  25     print("Initial data Riemann problem")
  26     dx = (b - a) / nx #space step
  27     x=[a+0.5*dx + i*dx for i in range(nx)]   # array of cell center (1D mesh)
  28     rho_initial = [ (xi<(a+b)/2)*rho0        + (xi>=(a+b)/2)*rho0*2       for xi in x]
  29     q_initial   = [ (xi<(a+b)/2)*rho0*(-300) + (xi>=(a+b)/2)*rho0*2*(-300)  for xi in x]
  30
  31     return rho_initial, q_initial
  32
  33 def matrix_coef(rho_l,q_l,rho_r,q_r):
  34     m_coef = -0.5*(q_l-q_r)/rho_r
  35     return m_coef
  36
  37 def jacobianMatricesm(coeff,rho_l,q_l,rho_r,q_r):
  38     RoeMat   = cdmath.Matrix(2,2);
  39     Dmac     = cdmath.Matrix(2,2);
  40
  41     u_l=q_l/rho_l
  42     u_r=q_r/rho_r
  43     Dmac_coef = matrix_coef(rho_l,q_l,rho_r,q_r)
  44     if rho_l<0 or rho_r<0 :
  45         print( "rho_l=",rho_l, " rho_r= ",rho_r)
  46         raise ValueError("Negative density")
  47     u = (u_l*sqrt(rho_l)+u_r*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));
  48     RoeMat[0,0]   = 0
  49     RoeMat[0,1]   = 1
  50     RoeMat[1,0]   = c0*c0 - u*u
  51     RoeMat[1,1]   = 2*u
  52
  53     Dmac[0,0]= abs(u)-u;
  54     Dmac[0,1]= 1;
  55     Dmac[1,0]= -c0*c0-u*u;
  56     Dmac[1,1]= abs(u)+u;
  57
  58     return (RoeMat-Dmac)*coeff*0.5
  59
  60 def jacobianMatricesp(coeff,rho_l,q_l,rho_r,q_r):
  61     RoeMat   = cdmath.Matrix(2,2);
  62     Dmac = cdmath.Matrix(2,2);
  63
  64     u_l=q_l/rho_l
  65     u_r=q_r/rho_r
  66     Dmac_coef = matrix_coef(rho_l,q_l,rho_r,q_r)
  67     if rho_l<0 or rho_r<0 :
  68         print( "rho_l=",rho_l, " rho_r= ",rho_r)
  69         raise ValueError("Negative density")
  70     u = (u_l*sqrt(rho_l)+u_r*sqrt(rho_r))/(sqrt(rho_l)+sqrt(rho_r));
  71     RoeMat[0,0]   = 0
  72     RoeMat[0,1]   = 1
  73     RoeMat[1,0]   = c0*c0 - u*u
  74     RoeMat[1,1]   = 2*u
  75
  76     Dmac[0,0]= abs(u)-u;
  77     Dmac[0,1]= 1;
  78     Dmac[1,0]= -c0*c0-u*u;
  79     Dmac[1,1]= abs(u)+u;
  80
  81     return (RoeMat+Dmac)*coeff*0.5
  82
  83 def EulerSystemStaggered(ntmax, tmax, cfl, a,b,nx, output_freq, meshName):
  84     dim=1
  85     nbComp=dim+1
  86     nbCells = nx
  87     dt = 0.
  88     dx=(b-a)/nx
  89     dt = cfl * dx / c0
  90     nbVoisinsMax=2
  91
  92     time = 0.
  93     it=0;
  94     iterGMRESMax=50
  95     newton_max = 50
  96     isStationary=False
  97
  98     #iteration vectors
  99     Un =cdmath.Vector(nbCells*(dim+1))
 100     dUn=cdmath.Vector(nbCells*(dim+1))
 101     dUk=cdmath.Vector(nbCells*(dim+1))
 102     Rhs=cdmath.Vector(nbCells*(dim+1))
 103
 104     # Initial conditions #
 105     print("Construction of the initial condition …")
 106     rho_field, q_field = initial_conditions_Riemann_problem(a,b,nx)
 107     v_field   = [   q_field[i]/rho_field[i]  for i in range(nx)]
 108     p_field   = [ rho_field[i]*(c0*c0)       for i in range(nx)]
 109
 110     max_initial_rho=max(rho_field)
 111     min_initial_rho=min(rho_field)
 112     max_initial_q=max(q_field)
 113     min_initial_q=min(q_field)
 114     max_initial_p=max(p_field)
 115     min_initial_p=min(p_field)
 116     max_initial_v=max(v_field)
 117     min_initial_v=min(v_field)
 118
 119     for k in range(nbCells):
 120         Un[k*nbComp+0] = rho_field[k]
 121         Un[k*nbComp+1] =   q_field[k]
 122
 123     print("Starting computation of the non linear Euler system with staggered scheme …")
 124     divMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
 125
 126     # Picture settings
 127     fig, ([axDensity, axMomentum],[axVelocity, axPressure]) = plt.subplots(2, 2,sharex=True, figsize=(10,10))
 128     fig.suptitle('Conservative staggered scheme')
 129     lineDensity, = axDensity.plot([a+0.5*dx + i*dx for i in range(nx)], rho_field, label='Density') #new picture for video # Returns a tuple of line objects, thus the comma
 130     axDensity.set(xlabel='x (m)', ylabel='Density')
 131     axDensity.set_xlim(a,b)
 132     axDensity.set_ylim(min_initial_rho - 0.1*(max_initial_rho-min_initial_rho), max_initial_rho +  0.1*(max_initial_rho-min_initial_rho) )
 133     axDensity.legend()
 134     lineMomentum, = axMomentum.plot([a+0.5*dx + i*dx for i in range(nx)], q_field, label='Momentum')
 135     axMomentum.set(xlabel='x (m)', ylabel='Momentum')
 136     axMomentum.set_xlim(a,b)
 137     axMomentum.set_ylim(min_initial_q - 0.1*(max_initial_q-min_initial_q), max_initial_q +  0.1*(max_initial_q-min_initial_q) )
 138     axMomentum.legend()
 139     lineVelocity, = axVelocity.plot([a+0.5*dx + i*dx for i in range(nx)], v_field, label='Velocity')
 140     axVelocity.set(xlabel='x (m)', ylabel='Velocity')
 141     axVelocity.set_xlim(a,b)
 142     axVelocity.set_ylim(min_initial_v - 0.4*abs(min_initial_v), max_initial_v +  0.05*abs(max_initial_v) )
 143     axVelocity.legend()
 144     linePressure, = axPressure.plot([a+0.5*dx + i*dx for i in range(nx)], p_field, label='Pressure')
 145     axPressure.set(xlabel='x (m)', ylabel='Pressure')
 146     axPressure.set_xlim(a,b)
 147     axPressure.set_ylim(min_initial_p - 0.05*abs(min_initial_p), max_initial_p +  0.05*abs(max_initial_p) )
 148     axPressure.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
 149     axPressure.legend()
 150
 151     # Video settings
 152     FFMpegWriter = manimation.writers['ffmpeg']
 153     metadata = dict(title="Conservative staggered scheme for the 1D isothermal Euler System", artist = "CEA Saclay", comment="Shock tube")
 154     writer=FFMpegWriter(fps=10, metadata=metadata, codec='h264')
 155     with writer.saving(fig, "1DEuler_System_ConservativeStaggered"+".mp4", ntmax):
 156         writer.grab_frame()
 157         plt.savefig("EulerSystem"+str(dim)+"D_ConservativeStaggered"+meshName+"_0"+".png")
 158
 159         # Starting time loop
 160         while (it<ntmax and time <= tmax and not isStationary):
 161             dUn = Un.deepCopy()
 162             Uk  = Un.deepCopy()
 163             residu = 1.
 164             k=0
 165             while (k<newton_max and residu > precision ):
 166                 #DEBUT BOUCLE NEWTON
 167                 divMat.zeroEntries()#sets the matrix coefficients to zero
 168                 for j in range(nbCells):#
 169                     if ( j==0) :
 170                         rho_r = Uk[j*nbComp+0]
 171                         q_r   = Uk[j*nbComp+1]
 172                         rho_l = rho_r # Conditions de Neumann
 173                         q_l   =   q_r
 174                         Am= jacobianMatricesm(dt/dx,rho_l,q_l,rho_r,q_r)
 175                         divMat.addValue(j*nbComp,(j+1)*nbComp,Am)
 176                         divMat.addValue(j*nbComp,    j*nbComp,Am*(-1.))
 177                         dUi=cdmath.Vector(2)
 178                         dUi[0] = Uk[(j+1)*nbComp+0]-Uk[j*nbComp+0]
 179                         dUi[1] = Uk[(j+1)*nbComp+1]-Uk[j*nbComp+1]
 180                         temp=cdmath.Vector(2)
 181                         temp = Am*dUi
 182                         #print("Bloc 0 matrix  :   ", Am)
 183                         Rhs[j*nbComp+0] = -temp[0]-(Uk[j*nbComp+0]-Un[j*nbComp+0])
 184                         Rhs[j*nbComp+1] = -temp[1]-(Uk[j*nbComp+1]-Un[j*nbComp+1])
 185                     elif ( j==nbCells-1) :
 186                         rho_l = Uk[j*nbComp+0]
 187                         q_l   = Uk[j*nbComp+1]
 188                         rho_r = rho_l # Conditions de Neumann
 189                         q_r   =   q_l
 190                         Ap= jacobianMatricesp(dt/dx,rho_l,q_l,rho_r,q_r)
 191                         divMat.addValue(j*nbComp,    j*nbComp,Ap)
 192                         divMat.addValue(j*nbComp,(j-1)*nbComp,Ap*(-1.))
 193                         dUi=cdmath.Vector(2)
 194                         dUi[0] = Uk[(j-1)*nbComp+0]-Uk[j*nbComp+0]
 195                         dUi[1] = Uk[(j-1)*nbComp+1]-Uk[j*nbComp+1]
 196                         temp=cdmath.Vector(2)
 197                         temp = Ap*dUi
 198                         Rhs[j*nbComp+0] = temp[0]-(Uk[j*nbComp+0]-Un[j*nbComp+0])
 199                         Rhs[j*nbComp+1] = temp[1]-(Uk[j*nbComp+1]-Un[j*nbComp+1])
 200                     else :
 201                         rho_l = Uk[(j-1)*nbComp+0]
 202                         q_l   = Uk[(j-1)*nbComp+1]
 203                         rho_r = Uk[j*nbComp+0]
 204                         q_r   = Uk[j*nbComp+1]
 205                         Ap = jacobianMatricesp(dt/dx,rho_l,q_l,rho_r,q_r)
 206                         ###############################################################
 207                         rho_l = Uk[j*nbComp+0]
 208                         q_l   = Uk[j*nbComp+1]
 209                         rho_r = Uk[(j+1)*nbComp+0]
 210                         q_r   = Uk[(j+1)*nbComp+1]
 211                         Am = jacobianMatricesm(dt/dx,rho_l,q_l,rho_r,q_r)
 212                         divMat.addValue(j*nbComp,(j+1)*nbComp,Am)
 213                         divMat.addValue(j*nbComp,    j*nbComp,Am*(-1.))
 214                         divMat.addValue(j*nbComp,    j*nbComp,Ap)
 215                         divMat.addValue(j*nbComp,(j-1)*nbComp,Ap*(-1.))
 216                         dUi1=cdmath.Vector(2)
 217                         dUi2=cdmath.Vector(2)
 218                         dUi1[0] = Uk[(j+1)*nbComp+0]-Uk[j*nbComp+0]
 219                         dUi1[1] = Uk[(j+1)*nbComp+1]-Uk[j*nbComp+1]
 220                         dUi2[0] = Uk[(j-1)*nbComp+0]-Uk[j*nbComp+0]
 221                         dUi2[1] = Uk[(j-1)*nbComp+1]-Uk[j*nbComp+1]
 222                         temp1 = cdmath.Vector(2)
 223                         temp2 = cdmath.Vector(2)
 224                         temp1 = Am*dUi1
 225                         temp2 = Ap*dUi2
 226                         Rhs[j*nbComp+0] = -temp1[0] + temp2[0] -(Uk[j*nbComp+0]-Un[j*nbComp+0])
 227                         Rhs[j*nbComp+1] = -temp1[1] + temp2[1] -(Uk[j*nbComp+1]-Un[j*nbComp+1])
 228                 divMat.diagonalShift(1)#only after  filling all coefficients
 229                 LS=cdmath.LinearSolver(divMat,Rhs,iterGMRESMax, precision, "GMRES","LU")
 230                 dUk=LS.solve();
 231                 residu = dUk.norm()
 232                 #stop
 233                 Uk+=dUk
 234                 if(not LS.getStatus()):
 235                     print("Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations")
 236                     raise ValueError("Pas de convergence du système linéaire");
 237                 k=k+1
 238             Un = Uk.deepCopy()
 239             dUn-=Un
 240             for k in range(nbCells):
 241                 rho_field[k] = Un[k*(dim+1)+0]
 242                 q_field[k]   = Un[k*(dim+1)+1]
 243             v_field   = [   q_field[i]/rho_field[i]  for i in range(nx)]
 244             p_field   = [ rho_field[i]*(c0*c0)       for i in range(nx)]
 245
 246             lineDensity.set_ydata(rho_field)
 247             lineMomentum.set_ydata(q_field)
 248             lineVelocity.set_ydata(v_field)
 249             linePressure.set_ydata(p_field)
 250             writer.grab_frame()
 251
 252             time=time+dt;
 253             it=it+1;
 254             #Sauvegardes
 255             if(it==1 or it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 256                 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 257                 #print("Last linear system converged in ", LS.getNumberOfIter(), " GMRES iterations", ", residu final:   ",residu)
 258
 259                 plt.savefig("EulerSystem"+str(dim)+"D_ConservativeStaggered"+meshName+"_"+str(it)+".png")
 260                 print
 261     print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 262     if(it>=ntmax):
 263         print("Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
 264         return
 265     elif(isStationary):
 266         print("Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
 267         print("------------------------------------------------------------------------------------")
 268         plt.savefig("EulerSystem"+str(dim)+"Staggered"+meshName+"_Stat.png")
 269         return
 270     else:
 271         print("Temps maximum Tmax= ", tmax, " atteint")
 272         return
 273
 274 def solve( a,b,nx, meshName, meshType, cfl):
 275     print("Resolution of the Euler system in dimension 1 on "+str(nx)+ " cells")
 276     print("Initial data : ", "Riemann problem")
 277     print("Boundary conditions : ", "Neumann")
 278     print("Mesh name : ",meshName , ", ", nx, " cells")
 279     # Problem data
 280     tmax = 10.
 281     ntmax = 50
 282     output_freq = 10
 283     EulerSystemStaggered(ntmax, tmax, cfl, a,b,nx, output_freq, meshName)
 284     return
 285
 286 if __name__ == """__main__""":
 287     a=0.
 288     b=1.
 289     nx=100
 290     cfl=0.99
 291     solve( a,b,nx,"SquareRegularSquares","RegularSquares",cfl)