DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
test_basic_eulerP.py
Go to the documentation of this file.
1from DNDSR import DNDS, Geom, CFV, EulerP
2import json
3from DNDSR.Geom.utils import *
4import numpy as np
5from DNDSR.DNDS.Debug_Py import MPIDebugHold
6
7
8def get_fv(mpi):
9 print(CFV.tUDof_1)
10 # CFV.VariationalReconstruction_2()
11 # meshFile = os.path.join(
12 # os.path.dirname(__file__), "..", "..", "data", "mesh", "NACA0012_H2.cgns"
13 # )
14 # meshFile = os.path.join(
15 # os.path.dirname(__file__), "..", "..", "data", "mesh", "Uniform_3x3.cgns"
16 # )
17 meshFile = os.path.join(
18 os.path.dirname(__file__), "..", "..", "data", "mesh", "Uniform32_Periodic.cgns"
19 )
20 mesh, reader, name2Id = create_mesh_from_CGNS(
21 meshFile,
22 mpi,
23 2,
24 periodic_geometry={
25 # "translation1": [3, 0, 0],
26 # "translation2": [0, 3, 0],
27 "translation1": [1, 0, 0],
28 "translation2": [0, 1, 0],
29 },
30 meshDirectBisect=4,
31 second_level_parts=16,
32 )
33 meshBnd, readerBnd = create_bnd_mesh(mesh)
34
35 fv = CFV.FiniteVolume(mpi, mesh)
36 settings = fv.GetSettings()
37 settings["intOrder"] = 3
38 settings["maxOrder"] = 3
39 fv.ParseSettings(settings)
40 if mpi.rank == 0:
41 print(fv.GetSettings())
42
43 # construction: volume
44 fv.SetCellAtrBasic()
45 fv.ConstructCellVolume()
46 fv.ConstructCellBary()
47 fv.ConstructCellCent()
48 fv.ConstructCellIntJacobiDet()
49 fv.ConstructCellIntPPhysics()
50 fv.ConstructCellAlignedHBox()
51 fv.ConstructCellMajorHBoxCoordInertia()
52 # construction: face
53 fv.SetFaceAtrBasic()
54 fv.ConstructFaceCent()
55 fv.ConstructFaceArea()
56 fv.ConstructFaceIntJacobiDet()
57 fv.ConstructFaceIntPPhysics()
58 fv.ConstructFaceUnitNorm()
59 fv.ConstructFaceMeanNorm()
60
61 fv.ConstructCellSmoothScale()
62
63 nB = fv.getArrayBytes() / 1024**2
64 nBMesh = mesh.getArrayBytes() / 1024**2
65 if mpi.rank == 0:
66 print(f"Bytes : {nB:.4g} MB")
67 print(f"Bytes mesh: {nBMesh:.4g} MB")
68
69 return mesh, reader, name2Id, meshBnd, readerBnd, fv
70
71
72def test_basic_eulerP(mpi: DNDS.MPIInfo, isCuda=False):
73
74 mesh, reader, name2Id, meshBnd, readerBnd, fv = get_fv(mpi)
75 print(name2Id.n2id_map)
76 print(name2Id.n2id_map is name2Id.n2id_map)
77 phys_dict = EulerP.Physics().to_dict()
78 phys_params = phys_dict["params"]
79 phys_params["TRef"] = 273.15
80 phys_params["cp"] = 1.0
81 phys_params["gamma"] = 1.4
82 phys_params["mu0"] = 1e-10
83 phys_dict["reference_values"] = [1.0, 1.0, 1.0, 1.0, 1.0]
84
85 phys = EulerP.Physics.from_dict(phys_dict)
86 print("Physics:")
87 print(phys.to_dict())
88
89 # bcInputs = [EulerP.BCInput()]
90 # bcInputs[0].name = "wall"
91 # bcInputs[0].type = EulerP.BCType.Wall
92 # bcInputs[0].value = [1]
93
94 bcInputs = []
95
96 bcHandler = EulerP.BCHandler(bcInputs, name2Id)
97 print(bcHandler.id2bc(1).type)
98
99 eval = EulerP.Evaluator(fv, bcHandler, phys)
100
101 u = CFV.tUDof_5()
102 fv.BuildUDof_5(u, 5)
103 gU = CFV.tUGrad_3x5()
104 fv.BuildUGrad_3x5(gU, 5)
105 s = CFV.tUDof_1()
106 fv.BuildUDof_1(s, 1)
107 gS = CFV.tUGrad_3x1()
108 fv.BuildUGrad_3x1(gS, 1)
109
110 uf = CFV.tUDof_5()
111 fv.BuildUDof_5(uf, 5, varloc=Geom.MeshLoc.Face)
112 gUf = CFV.tUGrad_3x5()
113 fv.BuildUGrad_3x5(gUf, 5, varloc=Geom.MeshLoc.Face)
114 sf = CFV.tUDof_1()
115 fv.BuildUDof_1(sf, 1, varloc=Geom.MeshLoc.Face)
116 gSf = CFV.tUGrad_3x1()
117 fv.BuildUGrad_3x1(gSf, 1, varloc=Geom.MeshLoc.Face)
118 if isCuda:
119 mesh.to_device("CUDA")
120 fv.to_device("CUDA")
121 eval.to_device("CUDA")
122 eval_device = eval.device()
123 print(eval_device)
124
125 u.to_device("CUDA")
126 gU.to_device("CUDA")
127 s.to_device("CUDA")
128 gS.to_device("CUDA")
129
130 uf.to_device("CUDA")
131 gUf.to_device("CUDA")
132 sf.to_device("CUDA")
133 gSf.to_device("CUDA")
134
135 data = {}
136 data["u"] = u
137 data["uGrad"] = gU
138
139 u.setConstant(np.array([1, 0, 0, 0, 2.5]).reshape(-1, 1))
140 for iCell in range(mesh.NumCell()):
141 x = fv.GetCellBary(iCell)
142 if np.linalg.norm(np.array([0.5, 0.5, 0]) - x, 2) < 0.25:
143 u[iCell] = np.array([1, 4, 0, 0, 10.5], dtype=np.float64).reshape(-1, 1)
144
145 uu = 1
146 vv = 1
147 r = 1 + 0.1 * np.cos(x[0] * 2 * np.pi) * np.cos(x[1] * 2 * np.pi)
148 u[iCell] = np.array(
149 [r, r * uu, r * vv, 0, 2.5 + 0.5 * r * (uu**2 + vv**2)], dtype=np.float64
150 ).reshape(-1, 1)
151 if isCuda:
152 u.to_device("CUDA")
153 print(u.father)
154
155 data["uGrad"].setConstant(100.0)
156 print(data["uGrad"].norm2())
157 RecGradient_Arg = eval.RecGradient_Arg()
158 RecGradient_Arg.u = data["u"]
159 RecGradient_Arg.uGrad = data["uGrad"]
160 RecGradient_Arg.uScalar = []
161 RecGradient_Arg.uScalarGrad = []
162 eval.RecGradient(RecGradient_Arg)
163
164 data["p"] = s.clone()
165 data["T"] = s.clone()
166 data["a"] = s.clone()
167 data["mu"] = s.clone()
168 data["gamma"] = s.clone()
169
170 data["uPrim"] = u.clone()
171 data["uGradPrim"] = gU.clone()
172
173 Cons2PrimMu_Arg = eval.Cons2PrimMu_Arg()
174 Cons2PrimMu_Arg.p = data["p"]
175 Cons2PrimMu_Arg.T = data["T"]
176 Cons2PrimMu_Arg.a = data["a"]
177 Cons2PrimMu_Arg.mu = data["mu"]
178 Cons2PrimMu_Arg.gamma = data["gamma"]
179
180 Cons2PrimMu_Arg.u = data["u"]
181 Cons2PrimMu_Arg.uGrad = data["uGrad"]
182 Cons2PrimMu_Arg.uPrim = data["uPrim"]
183 Cons2PrimMu_Arg.uGradPrim = data["uGradPrim"]
184
185 eval.Cons2PrimMu(Cons2PrimMu_Arg)
186
187 data["deltaLamCell"] = s.clone()
188 data["dt"] = s.clone()
189 data["uFL"] = uf
190 data["uGradFF"] = gUf
191 data["deltaLamFace"] = sf.clone()
192 data["faceLamEst"] = gSf.clone()
193 data["faceLamVisEst"] = sf.clone()
194
195 EstEigenDt_Arg = eval.EstEigenDt_Arg()
196 EstEigenDt_Arg.u = data["u"]
197 EstEigenDt_Arg.muCell = data["mu"]
198 EstEigenDt_Arg.aCell = data["a"]
199 EstEigenDt_Arg.deltaLamCell = data["deltaLamCell"]
200 EstEigenDt_Arg.deltaLamFace = data["deltaLamFace"]
201 EstEigenDt_Arg.faceLamEst = data["faceLamEst"]
202 EstEigenDt_Arg.faceLamVisEst = data["faceLamVisEst"]
203 EstEigenDt_Arg.dt = data["dt"]
204 eval.EstEigenDt(EstEigenDt_Arg)
205
206 data["uGradFF"] = gUf.clone()
207 data["uFL"] = uf.clone()
208 data["uFR"] = uf.clone()
209
210 RecFace2nd_Arg = eval.RecFace2nd_Arg()
211 RecFace2nd_Arg.u = data["u"]
212 RecFace2nd_Arg.uGrad = data["uGrad"]
213 RecFace2nd_Arg.uFL = data["uFL"]
214 RecFace2nd_Arg.uFR = data["uFR"]
215 RecFace2nd_Arg.uGradFF = data["uGradFF"]
216 # data["uGrad"].setConstant(0.0)
217 eval.RecFace2nd(RecFace2nd_Arg)
218
219 data["uPrimFL"] = uf.clone()
220 data["pFL"] = sf.clone()
221 data["TFL"] = sf.clone()
222 data["aFL"] = sf.clone()
223 data["gammaFL"] = sf.clone()
224 data["uPrimFR"] = uf.clone()
225 data["pFR"] = sf.clone()
226 data["TFR"] = sf.clone()
227 data["aFR"] = sf.clone()
228 data["gammaFR"] = sf.clone()
229
230 Cons2Prim_Arg_FL = eval.Cons2Prim_Arg()
231 Cons2Prim_Arg_FL.u = data["uFL"]
232 Cons2Prim_Arg_FL.uPrim = data["uPrimFL"]
233 Cons2Prim_Arg_FL.p = data["pFL"]
234 Cons2Prim_Arg_FL.T = data["TFL"]
235 Cons2Prim_Arg_FL.a = data["aFL"]
236 Cons2Prim_Arg_FL.gamma = data["gammaFL"]
237
238 Cons2Prim_Arg_FR = eval.Cons2Prim_Arg()
239 Cons2Prim_Arg_FR.u = data["uFR"]
240 Cons2Prim_Arg_FR.uPrim = data["uPrimFR"]
241 Cons2Prim_Arg_FR.p = data["pFR"]
242 Cons2Prim_Arg_FR.T = data["TFR"]
243 Cons2Prim_Arg_FR.a = data["aFR"]
244 Cons2Prim_Arg_FR.gamma = data["gammaFR"]
245 eval.Cons2Prim(Cons2Prim_Arg_FL)
246 eval.Cons2Prim(Cons2Prim_Arg_FR)
247
248 pFMin = min(data["pFR"].min(), data["pFL"].min())
249 print(f"pFMin, {pFMin}")
250
251 data["fluxFF"] = uf.clone()
252 data["rhs"] = u.clone()
253
254 Flux2nd_Arg = eval.Flux2nd_Arg()
255 Flux2nd_Arg.u = data["u"]
256 Flux2nd_Arg.uGrad = data["uGrad"]
257 Flux2nd_Arg.uPrim = data["uPrim"]
258 Flux2nd_Arg.uGradPrim = data["uGradPrim"]
259
260 Flux2nd_Arg.p = data["p"]
261 Flux2nd_Arg.T = data["T"]
262 Flux2nd_Arg.a = data["a"]
263 Flux2nd_Arg.mu = data["mu"]
264 Flux2nd_Arg.gamma = data["gamma"]
265
266 Flux2nd_Arg.uFL = data["uFL"]
267 Flux2nd_Arg.uFR = data["uFR"]
268 Flux2nd_Arg.pFL = data["pFL"]
269 Flux2nd_Arg.pFR = data["pFR"]
270
271 Flux2nd_Arg.uGradFF = data["uGradFF"]
272 Flux2nd_Arg.deltaLamFaceFF = data["deltaLamFace"]
273
274 Flux2nd_Arg.fluxFF = data["fluxFF"]
275 Flux2nd_Arg.rhs = data["rhs"]
276
277 data["rhs"].setConstant(0.0)
278 eval.Flux2nd(Flux2nd_Arg)
279
280 check_norm = data["p"].norm2()
281 print(f"p norm, {check_norm}")
282 check_norm = data["p"].min()
283 print(f"p min, {check_norm}")
284
285 data["p"].to_host()
286 check_norm = data["p"].min()
287 print(f"p min, {check_norm}")
288
289 for n, a in data.items():
290 a.to_host()
291 cell_out = s.clone()
292 cell_out1 = s.clone()
293 cell_out2 = s.clone()
294
295 for iCell in range(mesh.NumCell()):
296 np.array(cell_out[iCell], copy=False)[:] = np.array(
297 data["rhs"][iCell], copy=False
298 )[0]
299 np.array(cell_out1[iCell], copy=False)[:] = np.array(
300 data["u"][iCell], copy=False
301 )[1]
302 np.array(cell_out1[iCell], copy=False)[:] = iCell
303 np.array(cell_out2[iCell], copy=False)[:] = np.array(
304 data["uGrad"][iCell], copy=False
305 )[0, 1]
306
307 eval.PrintDataVTKHDF(
308 "test_0",
309 "test",
310 [cell_out, cell_out1, cell_out2, data["p"]],
311 ["cell_out", "cell_out1", "cell_out2", "p"],
312 )
313
314 # for iCell in range(mesh.NumCell()):
315 # print(iCell)
316 # print(f"{np.asarray(data['u'][iCell]).T}")
317 # for iFace in mesh.cell2face[iCell]:
318 # print(f"{iFace} -- {np.asarray(data['uFL'][iFace]).T}")
319
320
321if __name__ == "__main__":
322 mpi = DNDS.MPIInfo()
323 mpi.setWorld()
324
325 ## debug py
326
327 # MPIDebugHold()
328
329 # DNDS.Debug.MPIDebugHold(mpi)
330 test_basic_eulerP(mpi, isCuda=False)
test_basic_eulerP
Definition test_basic_eulerP.py:1
test_basic_eulerP.get_fv
get_fv(mpi)
Definition test_basic_eulerP.py:8
DNDS::MPIInfo
Lightweight bundle of an MPI communicator and the calling rank's coordinates.
Definition MPI.hpp:215
setConstant
dof1 setConstant(0.0)