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 MAC sur un maillage cartésien
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)
82 absA=cdmath.Matrix(dim+1,dim+1)
85 A[i+1,0]=normal[i]*coeff
86 absA[i+1,0]=-signun*A[i+1,0]
87 A[0,i+1]=c0*c0*normal[i]*coeff
88 absA[0,i+1]=signun*A[0,i+1]
93 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt):
94 nbCells = my_mesh.getNumberOfCells()
95 dim=my_mesh.getMeshDimension()
97 normal=cdmath.Vector(dim)
99 implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
101 idMoinsJacCL=cdmath.Matrix(nbComp)
103 v0=cdmath.Vector(dim)
104 for i in range(dim) :
107 for j in range(nbCells):#On parcourt les cellules
108 Cj = my_mesh.getCell(j)
109 nbFaces = Cj.getNumberOfFaces();
111 for k in range(nbFaces) :
112 indexFace = Cj.getFacesId()[k];
113 Fk = my_mesh.getFace(indexFace);
114 for i in range(dim) :
115 normal[i] = Cj.getNormalVector(k, i);#normale sortante
117 signun=sign(normal*v0)
118 Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure(),signun);
121 if ( not Fk.isBorder()) :
122 # hypothese: La cellule d'index indexC1 est la cellule courante index j */
123 if (Fk.getCellsId()[0] == j) :
125 cellAutre = Fk.getCellsId()[1];
126 elif(Fk.getCellsId()[1] == j) :
127 # hypothese non verifiée
128 cellAutre = Fk.getCellsId()[0];
130 raise ValueError("computeFluxes: problem with mesh, unknown cel number")
132 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
133 implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
135 if( Fk.getGroupName() != "Periodic" and Fk.getGroupName() != "Neumann"):#Wall boundary condition unless Periodic/Neumann specified explicitly
136 v=cdmath.Vector(dim+1)
137 for i in range(dim) :
139 idMoinsJacCL=v.tensProduct(v)*2
141 implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL)
143 elif( Fk.getGroupName() == "Periodic"):#Periodic boundary condition
144 indexFP=my_mesh.getIndexFacePeriodic(indexFace)
145 Fp = my_mesh.getFace(indexFP)
146 cellAutre = Fp.getCellsId()[0]
148 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
149 implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
150 elif(Fk.getGroupName() != "Neumann"):#Nothing to do for Neumann boundary condition
151 print Fk.getGroupName()
152 raise ValueError("computeFluxes: Unknown boundary condition name");
156 def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution):
157 dim=my_mesh.getMeshDimension()
158 nbCells = my_mesh.getNumberOfCells()
159 meshName=my_mesh.getName()
166 nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
170 Un=cdmath.Vector(nbCells*(dim+1))
171 dUn=cdmath.Vector(nbCells*(dim+1))
173 # Initial conditions #
174 print("Construction of the initial condition …")
175 if(dim==1 or filename.find("square")>-1 or filename.find("Square")>-1 or filename.find("cube")>-1 or filename.find("Cube")>-1):
176 pressure_field, velocity_field = initial_conditions_shock(my_mesh,False)
177 elif(filename.find("disk")>-1 or filename.find("Disk")>-1):
178 pressure_field, velocity_field = initial_conditions_shock(my_mesh,True)
180 print "Mesh name : ", filename
181 raise ValueError("Mesh name should contain substring square, cube or disk")
183 for k in range(nbCells):
184 Un[k*(dim+1)+0] = pressure_field[k]
185 Un[k*(dim+1)+1] = rho0*velocity_field[k,0]
187 Un[k + 2*nbCells] = rho0*velocity_field[k,1] # value on the bottom face
189 Un[k + 3*nbCells] = rho0*velocity_field[k,2]
191 #sauvegarde de la donnée initiale
192 pressure_field.setTime(time,it);
193 pressure_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure");
194 velocity_field.setTime(time,it);
195 velocity_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity");
196 #Postprocessing : save 2D picture
197 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure_initial")
198 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity_initial")
200 dx_min=my_mesh.minRatioVolSurf()
202 dt = cfl * dx_min / c0
204 divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
206 # Add the identity matrix on the diagonal
207 for j in range(nbCells*(dim+1)):
208 divMat.addValue(j,j,1.)
209 LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","LU")
211 print("Starting computation of the linear wave system with an pseudo staggered scheme …")
214 while (it<ntmax and time <= tmax and not isStationary):
218 cvgceLS=LS.getStatus();
219 iterGMRES=LS.getNumberOfIter();
221 print "Linear system did not converge ", iterGMRES, " GMRES iterations"
222 raise ValueError("Pas de convergence du système linéaire");
225 maxVector=dUn.maxVector(dim+1)
226 isStationary= maxVector[0]/p0<precision and maxVector[1]/rho0<precision and maxVector[2]/rho0<precision;
228 isStationary=isStationary and maxVector[3]/rho0<precision
233 if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
234 print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
235 print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0
236 print "Linear system converged in ", iterGMRES, " GMRES iterations"
238 for k in range(nbCells):
239 pressure_field[k]=Un[k*(dim+1)+0]
240 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
242 velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
244 velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
246 pressure_field.setTime(time,it);
247 pressure_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure",False);
248 velocity_field.setTime(time,it);
249 velocity_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity",False);
251 print"-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt)
252 print "Variation temporelle relative : pressure ", maxVector[0]/p0 ,", velocity x", maxVector[1]/rho0 ,", velocity y", maxVector[2]/rho0
256 print "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint"
258 print "Régime stationnaire atteint au pas de temps ", it, ", t= ", time
259 print "------------------------------------------------------------------------------------"
261 pressure_field.setTime(time,0);
262 pressure_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure_Stat");
263 velocity_field.setTime(time,0);
264 velocity_field.writeVTK("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity_Stat");
266 #Postprocessing : save 2D picture
267 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DPStag"+meshName+"_pressure_Stat")
268 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DPStag"+meshName+"_velocity_Stat")
271 print "Temps maximum Tmax= ", tmax, " atteint"
274 def solve(my_mesh,meshName,resolution):
275 print "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension()
276 print "Numerical method : implicit pseudo staggered"
277 print "Initial data : spherical wave"
278 print "Wall boundary conditions"
279 print "Mesh name : ",meshName , my_mesh.getNumberOfCells(), " cells"
284 cfl = 1./my_mesh.getSpaceDimension()
287 WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution)
289 def solve_file( filename,meshName, resolution):
290 my_mesh = cdmath.Mesh(filename+".med")
292 return solve(my_mesh, filename+str(my_mesh.getNumberOfCells()),resolution)
295 if __name__ == """__main__""":
296 if len(sys.argv) >1 :
298 my_mesh = cdmath.Mesh(filename)
299 solve(my_mesh,filename,100)
301 raise ValueError("WaveSystemUpwind.py expects a mesh file name")