diff --git a/examples/heat_1d_ap_direct.rs b/examples/heat_1d_ap_direct.rs
index 8e63af8..caefa8f 100644
--- a/examples/heat_1d_ap_direct.rs
+++ b/examples/heat_1d_ap_direct.rs
@@ -22,6 +22,8 @@ fn main() {
         i.add_line(0, input_domain.buffer());
         img = Some(i);
     }
+
+    let mut global_time = 0;
     for t in 1..args.lines as u32 {
         box_apply(
             &bc,
@@ -29,8 +31,10 @@ fn main() {
             &mut input_domain,
             &mut output_domain,
             args.steps_per_line,
+            global_time,
             args.chunk_size,
         );
+        global_time += args.steps_per_line;
         std::mem::swap(&mut input_domain, &mut output_domain);
         if let Some(i) = img.as_mut() {
             i.add_line(t, input_domain.buffer());
diff --git a/examples/heat_1d_ap_fft.rs b/examples/heat_1d_ap_fft.rs
index 9af8050..f06669b 100644
--- a/examples/heat_1d_ap_fft.rs
+++ b/examples/heat_1d_ap_fft.rs
@@ -6,9 +6,9 @@ fn main() {
     let (args, output_image_path) = Args::cli_parse("heat_1d_ap_fft");
 
     let stencil = nhls::standard_stencils::heat_1d(1.0, 1.0, 0.5);
+    let grid_bound = args.grid_bounds();
 
     // Create domains
-    let grid_bound = args.grid_bounds();
     let mut buffer_1 = OwnedDomain::new(grid_bound);
     let mut buffer_2 = OwnedDomain::new(grid_bound);
     let mut input_domain = buffer_1.as_slice_domain();
@@ -44,8 +44,11 @@ fn main() {
         i.add_line(0, input_domain.buffer());
         img = Some(i);
     }
+
+    let mut global_time = 0;
     for t in 1..args.lines as u32 {
-        solver.apply(&mut input_domain, &mut output_domain);
+        solver.apply(&mut input_domain, &mut output_domain, global_time);
+        global_time += args.steps_per_line;
         std::mem::swap(&mut input_domain, &mut output_domain);
         if let Some(i) = img.as_mut() {
             i.add_line(t, input_domain.buffer());
diff --git a/examples/heat_2d_ap_direct.rs b/examples/heat_2d_ap_direct.rs
index 93fcf18..5297432 100644
--- a/examples/heat_2d_ap_direct.rs
+++ b/examples/heat_2d_ap_direct.rs
@@ -1,7 +1,6 @@
 use nhls::domain::*;
 use nhls::image::*;
 use nhls::image_2d_example::*;
-use nhls::init;
 
 fn main() {
     let args = Args::cli_parse("heat_2d_ap_direct");
@@ -15,12 +14,6 @@ fn main() {
     let mut input_domain = OwnedDomain::new(grid_bound);
     let mut output_domain = OwnedDomain::new(grid_bound);
 
-    if args.rand_init {
-        init::rand(&mut input_domain, 1024, args.chunk_size);
-    } else {
-        init::normal_ic_2d(&mut input_domain, args.chunk_size);
-    }
-
     if args.write_images {
         image2d(&input_domain, &args.frame_name(0));
     }
@@ -29,6 +22,7 @@ fn main() {
     let bc = ConstantCheck::new(0.0, grid_bound);
 
     // Apply direct solver
+    let mut global_time = 0;
     for t in 1..args.images {
         nhls::solver::direct::box_apply(
             &bc,
@@ -36,9 +30,10 @@ fn main() {
             &mut input_domain,
             &mut output_domain,
             args.steps_per_image,
+            global_time,
             args.chunk_size,
         );
-
+        global_time += args.steps_per_image;
         std::mem::swap(&mut input_domain, &mut output_domain);
         if args.write_images {
             image2d(&input_domain, &args.frame_name(t));
diff --git a/examples/heat_2d_ap_fft.rs b/examples/heat_2d_ap_fft.rs
index 457006e..bd513e7 100644
--- a/examples/heat_2d_ap_fft.rs
+++ b/examples/heat_2d_ap_fft.rs
@@ -7,16 +7,7 @@ fn main() {
     let args = Args::cli_parse("heat_2d_ap_fft");
 
     let stencil = nhls::standard_stencils::heat_2d(1.0, 1.0, 1.0, 0.2, 0.2);
-
-    // Create domains
     let grid_bound = args.grid_bounds();
-    let mut buffer_1 = OwnedDomain::new(grid_bound);
-    let mut buffer_2 = OwnedDomain::new(grid_bound);
-    let mut input_domain = buffer_1.as_slice_domain();
-    let mut output_domain = buffer_2.as_slice_domain();
-    if args.write_images {
-        image2d(&input_domain, &args.frame_name(0));
-    }
 
     // Create BC
     let bc = ConstantCheck::new(1.0, grid_bound);
@@ -45,8 +36,19 @@ fn main() {
         solver.scratch_descriptor_file(&d_path);
     }
 
+    // Create domains
+    let mut buffer_1 = OwnedDomain::new(grid_bound);
+    let mut buffer_2 = OwnedDomain::new(grid_bound);
+    let mut input_domain = buffer_1.as_slice_domain();
+    let mut output_domain = buffer_2.as_slice_domain();
+    if args.write_images {
+        image2d(&input_domain, &args.frame_name(0));
+    }
+
+    let mut global_time = 0;
     for t in 1..args.images {
-        solver.apply(&mut input_domain, &mut output_domain);
+        solver.apply(&mut input_domain, &mut output_domain, global_time);
+        global_time += args.steps_per_image;
         std::mem::swap(&mut input_domain, &mut output_domain);
         if args.write_images {
             image2d(&input_domain, &args.frame_name(t));
diff --git a/examples/heat_3d_ap_fft.rs b/examples/heat_3d_ap_fft.rs
index dc89994..e20707d 100644
--- a/examples/heat_3d_ap_fft.rs
+++ b/examples/heat_3d_ap_fft.rs
@@ -43,8 +43,10 @@ fn main() {
         write_vtk3d(&input_domain, &args.frame_name(0));
     }
 
+    let mut global_time = 0;
     for t in 1..args.images {
-        solver.apply(&mut input_domain, &mut output_domain);
+        solver.apply(&mut input_domain, &mut output_domain, global_time);
+        global_time += args.steps_per_image;
         std::mem::swap(&mut input_domain, &mut output_domain);
         if args.write_images {
             write_vtk3d(&input_domain, &args.frame_name(t));
diff --git a/examples/time_varying_2d.rs b/examples/time_varying_2d.rs
new file mode 100644
index 0000000..9792b6b
--- /dev/null
+++ b/examples/time_varying_2d.rs
@@ -0,0 +1,117 @@
+use nhls::domain::*;
+use nhls::fft_solver::*;
+use nhls::image::*;
+use nhls::image_2d_example::*;
+use nhls::stencil::*;
+use nhls::util::*;
+
+pub struct PulseBC<const GRID_DIMENSION: usize> {
+    rate: f64,
+    n_f: f64,
+    sigma_sq: f64,
+    aabb: AABB<GRID_DIMENSION>,
+}
+
+impl<const GRID_DIMENSION: usize> PulseBC<GRID_DIMENSION> {
+    pub fn new(rate: f64, aabb: AABB<GRID_DIMENSION>) -> Self {
+        let n_f = aabb.exclusive_bounds()[1] as f64;
+        let sigma_sq: f64 = (n_f / 10.0) * (n_f / 10.0);
+        PulseBC {
+            rate,
+            aabb,
+            n_f,
+            sigma_sq,
+        }
+    }
+}
+
+impl<const GRID_DIMENSION: usize> BCCheck<GRID_DIMENSION>
+    for PulseBC<GRID_DIMENSION>
+{
+    fn check(
+        &self,
+        coord: &Coord<GRID_DIMENSION>,
+        global_time: usize,
+    ) -> Option<f64> {
+        if self.aabb.contains(coord) {
+            None
+        } else if coord[0] < self.aabb.bounds[(0, 0)] {
+            let x = (coord[1] as f64) - (self.n_f / 2.0);
+            let exp = -x * x / (2.0 * self.sigma_sq);
+            let normal_value = exp.exp();
+            let sin_mod = 0.5 * ((global_time as f64 * self.rate).sin() + 1.0);
+            Some(sin_mod * normal_value)
+        } else {
+            Some(0.0)
+        }
+    }
+}
+
+fn main() {
+    let args = Args::cli_parse("time_varying_2d");
+
+    let stencil = Stencil::new(
+        [[0, 0], [-1, 0], [1, 0], [0, -1], [0, 1]],
+        move |args: &[f64; 5]| {
+            let middle = args[0];
+            let left = args[1];
+            let right = args[2];
+            let bottom = args[3];
+            let top = args[4];
+            0.1 * middle + 0.4 * left + 0.1 * right + 0.2 * top + 0.2 * bottom
+        },
+    );
+
+    let w = stencil.extract_weights();
+    let s: f64 = w.iter().sum();
+    println!("{:?}, {}", w, s);
+    let grid_bound = args.grid_bounds();
+
+    // Create BC
+    let bc = PulseBC::new((2.0 * std::f64::consts::PI) / 1000.0, grid_bound);
+
+    // Create AP Solver
+    let planner_params = PlannerParameters {
+        plan_type: args.plan_type,
+        cutoff: 40,
+        ratio: 0.5,
+        chunk_size: args.chunk_size,
+    };
+    let solver = APSolver::new(
+        &bc,
+        &stencil,
+        grid_bound,
+        args.steps_per_image,
+        &planner_params,
+    );
+    if args.write_dot {
+        let mut dot_path = args.output_dir.clone();
+        dot_path.push("plan.dot");
+        solver.to_dot_file(&dot_path);
+
+        let mut d_path = args.output_dir.clone();
+        d_path.push("scratch.txt");
+        solver.scratch_descriptor_file(&d_path);
+    }
+
+    // Create domains
+    let mut buffer_1 = OwnedDomain::new(grid_bound);
+    let mut buffer_2 = OwnedDomain::new(grid_bound);
+    let mut input_domain = buffer_1.as_slice_domain();
+    let mut output_domain = buffer_2.as_slice_domain();
+    if args.write_images {
+        image2d(&input_domain, &args.frame_name(0));
+    }
+
+    let mut global_time = 0;
+    for t in 1..args.images {
+        solver.apply(&mut input_domain, &mut output_domain, global_time);
+        global_time += args.steps_per_image;
+        std::mem::swap(&mut input_domain, &mut output_domain);
+        if args.write_images {
+            image2d(&input_domain, &args.frame_name(t));
+        }
+    }
+
+    args.save_wisdom();
+}
diff --git a/src/domain/bc/constant.rs b/src/domain/bc/constant.rs
index 1a9ec79..2421f06 100644
--- a/src/domain/bc/constant.rs
+++ b/src/domain/bc/constant.rs
@@ -15,7 +15,11 @@ impl<const GRID_DIMENSION: usize> ConstantCheck<GRID_DIMENSION> {
 impl<const GRID_DIMENSION: usize> BCCheck<GRID_DIMENSION>
     for ConstantCheck<GRID_DIMENSION>
 {
-    fn check(&self, coord: &Coord<GRID_DIMENSION>) -> Option<f64> {
+    fn check(
+        &self,
+        coord: &Coord<GRID_DIMENSION>,
+        _global_time: usize,
+    ) -> Option<f64> {
         if self.bound.contains(coord) {
             return None;
         }
@@ -39,18 +43,18 @@ mod unit_tests {
         }
         let bc = ConstantCheck::new(-1.0, bound);
         for i in 0..n_r {
-            let v = bc.check(&vector![i as i32]);
+            let v = bc.check(&vector![i as i32], 0);
             assert_eq!(v, None);
         }
 
         {
-            let v = bc.check(&vector![-1]);
+            let v = bc.check(&vector![-1], 1);
             assert!(v.is_some());
             assert_approx_eq!(f64, v.unwrap(), -1.0);
         }
 
         {
-            let v = bc.check(&vector![11]);
+            let v = bc.check(&vector![11], 2);
             assert!(v.is_some());
             assert_approx_eq!(f64, v.unwrap(), -1.0);
         }
diff --git a/src/domain/bc/mod.rs b/src/domain/bc/mod.rs
index da9388f..e59c723 100644
--- a/src/domain/bc/mod.rs
+++ b/src/domain/bc/mod.rs
@@ -7,5 +7,9 @@ pub use periodic::*;
 use crate::util::*;
 
 pub trait BCCheck<const GRID_DIMENSION: usize>: Sync {
-    fn check(&self, world_coord: &Coord<GRID_DIMENSION>) -> Option<f64>;
+    fn check(
+        &self,
+        world_coord: &Coord<GRID_DIMENSION>,
+        global_time: usize,
+    ) -> Option<f64>;
 }
diff --git a/src/domain/bc/periodic.rs b/src/domain/bc/periodic.rs
index ca02895..508e17d 100644
--- a/src/domain/bc/periodic.rs
+++ b/src/domain/bc/periodic.rs
@@ -23,7 +23,11 @@ impl<
 impl<const GRID_DIMENSION: usize, DomainType: DomainView<GRID_DIMENSION>>
     BCCheck<GRID_DIMENSION> for PeriodicCheck<'_, GRID_DIMENSION, DomainType>
 {
-    fn check(&self, world_coord: &Coord<GRID_DIMENSION>) -> Option<f64> {
+    fn check(
+        &self,
+        world_coord: &Coord<GRID_DIMENSION>,
+        _global_time: usize,
+    ) -> Option<f64> {
         let p_coord = &self.domain.aabb().periodic_coord(world_coord);
         if p_coord != world_coord {
             return Some(self.domain.view(p_coord));
@@ -47,18 +51,18 @@ mod unit_tests {
             domain.par_set_values(|coord| coord[0] as f64, 1);
             let bc = PeriodicCheck::new(&domain);
             for (i, _) in domain.buffer().iter().enumerate() {
-                let v = bc.check(&vector![i as i32]);
+                let v = bc.check(&vector![i as i32], 0);
                 assert_eq!(v, None);
             }
 
             {
-                let v = bc.check(&vector![-1]);
+                let v = bc.check(&vector![-1], 1);
                 assert!(v.is_some());
                 assert_approx_eq!(f64, v.unwrap(), 10.0);
             }
 
             {
-                let v = bc.check(&vector![11]);
+                let v = bc.check(&vector![11], 2);
                 assert!(v.is_some());
                 assert_approx_eq!(f64, v.unwrap(), 0.0);
             }
diff --git a/src/domain/gather_args.rs b/src/domain/gather_args.rs
index 5b9be48..8522408 100644
--- a/src/domain/gather_args.rs
+++ b/src/domain/gather_args.rs
@@ -13,6 +13,7 @@ pub fn gather_args<
     bc: &BC,
     input: &DomainType,
     world_coord: &Coord<GRID_DIMENSION>,
+    global_time: usize,
 ) -> [f64; NEIGHBORHOOD_SIZE]
 where
     Operation: StencilOperation<f64, NEIGHBORHOOD_SIZE>,
@@ -22,7 +23,7 @@ where
     for (i, n_i) in stencil.offsets().iter().enumerate() {
         let n_world_coord = world_coord + n_i;
         result[i] = bc
-            .check(&n_world_coord)
+            .check(&n_world_coord, global_time)
             .unwrap_or_else(|| input.view(&n_world_coord));
     }
     result
@@ -47,7 +48,7 @@ mod unit_tests {
             [[0, -1], [0, 1], [1, 0], [-1, 0], [0, 0]],
             |_: &[f64; 5]| -1.0,
         );
-        let r = gather_args(&stencil, &bc, &domain, &vector![9, 9]);
+        let r = gather_args(&stencil, &bc, &domain, &vector![9, 9], 0);
         let e = [
             (9 + 3 * 8) as f64,
             -4.0,
@@ -76,7 +77,7 @@ mod unit_tests {
             [[0, -1], [0, 1], [1, 0], [-1, 0], [0, 0]],
             |_: &[f64; 5]| -1.0,
         );
-        let r = gather_args(&stencil, &bc, &domain, &vector![9, 9]);
+        let r = gather_args(&stencil, &bc, &domain, &vector![9, 9], 11);
         let e = [
             (9 + 3 * 8) as f64,
             9.0,
diff --git a/src/fft_solver/ap_solver.rs b/src/fft_solver/ap_solver.rs
index 0d4afad..5e6ef15 100644
--- a/src/fft_solver/ap_solver.rs
+++ b/src/fft_solver/ap_solver.rs
@@ -76,8 +76,9 @@ where
         &self,
         input_domain: &mut SliceDomain<'a, GRID_DIMENSION>,
         output_domain: &mut SliceDomain<'a, GRID_DIMENSION>,
+        global_time: usize,
     ) {
-        self.solve_root(input_domain, output_domain);
+        self.solve_root(input_domain, output_domain, global_time);
     }
 
     pub fn to_dot_file<P: AsRef<std::path::Path>>(&self, path: &P) {
@@ -126,15 +127,22 @@ where
         &self,
         input_domain: &mut SliceDomain<'a, GRID_DIMENSION>,
         output_domain: &mut SliceDomain<'a, GRID_DIMENSION>,
+        mut global_time: usize,
     ) {
         let repeat_solve = self.plan.unwrap_repeat_node(self.plan.root);
+        let repeat_periodic_solve =
+            self.plan.unwrap_periodic_node(repeat_solve.node);
+        let repeat_steps = repeat_periodic_solve.steps;
+
         for _ in 0..repeat_solve.n {
             self.periodic_solve_preallocated_io(
                 repeat_solve.node,
                 false,
                 input_domain,
                 output_domain,
+                global_time,
             );
+            global_time += repeat_steps;
             std::mem::swap(input_domain, output_domain);
         }
         if let Some(next) = repeat_solve.next {
@@ -143,6 +151,7 @@ where
                 false,
                 input_domain,
                 output_domain,
+                global_time,
             )
         } else {
             std::mem::swap(input_domain, output_domain);
@@ -154,13 +163,24 @@ where
         node_id: NodeId,
         input: &SliceDomain<'b, GRID_DIMENSION>,
         output: &mut SliceDomain<'b, GRID_DIMENSION>,
+        global_time: usize,
     ) {
         match self.plan.get_node(node_id) {
             PlanNode::DirectSolve(_) => {
-                self.direct_solve_allocate_io(node_id, input, output);
+                self.direct_solve_allocate_io(
+                    node_id,
+                    input,
+                    output,
+                    global_time,
+                );
             }
             PlanNode::PeriodicSolve(_) => {
-                self.periodic_solve_allocate_io(node_id, input, output);
+                self.periodic_solve_allocate_io(
+                    node_id,
+                    input,
+                    output,
+                    global_time,
+                );
             }
             PlanNode::Repeat(_) => {
                 panic!("ERROR: Not expecting repeat node");
@@ -173,14 +193,24 @@ where
         node_id: NodeId,
         input: &mut SliceDomain<'b, GRID_DIMENSION>,
         output: &mut SliceDomain<'b, GRID_DIMENSION>,
+        global_time: usize,
     ) {
         match self.plan.get_node(node_id) {
             PlanNode::DirectSolve(_) => {
-                self.direct_solve_preallocated_io(node_id, input, output);
+                self.direct_solve_preallocated_io(
+                    node_id,
+                    input,
+                    output,
+                    global_time,
+                );
             }
             PlanNode::PeriodicSolve(_) => {
                 self.periodic_solve_preallocated_io(
-                    node_id, true, input, output,
+                    node_id,
+                    true,
+                    input,
+                    output,
+                    global_time,
                 );
             }
             PlanNode::Repeat(_) => {
@@ -195,6 +225,7 @@ where
         resize: bool,
         input_domain: &mut SliceDomain<'b, GRID_DIMENSION>,
         output_domain: &mut SliceDomain<'b, GRID_DIMENSION>,
+        mut global_time: usize,
     ) {
         let periodic_solve = self.plan.unwrap_periodic_node(node_id);
         std::mem::swap(input_domain, output_domain);
@@ -234,6 +265,7 @@ where
                             node_id,
                             input_domain_const,
                             &mut node_output,
+                            global_time,
                         );
                     });
                 }
@@ -249,11 +281,13 @@ where
 
         // call time cut if needed
         if let Some(next_id) = periodic_solve.time_cut {
+            global_time += periodic_solve.steps;
             std::mem::swap(input_domain, output_domain);
             self.unknown_solve_preallocated_io(
                 next_id,
                 input_domain,
                 output_domain,
+                global_time,
             );
         }
     }
@@ -263,6 +297,7 @@ where
         node_id: NodeId,
         input: &SliceDomain<'b, GRID_DIMENSION>,
         output: &mut SliceDomain<'b, GRID_DIMENSION>,
+        global_time: usize,
     ) {
         let periodic_solve = self.plan.unwrap_periodic_node(node_id);
 
@@ -277,6 +312,7 @@ where
             true,
             &mut input_domain,
             &mut output_domain,
+            global_time,
         );
 
         // copy output to output
@@ -288,6 +324,7 @@ where
         node_id: NodeId,
         input: &SliceDomain<'b, GRID_DIMENSION>,
         output: &mut SliceDomain<'b, GRID_DIMENSION>,
+        global_time: usize,
     ) {
         let direct_solve = self.plan.unwrap_direct_node(node_id);
 
@@ -301,6 +338,7 @@ where
             node_id,
             &mut input_domain,
             &mut output_domain,
+            global_time,
         );
         debug_assert_eq!(*output_domain.aabb(), direct_solve.output_aabb);
 
@@ -313,6 +351,7 @@ where
         node_id: NodeId,
         input_domain: &mut SliceDomain<'b, GRID_DIMENSION>,
         output_domain: &mut SliceDomain<'b, GRID_DIMENSION>,
+        global_time: usize,
     ) {
         let direct_solve = self.plan.unwrap_direct_node(node_id);
 
@@ -336,6 +375,7 @@ where
             output_domain,
             &direct_solve.sloped_sides,
             direct_solve.steps,
+            global_time,
         );
 
         debug_assert_eq!(
diff --git a/src/fft_solver/direct_frustrum_solver.rs b/src/fft_solver/direct_frustrum_solver.rs
index 14e6d36..9d14fbc 100644
--- a/src/fft_solver/direct_frustrum_solver.rs
+++ b/src/fft_solver/direct_frustrum_solver.rs
@@ -37,6 +37,7 @@ where
         output_domain: &mut SliceDomain<'b, GRID_DIMENSION>,
         sloped_sides: &Bounds<GRID_DIMENSION>,
         steps: usize,
+        mut global_time: usize,
     ) {
         assert_eq!(input_domain.aabb(), output_domain.aabb());
 
@@ -47,6 +48,7 @@ where
 
         let mut output_box = *input_domain.aabb();
         for _ in 0..steps {
+            global_time += 1;
             output_box = output_box.add_bounds_diff(trapezoid_slopes);
             debug_assert!(
                 input_domain.aabb().buffer_size() >= output_box.buffer_size()
@@ -57,6 +59,7 @@ where
                 self.stencil,
                 input_domain,
                 output_domain,
+                global_time,
                 self.chunk_size,
             );
             std::mem::swap(input_domain, output_domain);
diff --git a/src/par_stencil.rs b/src/par_stencil.rs
index 29ae3a4..5b31c90 100644
--- a/src/par_stencil.rs
+++ b/src/par_stencil.rs
@@ -14,6 +14,7 @@ pub fn apply<
     stencil: &StencilF64<Operation, GRID_DIMENSION, NEIGHBORHOOD_SIZE>,
     input: &DomainType,
     output: &mut DomainType,
+    global_time: usize,
     chunk_size: usize,
 ) where
     Operation: StencilOperation<f64, NEIGHBORHOOD_SIZE>,
@@ -27,7 +28,13 @@ pub fn apply<
                     Coord<GRID_DIMENSION>,
                     &mut f64,
                 )| {
-                    let args = gather_args(stencil, bc, input, &world_coord);
+                    let args = gather_args(
+                        stencil,
+                        bc,
+                        input,
+                        &world_coord,
+                        global_time,
+                    );
                     let result = stencil.apply(&args);
                     *value_mut = result;
                 },
@@ -53,7 +60,14 @@ mod unit_test {
             input_domain.par_set_values(|_| 1.0, chunk_size);
 
             let bc = PeriodicCheck::new(&input_domain);
-            apply(&bc, &stencil, &input_domain, &mut output_domain, chunk_size);
+            apply(
+                &bc,
+                &stencil,
+                &input_domain,
+                &mut output_domain,
+                0,
+                chunk_size,
+            );
             for x in output_domain.buffer() {
                 assert_approx_eq!(f64, *x, 1.0);
             }
@@ -65,7 +79,7 @@ mod unit_test {
             input_domain.par_set_values(|_| 2.0, chunk_size);
 
             let bc = PeriodicCheck::new(&input_domain);
-            apply(&bc, &stencil, &input_domain, &mut output_domain, 1);
+            apply(&bc, &stencil, &input_domain, &mut output_domain, 0, 1);
             for x in output_domain.buffer() {
                 assert_approx_eq!(f64, *x, 2.0);
             }
@@ -77,7 +91,7 @@ mod unit_test {
         bounds: AABB<1>,
     }
     impl BCCheck<1> for ErrorCheck {
-        fn check(&self, c: &Coord<1>) -> Option<f64> {
+        fn check(&self, c: &Coord<1>, _global_time: usize) -> Option<f64> {
             assert!(self.bounds.contains(c));
             None
         }
@@ -105,7 +119,7 @@ mod unit_test {
             bounds: input_bound,
         };
 
-        apply(&bc, &stencil, &input_domain, &mut output_domain, 2);
+        apply(&bc, &stencil, &input_domain, &mut output_domain, 0, 2);
 
         for i in output_domain.buffer() {
             assert_approx_eq!(f64, *i, 1.0);
diff --git a/src/solver/direct.rs b/src/solver/direct.rs
index c6e4ccc..552fe1a 100644
--- a/src/solver/direct.rs
+++ b/src/solver/direct.rs
@@ -14,6 +14,7 @@ pub fn box_apply<
     input: &mut DomainType,
     output: &mut DomainType,
     steps: usize,
+    mut global_time: usize,
     chunk_size: usize,
 ) where
     Operation: StencilOperation<f64, NEIGHBORHOOD_SIZE>,
@@ -21,10 +22,12 @@ pub fn box_apply<
 {
     debug_assert_eq!(input.aabb(), output.aabb());
     for _ in 0..steps - 1 {
-        par_stencil::apply(bc, stencil, input, output, chunk_size);
+        global_time += 1;
+        par_stencil::apply(bc, stencil, input, output, global_time, chunk_size);
         std::mem::swap(input, output);
     }
-    par_stencil::apply(bc, stencil, input, output, chunk_size);
+    global_time += 1;
+    par_stencil::apply(bc, stencil, input, output, global_time, chunk_size);
 }
 
 #[cfg(test)]
@@ -63,6 +66,7 @@ mod unit_tests {
             &mut input_domain,
             &mut output_domain,
             steps,
+            0,
             chunk_size,
         );
 
@@ -177,6 +181,7 @@ mod unit_tests {
             &mut input_domain,
             &mut output_domain,
             steps,
+            0,
             chunk_size,
         );
         for i in 0..3 {
diff --git a/src/solver/periodic_direct.rs b/src/solver/periodic_direct.rs
index a109ba7..2b50f04 100644
--- a/src/solver/periodic_direct.rs
+++ b/src/solver/periodic_direct.rs
@@ -2,6 +2,10 @@ use crate::domain::*;
 use crate::par_stencil;
 use crate::stencil::*;
 
+/// Global time doesn't matter for periodic solves
+/// since its only used for boundary conditions
+const GLOBAL_TIME: usize = 0;
+
 pub fn direct_periodic_apply<
     Operation,
     const GRID_DIMENSION: usize,
@@ -20,13 +24,20 @@ pub fn direct_periodic_apply<
     for _ in 0..steps - 1 {
         {
             let bc = PeriodicCheck::new(input);
-            par_stencil::apply(&bc, stencil, input, output, chunk_size);
+            par_stencil::apply(
+                &bc,
+                stencil,
+                input,
+                output,
+                GLOBAL_TIME,
+                chunk_size,
+            );
         }
         std::mem::swap(input, output);
     }
 
     let bc = PeriodicCheck::new(input);
-    par_stencil::apply(&bc, stencil, input, output, chunk_size);
+    par_stencil::apply(&bc, stencil, input, output, GLOBAL_TIME, chunk_size);
 }
 
 #[cfg(test)]
diff --git a/tests/ap_solver_compare.rs b/tests/ap_solver_compare.rs
index 904fce2..7f69211 100644
--- a/tests/ap_solver_compare.rs
+++ b/tests/ap_solver_compare.rs
@@ -41,7 +41,7 @@ fn heat_1d_ap_compare() {
     };
     let fft_solver =
         APSolver::new(&bc, &stencil, grid_bound, n_steps, &planner_params);
-    fft_solver.apply(&mut fft_input_domain, &mut fft_output_domain);
+    fft_solver.apply(&mut fft_input_domain, &mut fft_output_domain, 0);
 
     box_apply(
         &bc,
@@ -49,6 +49,7 @@ fn heat_1d_ap_compare() {
         &mut direct_input_domain,
         &mut direct_output_domain,
         n_steps,
+        0,
         chunk_size,
     );
 
@@ -98,7 +99,7 @@ fn heat_2d_ap_compare() {
     };
     let fft_solver =
         APSolver::new(&bc, &stencil, grid_bound, n_steps, &planner_params);
-    fft_solver.apply(&mut fft_input_domain, &mut fft_output_domain);
+    fft_solver.apply(&mut fft_input_domain, &mut fft_output_domain, 0);
 
     box_apply(
         &bc,
@@ -106,6 +107,7 @@ fn heat_2d_ap_compare() {
         &mut direct_input_domain,
         &mut direct_output_domain,
         n_steps,
+        0,
         chunk_size,
     );