Note

Go to the end to download the full example code.

Testing Sod tube with Godunov#

CI test for Sod tube with Godunov

 8 import os
 9
10 import matplotlib.pyplot as plt
11 import numpy as np
12
13 import shamrock

Setup params

18 multx = 4
19 multy = 1
20 multz = 1
21
22 sz = 1 << 1
23 base = 32
24 gamma = 1.4

Init model

30 ctx = shamrock.Context()
31 ctx.pdata_layout_new()
32
33 model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")

Init config

38 cfg = model.gen_default_config()
39 scale_fact = 2 / (sz * base * multx)
40 cfg.set_scale_factor(scale_fact)
41
42 cfg.set_eos_gamma(gamma)
43 # cfg.set_riemann_solver_rusanov()
44 cfg.set_riemann_solver_hll()
45
46 # cfg.set_slope_lim_none()
47 # cfg.set_slope_lim_vanleer_f()
48 # cfg.set_slope_lim_vanleer_std()
49 # cfg.set_slope_lim_vanleer_sym()
50 cfg.set_slope_lim_minmod()
51 cfg.set_face_time_interpolation(True)
52 model.set_solver_config(cfg)

Init scheduler and grid

 58 model.init_scheduler(int(1e7), 1)
 59 model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz))
 60
 61
 62 # without face time interpolation
 63 # 0.07979993131348424 (0.17970690984930585, 0.0, 0.0) 0.12628776652228088
 64
 65 # with face time interpolation
 66 # 0.07894793711859852 (0.17754462339166546, 0.0, 0.0) 0.12498304725061045
 67
 68
 69 kx, ky, kz = 2 * np.pi, 0, 0
 70 delta_rho = 1e-2
 71
 72
 73 def rho_map(rmin, rmax):
 74
 75     x, y, z = rmin
 76     if x < 1:
 77         return 1
 78     else:
 79         return 0.125
 80
 81
 82 etot_L = 1.0 / (gamma - 1)
 83 etot_R = 0.1 / (gamma - 1)
 84
 85
 86 def rhoetot_map(rmin, rmax):
 87
 88     rho = rho_map(rmin, rmax)
 89
 90     x, y, z = rmin
 91     if x < 1:
 92         return etot_L
 93     else:
 94         return etot_R
 95
 96
 97 def rhovel_map(rmin, rmax):
 98     rho = rho_map(rmin, rmax)
 99
100     return (0, 0, 0)
101
102
103 model.set_field_value_lambda_f64("rho", rho_map)
104 model.set_field_value_lambda_f64("rhoetot", rhoetot_map)
105 model.set_field_value_lambda_f64_3("rhovel", rhovel_map)
106
107 t_target = 0.245
108
109 model.evolve_until(t_target)
110
111 # model.evolve_once()
112 xref = 1.0
113 xrange = 0.5
114 sod = shamrock.phys.SodTube(gamma=gamma, rho_1=1, P_1=1, rho_5=0.125, P_5=0.1)
115 sodanalysis = model.make_analysis_sodtube(sod, (1, 0, 0), t_target, xref, -xrange, xrange)
116
117
118 #################
119 ### Plot
120 #################
121 # do plot or not
122 if False:
123
124     def convert_to_cell_coords(dic):
125
126         cmin = dic["cell_min"]
127         cmax = dic["cell_max"]
128
129         xmin = []
130         ymin = []
131         zmin = []
132         xmax = []
133         ymax = []
134         zmax = []
135
136         for i in range(len(cmin)):
137
138             m, M = cmin[i], cmax[i]
139
140             mx, my, mz = m
141             Mx, My, Mz = M
142
143             for j in range(8):
144                 a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j)
145
146                 x, y, z = a
147                 xmin.append(x)
148                 ymin.append(y)
149                 zmin.append(z)
150
151                 x, y, z = b
152                 xmax.append(x)
153                 ymax.append(y)
154                 zmax.append(z)
155
156         dic["xmin"] = np.array(xmin)
157         dic["ymin"] = np.array(ymin)
158         dic["zmin"] = np.array(zmin)
159         dic["xmax"] = np.array(xmax)
160         dic["ymax"] = np.array(ymax)
161         dic["zmax"] = np.array(zmax)
162
163         return dic
164
165     dic = convert_to_cell_coords(ctx.collect_data())
166
167     X = []
168     rho = []
169     rhovelx = []
170     rhoetot = []
171
172     for i in range(len(dic["xmin"])):
173
174         X.append(dic["xmin"][i])
175         rho.append(dic["rho"][i])
176         rhovelx.append(dic["rhovel"][i][0])
177         rhoetot.append(dic["rhoetot"][i])
178
179     X = np.array(X)
180     rho = np.array(rho)
181     rhovelx = np.array(rhovelx)
182     rhoetot = np.array(rhoetot)
183
184     vx = rhovelx / rho
185
186     fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(9, 6), dpi=125)
187
188     plt.scatter(X, rho, rasterized=True, label="rho")
189     plt.scatter(X, vx, rasterized=True, label="v")
190     plt.scatter(X, (rhoetot - 0.5 * rho * (vx**2)) * (gamma - 1), rasterized=True, label="P")
191     # plt.scatter(X,rhoetot, rasterized=True,label="rhoetot")
192     plt.legend()
193     plt.grid()
194
195     #### add analytical soluce
196     arr_x = np.linspace(xref - xrange, xref + xrange, 1000)
197
198     arr_rho = []
199     arr_P = []
200     arr_vx = []
201
202     for i in range(len(arr_x)):
203         x_ = arr_x[i] - xref
204
205         _rho, _vx, _P = sod.get_value(t_target, x_)
206         arr_rho.append(_rho)
207         arr_vx.append(_vx)
208         arr_P.append(_P)
209
210     plt.plot(arr_x, arr_rho, color="black", label="analytic")
211     plt.plot(arr_x, arr_vx, color="black")
212     plt.plot(arr_x, arr_P, color="black")
213     plt.ylim(-0.1, 1.1)
214     plt.xlim(0.5, 1.5)
215     #######
216     plt.show()
217
218 #################
219 ### Test CD
220 #################
221 rho, v, P = sodanalysis.compute_L2_dist()
222 print(rho, v, P)
223 vx, vy, vz = v
224
225 # normally :
226 # rho 0.07979993131348424
227 # v (0.17970690984930585, 0.0, 0.0)
228 # P 0.12628776652228088
229
230 test_pass = True
231 pass_rho = 0.07979993131348424 + 1e-7
232 pass_vx = 0.17970690984930585 + 1e-7
233 pass_vy = 1e-09
234 pass_vz = 1e-09
235 pass_P = 0.12628776652228088 + 1e-7
236
237 err_log = ""
238
239 if rho > pass_rho:
240     err_log += ("error on rho is too high " + str(rho) + ">" + str(pass_rho)) + "\n"
241     test_pass = False
242 if vx > pass_vx:
243     err_log += ("error on vx is too high " + str(vx) + ">" + str(pass_vx)) + "\n"
244     test_pass = False
245 if vy > pass_vy:
246     err_log += ("error on vy is too high " + str(vy) + ">" + str(pass_vy)) + "\n"
247     test_pass = False
248 if vz > pass_vz:
249     err_log += ("error on vz is too high " + str(vz) + ">" + str(pass_vz)) + "\n"
250     test_pass = False
251 if P > pass_P:
252     err_log += ("error on P is too high " + str(P) + ">" + str(pass_P)) + "\n"
253     test_pass = False
254
255 if test_pass == False:
256     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery