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         cmin = dic["cell_min"]
112         cmax = dic["cell_max"]
113
114         xmin = []
115         ymin = []
116         zmin = []
117         xmax = []
118         ymax = []
119         zmax = []
120
121         for i in range(len(cmin)):
122             m, M = cmin[i], cmax[i]
123
124             mx, my, mz = m
125             Mx, My, Mz = M
126
127             for j in range(8):
128                 a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j)
129
130                 x, y, z = a
131                 xmin.append(x)
132                 ymin.append(y)
133                 zmin.append(z)
134
135                 x, y, z = b
136                 xmax.append(x)
137                 ymax.append(y)
138                 zmax.append(z)
139
140         dic["xmin"] = np.array(xmin)
141         dic["ymin"] = np.array(ymin)
142         dic["zmin"] = np.array(zmin)
143         dic["xmax"] = np.array(xmax)
144         dic["ymax"] = np.array(ymax)
145         dic["zmax"] = np.array(zmax)
146
147         return dic
148
149     dt = 0.005  # b_dt := 0.005 and c_dt := 0.05
150     t = 0
151     tend = 0.05  # b_tend := 0.05 and c_tend := 0.3
152     freq = 1
153
154     for i in range(13):
155         if i % freq == 0:
156             model.dump_vtk("colid_test" + str(i // freq) + ".vtk")
157
158         dic_i = convert_to_cell_coords(ctx.collect_data())
159         vg_num.append(dic_i["rhovel"][0][0] / dic_i["rho"][0])
160         vd1_num.append(dic_i["rhovel_dust"][0][0] / dic_i["rho_dust"][0])
161         vd2_num.append(dic_i["rhovel_dust"][1][0] / dic_i["rho_dust"][1])
162         model.evolve_once_override_time(dt * float(i), dt)
163         t = dt * i
164         times.append(t)
165
166         if t > tend:
167             break
168
169
170 # ============== post treatment ===================
171
172
173 ## ========= analytical function for velocity =====
174 def analytical_velocity(t, Vcom, c1, c2, lambda1, lambda2):
175     return Vcom + c1 * exp(lambda1 * t) + c2 * exp(lambda2 * t)
176
177
178 ## =========Test B setup=========
179 Vcom_B = 1.166666666666667
180 lambda1_B = -141.742430504416
181 lambda2_B = -1058.25756949558
182 cg1_B = -0.35610569612832
183 cg2_B = 0.18943902946166
184 cd11_B = 0.85310244713865
185 cd12_B = -0.01976911380532
186 cd21_B = -0.49699675101033
187 cd22_B = -0.16966991565634
188
189 ## ===== get numerical results ==================
190 times = []
191 vg_num = []
192 vd1_num = []
193 vd2_num = []
194
195 run_sim(times, vg_num, vd1_num, vd2_num)
196
197 ## =========== get analytical results =========
198 vg_anal = [analytical_velocity(t, Vcom_B, cg1_B, cg2_B, lambda1_B, lambda2_B) for t in times]
199 vd1_anal = [analytical_velocity(t, Vcom_B, cd11_B, cd12_B, lambda1_B, lambda2_B) for t in times]
200 vd2_anal = [analytical_velocity(t, Vcom_B, cd21_B, cd22_B, lambda1_B, lambda2_B) for t in times]
201
202
203 print(f"times = {times} with len = {len(times)}\n")
204 print(f" vg_num = {vg_num} with len = {len(vg_num)} \n")
205 print(f" vd1_num = {vd1_num} with len = {len(vd1_num)}  \n")
206 print(f" vd2_num = {vd2_num} with len = {len(vd2_num)} \n")
207
208
209 # ============ plots ======================
210 if False:
211     fig, axs = plt.subplots(1, 3)
212     axs[0].plot(times, vg_num, lw=2, label="vg_num")
213     axs[1].plot(times, vd1_num, lw=2, label="vd1_num")
214     axs[2].plot(times, vd2_num, lw=2, label="vd2_num")
215
216     axs[0].plot(times, vg_anal, "k--", lw=2, label="vg_ana")
217     axs[1].plot(times, vd1_anal, "k--", lw=2, label="vd1_ana")
218     axs[2].plot(times, vd2_anal, "k--", lw=2, label="vd2_ana")
219
220     axs[0].set_xlabel("Time")
221     axs[1].set_xlabel("Time")
222     axs[2].set_xlabel("Time")
223     axs[0].set_ylabel("Velocity")
224
225     axs[0].set_title("$V_{g}$")
226     axs[1].set_title("$V_{d,1}$")
227     axs[2].set_title("$V_{d,2}$")
228
229     axs[0].legend()
230     axs[1].legend()
231     axs[2].legend()
232     plt.savefig("dusty_collision_test_B.png")
233
234
235 # ============ CI test ======================
236 vg_ref = [
237     (1.0),
238     (0.9883720938694548),
239     (1.049486232042273),
240     (1.096048197829739),
241     (1.125013937681472),
242     (1.1422385019131143),
243     (1.1523623051572425),
244     (1.1582940234523524),
245     (1.1617665864783704),
246     (1.1637990418282576),
247     (1.1649885477526485),
248     (1.1656847059714033),
249 ]
250 vd1_ref = [
251     (2.0),
252     (1.6627907475115093),
253     (1.458355952484157),
254     (1.3375867407958744),
255     (1.2667291775092104),
256     (1.2252323221586126),
257     (1.2009423522823592),
258     (1.186726278039555),
259     (1.1784064159717405),
260     (1.1735373262313828),
261     (1.1706877682943513),
262     (1.169020115692882),
263 ]
264 vd2_ref = [
265     (0.5),
266     (0.8488372478527337),
267     (0.9921579951981785),
268     (1.066365330462393),
269     (1.1082572423046926),
270     (1.1325296212244393),
271     (1.1466958752948824),
272     (1.154980318467574),
273     (1.1598277046093364),
274     (1.1626644260265706),
275     (1.1643245650231309),
276     (1.1652961463646916),
277 ]
278
279 vg_diff = [abs(vg_num[i] - vg_ref[i]) for i in range(len(vg_num))]
280 vd1_diff = [abs(vd1_num[i] - vd1_ref[i]) for i in range(len(vd1_num))]
281 vd2_diff = [abs(vd2_num[i] - vd2_ref[i]) for i in range(len(vd2_num))]
282
283
284 print(f"vg_diff = {vg_diff} with len = {len(vg_diff)} \n")
285 print(f"vd1_diff = {vd1_diff} with len = {len(vd1_diff)} \n")
286 print(f"vd2_diff = {vd2_diff} with len = {len(vd2_diff)} \n")
287
288 test_pass = True
289 vg_max_pass = 1e-12
290 vd1_max_pass = 1e-12
291 vd2_max_pass = 1e-12
292
293 err_log = ""
294 if np.max(vg_diff) > vg_max_pass:
295     err_log += f"vg_diff = {np.max(vg_diff)} > {vg_max_pass} \n"
296     test_pass = False
297
298 if np.max(vd1_diff) > vd1_max_pass:
299     err_log += f"vd1_diff = {np.max(vd1_diff)} > {vd1_max_pass} \n"
300     test_pass = False
301
302 if np.max(vd2_diff) > vd2_max_pass:
303     err_log += f"vd2_diff = {np.max(vd2_diff)} > {vd2_max_pass} \n"
304     test_pass = False
305
306 if test_pass == False:
307     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery