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     x, y, z = rmin
 66     if x < 1:
 67         return 1
 68     else:
 69         return 0.125
 70
 71
 72 etot_L = 1.0 / (gamma - 1)
 73 etot_R = 0.1 / (gamma - 1)
 74
 75
 76 def rhoetot_map(rmin, rmax):
 77     rho = rho_map(rmin, rmax)
 78
 79     x, y, z = rmin
 80     if x < 1:
 81         return etot_L
 82     else:
 83         return etot_R
 84
 85
 86 def rhovel_map(rmin, rmax):
 87     rho = rho_map(rmin, rmax)
 88
 89     return (0, 0, 0)
 90
 91
 92 model.set_field_value_lambda_f64("rho", rho_map)
 93 model.set_field_value_lambda_f64("rhoetot", rhoetot_map)
 94 model.set_field_value_lambda_f64_3("rhovel", rhovel_map)
 95
 96 t_target = 0.245
 97
 98 # for i in range(1000):
 99 #    model.dump_vtk(f"test{i:04d}.vtk")
100 #    model.timestep()
101
102 model.evolve_until(t_target)
103
104 # model.evolve_once()
105 xref = 1.0
106 xrange = 0.5
107 sod = shamrock.phys.SodTube(gamma=gamma, rho_1=1, P_1=1, rho_5=0.125, P_5=0.1)
108 sodanalysis = model.make_analysis_sodtube(sod, (1, 0, 0), t_target, xref, -xrange, xrange)
109
110
111 #################
112 ### Plot
113 #################
114 # do plot or not
115 if True:
116
117     def convert_to_cell_coords(dic):
118         cmin = dic["cell_min"]
119         cmax = dic["cell_max"]
120
121         xmin = []
122         ymin = []
123         zmin = []
124         xmax = []
125         ymax = []
126         zmax = []
127
128         for i in range(len(cmin)):
129             m, M = cmin[i], cmax[i]
130
131             mx, my, mz = m
132             Mx, My, Mz = M
133
134             for j in range(8):
135                 a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j)
136
137                 x, y, z = a
138                 xmin.append(x)
139                 ymin.append(y)
140                 zmin.append(z)
141
142                 x, y, z = b
143                 xmax.append(x)
144                 ymax.append(y)
145                 zmax.append(z)
146
147         dic["xmin"] = np.array(xmin)
148         dic["ymin"] = np.array(ymin)
149         dic["zmin"] = np.array(zmin)
150         dic["xmax"] = np.array(xmax)
151         dic["ymax"] = np.array(ymax)
152         dic["zmax"] = np.array(zmax)
153
154         return dic
155
156     dic = convert_to_cell_coords(ctx.collect_data())
157
158     X = []
159     dX = []
160     rho = []
161     rhovelx = []
162     rhoetot = []
163
164     for i in range(len(dic["xmin"])):
165         X.append(dic["xmin"][i])
166         dX.append(dic["xmax"][i] - dic["xmin"][i])
167         rho.append(dic["rho"][i])
168         rhovelx.append(dic["rhovel"][i][0])
169         rhoetot.append(dic["rhoetot"][i])
170
171     X = np.array(X)
172     dX = np.array(dX)
173     rho = np.array(rho)
174     rhovelx = np.array(rhovelx)
175     rhoetot = np.array(rhoetot)
176
177     vx = rhovelx / rho
178
179     fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(9, 6), dpi=125)
180
181     ax1 = plt.gca()
182     ax2 = ax1.twinx()
183
184     l = -np.log2(dX / np.max(dX)) + 1
185
186     ax1.scatter(X, rho, rasterized=True, label="rho")
187     ax1.scatter(X, vx, rasterized=True, label="v")
188     ax1.scatter(X, (rhoetot - 0.5 * rho * (vx**2)) * (gamma - 1), rasterized=True, label="P")
189     ax2.scatter(X, l, rasterized=True, color="purple", label="AMR level")
190     # plt.scatter(X,rhoetot, rasterized=True,label="rhoetot")
191     ax1.legend(loc=0)
192     ax2.legend(loc=0)
193     ax1.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     ax1.plot(arr_x, arr_rho, color="black", label="analytic")
211     ax1.plot(arr_x, arr_vx, color="black")
212     ax1.plot(arr_x, arr_P, color="black")
213
214     ax1.set_ylim(-0.1, 1.1)
215     ax1.set_xlim(0.5, 1.5)
216     ax2.set_ylabel("AMR level")
217     plt.title(r"$m_{crit}=" + str(mass_crit) + "$")
218     plt.savefig("sod_tube.pdf")
219     plt.savefig("sod_tube.png")
220     #######
221     plt.show()
222
223 #################
224 ### Test CD
225 #################
226 rho, v, P = sodanalysis.compute_L2_dist()
227 print(rho, v, P)
228 vx, vy, vz = v
229
230 # normally :
231 # rho 0.07979993131348424
232 # v (0.17970690984930585, 0.0, 0.0)
233 # P 0.12628776652228088
234
235 test_pass = True
236 pass_rho = 0.07913442601255971 + 1e-7
237 pass_vx = 0.17762998971731672 + 1e-7
238 pass_vy = 0
239 pass_vz = 0
240 pass_P = 0.12516172582510562 + 1e-7
241
242 err_log = ""
243
244 if rho > pass_rho:
245     err_log += ("error on rho is too high " + str(rho) + ">" + str(pass_rho)) + "\n"
246     test_pass = False
247 if vx > pass_vx:
248     err_log += ("error on vx is too high " + str(vx) + ">" + str(pass_vx)) + "\n"
249     test_pass = False
250 if vy > pass_vy:
251     err_log += ("error on vy is too high " + str(vy) + ">" + str(pass_vy)) + "\n"
252     test_pass = False
253 if vz > pass_vz:
254     err_log += ("error on vz is too high " + str(vz) + ">" + str(pass_vz)) + "\n"
255     test_pass = False
256 if P > pass_P:
257     err_log += ("error on P is too high " + str(P) + ">" + str(pass_P)) + "\n"
258     test_pass = False
259
260 if test_pass == False:
261     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery