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 : Test de préservation d'un état stationnaire
11 # Utilisation du schéma upwind explicite ou implicite sur un maillage général
12 # Initialisation par une vortex stationnaire
13 # Conditions aux limites parois
14 # Création et sauvegarde du champ résultant et des figures
15 #================================================================================================================================
18 from math import sin, cos, pi, sqrt
24 p0=155e5#reference pressure in a pressurised nuclear vessel
25 c0=700.#reference sound speed for water at 155 bars, 600K
26 rho0=p0/c0*c0#reference density
29 def initial_conditions_disk_vortex(my_mesh):
30 print( "Disk vortex initial data")
31 dim = my_mesh.getMeshDimension()
32 nbCells = my_mesh.getNumberOfCells()
35 raise ValueError("Wave system on disk : mesh dimension should be 2")
37 pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
38 velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
40 for i in range(nbCells):
41 x = my_mesh.getCell(i).x()
42 y = my_mesh.getCell(i).y()
44 pressure_field[i] = p0
46 velocity_field[i,0] = -y
47 velocity_field[i,1] = x
48 velocity_field[i,2] = 0
50 return pressure_field, velocity_field
52 def initial_conditions_square_vortex(my_mesh):
53 print( "Initial data : Square vortex (Constant pressure, divergence free velocity)" )
54 dim = my_mesh.getMeshDimension()
55 nbCells = my_mesh.getNumberOfCells()
57 pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
58 velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
60 for i in range(nbCells):
61 x = my_mesh.getCell(i).x()
62 y = my_mesh.getCell(i).y()
63 z = my_mesh.getCell(i).z()
65 pressure_field[i] = p0
67 velocity_field[i,0] = 1
68 velocity_field[i,1] = 0
69 velocity_field[i,2] = 0
71 velocity_field[i,0] = sin(pi*x)*cos(pi*y)
72 velocity_field[i,1] = -sin(pi*y)*cos(pi*x)
73 velocity_field[i,2] = 0
75 velocity_field[i,0] = sin(pi*x)*cos(pi*y)*cos(pi*z)
76 velocity_field[i,1] = sin(pi*y)*cos(pi*x)*cos(pi*z)
77 velocity_field[i,2] = -2*sin(pi*z)*cos(pi*x)*cos(pi*y)
79 return pressure_field, velocity_field
81 def jacobianMatrices(normal,coeff):
83 A=cdmath.Matrix(dim+1,dim+1)
84 absA=cdmath.Matrix(dim+1,dim+1)
88 A[i+1,0]= normal[i]*coeff
89 A[0,i+1]=c0*c0*normal[i]*coeff
91 absA[i+1,j+1]=c0*normal[i]*normal[j]*coeff
95 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt):
96 nbCells = my_mesh.getNumberOfCells()
97 dim=my_mesh.getMeshDimension()
99 normal=cdmath.Vector(dim)
101 implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
103 idMoinsJacCL=cdmath.Matrix(nbComp)
105 for j in range(nbCells):#On parcourt les cellules
106 Cj = my_mesh.getCell(j)
107 nbFaces = Cj.getNumberOfFaces();
109 for k in range(nbFaces) :
110 indexFace = Cj.getFacesId()[k];
111 Fk = my_mesh.getFace(indexFace);
112 for i in range(dim) :
113 normal[i] = Cj.getNormalVector(k, i);#normale sortante
115 Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
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 cell 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, isImplicit):
154 dim=my_mesh.getMeshDimension()
155 nbCells = my_mesh.getNumberOfCells()
156 meshName=my_mesh.getName()
162 nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
164 # Initial conditions #
165 print("Construction of the initial condition …")
166 if(filename.find("square")>-1 or filename.find("Square")>-1 or filename.find("cube")>-1 or filename.find("Cube")>-1):
167 pressure_field, velocity_field = initial_conditions_square_vortex(my_mesh)
168 elif(filename.find("disk")>-1 or filename.find("Disk")>-1):
169 pressure_field, velocity_field = initial_conditions_disk_vortex(my_mesh)
171 print( "Mesh name : ", filename)
172 raise ValueError("Mesh name should contain substring square, cube or disk")
175 Un =cdmath.Vector(nbCells*(dim+1))
176 dUn=cdmath.Vector(nbCells*(dim+1))
178 for k in range(nbCells):
179 Un[k*(dim+1)+0] = pressure_field[k]
180 Un[k*(dim+1)+1] =rho0*velocity_field[k,0]
181 Un[k*(dim+1)+2] =rho0*velocity_field[k,1]
183 Un[k*(dim+1)+3] =rho0*velocity_field[k,2]
185 #sauvegarde de la donnée initiale
186 pressure_field.setTime(time,it);
187 pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure");
188 velocity_field.setTime(time,it);
189 velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity");
191 dx_min=my_mesh.minRatioVolSurf()
193 dt = cfl * dx_min / c0
195 divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
197 #Adding the identity matrix on the diagonal
198 divMat.diagonalShift(1)#only after filling all coefficients
201 LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","ILU")
203 LS.setComputeConditionNumber()
205 print("Starting computation of the linear wave system with an UPWIND scheme …")
208 while (it<ntmax and time <= tmax and not isStationary):
213 if(not LS.getStatus()):
214 print( "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations")
215 raise ValueError("Pas de convergence du système linéaire");
226 if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
227 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
229 for k in range(nbCells):
230 pressure_field[k] =Un[k*(dim+1)+0]
231 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
233 velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
235 velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
237 pressure_field.setTime(time,it);
238 pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure",False);
239 velocity_field.setTime(time,it);
240 velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity",False);
242 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
246 print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
248 print( "Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
250 pressure_field.setTime(time,0);
251 pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat");
252 velocity_field.setTime(time,0);
253 velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat");
255 #Postprocessing : Extraction of the diagonal data
256 diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,1,0],[1,0,0], resolution)
257 diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,1,0],[1,0,0], resolution)
258 #Postprocessing : save 2D picture
259 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat")
260 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat")
263 print( "Temps maximum Tmax= ", tmax, " atteint")
266 def solve(my_mesh,filename,resolution, isImplicit):
267 print( "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension())
268 print( "Numerical method : upwind")
269 print( "Initial data : stationary solution (constant pressure, divergence free velocity)")
270 print( "Wall boundary conditions")
271 print( "Mesh name : ",filename , my_mesh.getNumberOfCells(), " cells")
276 cfl = 1./my_mesh.getSpaceDimension()
279 WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution, isImplicit)
281 def solve_file( filename,resolution):
282 my_mesh = cdmath.Mesh(filename+".med")
283 solve(my_mesh, filename,resolution, isImplicit)
285 if __name__ == """__main__""":
286 if len(sys.argv) >2 :
288 isImplicit=bool(int(sys.argv[2]))
289 my_mesh = cdmath.Mesh(filename)
290 solve(my_mesh,filename,100, isImplicit)
292 raise ValueError("WaveSystemUpwind.py expects a mesh file name and a boolean (isImplicit)")