add verification script

src/BinaryOrbitCIC/particle_orbit_plot.py

@@ -0,0 +1,38 @@
+# note: the AMReX frontend is broken in yt 4.3.0
+import yt
+from yt.frontends.boxlib.data_structures import AMReXDataset
+
+def particle_dist(plotfiles):
+    t_arr = []
+    err_arr = []
+    d0 = 2.0*3.125e12
+
+    for pltfile in plotfiles:
+        ds = AMReXDataset(pltfile)
+        ad = ds.all_data()
+        x = ad["CIC_particles", "particle_position_x"]
+        y = ad["CIC_particles", "particle_position_y"]
+        z = ad["CIC_particles", "particle_position_z"]
+        dx = x[0] - x[1]
+        dy = y[0] - y[1]
+        dz = z[0] - z[1]
+        from math import sqrt
+        d = sqrt(dx*dx + dy*dy + dz*dz)
+        err = (d-d0)/d0
+        t_arr.append(float(ds.current_time) / 3.15e7)
+        err_arr.append(err)
+
+    return t_arr, err_arr
+
+import glob
+files = glob.glob("plt*")
+files = sorted(files)
+t, err = particle_dist(files)
+
+import matplotlib.pyplot as plt
+plt.figure(figsize=(6,4))
+plt.plot(t, err)
+plt.xlabel("time (yr)")
+plt.ylabel(r"$(d-d_0)/d_0$")
+plt.tight_layout()
+plt.savefig("orbit.png", dpi=150)

tests/BinaryOrbit.in

@@ -21,14 +21,14 @@ amr.n_error_buf     = 3     # minimum 3 cell buffer around tagged cells
 amr.grid_eff        = 1.0
 
 cfl = 0.3
-max_timesteps = 10000
-stop_time = 1.0e7   # s
+max_timesteps = 8000
+stop_time = 5.0e7   # s
 
 do_reflux = 1
 do_subcycle = 0
 do_cic_particles = 1    # turns on CIC particles
 
 ascent_interval = -1
 checkpoint_interval = -1
-plottime_interval = 3.0130519e5 # s
+plottime_interval = 2.0e5 # s
 derived_vars = gpot