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     x_min, y_min, z_min = rmin
 46     x_max, y_max, z_max = rmax
 47
 48     x = (x_min + x_max) / 2
 49
 50     if x < 0.6 and x > 0.4:
 51         return 2
 52
 53     return 1.0
 54
 55
 56 def rhoe_map(rmin, rmax):
 57     rho = rho_map(rmin, rmax)
 58     return 1.0 * rho
 59
 60
 61 def rhovel_map(rmin, rmax):
 62     rho = rho_map(rmin, rmax)
 63     return (1 * rho, 0 * rho, 0 * rho)
 64
 65
 66 model.set_field_value_lambda_f64("rho", rho_map)
 67 model.set_field_value_lambda_f64("rhoetot", rhoe_map)
 68 model.set_field_value_lambda_f64_3("rhovel", rhovel_map)
 69
 70 model.evolve_until(tmax)
 71
 72
 73 def convert_to_cell_coords(dic):
 74     cmin = dic["cell_min"]
 75     cmax = dic["cell_max"]
 76
 77     xmin = []
 78     ymin = []
 79     zmin = []
 80     xmax = []
 81     ymax = []
 82     zmax = []
 83
 84     for i in range(len(cmin)):
 85         m, M = cmin[i], cmax[i]
 86
 87         mx, my, mz = m
 88         Mx, My, Mz = M
 89
 90         for j in range(8):
 91             a, b = model.get_cell_coords(((mx, my, mz), (Mx, My, Mz)), j)
 92
 93             x, y, z = a
 94             xmin.append(x)
 95             ymin.append(y)
 96             zmin.append(z)
 97
 98             x, y, z = b
 99             xmax.append(x)
100             ymax.append(y)
101             zmax.append(z)
102
103     dic["xmin"] = np.array(xmin)
104     dic["ymin"] = np.array(ymin)
105     dic["zmin"] = np.array(zmin)
106     dic["xmax"] = np.array(xmax)
107     dic["ymax"] = np.array(ymax)
108     dic["zmax"] = np.array(zmax)
109
110     return dic
111
112
113 dic = convert_to_cell_coords(ctx.collect_data())
114
115
116 def get_l2_distance(a, b):
117     return np.sqrt(np.sum((np.array(a) - np.array(b)) ** 2)) / np.sqrt(len(a))
118
119
120 X = []
121
122 rho = []
123 rhovelx = []
124 rhovely = []
125 rhovelz = []
126 rhoetot = []
127
128 rho_ref = []
129 rhovelx_ref = []
130 rhovely_ref = []
131 rhovelz_ref = []
132 rhoetot_ref = []
133
134 for i in range(len(dic["xmin"])):
135     X.append(dic["xmin"][i])
136     rho.append(dic["rho"][i])
137     rhovelx.append(dic["rhovel"][i][0])
138     rhovely.append(dic["rhovel"][i][1])
139     rhovelz.append(dic["rhovel"][i][2])
140     rhoetot.append(dic["rhoetot"][i])
141
142     cell_min = (dic["xmin"][i], dic["ymin"][i], dic["zmin"][i])
143     cell_max = (dic["xmax"][i], dic["ymax"][i], dic["zmax"][i])
144
145     rho_ref.append(rho_map(cell_min, cell_max))
146     rhovelx_ref.append(rhovel_map(cell_min, cell_max)[0])
147     rhovely_ref.append(rhovel_map(cell_min, cell_max)[1])
148     rhovelz_ref.append(rhovel_map(cell_min, cell_max)[2])
149     rhoetot_ref.append(rhoe_map(cell_min, cell_max))
150
151 if do_plot:
152     plt.plot(X, rho, ".", label="rho")
153     plt.plot(X, rho_ref, ".", label="rho_ref")
154
155     plt.plot(X, rhovelx, ".", label="rhovelx")
156     plt.plot(X, rhovelx_ref, ".", label="rhovelx_ref")
157
158     plt.legend()
159     plt.grid()
160     # plt.ylim(0.9,2.5)
161     # plt.xlim(0,1)
162     plt.title("t=" + str(tmax))
163     plt.show()
164
165 l2_rho = get_l2_distance(rho, rho_ref)
166 l2_rhovelx = get_l2_distance(rhovelx, rhovelx_ref)
167 l2_rhovely = get_l2_distance(rhovely, rhovely_ref)
168 l2_rhovelz = get_l2_distance(rhovelz, rhovelz_ref)
169 l2_rhoetot = get_l2_distance(rhoetot, rhoetot_ref)
170
171 print(f"rho: {l2_rho}")
172 print(f"rhovelx: {l2_rhovelx}")
173 print(f"rhovely: {l2_rhovely}")
174 print(f"rhovelz: {l2_rhovelz}")
175 print(f"rhoetot: {l2_rhoetot}")
176
177 expected_l2_rho = 0.12026000336984567
178 expected_l2_rhovelx = 0.12026008009555872
179 expected_l2_rhovely = 0
180 expected_l2_rhovelz = 0
181 expected_l2_rhoetot = 0.12026008010012906
182
183 test_pass = True
184 pass_rho = 1.5e-15
185 pass_rhovelx = 1.5e-15
186 pass_rhovely = 0
187 pass_rhovelz = 0
188 pass_rhoetot = 1.5e-15
189
190 err_log = ""
191
192 if np.abs(l2_rho - expected_l2_rho) > pass_rho:
193     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"
194     test_pass = False
195 if np.abs(l2_rhovelx - expected_l2_rhovelx) > pass_rhovelx:
196     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"
197     test_pass = False
198 if np.abs(l2_rhovely - expected_l2_rhovely) > pass_rhovely:
199     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"
200     test_pass = False
201 if np.abs(l2_rhovelz - expected_l2_rhovelz) > pass_rhovelz:
202     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"
203     test_pass = False
204 if np.abs(l2_rhoetot - expected_l2_rhoetot) > pass_rhoetot:
205     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"
206     test_pass = False
207
208 if test_pass == False:
209     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery