From 7605e3b5e7fd40f4745c3391f3d64c900de00893 Mon Sep 17 00:00:00 2001
From: Ben Wibking <ben@wibking.com>
Date: Tue, 23 Jan 2024 21:49:15 -0500
Subject: [PATCH] add verification script

---
 src/BinaryOrbitCIC/particle_orbit_plot.py | 38 +++++++++++++++++++++++
 tests/BinaryOrbit.in                      |  6 ++--
 2 files changed, 41 insertions(+), 3 deletions(-)
 create mode 100644 src/BinaryOrbitCIC/particle_orbit_plot.py

diff --git a/src/BinaryOrbitCIC/particle_orbit_plot.py b/src/BinaryOrbitCIC/particle_orbit_plot.py
new file mode 100644
index 000000000..ca58d2d0f
--- /dev/null
+++ b/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)
diff --git a/tests/BinaryOrbit.in b/tests/BinaryOrbit.in
index be5c466c5..19bf49341 100644
--- a/tests/BinaryOrbit.in
+++ b/tests/BinaryOrbit.in
@@ -21,8 +21,8 @@ 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
@@ -30,5 +30,5 @@ 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