4 from math import sin, cos, pi, sqrt
12 def initial_conditions_transport_equation(my_mesh):
13 print( "Initial_data","Shock")
14 dim = my_mesh.getMeshDimension()
15 nbCells = my_mesh.getNumberOfCells()
17 initial_field = cdmath.Field("unknown", cdmath.CELLS, my_mesh, 1)
24 for i in range(nbCells):
25 x = my_mesh.getCell(i).x()
27 valX = (x - xcentre) * (x - xcentre)
32 y = my_mesh.getCell(i).y()
33 valY = (y - ycentre) * (y - ycentre)
34 val = sqrt(valX + valY)
36 y = my_mesh.getCell(i).y()
37 z = my_mesh.getCell(i).z()
38 valY = (y - ycentre) * (y - ycentre)
39 valZ = (z - zcentre) * (z - zcentre)
40 val = sqrt(valX + valY + valZ)
52 def upwinding_coeff(normal, coeff, velocity):
55 if(velocity*normal<0.):
56 return velocity*normal*coeff
61 def computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,test_bc,velocity):
62 nbCells = my_mesh.getNumberOfCells()
63 dim=my_mesh.getMeshDimension()
65 normal=cdmath.Vector(dim)
67 implMat=cdmath.SparseMatrixPetsc(nbCells*nbComp,nbCells*nbComp,(nbVoisinsMax+1)*nbComp)
69 for j in range(nbCells):#On parcourt les cellules
70 Cj = my_mesh.getCell(j)
71 nbFaces = Cj.getNumberOfFaces();
73 for k in range(nbFaces) :
74 indexFace = Cj.getFacesId()[k];
75 Fk = my_mesh.getFace(indexFace);
77 normal[i] = Cj.getNormalVector(k, i);#normale sortante
79 Am=upwinding_coeff( normal,dt*Fk.getMeasure()/Cj.getMeasure(),velocity);
82 if ( not Fk.isBorder()) :
83 # hypothese: La cellule d'index indexC1 est la cellule courante index j */
84 if (Fk.getCellsId()[0] == j) :
86 cellAutre = Fk.getCellsId()[1];
87 elif(Fk.getCellsId()[1] == j) :
88 # hypothese non verifiée
89 cellAutre = Fk.getCellsId()[0];
91 raise ValueError("computeFluxes: problem with mesh, unknown cel number")
93 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
94 implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
96 if( test_bc=="Periodic" and Fk.getGroupName() != "Neumann"):#Periodic boundary condition unless Wall/Neumann specified explicitly
97 indexFP = my_mesh.getIndexFacePeriodic(indexFace, my_mesh.getName()== "squareWithBrickWall", my_mesh.getName()== "squareWithHexagons")
98 Fp = my_mesh.getFace(indexFP)
99 cellAutre = Fp.getCellsId()[0]
101 implMat.addValue(j*nbComp,cellAutre*nbComp,Am)
102 implMat.addValue(j*nbComp, j*nbComp,Am*(-1.))
103 elif(test_bc!="Neumann" and Fk.getGroupName() != "Neumann"):#Nothing to do for Neumann boundary condition
104 print( Fk.getGroupName() )
105 raise ValueError("computeFluxes: Unknown boundary condition name");
109 def TransportEquationVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution,test_bc,velocity,isImplicit):
110 dim=my_mesh.getMeshDimension()
111 nbCells = my_mesh.getNumberOfCells()
118 nbVoisinsMax=my_mesh.getMaxNbNeighbours(cdmath.CELLS)
121 Un =cdmath.Vector(nbCells)
122 dUn=cdmath.Vector(nbCells)
124 # Initial conditions #
125 print("Construction of the initial condition …")
126 stat_field = initial_conditions_transport_equation(my_mesh)
127 for k in range(nbCells):
128 Un[k] = stat_field[k]
129 unknown_field = initial_conditions_transport_equation(my_mesh)
130 S = cdmath.Vector(nbCells)#source term is zero
132 total_mass_initial=unknown_field.integral()#For conservation test later
134 #sauvegarde de la donnée initiale
135 unknown_field.writeVTK("TransportEquation"+str(dim)+"DUpwind"+meshName);
137 dx_min=my_mesh.minRatioVolSurf()
139 dt = cfl * dx_min / velocity.norm()
141 divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,test_bc,velocity)
143 #Adding the identity matrix on the diagonal
144 divMat.diagonalShift(1)#only after filling all coefficients
145 #for j in range(nbCells*(dim+1)):
146 # divMat.addValue(j,j,1)
149 LS=cdmath.LinearSolver(divMat,Un+S*dt,iterGMRESMax, precision, "GMRES","ILU")
150 LS.setComputeConditionNumber()
152 print("Starting computation of the linear transport equation with an Upwind scheme …")
155 while (it<ntmax and time <= tmax and not isStationary):
158 LS.setSndMember(Un+S*dt)
160 if(not LS.getStatus()):
161 print("Linear system did not converge ", LS.getNumberOfIter(), " GMRES iterations")
162 raise ValueError("Pas de convergence du système linéaire");
169 maxVector=dUn.maxVector(1)
170 isStationary= maxVector[0]<precision
176 if(it%output_freq==0 or it>=ntmax or isStationary or time >=tmax):
177 print( "-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
178 print( "Variation temporelle relative : ", maxVector[0] )
180 print( "Linear system converged in ", LS.getNumberOfIter(), " GMRES iterations" )
182 for k in range(nbCells):
183 unknown_field[k] =Un[k]
185 unknown_field.setTime(time,it);
186 unknown_field.writeVTK("TransportEquation"+str(dim)+"DUpwind"+meshName,False);
188 print( "-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt) )
189 print( "Variation temporelle relative : ", maxVector[0] )
192 print( "Nombre de pas de temps maximum ntmax= ", ntmax, " atteint" )
193 raise ValueError("Maximum number of time steps reached : Stationary state not found !!!!!!!")
195 print( "Régime stationnaire atteint au pas de temps ", it, ", t= ", time)
196 if(test_bc=="Periodic"):
197 print( "Mass loss: ", (total_mass_initial-unknown_field.integral()).norm(), " precision required= ", precision )
198 assert (total_mass_initial-unknown_field.integral()).norm()<precision
199 print "------------------------------------------------------------------------------------"
201 unknown_field.setTime(time,0);
202 unknown_field.writeVTK("TransportEquation"+str(dim)+"DUpwind"+meshName+"_Stat");
204 #Postprocessing : Extraction of the diagonal data
206 diag_data_u=VTK_routines.Extract_field_data_over_line_to_numpyArray(unknown_field,[0,1,0],[1,0,0], resolution)
208 diag_data_u=VTK_routines.Extract_field_data_over_line_to_numpyArray(unknown_field,[0,0,0],[1,1,1], resolution)
209 #Postprocessing : save 2D picture
210 PV_routines.Save_PV_data_to_picture_file("TransportEquation"+str(dim)+"DUpwind"+meshName+"_Stat"+'_0.vtu',"unknown",'CELLS',"TransportEquation"+str(dim)+"DUpwind"+meshName+"_Stat")
212 return nbCells, time, it, unknown_field.getNormEuclidean().max(), diag_data_u
214 print( "Temps maximum Tmax= ", tmax, " atteint" )
215 raise ValueError("Maximum time reached : Stationary state not found !!!!!!!")
218 def solve(my_mesh, meshName, resolution, meshType, cfl, test_bc):
219 print( "Resolution of the Transport Equation in dimension ", my_mesh.getMeshDimension() )
220 print( "Numerical method : ", "Upwind" )
221 print( "Initial data : ", "Spherical shock" )
222 print( "Mesh name : ",meshName , ", ", my_mesh.getNumberOfCells(), " cells" )
231 dim=my_mesh.getMeshDimension()
232 velocity=cdmath.Vector(dim)
233 for i in range(dim) :
236 nbCells, t_final, ndt_final, max_unknown, diag_data_u = TransportEquationVF(ntmax, tmax, cfl, my_mesh, output_freq, meshName, resolution, test_bc, velocity,isImplicit)
238 return nbCells, t_final, ndt_final, max_unknown, diag_data_u
240 def solve_file( filename,meshName, resolution,meshType, cfl, test_bc):
241 my_mesh = cdmath.Mesh(filename+".med")
243 return solve(my_mesh, meshName+str(my_mesh.getNumberOfCells()),resolution, meshType, cfl, test_bc)
246 if __name__ == """__main__""":
247 M1=cdmath.Mesh(0.,1.,20,0.,1.,20,0)
249 solve(M1,"SquareRegularTriangles",100,"Regular triangles",cfl,"Periodic")