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
11 def initial_conditions_wave_system(my_mesh):
12 dim = my_mesh.getMeshDimension()
13 nbCells = my_mesh.getNumberOfCells()
16 raise ValueError("initial_conditions_wave_system: Mesh dimension should be 1")
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)
22 for i in range(nbCells):
23 x = my_mesh.getCell(i).x()
26 pressure_field[i] = p0
28 velocity_field[i,0] = 1
30 velocity_field[i,0] = -1
33 U[i,1] = rho0*velocity_field[i,0]
35 return U, pressure_field, velocity_field
37 def jacobianMatrices():
39 absA=cdmath.Matrix(2,2)
50 result=cdmath.Vector(2)
51 result[0] = c0*c0*U[1]
56 def numericalFlux(Uj,Ujp1,absA):
61 return Fj+Fjp1 +absA*(Uj-Ujp1)
63 def computeFluxes(U, SumFluxes):
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)
78 A, absA=jacobianMatrices()
80 for j in range(nbCells):#On parcourt les cellules
81 Cj = my_mesh.getCell(j)
83 for i in range(nbComp) :
88 for i in range(nbComp) :
91 elif ( j==nbCells-1) :
92 for i in range(nbComp) :
96 for i in range(nbComp) :
100 Fr=numericalFlux(Uj,Ujp1,absA)
101 Fl=numericalFlux(Ujm1,Uj,absA)
103 sumFluxCourant = (Fr - Fl)*0.5/Cj.getMeasure()
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];
110 def WaveSystem1DVF(ntmax, tmax, cfl, my_mesh, output_freq, resolution):
111 dim=my_mesh.getMeshDimension()
112 nbCells = my_mesh.getNumberOfCells()
119 SumFluxes = cdmath.Field("Fluxes", cdmath.CELLS, my_mesh, dim+1)
121 # Initial conditions #
122 print("Construction of the initial condition …")
123 U, pressure_field, velocity_field = initial_conditions_wave_system(my_mesh)
125 dx_min=my_mesh.minRatioVolSurf()
127 dt = cfl * dx_min / c0
129 print("Starting computation of the linear wave system with an UPWIND explicit scheme …")
132 while (it<ntmax and time <= tmax and not isStationary):
133 computeFluxes(U,SumFluxes);
136 maxVector=SumFluxes.normMax()
137 isStationary= maxVector[0]/p0<precision and maxVector[1]/rho0<precision
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 )
149 for k in range(nbCells):
150 pressure_field[k]=U[k,0]
151 velocity_field[k,0]=U[k,1]/rho0
153 pressure_field.setTime(time,it);
154 pressure_field.writeCSV("WaveSystem1DUpwind_pressure");
155 velocity_field.setTime(time,it);
156 velocity_field.writeCSV("WaveSystem1DUpwind_velocity");
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 )
163 print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint" )
164 raise ValueError("Maximum number of time steps reached : Stationary state not found !!!!!!!")
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
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");
177 print( "Temps maximum Tmax= ", tmax, " atteint" )
178 raise ValueError("Maximum time reached : Stationary state not found !!!!!!!")
181 def solve(my_mesh,resolution):
182 print("Resolution of the 1D Riemann problem for the Wave system with Upwind explicit scheme:")
190 WaveSystem1DVF(ntmax, tmax, cfl, my_mesh, output_freq,resolution)
192 if __name__ == """__main__""":
197 M=cdmath.Mesh(xinf,xsup,10)