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