diff --git a/src/include/fields.hxx b/src/include/fields.hxx
index 7696df04a..1194c7b7d 100644
--- a/src/include/fields.hxx
+++ b/src/include/fields.hxx
@@ -65,6 +65,16 @@ public:
     return type_traits<F>::get(e_(i, j, k, m));
   }
 
+  GT_INLINE const_reference operator()(int m, const Int3& i3) const
+  {
+    return (*this)(m, i3[0], i3[1], i3[2]);
+  }
+
+  GT_INLINE reference operator()(int m, const Int3& i3)
+  {
+    return (*this)(m, i3[0], i3[1], i3[2]);
+  }
+
 private:
   fields_t e_;
   Int3 ib_;
diff --git a/src/include/psc.hxx b/src/include/psc.hxx
index 93a8feed7..3adc7fd66 100644
--- a/src/include/psc.hxx
+++ b/src/include/psc.hxx
@@ -215,12 +215,12 @@ struct Psc
 
   void pre_first_step()
   {
-    bndf_.fill_ghosts_H(mflds_);
+    bndf.fill_ghosts_H(mflds_);
     bnd_.fill_ghosts(mflds_, HX, HX + 3);
 
     bnd_.fill_ghosts(mflds_, JXI, JXI + 3);
 
-    bndf_.fill_ghosts_E(mflds_);
+    bndf.fill_ghosts_E(mflds_);
     bnd_.fill_ghosts(mflds_, EX, EX + 3);
 
     // FIXME: do a half-step on p to bring it to its natural time,
@@ -399,7 +399,7 @@ struct Psc
 
     mpi_printf(comm, "***** Bnd fields J...\n");
     prof_start(pr_bndf);
-    bndf_.add_ghosts_J(mflds_);
+    bndf.add_ghosts_J(mflds_);
     bnd_.add_ghosts(mflds_, JXI, JXI + 3);
     bnd_.fill_ghosts(mflds_, JXI, JXI + 3);
     prof_stop(pr_bndf);
@@ -413,7 +413,7 @@ struct Psc
 
     mpi_printf(comm, "***** Bnd fields B (1 of 2)...\n");
     prof_restart(pr_bndf);
-    bndf_.fill_ghosts_H(mflds_);
+    bndf.fill_ghosts_H(mflds_);
     bnd_.fill_ghosts(mflds_, HX, HX + 3);
     prof_stop(pr_bndf);
     // state is now: x^{n+3/2}, p^{n+1}, E^{n+1/2}, B^{n+1}, j^{n+1}
@@ -426,7 +426,7 @@ struct Psc
 
     mpi_printf(comm, "***** Bnd fields E...\n");
     prof_restart(pr_bndf);
-    bndf_.fill_ghosts_E(mflds_);
+    bndf.fill_ghosts_E(mflds_);
     bnd_.fill_ghosts(mflds_, EX, EX + 3);
     prof_stop(pr_bndf);
     // state is now: x^{n+3/2}, p^{n+1}, E^{n+3/2}, B^{n+1}
@@ -439,7 +439,7 @@ struct Psc
 
     mpi_printf(comm, "***** Bnd fields B (2 of 2)...\n");
     prof_restart(pr_bndf);
-    bndf_.fill_ghosts_H(mflds_);
+    bndf.fill_ghosts_H(mflds_);
     bnd_.fill_ghosts(mflds_, HX, HX + 3);
     prof_stop(pr_bndf);
     // state is now: x^{n+3/2}, p^{n+1}, E^{n+3/2}, B^{n+3/2}
@@ -517,6 +517,7 @@ private:
 
 public:
   const Grid_t& grid() { return *grid_; }
+  BndFields bndf;
 
 private:
   double time_start_;
@@ -540,7 +541,6 @@ protected:
   PushParticles pushp_;
   PushFields pushf_;
   Bnd bnd_;
-  BndFields bndf_;
   BndParticles bndp_;
 
   Checkpointing checkpointing_;
diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx
index 2e8c18516..afe33fae5 100644
--- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx
+++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx
@@ -1,5 +1,6 @@
 
 #include "psc.h"
+#include "kg/VecRange.hxx"
 #include "fields.hxx"
 #include "bnd_fields.hxx"
 
@@ -7,22 +8,22 @@
 
 #include <limits>
 
-// FIXME, needs public access to Fields::ib, im
 // #define DEBUG
 
-template <typename MF, typename Dim>
+template <typename MFIELDS_STATE, typename Dim>
 struct BndFields_ : BndFieldsBase
 {
-  using Self = BndFields_<MF, Dim>;
-  using Mfields = MF;
-  using real_t = typename Mfields::real_t;
-  using fields_view_t = typename Mfields::fields_view_t;
+  using Self = BndFields_<MFIELDS_STATE, Dim>;
+  using MfieldsState = MFIELDS_STATE;
+  using real_t = typename MfieldsState::real_t;
+  using Real3 = Vec3<real_t>;
+  using fields_view_t = typename MfieldsState::fields_view_t;
   using dim_t = Dim;
 
   // ----------------------------------------------------------------------
   // fill_ghosts_E
 
-  void fill_ghosts_E(Mfields& mflds)
+  void fill_ghosts_E(MfieldsState& mflds)
   {
     const auto& grid = mflds.grid();
 
@@ -31,12 +32,20 @@ struct BndFields_ : BndFieldsBase
       for (int d = 0; d < 3; d++) {
         if (grid.atBoundaryLo(p, d)) {
           switch (grid.bc.fld_lo[d]) {
-            case BND_FLD_PERIODIC: break;
-            case BND_FLD_CONDUCTING_WALL:
+            case BND_FLD_PERIODIC: {
+              break;
+            }
+            case BND_FLD_CONDUCTING_WALL: {
               conducting_wall_E_lo(mflds, p, d);
               break;
-            case BND_FLD_OPEN: break;
-            default: assert(0);
+            }
+            case BND_FLD_OPEN: {
+              set_lower_ghosts(mflds, p, d, EX, background_e, false);
+              break;
+            }
+            default: {
+              assert(0);
+            }
           }
         }
       }
@@ -45,12 +54,20 @@ struct BndFields_ : BndFieldsBase
       for (int d = 0; d < 3; d++) {
         if (grid.atBoundaryHi(p, d)) {
           switch (grid.bc.fld_hi[d]) {
-            case BND_FLD_PERIODIC: break;
-            case BND_FLD_CONDUCTING_WALL:
+            case BND_FLD_PERIODIC: {
+              break;
+            }
+            case BND_FLD_CONDUCTING_WALL: {
               conducting_wall_E_hi(mflds, p, d);
               break;
-            case BND_FLD_OPEN: break;
-            default: assert(0);
+            }
+            case BND_FLD_OPEN: {
+              set_upper_ghosts(mflds, p, d, EX, background_e, false);
+              break;
+            }
+            default: {
+              assert(0);
+            }
           }
         }
       }
@@ -60,7 +77,7 @@ struct BndFields_ : BndFieldsBase
   // ----------------------------------------------------------------------
   // fill_ghosts_H
 
-  void fill_ghosts_H(Mfields& mflds)
+  void fill_ghosts_H(MfieldsState& mflds)
   {
     const auto& grid = mflds.grid();
 
@@ -69,12 +86,20 @@ struct BndFields_ : BndFieldsBase
       for (int d = 0; d < 3; d++) {
         if (grid.atBoundaryLo(p, d)) {
           switch (grid.bc.fld_lo[d]) {
-            case BND_FLD_PERIODIC: break;
-            case BND_FLD_CONDUCTING_WALL:
+            case BND_FLD_PERIODIC: {
+              break;
+            }
+            case BND_FLD_CONDUCTING_WALL: {
               conducting_wall_H_lo(mflds, p, d);
               break;
-            case BND_FLD_OPEN: open_H_lo(mflds, p, d); break;
-            default: assert(0);
+            }
+            case BND_FLD_OPEN: {
+              radiative_H_lo(mflds, p, d);
+              break;
+            }
+            default: {
+              assert(0);
+            }
           }
         }
       }
@@ -82,12 +107,20 @@ struct BndFields_ : BndFieldsBase
       for (int d = 0; d < 3; d++) {
         if (grid.atBoundaryHi(p, d)) {
           switch (grid.bc.fld_hi[d]) {
-            case BND_FLD_PERIODIC: break;
-            case BND_FLD_CONDUCTING_WALL:
+            case BND_FLD_PERIODIC: {
+              break;
+            }
+            case BND_FLD_CONDUCTING_WALL: {
               conducting_wall_H_hi(mflds, p, d);
               break;
-            case BND_FLD_OPEN: open_H_hi(mflds, p, d); break;
-            default: assert(0);
+            }
+            case BND_FLD_OPEN: {
+              radiative_H_hi(mflds, p, d);
+              break;
+            }
+            default: {
+              assert(0);
+            }
           }
         }
       }
@@ -97,7 +130,7 @@ struct BndFields_ : BndFieldsBase
   // ----------------------------------------------------------------------
   // add_ghosts_J
 
-  void add_ghosts_J(Mfields& mflds)
+  void add_ghosts_J(MfieldsState& mflds)
   {
     const auto& grid = mflds.grid();
 
@@ -106,12 +139,19 @@ struct BndFields_ : BndFieldsBase
       for (int d = 0; d < 3; d++) {
         if (grid.atBoundaryLo(p, d)) {
           switch (grid.bc.fld_lo[d]) {
-            case BND_FLD_PERIODIC:
-            case BND_FLD_OPEN: break;
-            case BND_FLD_CONDUCTING_WALL:
+            case BND_FLD_PERIODIC: {
+              break;
+            }
+            case BND_FLD_CONDUCTING_WALL: {
               conducting_wall_J_lo(mflds, p, d);
               break;
-            default: assert(0);
+            }
+            case BND_FLD_OPEN: {
+              break;
+            }
+            default: {
+              assert(0);
+            }
           }
         }
       }
@@ -119,43 +159,145 @@ struct BndFields_ : BndFieldsBase
       for (int d = 0; d < 3; d++) {
         if (grid.atBoundaryHi(p, d)) {
           switch (grid.bc.fld_hi[d]) {
-            case BND_FLD_PERIODIC:
-            case BND_FLD_OPEN: break;
-            case BND_FLD_CONDUCTING_WALL:
+            case BND_FLD_PERIODIC: {
+              break;
+            }
+            case BND_FLD_CONDUCTING_WALL: {
               conducting_wall_J_hi(mflds, p, d);
               break;
-            default: assert(0);
+            }
+            case BND_FLD_OPEN: {
+              break;
+            }
+            default: {
+              assert(0);
+            }
           }
         }
       }
     }
   }
 
-  static void fields_t_set_nan(real_t* f)
+  static void set_lower_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb,
+                                      bool include_edge)
   {
-    *f = std::numeric_limits<real_t>::quiet_NaN();
+#ifndef DEBUG
+    return;
+#endif
+    real_t nan = std::numeric_limits<real_t>::quiet_NaN();
+    set_lower_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge);
+  }
+
+  static void set_upper_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb,
+                                      bool include_edge)
+  {
+#ifndef DEBUG
+    return;
+#endif
+    real_t nan = std::numeric_limits<real_t>::quiet_NaN();
+    set_upper_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge);
+  }
+
+  /**
+   * @brief Set E or B lower ghosts to the given constants (each component has
+   * its own constant).
+   * @param mflds mflds
+   * @param p patch index
+   * @param d which dimension to set the ghosts of
+   * @param mb `EX` or `HX`; note that `mb+1` and `mb+2` are also set
+   * @param val the constants
+   * @param include_edge whether or not values located on exact domain edges
+   * should be considered "ghosts"
+   */
+  static void set_lower_ghosts(MfieldsState& mflds, int p, int d, int mb,
+                               Real3 val, bool include_edge)
+  {
+    auto F = make_Fields3d<dim_t>(mflds[p]);
+    Int3 start = mflds.ib();
+    Int3 stop = mflds.ib() + mflds.im();
+    stop[d] = 0;
+
+    // TODO use gtensor views instead of VecRange
+
+    for (int m = mb; m < mb + 3; m++) {
+      for (Int3 i3 : VecRange(start, stop)) {
+        F(m, i3) = val[m - mb];
+      }
+    }
+
+    if (!include_edge) {
+      return;
+    }
+
+    Int3 edge_start = mflds.ib();
+    Int3 edge_stop = mflds.ib() + mflds.im();
+    edge_start[d] = 0;
+    edge_stop[d] = 1;
+
+    for (int m = mb; m < mb + 3; m++) {
+      bool edge_ec = mb == EX && m - mb != d;
+      bool edge_fc = mb == HX && m - mb == d;
+
+      if (edge_ec || edge_fc) {
+        for (Int3 i3 : VecRange(edge_start, edge_stop)) {
+          F(m, i3) = val[m - mb];
+        }
+      }
+    }
+  }
+
+  /**
+   * @brief Set E or B upper ghosts to the given constants (each component has
+   * its own constant).
+   * @param mflds mflds
+   * @param p patch index
+   * @param d which dimension to set the ghosts of
+   * @param mb `EX` or `HX`; note that `mb+1` and `mb+2` are also set
+   * @param val the constants
+   * @param include_edge whether or not values located on exact domain edges
+   * should be considered "ghosts"
+   */
+  static void set_upper_ghosts(MfieldsState& mflds, int p, int d, int mb,
+                               Real3 val, bool include_edge)
+  {
+    auto F = make_Fields3d<dim_t>(mflds[p]);
+    Int3 start = mflds.ib();
+    Int3 stop = mflds.ib() + mflds.im();
+    start[d] = mflds.grid().ldims[d] + 1;
+
+    // TODO use gtensor views instead of VecRange
+
+    for (int m = mb; m < mb + 3; m++) {
+      for (Int3 i3 : VecRange(start, stop)) {
+        F(m, i3) = val[m - mb];
+      }
+    }
+
+    Int3 edge_start = mflds.ib();
+    Int3 edge_stop = mflds.ib() + mflds.im();
+    edge_start[d] = mflds.grid().ldims[d];
+    edge_stop[d] = mflds.grid().ldims[d] + 1;
+
+    for (int m = mb; m < mb + 3; m++) {
+      bool not_edge_ec = mb == EX && m - mb == d;
+      bool not_edge_fc = mb == HX && m - mb != d;
+
+      if (not_edge_ec || not_edge_fc || include_edge)
+        for (Int3 i3 : VecRange(edge_start, edge_stop)) {
+          F(m, i3) = val[m - mb];
+        }
+    }
   }
 
-  void conducting_wall_E_lo(Mfields& mflds, int p, int d)
+  void conducting_wall_E_lo(MfieldsState& mflds, int p, int d)
   {
+    set_lower_ghosts_to_nan(mflds, p, d, EX, true);
+
     auto F = make_Fields3d<dim_t>(mflds[p]);
     const int* ldims = mflds.grid().ldims;
     Int3 ib = mflds.ib(), im = mflds.im();
 
     if (d == 1) {
-#ifdef DEBUG
-      for (int iz = -2; iz < ldims[2] + 2; iz++) {
-        for (int ix = std::max(-2, F.ib_[0]);
-             ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) {
-          fields_t_set_nan(&F(EX, ix, -1, iz));
-          fields_t_set_nan(&F(EX, ix, -2, iz));
-          fields_t_set_nan(&F(EY, ix, -1, iz));
-          fields_t_set_nan(&F(EY, ix, -2, iz));
-          fields_t_set_nan(&F(EZ, ix, -1, iz));
-          fields_t_set_nan(&F(EZ, ix, -2, iz));
-        }
-      }
-#endif
       for (int iz = -2; iz < ldims[2] + 2; iz++) {
         // FIXME, needs to be for other dir, too, and it's ugly
         for (int ix = std::max(-2, ib[0]);
@@ -170,18 +312,6 @@ struct BndFields_ : BndFieldsBase
         }
       }
     } else if (d == 2) {
-#ifdef DEBUG
-      for (int iy = -2; iy < ldims[1] + 2; iy++) {
-        for (int ix = -2; ix < ldims[0] + 2; ix++) {
-          fields_t_set_nan(&F(EX, ix, iy, -1));
-          fields_t_set_nan(&F(EX, ix, iy, -2));
-          fields_t_set_nan(&F(EY, ix, iy, -1));
-          fields_t_set_nan(&F(EY, ix, iy, -2));
-          fields_t_set_nan(&F(EZ, ix, iy, -1));
-          fields_t_set_nan(&F(EZ, ix, iy, -2));
-        }
-      }
-#endif
       for (int iy = -2; iy < ldims[1] + 2; iy++) {
         for (int ix = std::max(-2, ib[0]);
              ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
@@ -199,27 +329,16 @@ struct BndFields_ : BndFieldsBase
     }
   }
 
-  void conducting_wall_E_hi(Mfields& mflds, int p, int d)
+  void conducting_wall_E_hi(MfieldsState& mflds, int p, int d)
   {
+    set_upper_ghosts_to_nan(mflds, p, d, EX, true);
+
     auto F = make_Fields3d<dim_t>(mflds[p]);
     const int* ldims = mflds.grid().ldims;
     Int3 ib = mflds.ib(), im = mflds.im();
 
     if (d == 1) {
       int my _mrc_unused = ldims[1];
-#ifdef DEBUG
-      for (int iz = -2; iz < ldims[2] + 2; iz++) {
-        for (int ix = std::max(-2, F.ib_[0]);
-             ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) {
-          fields_t_set_nan(&F(EX, ix, my, iz));
-          fields_t_set_nan(&F(EX, ix, my + 1, iz));
-          fields_t_set_nan(&F(EY, ix, my, iz));
-          fields_t_set_nan(&F(EY, ix, my + 1, iz));
-          fields_t_set_nan(&F(EZ, ix, my, iz));
-          fields_t_set_nan(&F(EZ, ix, my + 1, iz));
-        }
-      }
-#endif
       for (int iz = -2; iz < ldims[2] + 2; iz++) {
         for (int ix = std::max(-2, ib[0]);
              ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
@@ -234,18 +353,6 @@ struct BndFields_ : BndFieldsBase
       }
     } else if (d == 2) {
       int mz = ldims[2];
-#ifdef DEBUG
-      for (int iy = -2; iy < ldims[1] + 2; iy++) {
-        for (int ix = -2; ix < ldims[0] + 2; ix++) {
-          fields_t_set_nan(&F(EX, ix, iy, mz));
-          fields_t_set_nan(&F(EX, ix, iy, mz + 1));
-          fields_t_set_nan(&F(EY, ix, iy, mz));
-          fields_t_set_nan(&F(EY, ix, iy, mz + 1));
-          fields_t_set_nan(&F(EZ, ix, iy, mz));
-          fields_t_set_nan(&F(EZ, ix, iy, mz + 1));
-        }
-      }
-#endif
       for (int iy = -2; iy < ldims[1] + 2; iy++) {
         for (int ix = std::max(-2, ib[0]);
              ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
@@ -263,26 +370,15 @@ struct BndFields_ : BndFieldsBase
     }
   }
 
-  void conducting_wall_H_lo(Mfields& mflds, int p, int d)
+  void conducting_wall_H_lo(MfieldsState& mflds, int p, int d)
   {
+    set_lower_ghosts_to_nan(mflds, p, d, HX, false);
+
     auto F = make_Fields3d<dim_t>(mflds[p]);
     const int* ldims = mflds.grid().ldims;
     Int3 ib = mflds.ib(), im = mflds.im();
 
     if (d == 1) {
-#ifdef DEBUG
-      for (int iz = -2; iz < ldims[2] + 2; iz++) {
-        for (int ix = std::max(-2, F.ib_[0]);
-             ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) {
-          fields_t_set_nan(&F(HX, ix, -1, iz));
-          fields_t_set_nan(&F(HX, ix, -2, iz));
-          fields_t_set_nan(&F(HY, ix, -1, iz));
-          fields_t_set_nan(&F(HY, ix, -2, iz));
-          fields_t_set_nan(&F(HZ, ix, -1, iz));
-          fields_t_set_nan(&F(HZ, ix, -2, iz));
-        }
-      }
-#endif
       for (int iz = -1; iz < ldims[2] + 2; iz++) {
         for (int ix = std::max(-2, ib[0]);
              ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
@@ -294,18 +390,6 @@ struct BndFields_ : BndFieldsBase
         }
       }
     } else if (d == 2) {
-#ifdef DEBUG
-      for (int iy = -2; iy < ldims[1] + 2; iy++) {
-        for (int ix = -2; ix < ldims[0] + 2; ix++) {
-          fields_t_set_nan(&F(HX, ix, iy, -1));
-          fields_t_set_nan(&F(HX, ix, iy, -2));
-          fields_t_set_nan(&F(HY, ix, iy, -1));
-          fields_t_set_nan(&F(HY, ix, iy, -2));
-          fields_t_set_nan(&F(HZ, ix, iy, -1));
-          fields_t_set_nan(&F(HZ, ix, iy, -2));
-        }
-      }
-#endif
       for (int iy = -2; iy < ldims[1] + 2; iy++) {
         for (int ix = std::max(-2, ib[0]);
              ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
@@ -321,8 +405,10 @@ struct BndFields_ : BndFieldsBase
     }
   }
 
-  void conducting_wall_H_hi(Mfields& mflds, int p, int d)
+  void conducting_wall_H_hi(MfieldsState& mflds, int p, int d)
   {
+    set_upper_ghosts_to_nan(mflds, p, d, HX, false);
+
     auto F = make_Fields3d<dim_t>(mflds[p]);
 
     const int* ldims = mflds.grid().ldims;
@@ -330,18 +416,6 @@ struct BndFields_ : BndFieldsBase
 
     if (d == 1) {
       int my _mrc_unused = ldims[1];
-#ifdef DEBUG
-      for (int iz = -2; iz < ldims[2] + 2; iz++) {
-        for (int ix = std::max(-2, F.ib_[0]);
-             ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) {
-          fields_t_set_nan(&F(HX, ix, my, iz));
-          fields_t_set_nan(&F(HX, ix, my + 1, iz));
-          fields_t_set_nan(&F(HY, ix, my + 1, iz));
-          fields_t_set_nan(&F(HZ, ix, my, iz));
-          fields_t_set_nan(&F(HZ, ix, my + 1, iz));
-        }
-      }
-#endif
       for (int iz = -2; iz < ldims[2] + 2; iz++) {
         for (int ix = std::max(-2, ib[0]);
              ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
@@ -354,17 +428,6 @@ struct BndFields_ : BndFieldsBase
       }
     } else if (d == 2) {
       int mz = ldims[2];
-#ifdef DEBUG
-      for (int iy = -2; iy < ldims[1] + 2; iy++) {
-        for (int ix = -2; ix < ldims[0] + 2; ix++) {
-          fields_t_set_nan(&F(HX, ix, iy, mz));
-          fields_t_set_nan(&F(HX, ix, iy, mz + 1));
-          fields_t_set_nan(&F(HY, ix, iy, mz));
-          fields_t_set_nan(&F(HY, ix, iy, mz + 1));
-          fields_t_set_nan(&F(HZ, ix, iy, mz + 1));
-        }
-      }
-#endif
       for (int iy = -2; iy < ldims[1] + 2; iy++) {
         for (int ix = std::max(-2, ib[0]);
              ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
@@ -380,7 +443,7 @@ struct BndFields_ : BndFieldsBase
     }
   }
 
-  void conducting_wall_J_lo(Mfields& mflds, int p, int d)
+  void conducting_wall_J_lo(MfieldsState& mflds, int p, int d)
   {
     auto F = make_Fields3d<dim_t>(mflds[p]);
     const int* ldims = mflds.grid().ldims;
@@ -419,7 +482,7 @@ struct BndFields_ : BndFieldsBase
     }
   }
 
-  void conducting_wall_J_hi(Mfields& mflds, int p, int d)
+  void conducting_wall_J_hi(MfieldsState& mflds, int p, int d)
   {
     auto F = make_Fields3d<dim_t>(mflds[p]);
     const int* ldims = mflds.grid().ldims;
@@ -460,181 +523,101 @@ struct BndFields_ : BndFieldsBase
     }
   }
 
-  // ======================================================================
-  // open
-
-  // ----------------------------------------------------------------------
-  // open_H_lo
-
-  void open_H_lo(Mfields& mflds, int p, int d)
+  void radiative_H_lo(MfieldsState& mflds, int p, int d)
   {
+    set_lower_ghosts_to_nan(mflds, p, d, HX, false);
+
     auto F = make_Fields3d<dim_t>(mflds[p]);
     const Grid_t& grid = mflds.grid();
-    Int3 ldims = grid.ldims;
-    Int3 ib = mflds.ib();
-    Int3 im = mflds.im();
     real_t dt = grid.dt;
-    real_t dx = grid.domain.dx[0];
-    real_t dy = grid.domain.dx[1];
-    real_t dz = grid.domain.dx[2];
-
-#if 1
-
-    if (d == 1) {
-#ifdef DEBUG
-      for (int iz = -2; iz < ldims[2] + 2; iz++) {
-        for (int ix = std::max(-2, ib[0]);
-             ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
-          fields_t_set_nan(&F(HX, ix, -1, iz));
-          fields_t_set_nan(&F(HX, ix, -2, iz));
-          fields_t_set_nan(&F(HY, ix, -1, iz));
-          fields_t_set_nan(&F(HY, ix, -2, iz));
-          fields_t_set_nan(&F(HZ, ix, -1, iz));
-          fields_t_set_nan(&F(HZ, ix, -2, iz));
-        }
-      }
-#endif
-
-      for (int iz = -2; iz < ldims[2] + 2; iz++) {
-        for (int ix = std::max(-2, ib[0]);
-             ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
-          F(HX, ix, -1, iz) =
-            (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t) */
-             -2.f * F(EZ, ix, 0, iz)
-             /*- dt/dx * (F(HY, ix,0,iz) - F(HY, ix-1,0,iz)) */
-             - (1.f - dt / dy) * F(HX, ix, 0, iz) + dt * F(JZI, ix, 0, iz)) /
-            (1.f + dt / dy);
-          F(HZ, ix, -1, iz) =
-            (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t) */
-             +2.f * F(EX, ix, 0, iz) -
-             dt / dz * (F(HY, ix, 0, iz) - F(HY, ix, 0, iz - 1)) -
-             (1.f - dt / dy) * F(HZ, ix, 0, iz) + dt * F(JXI, ix, 0, iz)) /
-            (1.f + dt / dy);
-        }
-      }
-    } else if (d == 2) {
-#ifdef DEBUG
-      for (int iy = -2; iy < ldims[1] + 2; iy++) {
-        for (int ix = std::max(-2, ib[0]);
-             ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
-          fields_t_set_nan(&F(HX, ix, iy, -1));
-          fields_t_set_nan(&F(HX, ix, iy, -2));
-          fields_t_set_nan(&F(HY, ix, iy, -1));
-          fields_t_set_nan(&F(HY, ix, iy, -2));
-          fields_t_set_nan(&F(HZ, ix, iy, -1));
-          fields_t_set_nan(&F(HZ, ix, iy, -2));
-        }
-      }
-#endif
-      for (int iy = -2; iy < ldims[1] + 2; iy++) {
-        for (int ix = std::max(-2, ib[0]);
-             ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
-          F(HY, ix, iy, -1) =
-            (/* + 4.f * C_s_pulse_z1(x+0.5*dx,y,z,t) */
-             -2.f * F(EX, ix, iy, 0) -
-             dt / dy * (F(HZ, ix, iy, 0) - F(HZ, ix, iy - 1, 0)) -
-             (1.f - dt / dz) * F(HY, ix, iy, 0) + dt * F(JXI, ix, iy, 0)) /
-            (1.f + dt / dz);
-          F(HX, ix, iy, -1) =
-            (/* - 4.f * C_p_pulse_z1(x+0.5*dx,y,z,t) */
-             +2.f * F(EY, ix, iy, 0)
-             /*- dt/dx * (F(HZ, ix,iy,0) - F(HZ, ix-1,iy,0)) FIXME not in yz 2d
-              */
-             - (1.f - dt / dz) * F(HY, ix, iy, 0) - dt * F(JYI, ix, iy, 0)) /
-            (1.f + dt / dz);
-        }
-      }
-    } else {
-      assert(0);
+    Real3 dtdx = dt * Real3(grid.domain.dx_inv);
+
+    int d0 = (d + 0) % 3, d1 = (d + 1) % 3, d2 = (d + 2) % 3;
+    int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2;
+    int E1 = EX + d1, E2 = EX + d2;
+    int J1 = JXI + d1, J2 = JXI + d2;
+
+    Int3 start = mflds.ib();
+    Int3 stop = mflds.ib() + mflds.im();
+    start[d0] = -1;
+    stop[d0] = 0;
+
+    for (Int3 i3 : VecRange(start, stop)) {
+      Int3 edge_idx = i3;
+      edge_idx[d0] = 0;
+
+      Int3 neg1_d1{0, 0, 0};
+      neg1_d1[d1] = -1;
+
+      Int3 neg1_d2{0, 0, 0};
+      neg1_d2[d2] = -1;
+
+      F(H2, i3) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t), where d0=y */
+                   -2.f * (F(E1, edge_idx) - background_e[d1]) -
+                   dtdx[d2] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d2)) -
+                   (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) +
+                   dt * F(J1, edge_idx)) /
+                    (1.f + dtdx[d0]) +
+                  background_h[d2];
+      F(H1, i3) = (/* + 4.f * C_p_pulse_y1(x+.5*dx,y,z,t), where d0=y */
+                   +2.f * (F(E2, edge_idx) - background_e[d2]) -
+                   dtdx[d1] * (F(H0, edge_idx) - F(H0, edge_idx + neg1_d1)) -
+                   (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) +
+                   dt * F(J2, edge_idx)) /
+                    (1.f + dtdx[d0]) +
+                  background_h[d1];
     }
-#endif
   }
 
-  void open_H_hi(Mfields& mflds, int p, int d)
+  void radiative_H_hi(MfieldsState& mflds, int p, int d)
   {
+    set_upper_ghosts_to_nan(mflds, p, d, HX, false);
+
     auto F = make_Fields3d<dim_t>(mflds[p]);
     const Grid_t& grid = mflds.grid();
     Int3 ldims = grid.ldims;
-    Int3 ib = mflds.ib();
-    Int3 im = mflds.im();
     real_t dt = grid.dt;
-    real_t dx = grid.domain.dx[0];
-    real_t dy = grid.domain.dx[1];
-    real_t dz = grid.domain.dx[2];
-
-    // TODO
-    // assert(0);
-#if 1
-    if (d == 1) {
-      int my _mrc_unused = ldims[1];
-#ifdef DEBUG
-      for (int iz = -2; iz < ldims[2] + 2; iz++) {
-        for (int ix = std::max(-2, ib[0]);
-             ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
-          fields_t_set_nan(&F(HX, ix, my, iz));
-          fields_t_set_nan(&F(HX, ix, my + 1, iz));
-          fields_t_set_nan(&F(HY, ix, my, iz));
-          fields_t_set_nan(&F(HY, ix, my + 1, iz));
-          fields_t_set_nan(&F(HZ, ix, my + 1, iz));
-        }
-      }
-#endif
-      for (int iz = -2; iz < ldims[2] + 2; iz++) {
-        for (int ix = std::max(-2, ib[0]);
-             ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
-          F(HX, ix, my, iz) =
-            (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t) */
-             +2.f * F(EZ, ix, my, iz)
-             /*+ dt/dx * (F(HY, ix,my,iz) - F(HY, ix-1,my,iz)) */
-             - (1.f - dt / dy) * F(HX, ix, my - 1, iz) -
-             dt * F(JZI, ix, my, iz)) /
-            (1.f + dt / dy);
-          F(HZ, ix, my, iz) =
-            (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t) */
-             -2.f * F(EX, ix, my, iz) +
-             dt / dz * (F(HY, ix, my, iz) - F(HY, ix, my, iz - 1)) -
-             (1.f - dt / dy) * F(HZ, ix, my - 1, iz) -
-             dt * F(JXI, ix, my, iz)) /
-            (1.f + dt / dy);
-        }
-      }
-    } else if (d == 2) {
-      int mz = ldims[2];
-#ifdef DEBUG
-      for (int iy = -2; iy < ldims[1] + 2; iy++) {
-        for (int ix = -2; ix < ldims[0] + 2; ix++) {
-          fields_t_set_nan(&F(HX, ix, iy, mz));
-          fields_t_set_nan(&F(HX, ix, iy, mz + 1));
-          fields_t_set_nan(&F(HY, ix, iy, mz));
-          fields_t_set_nan(&F(HY, ix, iy, mz + 1));
-          fields_t_set_nan(&F(HZ, ix, iy, mz + 1));
-        }
-      }
-#endif
-      for (int iy = -2; iy < ldims[1] + 2; iy++) {
-        for (int ix = std::max(-2, ib[0]);
-             ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) {
-          F(HY, ix, iy, mz) =
-            (/* - 4.f * C_s_pulse_z2(x+0.5*dx,y,z,t) */
-             +2.f * F(EX, ix, iy, mz) +
-             dt / dy * (F(HZ, ix, iy, mz) - F(HZ, ix, iy - 1, mz)) -
-             (1.f - dt / dz) * F(HY, ix, iy, mz - 1) -
-             dt * F(JXI, ix, iy, mz)) /
-            (1.f + dt / dz);
-          F(HX, ix, iy, mz) = (/* + 4.f * C_p_pulse_z2(x+0.5*dx,y,z,t) */
-                               -2.f * F(EY, ix, iy, mz)
-                               /*+ dt/dx * (F(HZ, ix,iy,mz) - F(HZ, ix-1,iy,mz))
-                                  FIXME not in yz 2d*/
-                               - (1.f - dt / dz) * F(HX, ix, iy, mz - 1) +
-                               dt * F(JYI, ix, iy, mz)) /
-                              (1.f + dt / dz);
-        }
-      }
-    } else {
-      assert(0);
+    Real3 dtdx = dt * Real3(grid.domain.dx_inv);
+
+    int d0 = (d + 0) % 3, d1 = (d + 1) % 3, d2 = (d + 2) % 3;
+    int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2;
+    int E1 = EX + d1, E2 = EX + d2;
+    int J1 = JXI + d1, J2 = JXI + d2;
+
+    Int3 start = mflds.ib();
+    Int3 stop = mflds.ib() + mflds.im();
+    start[d0] = grid.ldims[d0];
+    stop[d0] = grid.ldims[d0] + 1;
+
+    for (Int3 i3 : VecRange(start, stop)) {
+      Int3 edge_idx = i3;
+      edge_idx[d0] -= 1;
+
+      Int3 neg1_d1{0, 0, 0};
+      neg1_d1[d1] = -1;
+
+      Int3 neg1_d2{0, 0, 0};
+      neg1_d2[d2] = -1;
+
+      F(H2, i3) = (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t), where d0=y */
+                   +2.f * (F(E1, i3) - background_e[d1]) +
+                   dtdx[d2] * (F(H0, i3) - F(H0, i3 + neg1_d2)) -
+                   (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) -
+                   dt * F(J1, i3)) /
+                    (1.f + dtdx[d0]) +
+                  background_h[d2];
+      F(H1, i3) = (/* + 4.f * C_p_pulse_y2(x+.5*dx,y,z,t), where d0=y */
+                   -2.f * (F(E2, i3) - background_e[d2]) +
+                   dtdx[d1] * (F(H0, i3) - F(H0, i3 + neg1_d1)) -
+                   (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) -
+                   dt * F(J2, i3)) /
+                    (1.f + dtdx[d0]) +
+                  background_h[d1];
     }
-#endif
   }
+
+  Vec3<real_t> background_e = {0.0, 0.0, 0.0};
+  Vec3<real_t> background_h = {0.0, 0.0, 0.0};
 };
 
 // ======================================================================
@@ -642,14 +625,14 @@ struct BndFields_ : BndFieldsBase
 
 // used by CUDA
 
-template <class MF>
+template <class MFIELDS_STATE>
 struct BndFieldsNone : BndFieldsBase
 {
-  using Mfields = MF;
+  using MfieldsState = MFIELDS_STATE;
 
   // clang-format off
-  void fill_ghosts_E(Mfields& mflds) {};
-  void fill_ghosts_H(Mfields& mflds) {};
-  void add_ghosts_J(Mfields& mflds) {};
+  void fill_ghosts_E(MfieldsState& mflds) {};
+  void fill_ghosts_H(MfieldsState& mflds) {};
+  void add_ghosts_J(MfieldsState& mflds) {};
   // clang-format on
 };
diff --git a/src/libpsc/psc_checks/checks_impl.hxx b/src/libpsc/psc_checks/checks_impl.hxx
index 1b098fec0..d9d2fea75 100644
--- a/src/libpsc/psc_checks/checks_impl.hxx
+++ b/src/libpsc/psc_checks/checks_impl.hxx
@@ -146,12 +146,15 @@ public:
 
     for (int p = 0; p < grid.n_patches(); p++) {
       for (int d = 0; d < 3; d++) {
-        Int3 r = grid.ldims;
-        r[d] = 1;
+        if (grid.atBoundaryLo(p, d) &&
+            (grid.bc.fld_lo[d] == BND_FLD_CONDUCTING_WALL ||
+             grid.bc.fld_lo[d] == BND_FLD_OPEN)) {
+
+          // account for implicit surface charges
+
+          Int3 r = grid.ldims;
+          r[d] = 1;
 
-        // account for implicit surface charges
-        if (grid.bc.fld_lo[d] == BND_FLD_CONDUCTING_WALL &&
-            grid.atBoundaryLo(p, d)) {
           rho.view(_s(0, r[0]), _s(0, r[1]), _s(0, r[2]), 0, p) =
             dive.view(_s(0, r[0]), _s(0, r[1]), _s(0, r[2]), 0, p);
         }