9 def initial_conditions(my_mesh):
13 y_field = cdmath.Field("Y field", cdmath.CELLS, my_mesh, 1)
14 nbCells = my_mesh.getNumberOfCells()
15 for j in range(nbCells):
16 x = my_mesh.getCell(j).x()
17 y = my_mesh.getCell(j).y()
18 valX = (x - xcentre) * (x - xcentre)
19 valY = (y - ycentre) * (y - ycentre)
20 val = math.sqrt(valX + valY)
31 def sigma_flux(VitesseX, VitesseY, cfl, y_field, indexFacesPerio):
32 # Calculation of fluxes #
33 SumFlux = cdmath.Field("Fluxes", cdmath.CELLS, y_field.getMesh(), 1)
34 my_mesh = y_field.getMesh()
35 nbCells = my_mesh.getNumberOfCells()
36 normU = math.sqrt(VitesseX * VitesseX + VitesseY * VitesseY)
37 for j in range(nbCells):
38 Cj = my_mesh.getCell(j)
39 nbFace = Cj.getNumberOfFaces()
42 for k in range(nbFace):
43 indexFace = Cj.getFacesId()[k]
44 Fk = my_mesh.getFace(indexFace)
45 NormalX = Cj.getNormalVector(k, 0)
46 NormalY = Cj.getNormalVector(k, 1)
47 LengthFk = Fk.getMeasure()
48 UN = VitesseX * NormalX + VitesseY * NormalY
49 minlengthFk = min(minlengthFk, LengthFk / abs(UN))
50 minlengthFk = min(minlengthFk, LengthFk / abs(VitesseX))
51 minlengthFk = min(minlengthFk, LengthFk / abs(VitesseY))
55 if (not Fk.isBorder()):
56 indexC1 = Fk.getCellsId()[0]
57 indexC2 = Fk.getCellsId()[1]
58 # hypothesis: the cell of index indexC1 is the current cell of index j #
60 # hypothese is verified #
61 cellCourante = indexC1
65 # hypothesis is not verified #
66 cellCourante = indexC2
69 # define left and right cell with the product of velocity * outward normal vector
70 # if u*n>0: nothing to do, else invert left and right
72 conc = y_field[cellCourante]
75 conc = y_field[cellAutre]
79 # homogeneous Neumann boundary conditions #
80 if (Fk.getGroupName() == "GAUCHE" or Fk.getGroupName() == "DROITE"):
82 conc = y_field[cellCourante]
88 # periodical boundary conditions #
89 if (Fk.getGroupName() == "BAS" or Fk.getGroupName() == "HAUT"):
90 indexFP = indexFacesPerio[indexFace]
91 # another way to get the index of the periodical face #
92 # int indexFP=my_mesh.getIndexFacePeriodic(indexFace);
93 Fp = my_mesh.getFace(indexFP)
94 indexCp = Fp.getCellsId()[0]
96 conc = y_field[cellCourante]
99 conc = y_field[indexCp]
103 SumF = SumF + UN * LengthFk * conc
105 dt = cfl * minlengthFk / normU
106 SumFlux[j] = dt * SumF / Cj.getMeasure()
111 def EquationTransport2D(tmax, VitesseX, VitesseY, cfl, freqSortie, my_mesh, output_filename):
113 # Initial conditions #
114 print("Construction of the initial condition …")
115 y_field = initial_conditions(my_mesh)
117 # MED output of the initial condition at t=0 and iter = 0
122 print("Saving the solution at T=" + str(time) + "…")
123 y_field.setTime(time, it)
124 y_field.writeMED(output_filename)
125 y_field.writeVTK(output_filename)
126 y_field.writeCSV(output_filename)
129 print("Resolution of the transport equation with an UPWIND scheme…")
131 indexFacesPerio = my_mesh.getIndexFacePeriodic()
133 while (it < ntmax and time <= tmax):
134 dt, SumFlux = sigma_flux(VitesseX, VitesseY, cfl, y_field, indexFacesPerio)
135 print("-- Iter: " + str(it) + ", Time: " + str(time) + ", dt: " + str(dt))
137 # Advancing one time step #
141 # Output every freq times
142 if (it % freqSortie == 0):
143 y_field.setTime(time, it)
144 y_field.writeMED(output_filename, False)
145 y_field.writeVTK(output_filename, False)
146 y_field.writeCSV(output_filename)
153 print("Resolution of the 2D transport equation:")
154 print("- DOMAIN: SQUARE")
155 print("- MESH: TRIANGULAR, GENERATED WITH SALOME")
156 print("- PERIODICAL BC UP AND DOWN")
157 print("- HOMOGENEOUS NEUMANN BC LEFT AND RIGHT")
166 print("Loading triangular mesh …")
167 my_mesh = cdmath.Mesh("../../tests/ressources/meshSquare.med")
168 output_filename = "Exercice2PyTest"
169 EquationTransport2D(tmax, VitesseX, VitesseY, cfl, freqSortie, my_mesh, output_filename)
170 print("CDMATH calculation done.")
174 if __name__ == """__main__""":