Note

Go to the end to download the full example code.

Testing dusty box with Godunov#

CI test for dusty box with Godunov

  8 from math import *
  9
 10 import matplotlib as mpl
 11 import matplotlib.pyplot as plt
 12 import numpy as np
 13
 14 import shamrock
 15
 16
 17 def run_sim(times, vg_num, vd1_num, vd2_num):
 18     ctx = shamrock.Context()
 19     ctx.pdata_layout_new()
 20
 21     model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")
 22
 23     multx = 1
 24     multy = 1
 25     multz = 1
 26
 27     sz = 1 << 1
 28     base = 2
 29
 30     cfg = model.gen_default_config()
 31     scale_fact = 1 / (sz * base * multx)
 32     cfg.set_scale_factor(scale_fact)
 33     cfg.set_Csafe(
 34         0.44
 35     )  #  TODO : remember to add possibility of different CFL for fluids(e.g Csafe_gas and Csafe_dust ...)
 36     cfg.set_eos_gamma(1.4)  # set adiabatic index gamma , here adiabatic EOS
 37     cfg.set_dust_mode_dhll(2)  # enable dust config
 38     cfg.set_drag_mode_irk1(True)  # enable drag config
 39     cfg.set_face_time_interpolation(False)
 40
 41     # ======= set drag coefficients for test B ========
 42     cfg.set_alpha_values(100)  # ts := 0.01
 43     cfg.set_alpha_values(500)  # ts := 0.002
 44
 45     """
 46     #======= set drag coefficients for test C ========
 47     cfg.set_alpha_values(0.5)          # ts  := 2
 48     cfg.set_alpha_values(1)            # ts  := 1
 49     """
 50
 51     model.set_solver_config(cfg)
 52     model.init_scheduler(int(1e7), 1)
 53     model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz))
 54
 55     # ============= Fileds maps for gas ==============
 56
 57     def rho_map(rmin, rmax):
 58         return 1  # 1 is the initial density
 59
 60     def rhovel_map(rmin, rmax):
 61         return (1, 0, 0)  # vg_x:=1, vg_y:=0, vg_z:=0
 62
 63     def rhoe_map(rmin, rmax):
 64         cs_0 = 1.4
 65         gamma = 1.4
 66         rho_0 = 1
 67         press = (cs_0 * rho_0) / gamma
 68         rhoeint = press / (gamma - 1.0)
 69         rhoekin = 0.5 * rho_0  # vg_x:=1, vg_y:=0, vg_z:=0
 70         return rhoeint + rhoekin
 71
 72     # =========== Fields maps for dust in test B ============
 73     def b_rho_d_1_map(rmin, rmax):
 74         return 1  # rho_d_1 := 1
 75
 76     def b_rho_d_2_map(rmin, rmax):
 77         return 1  # rho_d_2 := 1
 78
 79     def b_rhovel_d_1_map(rmin, rmax):
 80         return (2, 0, 0)  # vd_1_x:=2, vd_1_y:=0, vd_1_z:=0
 81
 82     def b_rhovel_d_2_map(rmin, rmax):
 83         return (0.5, 0, 0)  # vd_2_x:=0.5, vd_2_y:=0, vd_2_z:=0
 84
 85     # =========== Fields maps for dust in test C ============
 86     def c_rho_d_1_map(rmin, rmax):
 87         return 10  # rho_d_1 := 10
 88
 89     def c_rho_d_2_map(rmin, rmax):
 90         return 100  # rho_d_2 := 100
 91
 92     def c_rhovel_d_1_map(rmin, rmax):
 93         return (20, 0, 0)  # vd_1_x:=2, vd_1_y:=0, vd_1_z:=0
 94
 95     def c_rhovel_d_2_map(rmin, rmax):
 96         return (50, 0, 0)  # vd_2_x:=0.5, vd_2_y:=0, vd_2_z:=0
 97
 98     # ============ set init fields values for gas =============
 99
100     model.set_field_value_lambda_f64("rho", rho_map)
101     model.set_field_value_lambda_f64("rhoetot", rhoe_map)
102     model.set_field_value_lambda_f64_3("rhovel", rhovel_map)
103
104     # ============ set init fields values for dusts in test B ==========
105     model.set_field_value_lambda_f64("rho_dust", b_rho_d_1_map, 0)
106     model.set_field_value_lambda_f64_3("rhovel_dust", b_rhovel_d_1_map, 0)
107     model.set_field_value_lambda_f64("rho_dust", b_rho_d_2_map, 1)
108     model.set_field_value_lambda_f64_3("rhovel_dust", b_rhovel_d_2_map, 1)
109
110     def convert_to_cell_coords(dic):
111
112         cmin = dic["cell_min"]
113         cmax = dic["cell_max"]
114
115         xmin = []
116         ymin = []
117         zmin = []
118         xmax = []
119         ymax = []
120         zmax = []
121
122         for i in range(len(cmin)):
123
124             m, M = cmin[i], cmax[i]
125
126             mx, my, mz = m
127             Mx, My, Mz = M
128
129             for j in range(8):
130                 a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j)
131
132                 x, y, z = a
133                 xmin.append(x)
134                 ymin.append(y)
135                 zmin.append(z)
136
137                 x, y, z = b
138                 xmax.append(x)
139                 ymax.append(y)
140                 zmax.append(z)
141
142         dic["xmin"] = np.array(xmin)
143         dic["ymin"] = np.array(ymin)
144         dic["zmin"] = np.array(zmin)
145         dic["xmax"] = np.array(xmax)
146         dic["ymax"] = np.array(ymax)
147         dic["zmax"] = np.array(zmax)
148
149         return dic
150
151     dt = 0.005  # b_dt := 0.005 and c_dt := 0.05
152     t = 0
153     tend = 0.05  # b_tend := 0.05 and c_tend := 0.3
154     freq = 1
155
156     for i in range(13):
157
158         if i % freq == 0:
159             model.dump_vtk("colid_test" + str(i // freq) + ".vtk")
160
161         dic_i = convert_to_cell_coords(ctx.collect_data())
162         vg_num.append(dic_i["rhovel"][0][0] / dic_i["rho"][0])
163         vd1_num.append(dic_i["rhovel_dust"][0][0] / dic_i["rho_dust"][0])
164         vd2_num.append(dic_i["rhovel_dust"][1][0] / dic_i["rho_dust"][1])
165         model.evolve_once_override_time(dt * float(i), dt)
166         t = dt * i
167         times.append(t)
168
169         if t > tend:
170             break
171
172
173 # ============== post treatment ===================
174
175
176 ## ========= analytical function for velocity =====
177 def analytical_velocity(t, Vcom, c1, c2, lambda1, lambda2):
178     return Vcom + c1 * exp(lambda1 * t) + c2 * exp(lambda2 * t)
179
180
181 ## =========Test B setup=========
182 Vcom_B = 1.166666666666667
183 lambda1_B = -141.742430504416
184 lambda2_B = -1058.25756949558
185 cg1_B = -0.35610569612832
186 cg2_B = 0.18943902946166
187 cd11_B = 0.85310244713865
188 cd12_B = -0.01976911380532
189 cd21_B = -0.49699675101033
190 cd22_B = -0.16966991565634
191
192 ## ===== get numerical results ==================
193 times = []
194 vg_num = []
195 vd1_num = []
196 vd2_num = []
197
198 run_sim(times, vg_num, vd1_num, vd2_num)
199
200 ## =========== get analytical results =========
201 vg_anal = [analytical_velocity(t, Vcom_B, cg1_B, cg2_B, lambda1_B, lambda2_B) for t in times]
202 vd1_anal = [analytical_velocity(t, Vcom_B, cd11_B, cd12_B, lambda1_B, lambda2_B) for t in times]
203 vd2_anal = [analytical_velocity(t, Vcom_B, cd21_B, cd22_B, lambda1_B, lambda2_B) for t in times]
204
205
206 print(f"times = {times} with len = {len(times)}\n")
207 print(f" vg_num = {vg_num} with len = {len(vg_num)} \n")
208 print(f" vd1_num = {vd1_num} with len = {len(vd1_num)}  \n")
209 print(f" vd2_num = {vd2_num} with len = {len(vd2_num)} \n")
210
211
212 # ============ plots ======================
213 if False:
214     fig, axs = plt.subplots(1, 3)
215     axs[0].plot(times, vg_num, lw=2, label="vg_num")
216     axs[1].plot(times, vd1_num, lw=2, label="vd1_num")
217     axs[2].plot(times, vd2_num, lw=2, label="vd2_num")
218
219     axs[0].plot(times, vg_anal, "k--", lw=2, label="vg_ana")
220     axs[1].plot(times, vd1_anal, "k--", lw=2, label="vd1_ana")
221     axs[2].plot(times, vd2_anal, "k--", lw=2, label="vd2_ana")
222
223     axs[0].set_xlabel("Time")
224     axs[1].set_xlabel("Time")
225     axs[2].set_xlabel("Time")
226     axs[0].set_ylabel("Velocity")
227
228     axs[0].set_title("$V_{g}$")
229     axs[1].set_title("$V_{d,1}$")
230     axs[2].set_title("$V_{d,2}$")
231
232     axs[0].legend()
233     axs[1].legend()
234     axs[2].legend()
235     plt.savefig("dusty_collision_test_B.png")
236
237
238 # ============ CI test ======================
239 vg_ref = [
240     (1.0),
241     (0.9883720938694548),
242     (1.049486232042273),
243     (1.096048197829739),
244     (1.125013937681472),
245     (1.1422385019131143),
246     (1.1523623051572425),
247     (1.1582940234523524),
248     (1.1617665864783704),
249     (1.1637990418282576),
250     (1.1649885477526485),
251     (1.1656847059714033),
252 ]
253 vd1_ref = [
254     (2.0),
255     (1.6627907475115093),
256     (1.458355952484157),
257     (1.3375867407958744),
258     (1.2667291775092104),
259     (1.2252323221586126),
260     (1.2009423522823592),
261     (1.186726278039555),
262     (1.1784064159717405),
263     (1.1735373262313828),
264     (1.1706877682943513),
265     (1.169020115692882),
266 ]
267 vd2_ref = [
268     (0.5),
269     (0.8488372478527337),
270     (0.9921579951981785),
271     (1.066365330462393),
272     (1.1082572423046926),
273     (1.1325296212244393),
274     (1.1466958752948824),
275     (1.154980318467574),
276     (1.1598277046093364),
277     (1.1626644260265706),
278     (1.1643245650231309),
279     (1.1652961463646916),
280 ]
281
282 vg_diff = [abs(vg_num[i] - vg_ref[i]) for i in range(len(vg_num))]
283 vd1_diff = [abs(vd1_num[i] - vd1_ref[i]) for i in range(len(vd1_num))]
284 vd2_diff = [abs(vd2_num[i] - vd2_ref[i]) for i in range(len(vd2_num))]
285
286
287 print(f"vg_diff = {vg_diff} with len = {len(vg_diff)} \n")
288 print(f"vd1_diff = {vd1_diff} with len = {len(vd1_diff)} \n")
289 print(f"vd2_diff = {vd2_diff} with len = {len(vd2_diff)} \n")
290
291 test_pass = True
292 vg_max_pass = 1e-12
293 vd1_max_pass = 1e-12
294 vd2_max_pass = 1e-12
295
296 err_log = ""
297 if np.max(vg_diff) > vg_max_pass:
298     err_log += f"vg_diff = {np.max(vg_diff)} > {vg_max_pass} \n"
299     test_pass = False
300
301 if np.max(vd1_diff) > vd1_max_pass:
302     err_log += f"vd1_diff = {np.max(vd1_diff)} > {vd1_max_pass} \n"
303     test_pass = False
304
305 if np.max(vd2_diff) > vd2_max_pass:
306     err_log += f"vd2_diff = {np.max(vd2_diff)} > {vd2_max_pass} \n"
307     test_pass = False
308
309 if test_pass == False:
310     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery