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 upwind explicite ou 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 #================================================================================================================================
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_shock(my_mesh, isCircle):
30 print( "Initial data : Spherical wave" )
31 dim = my_mesh.getMeshDimension()
32 nbCells = my_mesh.getNumberOfCells()
44 pressure_field = cdmath.Field("Pressure", cdmath.CELLS, my_mesh, 1)
45 velocity_field = cdmath.Field("Velocity", cdmath.CELLS, my_mesh, 3)
47 for i in range(nbCells):
48 velocity_field[i,0] = 0
49 velocity_field[i,1] = 0
50 velocity_field[i,2] = 0
52 x = my_mesh.getCell(i).x()
53 valX = (x - xcentre) * (x - xcentre)
58 y = my_mesh.getCell(i).y()
59 valY = (y - ycentre) * (y - ycentre)
60 val = sqrt(valX + valY)
62 y = my_mesh.getCell(i).y()
63 z = my_mesh.getCell(i).z()
64 valY = (y - ycentre) * (y - ycentre)
65 valZ = (z - zcentre) * (z - zcentre)
66 val = sqrt(valX + valY + valZ)
69 pressure_field[i] = p0
72 pressure_field[i] = p0/2
76 return pressure_field, velocity_field
78 def jacobianMatrices(normal,coeff):
80 A=cdmath.Matrix(dim+1,dim+1)
81 absA=cdmath.Matrix(dim+1,dim+1)
85 A[i+1,0]= normal[i]*coeff
86 A[0,i+1]=c0*c0*normal[i]*coeff
88 absA[i+1,j+1]=c0*normal[i]*normal[j]*coeff
92 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt):
93 nbCells = my_mesh.getNumberOfCells()
94 dim=my_mesh.getMeshDimension()
96 normal=cdmath.Vector(dim)
98 implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
100 idMoinsJacCL=cdmath.Matrix(nbComp)
102 for j in range(nbCells):#On parcourt les cellules
103 Cj = my_mesh.getCell(j)
104 nbFaces = Cj.getNumberOfFaces();
106 for k in range(nbFaces) :
107 indexFace = Cj.getFacesId()[k];
108 Fk = my_mesh.getFace(indexFace);
109 for i in range(dim) :
110 normal[i] = Cj.getNormalVector(k, i);#normale sortante
112 Am=jacobianMatrices( normal,dt*Fk.getMeasure()/Cj.getMeasure());
115 if ( not Fk.isBorder()) :
116 # hypothese: La cellule d'index indexC1 est la cellule courante index j */
117 if (Fk.getCellsId()[0] == j) :
119 cellAutre = Fk.getCellsId()[1];
120 elif(Fk.getCellsId()[1] == j) :
121 # hypothese non verifiée
122 cellAutre = Fk.getCellsId()[0];
124 raise ValueError("computeFluxes: problem with mesh, unknown cell number")
126 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
127 implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
129 if( Fk.getGroupName() != "Periodic" and Fk.getGroupName() != "Neumann"):#Wall boundary condition unless Periodic/Neumann specified explicitly
130 v=cdmath.Vector(dim+1)
131 for i in range(dim) :
133 idMoinsJacCL=v.tensProduct(v)*2
135 implMat.addValue(j*nbComp,j*nbComp,Am*(-1.)*idMoinsJacCL)
137 elif( Fk.getGroupName() == "Periodic"):#Periodic boundary condition
138 indexFP=my_mesh.getIndexFacePeriodic(indexFace)
139 Fp = my_mesh.getFace(indexFP)
140 cellAutre = Fp.getCellsId()[0]
142 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
143 implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
144 elif(Fk.getGroupName() != "Neumann"):#Nothing to do for Neumann boundary condition
145 print( Fk.getGroupName() )
146 raise ValueError("computeFluxes: Unknown boundary condition name");
150 def WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution, isImplicit):
151 dim=my_mesh.getMeshDimension()
152 nbCells = my_mesh.getNumberOfCells()
153 meshName=my_mesh.getName()
159 nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
161 # Initial conditions #
162 print("Construction of the initial condition …")
163 if(dim==1 or filename.find("square")>-1 or filename.find("Square")>-1 or filename.find("cube")>-1 or filename.find("Cube")>-1):
164 pressure_field, velocity_field = initial_conditions_shock(my_mesh,False)
165 elif(filename.find("disk")>-1 or filename.find("Disk")>-1):
166 pressure_field, velocity_field = initial_conditions_shock(my_mesh,True)
168 print( "Mesh name : ", filename )
169 raise ValueError("Mesh name should contain substring square, cube or disk")
172 Un =cdmath.Vector(nbCells*(dim+1))
173 dUn=cdmath.Vector(nbCells*(dim+1))
175 for k in range(nbCells):
176 Un[k*(dim+1)+0] = pressure_field[k]
177 Un[k*(dim+1)+1] =rho0*velocity_field[k,0]
179 Un[k + 2*nbCells] = rho0*velocity_field[k,1] # value on the bottom face
181 Un[k + 3*nbCells] = rho0*velocity_field[k,2]
183 #sauvegarde de la donnée initiale
184 pressure_field.setTime(time,it);
185 pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_pressure");
186 velocity_field.setTime(time,it);
187 velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_velocity");
189 dx_min=my_mesh.minRatioVolSurf()
191 dt = cfl * dx_min / c0
193 divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt)
195 #Adding the identity matrix on the diagonal
196 divMat.diagonalShift(1)#only after filling all coefficients
199 LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, "GMRES","ILU")
201 LS.setComputeConditionNumber()
203 print("Starting computation of the linear wave system with an explicit UPWIND scheme …")
206 while (it<ntmax and time <= tmax and not isStationary):
211 if(not LS.getStatus()):
212 print( "Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations" )
213 raise ValueError("Pas de convergence du système linéaire");
224 if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
225 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
227 for k in range(nbCells):
228 pressure_field[k] =Un[k*(dim+1)+0]
229 velocity_field[k,0]=Un[k*(dim+1)+1]/rho0
231 velocity_field[k,1]=Un[k*(dim+1)+2]/rho0
233 velocity_field[k,2]=Un[k*(dim+1)+3]/rho0
235 pressure_field.setTime(time,it);
236 pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_pressure",False);
237 velocity_field.setTime(time,it);
238 velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_velocity",False);
240 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
244 print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint")
246 print( "Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
248 pressure_field.setTime(time,0);
249 pressure_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_pressure_Stat");
250 velocity_field.setTime(time,0);
251 velocity_field.writeVTK("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_velocity_Stat");
253 #Postprocessing : Extraction of the diagonal data
254 diag_data_press=VTK_routines.Extract_field_data_over_line_to_numpyArray(pressure_field,[0,1,0],[1,0,0], resolution)
255 diag_data_vel =VTK_routines.Extract_field_data_over_line_to_numpyArray(velocity_field,[0,1,0],[1,0,0], resolution)
256 #Postprocessing : save 2D picture
257 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_pressure_Stat"+'_0.vtu',"Pressure",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_pressure_Stat")
258 PV_routines.Save_PV_data_to_picture_file("WaveSystem"+str(dim)+"DUpwind"+"_isImplicit"+str(isImplicit)+meshName+"_velocity_Stat"+'_0.vtu',"Velocity",'CELLS',"WaveSystem"+str(dim)+"DUpwind"+meshName+"_velocity_Stat")
261 print( "Temps maximum Tmax= ", tmax, " atteint")
264 def solve(my_mesh,filename,resolution, isImplicit):
265 print( "Resolution of the Wave system in dimension ", my_mesh.getSpaceDimension() )
266 print( "Numerical method : upwind" )
267 print( "Initial data : spherical wave" )
268 print( "Wall boundary conditions" )
269 print( "Mesh name : ",filename , my_mesh.getNumberOfCells(), " cells" )
274 cfl = 1./my_mesh.getSpaceDimension()
277 WaveSystemVF(ntmax, tmax, cfl, my_mesh, output_freq, filename,resolution, isImplicit)
279 def solve_file( filename,resolution, isImplicit):
280 my_mesh = cdmath.Mesh(filename+".med")
281 solve(my_mesh, filename,resolution, isImplicit)
283 if __name__ == """__main__""":
284 if len(sys.argv) >2 :
286 isImplicit=bool(int(sys.argv[2]))
287 my_mesh = cdmath.Mesh(filename)
288 solve(my_mesh,filename,100, isImplicit)
290 raise ValueError("WaveSystemUpwind.py expects a mesh file name and a boolean (isImplicit)")