Note

Go to the end to download the full example code.

Testing Sod tube with SPH#

CI test for Sod tube with SPH

  8 import matplotlib.pyplot as plt
  9
 10 import shamrock
 11
 12 gamma = 1.4
 13
 14 rho_g = 1
 15 rho_d = 0.125
 16
 17 fact = (rho_g / rho_d) ** (1.0 / 3.0)
 18
 19 P_g = 1
 20 P_d = 0.1
 21
 22 u_g = P_g / ((gamma - 1) * rho_g)
 23 u_d = P_d / ((gamma - 1) * rho_d)
 24
 25 resol = 128
 26
 27 ctx = shamrock.Context()
 28 ctx.pdata_layout_new()
 29
 30 model = shamrock.get_Model_SPH(context=ctx, vector_type="f64_3", sph_kernel="M6")
 31
 32 cfg = model.gen_default_config()
 33 # cfg.set_artif_viscosity_Constant(alpha_u = 1, alpha_AV = 1, beta_AV = 2)
 34 # cfg.set_artif_viscosity_VaryingMM97(alpha_min = 0.1,alpha_max = 1,sigma_decay = 0.1, alpha_u = 1, beta_AV = 2)
 35 cfg.set_artif_viscosity_VaryingCD10(
 36     alpha_min=0.0, alpha_max=1, sigma_decay=0.1, alpha_u=1, beta_AV=2
 37 )
 38 cfg.set_boundary_periodic()
 39 cfg.set_eos_adiabatic(gamma)
 40 cfg.print_status()
 41 model.set_solver_config(cfg)
 42
 43 model.init_scheduler(int(1e8), 1)
 44
 45
 46 (xs, ys, zs) = model.get_box_dim_fcc_3d(1, resol, 24, 24)
 47 dr = 1 / xs
 48 (xs, ys, zs) = model.get_box_dim_fcc_3d(dr, resol, 24, 24)
 49
 50 model.resize_simulation_box((-xs, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
 51
 52
 53 setup = model.get_setup()
 54 gen1 = setup.make_generator_lattice_hcp(dr, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
 55 gen2 = setup.make_generator_lattice_hcp(dr * fact, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
 56 comb = setup.make_combiner_add(gen1, gen2)
 57 # print(comb.get_dot())
 58 setup.apply_setup(comb)
 59
 60 # model.add_cube_fcc_3d(dr, (-xs,-ys/2,-zs/2),(0,ys/2,zs/2))
 61 # model.add_cube_fcc_3d(dr*fact, (0,-ys/2,-zs/2),(xs,ys/2,zs/2))
 62
 63 model.set_value_in_a_box("uint", "f64", u_g, (-xs, -ys / 2, -zs / 2), (0, ys / 2, zs / 2))
 64 model.set_value_in_a_box("uint", "f64", u_d, (0, -ys / 2, -zs / 2), (xs, ys / 2, zs / 2))
 65
 66
 67 vol_b = xs * ys * zs
 68
 69 totmass = (rho_d * vol_b) + (rho_g * vol_b)
 70
 71 print("Total mass :", totmass)
 72
 73
 74 pmass = model.total_mass_to_part_mass(totmass)
 75 model.set_particle_mass(pmass)
 76 print("Current part mass :", pmass)
 77
 78
 79 model.set_cfl_cour(0.1)
 80 model.set_cfl_force(0.1)
 81
 82 t_target = 0.245
 83
 84 print("solver tex :")
 85 print("--------------------------------")
 86 print(model.get_solver_tex())
 87 print("--------------------------------")
 88 print("solver dot graph :")
 89 print("--------------------------------")
 90 print("digraph G {")
 91 print(model.get_solver_dot_graph())
 92 print("}")
 93 print("--------------------------------")
 94
 95 model.evolve_until(t_target)
 96
 97 # model.evolve_once()
 98
 99 sod = shamrock.phys.SodTube(gamma=gamma, rho_1=1, P_1=1, rho_5=0.125, P_5=0.1)
100 sodanalysis = model.make_analysis_sodtube(sod, (1, 0, 0), t_target, 0.0, -0.5, 0.5)
101
102 #################
103 ### Test CD
104 #################
105 rho, v, P = sodanalysis.compute_L2_dist()
106 vx, vy, vz = v
107
108 print("current errors :")
109 print(f"err_rho = {rho}")
110 print(f"err_vx = {vx}")
111 print(f"err_vy = {vy}")
112 print(f"err_vz = {vz}")
113 print(f"err_P = {P}")
114
115 # normally :
116 # rho 0.0001615491818848632
117 # v (0.0011627047434807855, 2.9881306160215856e-05, 1.7413547093275864e-07)
118 # P0.0001248364612976704
119
120 test_pass = True
121
122 expect_rho = 0.00016154918188486815
123 expect_vx = 0.001162704743480841
124 expect_vy = 2.988130616021184e-05
125 expect_vz = 1.7413547093230376e-07
126 expect_P = 0.00012483646129766217
127
128 tol = 1e-11
129
130
131 def float_equal(val1, val2, prec):
132     return abs(val1 - val2) < prec
133
134
135 err_log = ""
136
137 if not float_equal(rho, expect_rho, tol * expect_rho):
138     err_log += "error on rho is outside of tolerances:\n"
139     err_log += f"  expected error = {expect_rho} +- {tol*expect_rho}\n"
140     err_log += f"  obtained error = {rho} (relative error = {(rho - expect_rho)/expect_rho})\n"
141     test_pass = False
142
143 if not float_equal(vx, expect_vx, tol * expect_vx):
144     err_log += "error on vx is outside of tolerances:\n"
145     err_log += f"  expected error = {expect_vx} +- {tol*expect_vx}\n"
146     err_log += f"  obtained error = {vx} (relative error = {(vx - expect_vx)/expect_vx})\n"
147     test_pass = False
148
149 if not float_equal(vy, expect_vy, tol * expect_vy):
150     err_log += "error on vy is outside of tolerances:\n"
151     err_log += f"  expected error = {expect_vy} +- {tol*expect_vy}\n"
152     err_log += f"  obtained error = {vy} (relative error = {(vy - expect_vy)/expect_vy})\n"
153     test_pass = False
154
155 if not float_equal(vz, expect_vz, tol * expect_vz):
156     err_log += "error on vz is outside of tolerances:\n"
157     err_log += f"  expected error = {expect_vz} +- {tol*expect_vz}\n"
158     err_log += f"  obtained error = {vz} (relative error = {(vz - expect_vz)/expect_vz})\n"
159     test_pass = False
160
161 if not float_equal(P, expect_P, tol * expect_P):
162     err_log += "error on P is outside of tolerances:\n"
163     err_log += f"  expected error = {expect_P} +- {tol*expect_P}\n"
164     err_log += f"  obtained error = {P} (relative error = {(P - expect_P)/expect_P})\n"
165     test_pass = False
166
167 if test_pass == False:
168     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery