Note

Go to the end to download the full example code.

Godunov advection test#

This test is used to test the Godunov advection setup

  8 import os
  9
 10 import matplotlib.pyplot as plt
 11 import numpy as np
 12
 13 import shamrock
 14
 15 tmax = 1.0
 16 do_plot = True
 17
 18 multx = 1
 19 multy = 1
 20 multz = 1
 21
 22 sz = 1 << 1
 23 base = 32
 24
 25 ctx = shamrock.Context()
 26 ctx.pdata_layout_new()
 27
 28 model = shamrock.get_Model_Ramses(context=ctx, vector_type="f64_3", grid_repr="i64_3")
 29
 30 cfg = model.gen_default_config()
 31 scale_fact = 1 / (sz * base * multx)
 32 cfg.set_scale_factor(scale_fact)
 33 cfg.set_eos_gamma(1.000001)
 34 model.set_solver_config(cfg)
 35
 36 model.init_scheduler(int(1e7), 1)
 37 model.make_base_grid((0, 0, 0), (sz, sz, sz), (base * multx, base * multy, base * multz))
 38
 39 kx, ky, kz = 2 * np.pi, 0, 0
 40 delta_rho = 0
 41 delta_v = 1e-5
 42
 43
 44 def rho_map(rmin, rmax):
 45
 46     x_min, y_min, z_min = rmin
 47     x_max, y_max, z_max = rmax
 48
 49     x = (x_min + x_max) / 2
 50
 51     if x < 0.6 and x > 0.4:
 52         return 2
 53
 54     return 1.0
 55
 56
 57 def rhoe_map(rmin, rmax):
 58     rho = rho_map(rmin, rmax)
 59     return 1.0 * rho
 60
 61
 62 def rhovel_map(rmin, rmax):
 63     rho = rho_map(rmin, rmax)
 64     return (1 * rho, 0 * rho, 0 * rho)
 65
 66
 67 model.set_field_value_lambda_f64("rho", rho_map)
 68 model.set_field_value_lambda_f64("rhoetot", rhoe_map)
 69 model.set_field_value_lambda_f64_3("rhovel", rhovel_map)
 70
 71 model.evolve_until(tmax)
 72
 73
 74 def convert_to_cell_coords(dic):
 75
 76     cmin = dic["cell_min"]
 77     cmax = dic["cell_max"]
 78
 79     xmin = []
 80     ymin = []
 81     zmin = []
 82     xmax = []
 83     ymax = []
 84     zmax = []
 85
 86     for i in range(len(cmin)):
 87
 88         m, M = cmin[i], cmax[i]
 89
 90         mx, my, mz = m
 91         Mx, My, Mz = M
 92
 93         for j in range(8):
 94             a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j)
 95
 96             x, y, z = a
 97             xmin.append(x)
 98             ymin.append(y)
 99             zmin.append(z)
