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

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery