diff --git a/pgens/piston_shock/pgen.hpp b/pgens/piston_shock/pgen.hpp
index 745f751d..22f373b4 100644
--- a/pgens/piston_shock/pgen.hpp
+++ b/pgens/piston_shock/pgen.hpp
@@ -4,25 +4,72 @@
 #include "enums.h"
 #include "global.h"
 
-#include "arch/traits.h"
+#include "utils/error.h"
+#include "utils/numeric.h"
 
-#include "archetypes/energy_dist.h"
-#include "archetypes/particle_injector.h"
-#include "archetypes/piston.h"
+#include "archetypes/field_setter.h"
 #include "archetypes/problem_generator.h"
-#include "archetypes/spatial_dist.h"
+#include "archetypes/traits.h"
 #include "archetypes/utils.h"
-#include "framework/domain/domain.h"
+#include "archetypes/piston.h"
+#include "archetypes/moving_window.h"
 #include "framework/domain/metadomain.h"
 
-#include <vector>
+#include <algorithm>
+#include <utility>
 
 namespace user {
   using namespace ntt;
 
+  template <Dimension D>
+  struct InitFields {
+    /*
+      Sets up magnetic and electric field components for the simulation.
+      Must satisfy E = -v x B for Lorentz Force to be zero.
+
+      @param bmag: magnetic field scaling
+      @param btheta: magnetic field polar angle
+      @param bphi: magnetic field azimuthal angle
+      @param drift_ux: drift velocity in the x direction
+    */
+    InitFields(real_t bmag, real_t btheta, real_t bphi, real_t drift_ux)
+      : Bmag { bmag }
+      , Btheta { btheta * static_cast<real_t>(convert::deg2rad) }
+      , Bphi { bphi * static_cast<real_t>(convert::deg2rad) }
+      , Vx { drift_ux } {}
+
+    // magnetic field components
+    Inline auto bx1(const coord_t<D>&) const -> real_t {
+      return Bmag * math::cos(Btheta);
+    }
+
+    Inline auto bx2(const coord_t<D>&) const -> real_t {
+      return Bmag * math::sin(Btheta) * math::sin(Bphi);
+    }
+
+    Inline auto bx3(const coord_t<D>&) const -> real_t {
+      return Bmag * math::sin(Btheta) * math::cos(Bphi);
+    }
+
+    // electric field components
+    Inline auto ex1(const coord_t<D>&) const -> real_t {
+      return ZERO;
+    }
+
+    Inline auto ex2(const coord_t<D>&) const -> real_t {
+      return -Vx * Bmag * math::sin(Btheta) * math::cos(Bphi);
+    }
+
+    Inline auto ex3(const coord_t<D>&) const -> real_t {
+      return Vx * Bmag * math::sin(Btheta) * math::sin(Bphi);
+    }
+
+  private:
+    const real_t Btheta, Bphi, Vx, Bmag;
+  };
+
   template <SimEngine::type S, class M>
   struct PGen : public arch::ProblemGenerator<S, M> {
-
     // compatibility traits for the problem generator
     static constexpr auto engines {
       arch::traits::pgen::compatible_with<SimEngine::SRPIC>::value
@@ -39,46 +86,59 @@ namespace user {
     using arch::ProblemGenerator<S, M>::C;
     using arch::ProblemGenerator<S, M>::params;
 
-    const Metadomain<S, M>& global_domain;
+    Metadomain<S, M>& global_domain;
 
     // domain properties
-    const real_t global_xmin, global_xmax;
+    const real_t  global_xmin, global_xmax;
+    const real_t  cell_size;
+    // gas properties
+    const real_t  temperature, temperature_ratio;
+    // injector properties
+    const real_t  dt;
+    // magnetic field properties
+    real_t        Btheta, Bphi, Bmag;
+    InitFields<D> init_flds;
 
-    // plasma
-    const real_t temperature, temperature_ratio;
     // piston properties
-    const real_t piston_velocity, piston_initial_position;
+    const real_t _piston_velocity;
+    int i_piston;
+    real_t di_piston, piston_position, piston_velocity;
 
-    inline PGen(const SimulationParams& p, const Metadomain<S, M>& global_domain)
+    // window properties
+    const int window_update_frequency;
+
+    inline PGen(const SimulationParams& p, Metadomain<S, M>& global_domain)
       : arch::ProblemGenerator<S, M> { p }
       , global_domain { global_domain }
       , global_xmin { global_domain.mesh().extent(in::x1).first }
       , global_xmax { global_domain.mesh().extent(in::x1).second }
-      , piston_velocity { p.template get<real_t>("setup.piston_velocity", 0.0) }
-      , piston_initial_position { p.template get<real_t>(
-          "setup.piston_initial_position",
-          0.0) }
-      , temperature { p.template get<real_t>("setup.temperature", 0.0) }
-      , temperature_ratio { p.template get<real_t>(
-          "setup.temperature_ratio") } {}
-
-    inline void InitPrtls(Domain<S, M>& local_domain) {
-      real_t xg_min = global_xmin + piston_initial_position;
-      real_t xg_max = global_xmax;
-
-      // define box to inject into
-      boundaries_t<real_t> box;
-      // loop over all dimensions
-      for (auto d { 0u }; d < (unsigned int)M::Dim; ++d) {
-        // compute the range for the x-direction
-        if (d == static_cast<decltype(d)>(in::x1)) {
-          box.push_back({ xg_min, xg_max });
-        } else {
-          // inject into full range in other directions
-          box.push_back(Range::All);
-        }
-      }
+      , cell_size { (global_xmax - global_xmin) / global_domain.mesh().n_all(in::x1) }
+      , temperature { p.template get<real_t>("setup.temperature") }
+      , temperature_ratio { p.template get<real_t>("setup.temperature_ratio") }
+      , Bmag { p.template get<real_t>("setup.Bmag", ZERO) }
+      , Btheta { p.template get<real_t>("setup.Btheta", ZERO) }
+      , Bphi { p.template get<real_t>("setup.Bphi", ZERO) }
+      , init_flds { Bmag, Btheta, Bphi, ZERO }
+      , dt { p.template get<real_t>("algorithms.timestep.dt") }
+      , window_update_frequency { p.template get<int>("setup.window_update_frequency", N_GHOSTS) }
+      , _piston_velocity { p.template get<real_t>("setup.piston_velocity", ZERO)} {}
+
+    inline PGen() {}
+
+    auto MatchFields(simtime_t) const -> InitFields<D> {
+      return init_flds;
+    }
 
+    inline void InitPrtls(Domain<S, M>& domain) {
+
+      // set initial position of piston
+      i_piston = 0;
+      di_piston = ZERO;
+      piston_position = ZERO;
+      
+      coord_t<M::PrtlDim> xp_Cd { ZERO };
+      // piston velocity in code units
+      piston_velocity = domain.mesh.metric.template transform<1, Idx::XYZ, Idx::U>(xp_Cd, _piston_velocity);
       // define temperatures of species
       const auto temperatures = std::make_pair(temperature,
                                                temperature_ratio * temperature);
@@ -88,59 +148,101 @@ namespace user {
 
       // inject particles
       arch::InjectUniformMaxwellians<S, M>(params,
-                                           local_domain,
+                                           domain,
                                            ONE,
                                            temperatures,
                                            { 1, 2 },
                                            drifts,
-                                           false,
-                                           box);
+                                           false);
     }
 
+    void CustomPostStep(timestep_t step, simtime_t time, Domain<S, M>& domain) {
+
+      /* 
+        update piston position
+      */
+      // piston movement over timestep
+      di_piston += piston_velocity * dt;
+      // check if the piston has moved to the next cell
+      i_piston += static_cast<int>(di_piston >= ONE);
+      // keep track of how much the piston has moved into the next cell
+      di_piston -= (di_piston >= ONE);
+
+      // check if the window should be moved
+      if ((i_piston >= window_update_frequency) && (window_update_frequency > 0)) {
+        // move the window and all fields and particles in it
+        arch::MoveWindow<S, M, in::x1>(domain, window_update_frequency);
+        // synch ghost zones after moving the window
+        global_domain.CommunicateFields(domain, Comm::E | Comm::B);
+        // communicate particles after moving
+        global_domain.CommunicateParticles(domain);
+
+        /*
+            Inject slab of fresh plasma
+        */
+        const real_t xmax = global_xmax;
+        const real_t xmin = xmax - window_update_frequency * cell_size;
+        // define box to inject into
+        boundaries_t<real_t> inj_box;
+        // loop over all dimension
+        for (auto d = 0u; d < M::Dim; ++d) {
+          if (d == 0) {
+            inj_box.push_back({ xmin, xmax });
+          } else {
+            inj_box.push_back(Range::All);
+          }
+        }
+
+        // same maxwell distribution as above
+        const auto temperatures = std::make_pair(temperature,
+                                                temperature_ratio * temperature);
+        const auto drifts       = std::make_pair(
+          std::vector<real_t> { ZERO, ZERO, ZERO },
+          std::vector<real_t> { ZERO, ZERO, ZERO });
+        arch::InjectUniformMaxwellians<S, M>(params,
+                                            domain,
+                                            ONE,
+                                            temperatures,
+                                            { 1, 2 },
+                                            drifts,
+                                            false,
+                                            inj_box);
+
+        i_piston -= window_update_frequency;
+      }
+
+      // compute current position of piston
+      piston_position = static_cast<real_t>(i_piston) + static_cast<real_t>(di_piston);
+    }
+    
+
     struct CustomPrtlUpdate {
-      real_t piston_position_current; // current position of piston
-      real_t piston_velocity_current;
+      real_t x_piston;
+      real_t v_piston;
       real_t xg_max;
-
       bool is_left;
-
       bool massive;
-
+      
       template <class Coord, class PusherKernel>
       Inline void operator()(index_t p, Coord& xp, PusherKernel& pusher) const {
 
         real_t piston_position_use;
 
-        // make sure piston has not reached the right wall
-        if (piston_position_current < xg_max) {
-          piston_position_use = piston_position_current;
+        if (x_piston < xg_max){ //make sure piston has not reached the right wall
+            piston_position_use = x_piston;
         } else {
           piston_position_use = xg_max;
         }
-
-        if (arch::CrossesPiston<S, M, PusherKernel>(p,
-                                                    pusher,
-                                                    piston_position_use,
-                                                    piston_velocity_current,
-                                                    is_left)) {
-          arch::PistonUpdate<S, M, PusherKernel>(p,
-                                                 pusher,
-                                                 piston_position_use,
-                                                 piston_velocity_current,
-                                                 massive);
+        if (arch::CrossesPiston<S, M, PusherKernel>(p, pusher, piston_position_use, v_piston, is_left)){
+            arch::PistonUpdate<S, M, PusherKernel>(p, pusher, piston_position_use, v_piston, massive);
         }
       }
     };
 
     template <class D>
     auto CustomParticleUpdate(simtime_t time, spidx_t sp, D& domain) const {
-      return CustomPrtlUpdate { piston_initial_position +
-                                  static_cast<real_t>(time) * piston_velocity,
-                                piston_velocity,
-                                global_domain.mesh().extent(in::x1).second,
-                                true,
-                                true };
-    };
+      return CustomPrtlUpdate { piston_position, piston_velocity, global_xmax, true, true};
+      };
   };
 
 } // namespace user
diff --git a/src/archetypes/moving_window.h b/src/archetypes/moving_window.h
new file mode 100644
index 00000000..49865b26
--- /dev/null
+++ b/src/archetypes/moving_window.h
@@ -0,0 +1,288 @@
+/**
+ * @file archetypes/moving_window.h
+ * @brief Moving window functions for implementing the moving window in the CustomPostStep
+ * @implements
+ *   - arch::MoveWindow<> -> void
+ * @namespaces:
+ *   - arch::
+ */
+
+#ifndef ARCHETYPES_WINDOW_H
+#define ARCHETYPES_WINDOW_H
+
+#include "enums.h"
+#include "global.h"
+
+#include "archetypes/energy_dist.h"
+#include "framework/domain/domain.h"
+#include "framework/parameters/parameters.h"
+#include "framework/domain/metadomain.h"
+
+#include <utility>
+
+
+namespace arch {
+
+  template <Dimension D, in o>
+  struct FieldShift_kernel {
+    static_assert(static_cast<dim_t>(o) < static_cast<dim_t>(D),
+                  "Invalid component index");
+
+    ndfield_t<D, 6>   Fld;
+    ndfield_t<D, 6>   backup_Fld;
+    const BCTags      tags;
+    const int         window_shift;
+
+    FieldShift_kernel(ndfield_t<D, 6> Fld, ndfield_t<D, 6> backup_Fld, const int window_shift, BCTags tags)
+      : Fld { Fld }
+      , backup_Fld { backup_Fld }
+      , window_shift { window_shift }
+      , tags { tags } {}
+
+    Inline void operator()(index_t i1) const {
+      if constexpr (D == Dim::_1D) {
+        if (tags & BC::E) {
+          Fld(i1, em::ex1) = backup_Fld(i1 + window_shift, em::ex1);
+          Fld(i1, em::ex2) = backup_Fld(i1 + window_shift, em::ex2);
+          Fld(i1, em::ex3) = backup_Fld(i1 + window_shift, em::ex3);
+        }
+        if (tags & BC::B) {
+          Fld(i1, em::bx1) = backup_Fld(i1 + window_shift, em::bx1);
+          Fld(i1, em::bx2) = backup_Fld(i1 + window_shift, em::bx2);
+          Fld(i1, em::bx3) = backup_Fld(i1 + window_shift, em::bx3);
+        }
+      } else {
+        raise::KernelError(
+          HERE,
+          "FieldShift_kernel: 1D implementation called for D != 1");
+      }
+    }
+
+    Inline void operator()(index_t i1, index_t i2) const {
+      if constexpr (D == Dim::_2D) {
+        if constexpr (o == in::x1) {
+          if (tags & BC::E) {
+            Fld(i1, i2, em::ex1) = backup_Fld(i1 + window_shift, i2, em::ex1);
+            Fld(i1, i2, em::ex2) = backup_Fld(i1 + window_shift, i2, em::ex2);
+            Fld(i1, i2, em::ex3) = backup_Fld(i1 + window_shift, i2, em::ex3);
+          }
+          if (tags & BC::B) {
+            Fld(i1, i2, em::bx1) = backup_Fld(i1 + window_shift, i2, em::bx1);
+            Fld(i1, i2, em::bx2) = backup_Fld(i1 + window_shift, i2, em::bx2);
+            Fld(i1, i2, em::bx3) = backup_Fld(i1 + window_shift, i2, em::bx3);
+          }
+        } else if constexpr (o == in::x2) {
+          if (tags & BC::E) {
+            Fld(i1, i2, em::ex1) = backup_Fld(i1, i2 + window_shift, em::ex1);
+            Fld(i1, i2, em::ex2) = backup_Fld(i1, i2 + window_shift, em::ex2);
+            Fld(i1, i2, em::ex3) = backup_Fld(i1, i2 + window_shift, em::ex3);
+          }
+          if (tags & BC::B) {
+            Fld(i1, i2, em::bx1) = backup_Fld(i1, i2 + window_shift, em::bx1);
+            Fld(i1, i2, em::bx2) = backup_Fld(i1, i2 + window_shift, em::bx2);
+            Fld(i1, i2, em::bx3) = backup_Fld(i1, i2 + window_shift, em::bx3);
+          }
+        }
+      } else {
+        raise::KernelError(
+          HERE,
+          "FieldShift_kernel: 2D implementation called for D != 2");
+      }
+    }
+
+    Inline void operator()(index_t i1, index_t i2, index_t i3) const {
+      if constexpr (D == Dim::_3D) {
+        if constexpr (o == in::x1) {
+          if (tags & BC::E) {
+            Fld(i1, i2, i3, em::ex1) = backup_Fld(i1 + window_shift, i2, i3, em::ex1);
+            Fld(i1, i2, i3, em::ex2) = backup_Fld(i1 + window_shift, i2, i3, em::ex2);
+            Fld(i1, i2, i3, em::ex3) = backup_Fld(i1 + window_shift, i2, i3, em::ex3);
+          }
+          if (tags & BC::B) {
+            Fld(i1, i2, i3, em::bx1) = backup_Fld(i1 + window_shift, i2, i3, em::bx1);
+            Fld(i1, i2, i3, em::bx2) = backup_Fld(i1 + window_shift, i2, i3, em::bx2);
+            Fld(i1, i2, i3, em::bx3) = backup_Fld(i1 + window_shift, i2, i3, em::bx3);
+          }
+        } else if constexpr (o == in::x2) {
+          if (tags & BC::E) {
+            Fld(i1, i2, i3, em::ex1) = backup_Fld(i1, i2 + window_shift, i3, em::ex1);
+            Fld(i1, i2, i3, em::ex2) = backup_Fld(i1, i2 + window_shift, i3, em::ex2);
+            Fld(i1, i2, i3, em::ex3) = backup_Fld(i1, i2 + window_shift, i3, em::ex3);
+          }
+          if (tags & BC::B) {
+            Fld(i1, i2, i3, em::bx1) = backup_Fld(i1, i2 + window_shift, i3, em::bx1);
+            Fld(i1, i2, i3, em::bx2) = backup_Fld(i1, i2 + window_shift, i3, em::bx2);
+            Fld(i1, i2, i3, em::bx3) = backup_Fld(i1, i2 + window_shift, i3, em::bx3);
+          }
+        } else {
+          if (tags & BC::E) {
+            Fld(i1, i2, i3, em::ex1) = backup_Fld(i1, i2, i3 + window_shift, em::ex1);
+            Fld(i1, i2, i3, em::ex2) = backup_Fld(i1, i2, i3 + window_shift, em::ex2);
+            Fld(i1, i2, i3, em::ex3) = backup_Fld(i1, i2, i3 + window_shift, em::ex3);
+          }
+          if (tags & BC::B) {
+            Fld(i1, i2, i3, em::bx1) = backup_Fld(i1, i2, i3 + window_shift, em::bx1);
+            Fld(i1, i2, i3, em::bx2) = backup_Fld(i1, i2, i3 + window_shift, em::bx2);
+            Fld(i1, i2, i3, em::bx3) = backup_Fld(i1, i2, i3 + window_shift, em::bx3);
+          }
+        }
+      } else {
+        raise::KernelError(
+          HERE,
+          "FieldShift_kernel: 3D implementation called for D != 3");
+      }
+    }
+  };
+
+  /**
+   * @brief Updates particle position and fields in the moving window.
+
+   */
+  template <SimEngine::type S, class M, in o>
+  Inline void MoveWindow(Domain<S, M>& domain, 
+                         const int window_shift) {
+    
+    if constexpr (M::CoordType != Coord::Cart) {
+      raise::Error(
+        "Moving window only applicable to cartesian coordinates",
+        HERE);
+    } else {
+
+      /*
+        move particles in the window back by the window size
+      */
+      const auto nspec = domain.species.size();
+      const auto& mesh = domain.mesh;
+      if (o == in::x1) {
+        // loop over particle species
+        for (auto s { 0u }; s < nspec; ++s) {
+          // get particle properties
+          auto& species = domain.species[s];
+          auto  i1      = species.i1;
+          Kokkos::parallel_for(
+            "MoveParticles",
+            species.rangeActiveParticles(),
+            Lambda(index_t p) {
+              // shift particle position back by window update frequency
+              i1(p) -= window_shift;
+            }
+          );
+        }
+      } else if (o == in::x2) {
+        if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) {
+          // loop over particle species
+          for (auto s { 0u }; s < nspec; ++s) {
+            // get particle properties
+            auto& species = domain.species[s];
+            auto  i2      = species.i2;
+            Kokkos::parallel_for(
+              "MoveParticles",
+              species.rangeActiveParticles(),
+              Lambda(index_t p) {
+                // shift particle position back by window update frequency
+                i2(p) -= window_shift;
+              }
+            );
+          }
+        } else {
+          raise::Error("Invalid dimension", HERE);
+        }
+      } else if (o == in::x3) {
+        if constexpr (M::Dim == Dim::_3D) {
+          // loop over particle species
+          for (auto s { 0u }; s < nspec; ++s) {
+            // get particle properties
+            auto& species = domain.species[s];
+            auto  i3      = species.i3;
+            Kokkos::parallel_for(
+              "MoveParticles",
+              species.rangeActiveParticles(),
+              Lambda(index_t p) {
+                // shift particle position back by window update frequency
+                i3(p) -= window_shift;
+              }
+            );
+          }
+        } else {
+          raise::Error("Invalid direction", HERE);
+        }
+      } else {
+        raise::Error("Invalid direction", HERE);
+      }
+
+      // shift fields in the window back by the window size
+      std::vector<std::size_t> xi_min, xi_max;
+      const std::vector<in> all_dirs { in::x1, in::x2, in::x3 };
+      for (auto d { 0u }; d < M::Dim; ++d) {
+        const auto dd = all_dirs[d];
+        if (o == dd) {
+          xi_min.push_back(0);
+          xi_max.push_back(domain.mesh.n_all(dd) - window_shift);
+        } else {
+          xi_min.push_back(0);
+          xi_max.push_back(domain.mesh.n_all(dd));
+        }
+      }
+      raise::ErrorIf(xi_min.size() != xi_max.size() or
+                       xi_min.size() != static_cast<std::size_t>(M::Dim),
+                     "Invalid range size",
+                     HERE);
+
+      // loop range for shifting fields
+      range_t<M::Dim> range;
+      if constexpr (M::Dim == Dim::_1D) {
+        range = CreateRangePolicy<M::Dim>({ xi_min[0] }, { xi_max[0] });
+      } else if constexpr (M::Dim == Dim::_2D) {
+        range = CreateRangePolicy<M::Dim>({ xi_min[0], xi_min[1] },
+                                          { xi_max[0], xi_max[1] });
+      } else if constexpr (M::Dim == Dim::_3D) {
+        range = CreateRangePolicy<M::Dim>({ xi_min[0], xi_min[1], xi_min[2] },
+                                          { xi_max[0], xi_max[1], xi_max[2] });
+      } else {
+        raise::Error("Invalid dimension", HERE);
+      }
+
+      // copy fields to backup before shifting
+      Kokkos::deep_copy(domain.fields.bckp, domain.fields.em); 
+
+      if (o == in::x1) {
+        Kokkos::parallel_for(
+            "ShiftFields",
+            range,
+            FieldShift_kernel<M::Dim, in::x1>(
+              domain.fields.em,
+              domain.fields.bckp,
+              window_shift,
+              BC::B | BC::E));
+      } else if (o == in::x2) {
+        if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) {
+          Kokkos::parallel_for(
+            "ShiftFields",
+            range,
+            FieldShift_kernel<M::Dim, in::x2>(
+              domain.fields.em,
+              domain.fields.bckp,
+              window_shift,
+              BC::B | BC::E));
+        } else {
+          raise::Error("Invalid dimension", HERE);
+        }
+      } else {
+        if constexpr (M::Dim == Dim::_3D) {
+          Kokkos::parallel_for(
+            "ShiftFields",
+            range,
+            FieldShift_kernel<M::Dim, in::x3>(
+              domain.fields.em,
+              domain.fields.bckp,
+              window_shift,
+              BC::B | BC::E));
+        } else {
+          raise::Error("Invalid dimension", HERE);
+        }
+      }
+    }
+  }
+} // namespace arch
+
+#endif // ARCHETYPES_UTILS_H