100
101             x, y, z = b
102             xmax.append(x)
103             ymax.append(y)
104             zmax.append(z)
105
106     dic["xmin"] = np.array(xmin)
107     dic["ymin"] = np.array(ymin)
108     dic["zmin"] = np.array(zmin)
109     dic["xmax"] = np.array(xmax)
110     dic["ymax"] = np.array(ymax)
111     dic["zmax"] = np.array(zmax)
112
113     return dic
114
115
116 dic = convert_to_cell_coords(ctx.collect_data())
117
118
119 def get_l2_distance(a, b):
120     return np.sqrt(np.sum((np.array(a) - np.array(b)) ** 2)) / np.sqrt(len(a))
121
122
123 X = []
124
125 rho = []
126 rhovelx = []
127 rhovely = []
128 rhovelz = []
129 rhoetot = []
130
131 rho_ref = []
132 rhovelx_ref = []
133 rhovely_ref = []
134 rhovelz_ref = []
135 rhoetot_ref = []
136
137 for i in range(len(dic["xmin"])):
138     X.append(dic["xmin"][i])
139     rho.append(dic["rho"][i])
140     rhovelx.append(dic["rhovel"][i][0])
141     rhovely.append(dic["rhovel"][i][1])
142     rhovelz.append(dic["rhovel"][i][2])
143     rhoetot.append(dic["rhoetot"][i])
144
145     cell_min = (dic["xmin"][i], dic["ymin"][i], dic["zmin"][i])
146     cell_max = (dic["xmax"][i], dic["ymax"][i], dic["zmax"][i])
147
148     rho_ref.append(rho_map(cell_min, cell_max))
149     rhovelx_ref.append(rhovel_map(cell_min, cell_max)[0])
150     rhovely_ref.append(rhovel_map(cell_min, cell_max)[1])
151     rhovelz_ref.append(rhovel_map(cell_min, cell_max)[2])
152     rhoetot_ref.append(rhoe_map(cell_min, cell_max))
153
154 if do_plot:
155     plt.plot(X, rho, ".", label="rho")
156     plt.plot(X, rho_ref, ".", label="rho_ref")
157
158     plt.plot(X, rhovelx, ".", label="rhovelx")
159     plt.plot(X, rhovelx_ref, ".", label="rhovelx_ref")
160
161     plt.legend()
162     plt.grid()
163     # plt.ylim(0.9,2.5)
164     # plt.xlim(0,1)
165     plt.title("t=" + str(tmax))
166     plt.show()
167
168 l2_rho = get_l2_distance(rho, rho_ref)
169 l2_rhovelx = get_l2_distance(rhovelx, rhovelx_ref)
170 l2_rhovely = get_l2_distance(rhovely, rhovely_ref)
171 l2_rhovelz = get_l2_distance(rhovelz, rhovelz_ref)
172 l2_rhoetot = get_l2_distance(rhoetot, rhoetot_ref)
173
174 print(f"rho: {l2_rho}")
175 print(f"rhovelx: {l2_rhovelx}")
176 print(f"rhovely: {l2_rhovely}")
177 print(f"rhovelz: {l2_rhovelz}")
178 print(f"rhoetot: {l2_rhoetot}")
179
180 expected_l2_rho = 0.12026000336984567
181 expected_l2_rhovelx = 0.12026008009555872
182 expected_l2_rhovely = 0
183 expected_l2_rhovelz = 0
184 expected_l2_rhoetot = 0.12026008010012906
185
186 test_pass = True
187 pass_rho = 1.5e-15
188 pass_rhovelx = 1.5e-15
189 pass_rhovely = 0
190 pass_rhovelz = 0
191 pass_rhoetot = 1.5e-15
192
193 err_log = ""
194
195 if np.abs(l2_rho - expected_l2_rho) > pass_rho:
196     err_log += f"error on rho is too far from expected value: {l2_rho} != {expected_l2_rho}, delta: {np.abs(l2_rho - expected_l2_rho)}\n"
197     test_pass = False
198 if np.abs(l2_rhovelx - expected_l2_rhovelx) > pass_rhovelx:
199     err_log += f"error on rhovelx is too far from expected value: {l2_rhovelx} != {expected_l2_rhovelx}, delta: {np.abs(l2_rhovelx - expected_l2_rhovelx)}\n"
200     test_pass = False
201 if np.abs(l2_rhovely - expected_l2_rhovely) > pass_rhovely:
202     err_log += f"error on rhovely is too far from expected value: {l2_rhovely} != {expected_l2_rhovely}, delta: {np.abs(l2_rhovely - expected_l2_rhovely)}\n"
203     test_pass = False
204 if np.abs(l2_rhovelz - expected_l2_rhovelz) > pass_rhovelz:
205     err_log += f"error on rhovelz is too far from expected value: {l2_rhovelz} != {expected_l2_rhovelz}, delta: {np.abs(l2_rhovelz - expected_l2_rhovelz)}\n"
206     test_pass = False
207 if np.abs(l2_rhoetot - expected_l2_rhoetot) > pass_rhoetot:
208     err_log += f"error on rhoetot is too far from expected value: {l2_rhoetot} != {expected_l2_rhoetot}, delta: {np.abs(l2_rhoetot - expected_l2_rhoetot)}\n"
209     test_pass = False
210
211 if test_pass == False:
212     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery