CDMATH/tests/examples/WaveSystem1DUpwind/WaveSystem1DUpwind.py

   1 #!/usr/bin/env python
   2 # -*-coding:utf-8 -*
   3
   4 import cdmath
   5
   6 p0=155e5#reference pressure in a pressurised nuclear vessel
   7 c0=700.#reference sound speed for water at 155 bars, 600K
   8 rho0=p0/c0*c0#reference density
   9 precision=1e-5
  10
  11 def initial_conditions_wave_system(my_mesh):
  12     dim     = my_mesh.getMeshDimension()
  13     nbCells = my_mesh.getNumberOfCells()
  14
  15     if(dim!=1):
  16         raise ValueError("initial_conditions_wave_system: Mesh dimension should be 1")
  17
  18     pressure_field = cdmath.Field("Pressure",            cdmath.CELLS, my_mesh, 1)
  19     velocity_field = cdmath.Field("Velocity",            cdmath.CELLS, my_mesh, dim)
  20     U              = cdmath.Field("Conservative vector", cdmath.CELLS, my_mesh, dim+1)
  21
  22     for i in range(nbCells):
  23         x = my_mesh.getCell(i).x()
  24
  25
  26         pressure_field[i] = p0
  27         if(x>0.5):
  28             velocity_field[i,0] =   1
  29         else:
  30             velocity_field[i,0] =  -1
  31
  32         U[i,0] =   p0
  33         U[i,1] =  rho0*velocity_field[i,0]
  34
  35     return U, pressure_field, velocity_field
  36
  37 def jacobianMatrices():
  38     A=cdmath.Matrix(2,2)
  39     absA=cdmath.Matrix(2,2)
  40
  41     absA[0,0]=c0
  42     absA[1,1]=c0
  43     A[1,0]=1
  44     A[0,1]=c0*c0
  45
  46     return A, absA
  47
  48 def Flux(U):
  49
  50     result=cdmath.Vector(2)
  51     result[0] = c0*c0*U[1]
  52     result[1] = U[0]
  53
  54     return result
  55
  56 def numericalFlux(Uj,Ujp1,absA):
  57
  58     Fj   = Flux(Uj)
  59     Fjp1 = Flux(Ujp1)
  60
  61     return Fj+Fjp1 +absA*(Uj-Ujp1)
  62
  63 def computeFluxes(U, SumFluxes):
  64     my_mesh =U.getMesh();
  65     nbCells = my_mesh.getNumberOfCells();
  66     dim=my_mesh.getMeshDimension();
  67     nbComp=U.getNumberOfComponents();
  68     Fj=cdmath.Vector(nbComp)
  69     Fjp1=cdmath.Vector(nbComp)
  70     Fjm1=cdmath.Vector(nbComp)
  71     Uj=cdmath.Vector(nbComp)
  72     Ujp1=cdmath.Vector(nbComp)
  73     Ujm1=cdmath.Vector(nbComp)
  74     normal=cdmath.Vector(dim)
  75     sumFluxCourant=cdmath.Vector(nbComp)
  76     sumFluxCourant2=cdmath.Vector(nbComp)
  77
  78     A, absA=jacobianMatrices()
  79
  80     for j in range(nbCells):#On parcourt les cellules
  81         Cj = my_mesh.getCell(j)
  82
  83         for i in range(nbComp) :
  84             Uj[i]=U[j,i];
  85             sumFluxCourant[i]=0;
  86
  87         if ( j==0) :
  88             for i in range(nbComp) :
  89                 Ujp1[i]=U[j+1,i];
  90                 Ujm1[i]=U[j  ,i];
  91         elif ( j==nbCells-1) :
  92             for i in range(nbComp) :
  93                 Ujp1[i]=U[j  ,i];
  94                 Ujm1[i]=U[j-1,i];
  95         else :
  96             for i in range(nbComp) :
  97                 Ujp1[i]=U[j+1,i];
  98                 Ujm1[i]=U[j-1,i];
  99
 100         Fr=numericalFlux(Uj,Ujp1,absA)
 101         Fl=numericalFlux(Ujm1,Uj,absA)
 102
 103         sumFluxCourant = (Fr - Fl)*0.5/Cj.getMeasure()
 104
 105         #On divise par le volume de la cellule la contribution des flux au snd membre
 106         for i in range(nbComp):
 107             SumFluxes[j,i]=sumFluxCourant[i];
 108
 109
 110 def WaveSystem1DVF(ntmax, tmax, cfl, my_mesh, output_freq, resolution):
 111     dim=my_mesh.getMeshDimension()
 112     nbCells = my_mesh.getNumberOfCells()
 113
 114     dt = 0.
 115     time = 0.
 116     it=0;
 117     isStationary=False;
 118
 119     SumFluxes = cdmath.Field("Fluxes", cdmath.CELLS, my_mesh, dim+1)
 120
 121     # Initial conditions #
 122     print("Construction of the initial condition …")
 123     U, pressure_field, velocity_field = initial_conditions_wave_system(my_mesh)
 124
 125     dx_min=my_mesh.minRatioVolSurf()
 126
 127     dt = cfl * dx_min / c0
 128
 129     print("Starting computation of the linear wave system with an UPWIND explicit scheme …")
 130
 131     # Starting time loop
 132     while (it<ntmax and time <= tmax and not isStationary):
 133         computeFluxes(U,SumFluxes);
 134
 135         SumFluxes*=dt;
 136         maxVector=SumFluxes.normMax()
 137         isStationary= maxVector[0]/p0<precision and maxVector[1]/rho0<precision
 138         U-=SumFluxes;
 139
 140         time=time+dt;
 141         it=it+1;
 142
 143         #Sauvegardes
 144         if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
 145             print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 146             print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0
 147             print
 148
 149             for k in range(nbCells):
 150                 pressure_field[k]=U[k,0]
 151                 velocity_field[k,0]=U[k,1]/rho0
 152
 153             pressure_field.setTime(time,it);
 154             pressure_field.writeCSV("WaveSystem1DUpwind_pressure");
 155             velocity_field.setTime(time,it);
 156             velocity_field.writeCSV("WaveSystem1DUpwind_velocity");
 157
 158     print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
 159     print "|| Un+1 - Un || : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0
 160     print
 161
 162     if(it>=ntmax):
 163         print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
 164         raise ValueError("Maximum number of time steps reached : Stationary state not found !!!!!!!")
 165     elif(isStationary):
 166         print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
 167         for k in range(nbCells):
 168             pressure_field[k]=U[k,0]
 169             velocity_field[k,0]=U[k,1]/rho0
 170
 171         pressure_field.setTime(time,0);
 172         pressure_field.writeCSV("WaveSystem1DUpwind_pressure_Stat");
 173         velocity_field.setTime(time,0);
 174         velocity_field.writeCSV("WaveSystem1DUpwind_velocity_Stat");
 175
 176     else:
 177         print "Temps maximum Tmax= ", tmax, " atteint"
 178         raise ValueError("Maximum time reached : Stationary state not found !!!!!!!")
 179
 180
 181 def solve(my_mesh,resolution):
 182     print("Resolution of the 1D Riemann problem for the Wave system with Upwind explicit scheme:")
 183
 184     # Problem data
 185     tmax = 1.
 186     ntmax = 100
 187     cfl = 0.95
 188     output_freq = 10
 189
 190     WaveSystem1DVF(ntmax, tmax, cfl, my_mesh, output_freq,resolution)
 191
 192 if __name__ == """__main__""":
 193
 194     xinf=0
 195     xsup=1
 196
 197     M=cdmath.Mesh(xinf,xsup,10)
 198
 199     solve(M,100)