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 model.evolve_until(t_target)
 85
 86 # model.evolve_once()
 87
 88 sod = shamrock.phys.SodTube(gamma=gamma, rho_1=1, P_1=1, rho_5=0.125, P_5=0.1)
 89 sodanalysis = model.make_analysis_sodtube(sod, (1, 0, 0), t_target, 0.0, -0.5, 0.5)
 90
 91 #################
 92 ### Test CD
 93 #################
 94 rho, v, P = sodanalysis.compute_L2_dist()
 95 vx, vy, vz = v
 96
 97 print("current errors :")
 98 print(f"err_rho = {rho}")
 99 print(f"err_vx = {vx}")
100 print(f"err_vy = {vy}")
101 print(f"err_vz = {vz}")
102 print(f"err_P = {P}")
103
104 # normally :
105 # rho 0.0001615491818848632
106 # v (0.0011627047434807855, 2.9881306160215856e-05, 1.7413547093275864e-07)
107 # P0.0001248364612976704
108
109 test_pass = True
110
111 expect_rho = 0.00016154918188486815
112 expect_vx = 0.001162704743480841
113 expect_vy = 2.988130616021184e-05
114 expect_vz = 1.7413547093230376e-07
115 expect_P = 0.00012483646129766217
116
117 tol = 1e-11
118
119
120 def float_equal(val1, val2, prec):
121     return abs(val1 - val2) < prec
122
123
124 err_log = ""
125
126 if not float_equal(rho, expect_rho, tol * expect_rho):
127     err_log += "error on rho is outside of tolerances:\n"
128     err_log += f"  expected error = {expect_rho} +- {tol*expect_rho}\n"
129     err_log += f"  obtained error = {rho} (relative error = {(rho - expect_rho)/expect_rho})\n"
130     test_pass = False
131
132 if not float_equal(vx, expect_vx, tol * expect_vx):
133     err_log += "error on vx is outside of tolerances:\n"
134     err_log += f"  expected error = {expect_vx} +- {tol*expect_vx}\n"
135     err_log += f"  obtained error = {vx} (relative error = {(vx - expect_vx)/expect_vx})\n"
136     test_pass = False
137
138 if not float_equal(vy, expect_vy, tol * expect_vy):
139     err_log += "error on vy is outside of tolerances:\n"
140     err_log += f"  expected error = {expect_vy} +- {tol*expect_vy}\n"
141     err_log += f"  obtained error = {vy} (relative error = {(vy - expect_vy)/expect_vy})\n"
142     test_pass = False
143
144 if not float_equal(vz, expect_vz, tol * expect_vz):
145     err_log += "error on vz is outside of tolerances:\n"
146     err_log += f"  expected error = {expect_vz} +- {tol*expect_vz}\n"
147     err_log += f"  obtained error = {vz} (relative error = {(vz - expect_vz)/expect_vz})\n"
148     test_pass = False
149
150 if not float_equal(P, expect_P, tol * expect_P):
151     err_log += "error on P is outside of tolerances:\n"
152     err_log += f"  expected error = {expect_P} +- {tol*expect_P}\n"
153     err_log += f"  obtained error = {P} (relative error = {(P - expect_P)/expect_P})\n"
154     test_pass = False
155
156 if test_pass == False:
157     exit("Test did not pass L2 margins : \n" + err_log)

Estimated memory usage: 0 MB

Gallery generated by Sphinx-Gallery