4 #===============================================================================================================================
5 # Name : Résolution VF du système des ondes 2D sans terme source
6 # \partial_t p + c^2 \div q = 0
7 # \partial_t q + \grad p = 0
8 # Author : Michaël Ndjinga
9 # Copyright : CEA Saclay 2019
10 # Description : Propagation d'une onde de choc sphérique
11 # Utilisation du schéma centré (implicite) sur un maillage général
12 # Initialisation par une surpression sphérique
13 # Conditions aux limites parois
14 # Création et sauvegarde du champ résultant et des figures
15 #================================================================================================================================
19 from numpy import sign
25 p0=155e5#reference pressure in a pressurised nuclear vessel
26 c0=700.#reference sound speed for water at 155 bars, 600K
27 rho0=p0/c0*c0#reference density
30 def initial_conditions_shock(my_mesh, isCircle):
31 print "Initial data : Spherical wave"
32 dim = my_mesh.getMeshDimension()
33 nbCells = my_mesh.getNumberOfCells()
45 pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
46 velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
48 for i in range(nbCells):
49 velocity_field[i,0] = 0
50 velocity_field[i,1] = 0
51 velocity_field[i,2] = 0
53 x = my_mesh.getCell(i).x()
54 valX = (x - xcentre) * (x - xcentre)
59 y = my_mesh.getCell(i).y()
60 valY = (y - ycentre) * (y - ycentre)
61 val = sqrt(valX + valY)
63 y = my_mesh.getCell(i).y()
64 z = my_mesh.getCell(i).z()
65 valY = (y - ycentre) * (y - ycentre)
66 valZ = (z - zcentre) * (z - zcentre)
67 val = sqrt(valX + valY + valZ)
70 pressure_field[i] = p0
73 pressure_field[i] = p0/2
77 return pressure_field, velocity_field
79 def jacobianMatrices(normal, coeff, signun):
81 A=cdmath.Matrix(dim+1,dim+1)
84 A[i+1,0]=normal[i]*coeff
85 A[0,i+1]=c0*c0*normal[i]*coeff
90 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt):
91 nbCells = my_mesh.getNumberOfCells()
92 dim=my_mesh.getMeshDimension()
94 normal=cdmath.Vector(dim)
96 implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
98 idMoinsJacCL=cdmath.Matrix(nbComp)
100 v0=cdmath.Vector(dim)
101 for i in range(dim) :
104 for j in range(nbCells):#On parcourt les cellules
105 Cj = my_mesh.getCell(j)
106 nbFaces = Cj.getNumberOfFaces();
108 for k in range(nbFaces) :
109 indexFace = Cj.getFacesId()[k];
110 Fk = my_mesh.getFace(indexFace);
111 for i in range(dim) :
112 normal[i] = Cj.getNormalVector(k, i);#normale sortante
114 signun=sign(normal*v0)
115 Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure(),signun);
118 if ( not Fk.isBorder()) :
119 # hypothese: La cellule d'index indexC1 est la cellule courante index j */
120 if (Fk.getCellsId()[0] == j) :
122 cellAutre = Fk.getCellsId()[1];
123 elif(Fk.getCellsId()[1] == j) :
124 # hypothese non verifiée
125 cellAutre = Fk.getCellsId()[0];
127 raise ValueError("computeFluxes: problem with mesh, unknown cel number")
129 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
130 implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
132 if( Fk.getGroupName() != "Periodic" and Fk.getGroupName() != "Neumann"):#Wall boundary condition unless Periodic/Neumann specified explicitly
133 v=cdmath.Vector(dim+1)
134 for i in range(dim) :
136 idMoinsJacCL=v.tensProduct(v)*2
138 implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL)
140 elif( Fk.getGroupName() == "Periodic"):#Periodic boundary condition
141 indexFP=my_mesh.getIndexFacePeriodic(indexFace)
142 Fp = my_mesh.getFace(indexFP)
143 cellAutre = Fp.getCellsId()[0]
145 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
146 implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
147 elif(Fk.getGroupName() != "Neumann"):#Nothing to do for Neumann boundary condition
148 print Fk.getGroupName()
149 raise ValueError("computeFluxes: Unknown boundary condition name");
153 def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution):
154 dim=my_mesh.getMeshDimension()
155 nbCells = my_mesh.getNumberOfCells()
156 meshName=my_mesh.getName()
163 nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
167 Un=cdmath.Vector(nbCells*(dim+1))
168 dUn=cdmath.Vector(nbCells*(dim+1))
170 # Initial conditions #
171 print("Construction of the initial condition …")
172 if(dim==1 or filename.find("square")>-1 or filename.find("Square")>-1 or filename.find("cube")>-1 or filename.find("Cube")>-1):
173 pressure_field, velocity_field = initial_conditions_shock(my_mesh,False)
174 elif(filename.find("disk")>-1 or filename.find("Disk")>-1):
175 pressure_field, velocity_field = initial_conditions_shock(my_mesh,True)
177 print "Mesh name : ", filename
178 raise ValueError("Mesh name should contain substring square, cube or disk")
180 for k in range(nbCells):
181 Un[k*(dim+1)+0] = pressure_field[k]
182 Un[k*(dim+1)+1] = rho0*velocity_field[k,0]
184 Un[k + 2*nbCells] = rho0*velocity_field[k,1] # value on the bottom face
186 Un[k + 3*nbCells] = rho0*velocity_field[k,2]
188 #sauvegarde de la donnée initiale
189 pressure_field.setTime(time,it);
190 pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure");
191 velocity_field.setTime(time,it);
192 velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity");
193 #Postprocessing : save 2D picture
194 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_initial")
195 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_initial")
197 dx_min=my_mesh.minRatioVolSurf()
199 dt = cfl * dx_min / c0
201 divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
203 # Add the identity matrix on the diagonal
204 for j in range(nbCells*(dim+1)):
205 divMat.addValue(j,j,1.)
206 LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","LU")
208 print("Starting computation of the linear wave system with a centered scheme …")
211 while (it<ntmax and time <= tmax and not isStationary):
215 cvgceLS=LS.getStatus();
216 iterGMRES=LS.getNumberOfIter();
218 print "Linear system did not converge ", iterGMRES, " GMRES iterations"
219 raise ValueError("Pas de convergence du système linéaire");
222 maxVector=dUn.maxVector(dim+1)
223 isStationary= maxVector[0]/p0<precision and maxVector[1]/rho0<precision and maxVector[2]/rho0<precision;
225 isStationary=isStationary and maxVector[3]/rho0<precision
230 if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
231 print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
232 print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0
233 print "Linear system converged in ", iterGMRES, " GMRES iterations"
235 for k in range(nbCells):
236 pressure_field[k]=Un[k*(dim+1)+0]
237 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
239 velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
241 velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
243 pressure_field.setTime(time,it);
244 pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure",False);
245 velocity_field.setTime(time,it);
246 velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity",False);
248 print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
249 print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0
253 print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
255 print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
256 print "------------------------------------------------------------------------------------"
258 pressure_field.setTime(time,0);
259 pressure_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_Stat");
260 velocity_field.setTime(time,0);
261 velocity_field.writeVTK("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_Stat");
263 #Postprocessing : save 2D picture
264 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DCentered"+meshName+"_pressure_Stat")
265 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DCentered"+meshName+"_velocity_Stat")
268 print "Temps maximum Tmax= ", tmax, " atteint"
271 def solve(my_mesh,meshName,resolution):
272 print "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension()
273 print "Numerical method : implicit centered"
274 print "Initial data : spherical wave"
275 print "Wall boundary conditions"
276 print "Mesh name : ",meshName , my_mesh.getNumberOfCells(), " cells"
281 cfl = 1./my_mesh.getSpaceDimension()
284 WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution)
286 def solve_file( filename,meshName, resolution):
287 my_mesh = cdmath.Mesh(filename+".med")
289 return solve(my_mesh, filename+str(my_mesh.getNumberOfCells()),resolution)
292 if __name__ == """__main__""":
293 if len(sys.argv) >1 :
295 my_mesh = cdmath.Mesh(filename)
296 solve(my_mesh,filename,100)
298 raise ValueError("WaveSystemUpwind.py expects a mesh file name")