Note

Go to the end to download the full example code.

Testing Sod tube with Godunov with AMR enabled#

CI test for Sod tube with Godunov with AMR enabled

  8 import os
  9
 10 import matplotlib.pyplot as plt
 11 import numpy as np
 12
 13 import shamrock
 14
 15 ctx = shamrock.Context()
 16 ctx.pdata_layout_new()
 17
 18 model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")
 19
 20
 21 multx = 4
 22 multy = 1
 23 multz = 1
 24
 25 cell_size = 1 << 2  # refinement is limited to cell_size = 2
 26 base = 16
 27
 28 cfg = model.gen_default_config()
 29 scale_fact = 2 / (cell_size * base * multx)
 30 cfg.set_scale_factor(scale_fact)
 31
 32 gamma = 1.4
 33 cfg.set_eos_gamma(gamma)
 34 # cfg.set_riemann_solver_rusanov()
 35 cfg.set_riemann_solver_hll()
 36
 37 # cfg.set_slope_lim_none()
 38 # cfg.set_slope_lim_vanleer_f()
 39 # cfg.set_slope_lim_vanleer_std()
 40 # cfg.set_slope_lim_vanleer_sym()
 41 cfg.set_slope_lim_minmod()
 42 cfg.set_face_time_interpolation(True)
 43 mass_crit = 0.0000001 * 5 * 2 * 1.2
 44 cfg.set_amr_mode_density_based(crit_mass=mass_crit)
 45 model.set_solver_config(cfg)
 46
 47
 48 model.init_scheduler(int(1e7), 1)
 49 model.make_base_grid(
 50     (0, 0, 0), (cell_size, cell_size, cell_size), (base * multx, base * multy, base * multz)
 51 )
 52
 53 # without face time interpolation
 54 # 0.07979993131348424 (0.17970690984930585, 0.0, 0.0) 0.12628776652228088
 55
 56 # with face time interpolation
 57 # 0.07894793711859852 (0.17754462339166546, 0.0, 0.0) 0.12498304725061045
 58
 59
 60 kx, ky, kz = 2 * np.pi, 0, 0
 61 delta_rho = 1e-2
 62
 63
 64 def rho_map(rmin, rmax):
 65
 66     x, y, z = rmin
 67     if x < 1:
 68         return 1
 69     else:
 70         return 0.125
 71
 72
 73 etot_L = 1.0 / (gamma - 1)
 74 etot_R = 0.1 / (gamma - 1)
 75
 76
 77 def rhoetot_map(rmin, rmax):
 78
 79     rho = rho_map(rmin, rmax)
 80
 81     x, y, z = rmin
 82     if x < 1:
 83         return etot_L
 84     else:
 85         return etot_R
 86
 87
 88 def rhovel_map(rmin, rmax):
 89     rho = rho_map(rmin, rmax)
 90
 91     return (0, 0, 0)
 92
 93
 94 model.set_field_value_lambda_f64("rho", rho_map)
 95 model.set_field_value_lambda_f64("rhoetot", rhoetot_map)
 96 model.set_field_value_lambda_f64_3("rhovel", rhovel_map)
 97
 98 t_target = 0.245
 99
100 # for i in range(1000):
101 #    model.dump_vtk(f"test{i:04d}.vtk")
102 #    model.timestep()
103
104 model.evolve_until(t_target)
105
106 # model.evolve_once()
107 xref = 1.0
108 xrange = 0.5
109 sod = shamrock.phys.SodTube(gamma=gamma, rho_1=1, P_1=1, rho_5=0.125, P_5=0.1)
110 sodanalysis = model.make_analysis_sodtube(sod, (1, 0, 0), t_target, xref, -xrange, xrange)
111
112
113 #################
114 ### Plot
115 #################
116 # do plot or not
117 if True:
118
119     def convert_to_cell_coords(dic):
120
121         cmin = dic["cell_min"]
122         cmax = dic["cell_max"]
123
124         xmin = []
125         ymin = []
126         zmin = []
127         xmax = []
128         ymax = []
129         zmax = []
130
131         for i in range(len(cmin)):
132
133             m, M = cmin[i], cmax[i]
134
135             mx, my, mz = m
136             Mx, My, Mz = M
137
138             for j in range(8):
139                 a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j)
140
141                 x, y, z = a
142                 xmin.append(x)
143                 ymin.append(y)
144                 zmin.append(z)
145
146                 x, y, z = b
147                 xmax.append(x)
148                 ymax.append(y)
149                 zmax.append(z)
150
151         dic["xmin"] = np.array(xmin)
152         dic["ymin"] = np.array(ymin)
153         dic["zmin"] = np.array(zmin)
154         dic["xmax"] = np.array(xmax)
155         dic["ymax"] = np.array(ymax)
156         dic["zmax"] = np.array(zmax)
157
158         return dic
159
160     dic = convert_to_cell_coords(ctx.collect_data())
161
162     X = []
163     dX = []
164     rho = []
165     rhovelx = []
166     rhoetot = []
167
168     for i in range(len(dic["xmin"])):
169
170         X.append(dic["xmin"][i])
171         dX.append(dic["xmax"][i] - dic["xmin"][i])
172         rho.append(dic["rho"][i])
173         rhovelx.append(dic["rhovel"][i][0])
174         rhoetot.append(dic["rhoetot"][i])
175
176     X = np.array(X)
177     dX = np.array(dX)
178     rho = np.array(rho)
179     rhovelx = np.array(rhovelx)
180     rhoetot = np.array(rhoetot)
181
182     vx = rhovelx / rho
183
184     fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(9, 6), dpi=125)
185
186     ax1 = plt.gca()
187     ax2 = ax1.twinx()
188
189     l = -np.log2(dX / np.max(dX)) + 1
190
191     ax1.scatter(X, rho, rasterized=True, label="rho")
192     ax1.scatter(X, vx, rasterized=True, label="v")
193     ax1.scatter(X, (rhoetot - 0.5 * rho * (vx**2)) * (gamma - 1), rasterized=True, label="P")
194     ax2.scatter(X, l, rasterized=True, color="purple", label="AMR level")
195     # plt.scatter(X,rhoetot, rasterized=True,label="rhoetot")
196     ax1.legend(loc=0)
197     ax2.legend(loc=0)
198     ax1.grid()
199
200     #### add analytical soluce
201     arr_x = np.linspace(xref - xrange, xref + xrange, 1000)
202
203     arr_rho = []
204     arr_P = []
205     arr_vx = []
206
207     for i in range(len(arr_x)):
208         x_ = arr_x[i] - xref
209
210         _rho, _vx, _P = sod.get_value(t_target, x_)
211         arr_rho.append(_rho)
212         arr_vx.append(_vx)
213         arr_P.append(_P)
214
215     ax1.plot(arr_x, arr_rho, color="black", label="analytic")
216     ax1.plot(arr_x, arr_vx, color="black")
217     ax1.plot(arr_x, arr_P, color="black")
218
219     ax1.set_ylim(-0.1, 1.1)
220     ax1.set_xlim(0.5, 1.5)
221     ax2.set_ylabel("AMR level")
222     plt.title(r"$m_{crit}=" + str(mass_crit) + "$")
223     plt.savefig("sod_tube.pdf")
224     plt.savefig("sod_tube.png")
225     #######
226     plt.show()
227
228 #################
229 ### Test CD
230 #################
231 rho, v, P = sodanalysis.compute_L2_dist()
232 print(rho, v, P)
233 vx, vy, vz = v
234
235 # normally :
236 # rho 0.07979993131348424
237 # v (0.17970690984930585, 0.0, 0.0)
238 # P 0.12628776652228088
239
240 test_pass = True
241 pass_rho = 0.07913442601255971 + 1e-7
242 pass_vx = 0.17762998971731672 + 1e-7
243 pass_vy = 0
244 pass_vz = 0
245 pass_P = 0.12516172582510562 + 1e-7
246
247 err_log = ""
248
249 if rho > pass_rho:
250     err_log += ("error on rho is too high " + str(rho) + ">" + str(pass_rho)) + "\n"
251     test_pass = False
252 if vx > pass_vx:
253     err_log += ("error on vx is too high " + str(vx) + ">" + str(pass_vx)) + "\n"
254     test_pass = False
255 if vy > pass_vy:
256     err_log += ("error on vy is too high " + str(vy) + ">" + str(pass_vy)) + "\n"
257     test_pass = False
258 if vz > pass_vz:
259     err_log += ("error on vz is too high " + str(vz) + ">" + str(pass_vz)) + "\n"
260     test_pass = False
261 if P > pass_P:
262     err_log += ("error on P is too high " + str(P) + ">" + str(pass_P)) + "\n"
263     test_pass = False
264
265 if test_pass == False:
266     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery