diff --git a/dev/nix/shell.nix b/dev/nix/shell.nix
index 13509dac2..1f21e82b0 100644
--- a/dev/nix/shell.nix
+++ b/dev/nix/shell.nix
@@ -26,10 +26,7 @@ pkgs.mkShell {
     python312Packages.jupyter
 
     cmake-format
-<<<<<<< HEAD
     cmake-lint
-=======
->>>>>>> 8b0f205a866f7f7534d8190e3ebee580ea09f7d8
     neocmakelsp
     black
     pyright
diff --git a/input.example.toml b/input.example.toml
index c5622d65a..788c30685 100644
--- a/input.example.toml
+++ b/input.example.toml
@@ -90,11 +90,11 @@
     # Boundary conditions for fields:
     #   @required
     #   @type: 1/2/3-size array of string tuples, each of size 1 or 2
-    #   @valid: "PERIODIC", "ABSORB", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON"
-    #   @example: [["CUSTOM", "ABSORB"]] (for 2D spherical [[rmin, rmax]])
+    #   @valid: "PERIODIC", "MATCH", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON"
+    #   @example: [["CUSTOM", "MATCH"]] (for 2D spherical [[rmin, rmax]])
     #   @note: When periodic in any of the directions, you should only set one value: [..., ["PERIODIC"], ...]
-    #   @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]): [["ATMOSPHERE", "ABSORB"]]
-    #   @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["ABSORB"]]
+    #   @note: In spherical, bondaries in theta/phi are set automatically (only specify bc @ [rmin, rmax]): [["ATMOSPHERE", "MATCH"]]
+    #   @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["MATCH"]]
     fields = ""
     # Boundary conditions for fields:
     #   @required
@@ -106,8 +106,8 @@
     #   @note: In GR, the horizon boundary is set automatically (only specify bc @ rmax): [["ABSORB"]]
     particles = ""
 
-    [grid.boundaries.absorb]
-      # Size of the absorption layer in physical (code) units:
+    [grid.boundaries.match]
+      # Size of the matching layer for fields in physical (code) units:
       #   @type: float
       #   @default: 1% of the domain size (in shortest dimension)
       #   @note: In spherical, this is the size of the layer in r from the outer wall
@@ -118,6 +118,14 @@
       #   @default: 1.0
       coeff = ""
 
+    [grid.boundaries.absorb]
+      # Size of the absorption layer for particles in physical (code) units:
+      #   @type: float
+      #   @default: 1% of the domain size (in shortest dimension)
+      #   @note: In spherical, this is the size of the layer in r from the outer wall
+      #   @note: In cartesian, this is the same for all dimensions where applicable
+      ds = ""
+
     [grid.boundaries.atmosphere]
       # @required: if ATMOSPHERE is one of the boundaries
       # Temperature of the atmosphere in units of m0 c^2
diff --git a/setups/srpic/monopole/monopole.toml b/legacy/_monopole/monopole.toml
similarity index 100%
rename from setups/srpic/monopole/monopole.toml
rename to legacy/_monopole/monopole.toml
diff --git a/setups/srpic/monopole/pgen.hpp b/legacy/_monopole/pgen.hpp
similarity index 97%
rename from setups/srpic/monopole/pgen.hpp
rename to legacy/_monopole/pgen.hpp
index 389a6c6f7..ed8877b71 100644
--- a/setups/srpic/monopole/pgen.hpp
+++ b/legacy/_monopole/pgen.hpp
@@ -86,7 +86,7 @@ namespace user {
 
     inline PGen() {}
 
-    auto FieldDriver(real_t time) const -> DriveFields<D> {
+    auto AtmFields(real_t time) const -> DriveFields<D> {
       return DriveFields<D> { time, Bsurf, Rstar, Omega };
     }
   };
diff --git a/setups/srpic/magnetar/pgen.hpp b/setups/srpic/magnetar/pgen.hpp
index cacbb7c9a..10f98ea5d 100644
--- a/setups/srpic/magnetar/pgen.hpp
+++ b/setups/srpic/magnetar/pgen.hpp
@@ -85,7 +85,7 @@ namespace user {
 
     const real_t  Bsurf, Rstar, Omega, gamma_pairs, pp_thres;
     InitFields<D> init_flds;
-    
+
     inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
       : arch::ProblemGenerator<S, M>(p)
       , global_domain { m }
@@ -94,12 +94,11 @@ namespace user {
       , Omega { p.template get<real_t>("setup.omega") }
       , pp_thres { p.template get<real_t>("setup.pp_thres") }
       , gamma_pairs { p.template get<real_t>("setup.gamma_pairs") }
-      , init_flds { Bsurf, Rstar } {
-      }
+      , init_flds { Bsurf, Rstar } {}
 
     inline PGen() {}
 
-    auto FieldDriver(real_t time) const -> DriveFields<D> {
+    auto AtmFields(real_t time) const -> DriveFields<D> {
       const real_t omega_t =
         Omega *
         ((ONE - math::tanh((static_cast<real_t>(5.0) - time) * HALF)) *
@@ -109,170 +108,172 @@ namespace user {
       return DriveFields<D> { time, Bsurf, Rstar, omega_t };
     }
 
-        void CustomPostStep(std::size_t , long double, Domain<S, M>& domain) {
-
-    // Ad-hoc PP kernel
-    {
-
-        auto& species2_e = domain.species[2];
-        auto& species2_p = domain.species[3];
-        auto& species3_e = domain.species[4];
-        auto& species3_p = domain.species[5];
-        auto metric = domain.mesh.metric;
-        auto pp_thres_ = this->pp_thres;
-        auto gamma_pairs_ = this->gamma_pairs;
-
-         for (std::size_t s { 0 }; s < 6; ++s) {
-            if (s == 1) {
-              continue;
-            }
-
-            array_t<std::size_t> elec_ind("elec_ind");
-            array_t<std::size_t> pos_ind("pos_ind");
-              
-            auto offset_e = species3_e.npart();
-            auto offset_p = species3_p.npart();
-
-            auto ux1_e    = species3_e.ux1;
-            auto ux2_e    = species3_e.ux2;
-            auto ux3_e    = species3_e.ux3;
-            auto i1_e     = species3_e.i1;
-            auto i2_e     = species3_e.i2;
-            auto dx1_e    = species3_e.dx1;
-            auto dx2_e    = species3_e.dx2;
-            auto phi_e    = species3_e.phi;
-            auto weight_e = species3_e.weight;
-            auto tag_e    = species3_e.tag;
-
-            auto ux1_p    = species3_p.ux1;
-            auto ux2_p    = species3_p.ux2;
-            auto ux3_p    = species3_p.ux3;
-            auto i1_p     = species3_p.i1;
-            auto i2_p     = species3_p.i2;
-            auto dx1_p    = species3_p.dx1;
-            auto dx2_p    = species3_p.dx2;
-            auto phi_p    = species3_p.phi;
-            auto weight_p = species3_p.weight;
-            auto tag_p    = species3_p.tag;
-
-            if (s == 0) {
-
-              offset_e = species2_e.npart();
-              offset_p = species2_p.npart();
-
-              ux1_e    = species2_e.ux1;
-              ux2_e    = species2_e.ux2;
-              ux3_e    = species2_e.ux3;
-              i1_e     = species2_e.i1;
-              i2_e     = species2_e.i2;
-              dx1_e    = species2_e.dx1;
-              dx2_e    = species2_e.dx2;
-              phi_e    = species2_e.phi;
-              weight_e = species2_e.weight;
-              tag_e    = species2_e.tag;
-
-              ux1_p    = species2_p.ux1;
-              ux2_p    = species2_p.ux2;
-              ux3_p    = species2_p.ux3;
-              i1_p     = species2_p.i1;
-              i2_p     = species2_p.i2;
-              dx1_p    = species2_p.dx1;
-              dx2_p    = species2_p.dx2;
-              phi_p    = species2_p.phi;
-              weight_p = species2_p.weight;
-              tag_p    = species2_p.tag;
-
-            }
-
-            auto& species = domain.species[s];
-            auto ux1    = species.ux1;
-            auto ux2    = species.ux2;
-            auto ux3    = species.ux3;
-            auto i1     = species.i1;
-            auto i2     = species.i2;
-            auto dx1    = species.dx1;
-            auto dx2    = species.dx2;
-            auto phi    = species.phi;
-            auto weight = species.weight;
-            auto tag    = species.tag;
-
-    Kokkos::parallel_for(
-        "InjectPairs", species.rangeActiveParticles(), Lambda(index_t p) {
-          if (tag(p) == ParticleTag::dead) {
-            return;
-          }
+    auto MatchFields(real_t) const -> InitFields<D> {
+      return InitFields<D> { Bsurf, Rstar };
+    }
 
-            auto px      = ux1(p);
-            auto py      = ux2(p);
-            auto pz      = ux3(p);
-            auto gamma   = math::sqrt(ONE + SQR(px) + SQR(py) + SQR(pz));
-
-            const coord_t<D> xCd{
-                static_cast<real_t>(i1(p)) + dx1(p),
-                static_cast<real_t>(i2(p)) + dx2(p)};
-
-            coord_t<D> xPh { ZERO };
-            metric.template convert<Crd::Cd, Crd::Ph>(xCd, xPh);
-
-          if ((gamma > pp_thres_) && (math::sin(xPh[1]) > 0.1)) {
-
-            auto new_gamma = gamma - 2.0 * gamma_pairs_;
-            auto new_fac = math::sqrt(SQR(new_gamma) - 1.0) / math::sqrt(SQR(gamma) - 1.0);
-            auto pair_fac = math::sqrt(SQR(gamma_pairs_) - 1.0) / math::sqrt(SQR(gamma) - 1.0);
-
-            auto elec_p = Kokkos::atomic_fetch_add(&elec_ind(), 1);
-            auto pos_p  = Kokkos::atomic_fetch_add(&pos_ind(), 1);
-
-              i1_e(elec_p + offset_e) = i1(p);
-              dx1_e(elec_p + offset_e) = dx1(p);
-              i2_e(elec_p + offset_e) = i2(p);
-              dx2_e(elec_p + offset_e) = dx2(p);
-              phi_e(elec_p + offset_e) = phi(p);
-              ux1_e(elec_p + offset_e) = px * pair_fac;
-              ux2_e(elec_p + offset_e) = py * pair_fac;
-              ux3_e(elec_p + offset_e) = pz * pair_fac;
-              weight_e(elec_p + offset_e) = weight(p);
-              tag_e(elec_p + offset_e) = ParticleTag::alive;
-
-              i1_p(pos_p + offset_p) = i1(p);
-              dx1_p(pos_p + offset_p) = dx1(p);
-              i2_p(pos_p + offset_p) = i2(p);
-              dx2_p(pos_p + offset_p) = dx2(p);
-              phi_p(pos_p + offset_p) = phi(p);
-              ux1_p(pos_p + offset_p) = px * pair_fac;
-              ux2_p(pos_p + offset_p) = py * pair_fac;
-              ux3_p(pos_p + offset_p) = pz * pair_fac;
-              weight_p(pos_p + offset_p) = weight(p);
-              tag_p(pos_p + offset_p) = ParticleTag::alive;
-
-              ux1(p) *= new_fac;
-              ux2(p) *= new_fac;
-              ux3(p) *= new_fac;
-          }
+    void CustomPostStep(std::size_t, long double, Domain<S, M>& domain) {
 
-        });
+      // Ad-hoc PP kernel
+      {
 
-            auto elec_ind_h = Kokkos::create_mirror(elec_ind);
-            Kokkos::deep_copy(elec_ind_h, elec_ind);
-            if (s == 0) {
-              species2_e.set_npart(offset_e + elec_ind_h());
-            } else {
-              species3_e.set_npart(offset_e + elec_ind_h());
-            }
+        auto& species2_e   = domain.species[2];
+        auto& species2_p   = domain.species[3];
+        auto& species3_e   = domain.species[4];
+        auto& species3_p   = domain.species[5];
+        auto  metric       = domain.mesh.metric;
+        auto  pp_thres_    = this->pp_thres;
+        auto  gamma_pairs_ = this->gamma_pairs;
 
-            auto pos_ind_h = Kokkos::create_mirror(pos_ind);
-            Kokkos::deep_copy(pos_ind_h, pos_ind);
-            if (s == 0) {
-              species2_p.set_npart(offset_p + pos_ind_h());
-            } else {
-              species3_p.set_npart(offset_p + pos_ind_h());
-            }
+        for (std::size_t s { 0 }; s < 6; ++s) {
+          if (s == 1) {
+            continue;
+          }
 
+          array_t<std::size_t> elec_ind("elec_ind");
+          array_t<std::size_t> pos_ind("pos_ind");
+
+          auto offset_e = species3_e.npart();
+          auto offset_p = species3_p.npart();
+
+          auto ux1_e    = species3_e.ux1;
+          auto ux2_e    = species3_e.ux2;
+          auto ux3_e    = species3_e.ux3;
+          auto i1_e     = species3_e.i1;
+          auto i2_e     = species3_e.i2;
+          auto dx1_e    = species3_e.dx1;
+          auto dx2_e    = species3_e.dx2;
+          auto phi_e    = species3_e.phi;
+          auto weight_e = species3_e.weight;
+          auto tag_e    = species3_e.tag;
+
+          auto ux1_p    = species3_p.ux1;
+          auto ux2_p    = species3_p.ux2;
+          auto ux3_p    = species3_p.ux3;
+          auto i1_p     = species3_p.i1;
+          auto i2_p     = species3_p.i2;
+          auto dx1_p    = species3_p.dx1;
+          auto dx2_p    = species3_p.dx2;
+          auto phi_p    = species3_p.phi;
+          auto weight_p = species3_p.weight;
+          auto tag_p    = species3_p.tag;
+
+          if (s == 0) {
+
+            offset_e = species2_e.npart();
+            offset_p = species2_p.npart();
+
+            ux1_e    = species2_e.ux1;
+            ux2_e    = species2_e.ux2;
+            ux3_e    = species2_e.ux3;
+            i1_e     = species2_e.i1;
+            i2_e     = species2_e.i2;
+            dx1_e    = species2_e.dx1;
+            dx2_e    = species2_e.dx2;
+            phi_e    = species2_e.phi;
+            weight_e = species2_e.weight;
+            tag_e    = species2_e.tag;
+
+            ux1_p    = species2_p.ux1;
+            ux2_p    = species2_p.ux2;
+            ux3_p    = species2_p.ux3;
+            i1_p     = species2_p.i1;
+            i2_p     = species2_p.i2;
+            dx1_p    = species2_p.dx1;
+            dx2_p    = species2_p.dx2;
+            phi_p    = species2_p.phi;
+            weight_p = species2_p.weight;
+            tag_p    = species2_p.tag;
+          }
+
+          auto& species = domain.species[s];
+          auto  ux1     = species.ux1;
+          auto  ux2     = species.ux2;
+          auto  ux3     = species.ux3;
+          auto  i1      = species.i1;
+          auto  i2      = species.i2;
+          auto  dx1     = species.dx1;
+          auto  dx2     = species.dx2;
+          auto  phi     = species.phi;
+          auto  weight  = species.weight;
+          auto  tag     = species.tag;
+
+          Kokkos::parallel_for(
+            "InjectPairs",
+            species.rangeActiveParticles(),
+            Lambda(index_t p) {
+              if (tag(p) == ParticleTag::dead) {
+                return;
+              }
+
+              auto px    = ux1(p);
+              auto py    = ux2(p);
+              auto pz    = ux3(p);
+              auto gamma = math::sqrt(ONE + SQR(px) + SQR(py) + SQR(pz));
+
+              const coord_t<D> xCd { static_cast<real_t>(i1(p)) + dx1(p),
+                                     static_cast<real_t>(i2(p)) + dx2(p) };
+
+              coord_t<D> xPh { ZERO };
+              metric.template convert<Crd::Cd, Crd::Ph>(xCd, xPh);
+
+              if ((gamma > pp_thres_) && (math::sin(xPh[1]) > 0.1)) {
+
+                auto new_gamma = gamma - 2.0 * gamma_pairs_;
+                auto new_fac   = math::sqrt(SQR(new_gamma) - 1.0) /
+                               math::sqrt(SQR(gamma) - 1.0);
+                auto pair_fac = math::sqrt(SQR(gamma_pairs_) - 1.0) /
+                                math::sqrt(SQR(gamma) - 1.0);
+
+                auto elec_p = Kokkos::atomic_fetch_add(&elec_ind(), 1);
+                auto pos_p  = Kokkos::atomic_fetch_add(&pos_ind(), 1);
+
+                i1_e(elec_p + offset_e)     = i1(p);
+                dx1_e(elec_p + offset_e)    = dx1(p);
+                i2_e(elec_p + offset_e)     = i2(p);
+                dx2_e(elec_p + offset_e)    = dx2(p);
+                phi_e(elec_p + offset_e)    = phi(p);
+                ux1_e(elec_p + offset_e)    = px * pair_fac;
+                ux2_e(elec_p + offset_e)    = py * pair_fac;
+                ux3_e(elec_p + offset_e)    = pz * pair_fac;
+                weight_e(elec_p + offset_e) = weight(p);
+                tag_e(elec_p + offset_e)    = ParticleTag::alive;
+
+                i1_p(pos_p + offset_p)     = i1(p);
+                dx1_p(pos_p + offset_p)    = dx1(p);
+                i2_p(pos_p + offset_p)     = i2(p);
+                dx2_p(pos_p + offset_p)    = dx2(p);
+                phi_p(pos_p + offset_p)    = phi(p);
+                ux1_p(pos_p + offset_p)    = px * pair_fac;
+                ux2_p(pos_p + offset_p)    = py * pair_fac;
+                ux3_p(pos_p + offset_p)    = pz * pair_fac;
+                weight_p(pos_p + offset_p) = weight(p);
+                tag_p(pos_p + offset_p)    = ParticleTag::alive;
+
+                ux1(p) *= new_fac;
+                ux2(p) *= new_fac;
+                ux3(p) *= new_fac;
+              }
+            });
+
+          auto elec_ind_h = Kokkos::create_mirror(elec_ind);
+          Kokkos::deep_copy(elec_ind_h, elec_ind);
+          if (s == 0) {
+            species2_e.set_npart(offset_e + elec_ind_h());
+          } else {
+            species3_e.set_npart(offset_e + elec_ind_h());
+          }
 
+          auto pos_ind_h = Kokkos::create_mirror(pos_ind);
+          Kokkos::deep_copy(pos_ind_h, pos_ind);
+          if (s == 0) {
+            species2_p.set_npart(offset_p + pos_ind_h());
+          } else {
+            species3_p.set_npart(offset_p + pos_ind_h());
           }
-        } // Ad-hoc PP kernel
         }
-  
+      } // Ad-hoc PP kernel
+    }
   };
 
 } // namespace user
diff --git a/setups/srpic/magnetosphere/pgen.hpp b/setups/srpic/magnetosphere/pgen.hpp
index 681c4d6d1..64fe13cfe 100644
--- a/setups/srpic/magnetosphere/pgen.hpp
+++ b/setups/srpic/magnetosphere/pgen.hpp
@@ -86,9 +86,13 @@ namespace user {
 
     inline PGen() {}
 
-    auto FieldDriver(real_t time) const -> DriveFields<D> {
+    auto AtmFields(real_t time) const -> DriveFields<D> {
       return DriveFields<D> { time, Bsurf, Rstar, Omega };
     }
+
+    auto MatchFields(real_t) const -> InitFields<D> {
+      return InitFields<D> { Bsurf, Rstar };
+    }
   };
 
 } // namespace user
diff --git a/setups/srpic/shock/pgen.hpp b/setups/srpic/shock/pgen.hpp
index 1eedb3a01..b8f169521 100644
--- a/setups/srpic/shock/pgen.hpp
+++ b/setups/srpic/shock/pgen.hpp
@@ -13,6 +13,8 @@
 #include "archetypes/problem_generator.h"
 #include "framework/domain/metadomain.h"
 
+#include <utility>
+
 namespace user {
   using namespace ntt;
 
@@ -28,7 +30,7 @@ namespace user {
       @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 * static_cast<real_t>(convert::deg2rad) }
+      : Bmag { bmag }
       , Btheta { btheta * static_cast<real_t>(convert::deg2rad) }
       , Bphi { bphi * static_cast<real_t>(convert::deg2rad) }
       , Vx { drift_ux } {}
@@ -93,20 +95,21 @@ namespace user {
 
     inline PGen() {}
 
-    auto FixField(const em& comp) const -> real_t {
+    auto FixFieldsConst(const bc_in&, const em& comp) const
+      -> std::pair<real_t, bool> {
       if (comp == em::ex2) {
-        return init_flds.ex2({ ZERO });
+        return { init_flds.ex2({ ZERO }), true };
       } else if (comp == em::ex3) {
-        return init_flds.ex3({ ZERO });
-      } else if (comp == em::bx1) {
-        return init_flds.bx1({ ZERO });
+        return { init_flds.ex3({ ZERO }), true };
       } else {
-        raise::Error("Other components should not be requested when BC is in X",
-                     HERE);
-        return ZERO;
+        return { ZERO, false };
       }
     }
 
+    auto MatchFields(real_t time) const -> InitFields<D> {
+      return init_flds;
+    }
+
     inline void InitPrtls(Domain<S, M>& local_domain) {
       const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
                                                       local_domain.random_pool,
diff --git a/setups/srpic/shock/shock.toml b/setups/srpic/shock/shock.toml
index c8a6f7c9d..4ed3a2b9e 100644
--- a/setups/srpic/shock/shock.toml
+++ b/setups/srpic/shock/shock.toml
@@ -11,7 +11,7 @@
     metric = "minkowski"
 
   [grid.boundaries]
-    fields = [["FIXED", "ABSORB"], ["PERIODIC"]]
+    fields = [["ABSORB", "FIXED"], ["PERIODIC"]]
     particles = [["REFLECT", "ABSORB"], ["PERIODIC"]]
 
 [scales]
diff --git a/setups/wip/magpump/pgen.hpp b/setups/wip/magpump/pgen.hpp
new file mode 100644
index 000000000..21d4c8882
--- /dev/null
+++ b/setups/wip/magpump/pgen.hpp
@@ -0,0 +1,170 @@
+#ifndef PROBLEM_GENERATOR_H
+#define PROBLEM_GENERATOR_H
+
+#include "enums.h"
+#include "global.h"
+
+#include "arch/traits.h"
+
+#include "archetypes/particle_injector.h"
+#include "archetypes/problem_generator.h"
+#include "framework/domain/metadomain.h"
+
+#include <plog/Log.h>
+
+namespace user {
+  using namespace ntt;
+
+  template <Dimension D>
+  struct InitFields {
+    InitFields(real_t bsurf, real_t rstar) : Bsurf { bsurf }, Rstar { rstar } {}
+
+    Inline auto bx1(const coord_t<D>& x_Ph) const -> real_t {
+      return Bsurf * math::cos(x_Ph[1]) / CUBE(x_Ph[0] / Rstar);
+    }
+
+    Inline auto bx2(const coord_t<D>& x_Ph) const -> real_t {
+      return Bsurf * HALF * math::sin(x_Ph[1]) / CUBE(x_Ph[0] / Rstar);
+    }
+
+  private:
+    const real_t Bsurf, Rstar;
+  };
+
+  template <Dimension D>
+  struct DriveFields : public InitFields<D> {
+    DriveFields(real_t time, real_t bsurf, real_t rstar)
+      : InitFields<D> { bsurf, rstar }
+      , time { time } {}
+
+    using InitFields<D>::bx1;
+    using InitFields<D>::bx2;
+
+    Inline auto bx3(const coord_t<D>&) const -> real_t {
+      return ZERO;
+    }
+
+    Inline auto ex1(const coord_t<D>&) const -> real_t {
+      return ZERO;
+    }
+
+    Inline auto ex2(const coord_t<D>& x_Ph) const -> real_t {
+      return ZERO;
+    }
+
+    Inline auto ex3(const coord_t<D>&) const -> real_t {
+      return ZERO;
+    }
+
+  private:
+    const real_t time;
+  };
+
+  template <SimEngine::type S, class M>
+  struct Inflow : public arch::EnergyDistribution<S, M> {
+    Inflow(const M& metric, real_t vin)
+      : arch::EnergyDistribution<S, M> { metric }
+      , vin { vin } {}
+
+    Inline void operator()(const coord_t<M::Dim>&,
+                           vec_t<Dim::_3D>& v_Ph,
+                           unsigned short) const override {
+      v_Ph[0] = -vin;
+    }
+
+  private:
+    const real_t vin;
+  };
+
+  template <SimEngine::type S, class M>
+  struct Sphere : public arch::SpatialDistribution<S, M> {
+    Sphere(const M& metric, real_t r0, real_t dr)
+      : arch::SpatialDistribution<S, M> { metric }
+      , r0 { r0 }
+      , dr { dr } {}
+
+    Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t override {
+      return math::exp(-SQR((x_Ph[0] - r0) / dr)) *
+             (x_Ph[1] > 0.25 && x_Ph[1] < constant::PI - 0.25);
+    }
+
+  private:
+    const real_t r0, dr;
+  };
+
+  template <SimEngine::type S, class M>
+  struct PGen : public arch::ProblemGenerator<S, M> {
+    static constexpr auto engines { traits::compatible_with<SimEngine::SRPIC>::value };
+    static constexpr auto metrics {
+      traits::compatible_with<Metric::Spherical, Metric::QSpherical>::value
+    };
+    static constexpr auto dimensions { traits::compatible_with<Dim::_2D>::value };
+
+    using arch::ProblemGenerator<S, M>::D;
+    using arch::ProblemGenerator<S, M>::C;
+    using arch::ProblemGenerator<S, M>::params;
+
+    const real_t  Bsurf, pump_period, pump_ampl, pump_radius, Rstar;
+    const real_t  vin, drinj;
+    InitFields<D> init_flds;
+
+    inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
+      : arch::ProblemGenerator<S, M> { p }
+      , Bsurf { p.template get<real_t>("setup.Bsurf", ONE) }
+      , pump_period { p.template get<real_t>("setup.pump_period") }
+      , pump_ampl { p.template get<real_t>("setup.pump_ampl") }
+      , pump_radius { p.template get<real_t>("setup.pump_radius") }
+      , Rstar { m.mesh().extent(in::x1).first }
+      , vin { p.template get<real_t>("setup.vin") }
+      , drinj { p.template get<real_t>("setup.drinj") }
+      , init_flds { Bsurf, Rstar } {}
+
+    auto FieldDriver(real_t time) const -> DriveFields<D> {
+      return DriveFields<D> { time, Bsurf, Rstar };
+    }
+
+    void CustomPostStep(std::size_t, long double time, Domain<S, M>& domain) {
+      const real_t radius = pump_radius +
+                            pump_ampl *
+                              math::sin(time * constant::TWO_PI / pump_period);
+      const real_t dr     = 1.0;
+      const auto&  metric = domain.mesh.metric;
+      auto         EM     = domain.fields.em;
+      Kokkos::parallel_for(
+        "outerBC",
+        domain.mesh.rangeActiveCells(),
+        Lambda(index_t i1, index_t i2) {
+          const auto i1_ = COORD(i1), i2_ = COORD(i2);
+          const auto r = metric.template convert<1, Crd::Cd, Crd::Ph>(i1_);
+          if (r > radius - 5 * dr) {
+            const auto smooth   = HALF * (ONE - math::tanh((r - radius) / dr));
+            EM(i1, i2, em::ex1) = smooth * EM(i1, i2, em::ex1);
+            EM(i1, i2, em::ex2) = smooth * EM(i1, i2, em::ex2);
+            EM(i1, i2, em::ex3) = smooth * EM(i1, i2, em::ex3);
+            EM(i1, i2, em::bx1) = smooth * EM(i1, i2, em::bx1);
+            EM(i1, i2, em::bx2) = smooth * EM(i1, i2, em::bx2);
+            EM(i1, i2, em::bx3) = smooth * EM(i1, i2, em::bx3);
+          }
+        });
+
+      if (time < pump_period * 0.25) {
+        const auto energy_dist = Inflow<S, M>(domain.mesh.metric, vin);
+        const auto spatial_dist = Sphere<S, M>(domain.mesh.metric, radius, drinj);
+        const auto injector = arch::NonUniformInjector<S, M, Inflow, Sphere>(
+          energy_dist,
+          spatial_dist,
+          { 1, 2 });
+
+        arch::InjectNonUniform<S, M, arch::NonUniformInjector<S, M, Inflow, Sphere>>(
+          params,
+          domain,
+          injector,
+          ONE,
+          true);
+      }
+    }
+  };
+
+} // namespace user
+
+#endif
diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp
index d8cb1e64c..9f5e4551f 100644
--- a/src/engines/srpic.hpp
+++ b/src/engines/srpic.hpp
@@ -584,8 +584,8 @@ namespace ntt {
 
     void FieldBoundaries(domain_t& domain, BCTags tags) {
       for (auto& direction : dir::Directions<M::Dim>::orth) {
-        if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::ABSORB) {
-          AbsorbFieldsIn(direction, domain, tags);
+        if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::MATCH) {
+          MatchFieldsIn(direction, domain, tags);
         } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::AXIS) {
           if (domain.mesh.flds_bc_in(direction) == FldsBC::AXIS) {
             AxisFieldsIn(direction, domain, tags);
@@ -606,14 +606,13 @@ namespace ntt {
       } // loop over directions
     }
 
-    void AbsorbFieldsIn(dir::direction_t<M::Dim> direction,
-                        domain_t&                domain,
-                        BCTags                   tags) {
+    void MatchFieldsIn(dir::direction_t<M::Dim> direction,
+                       domain_t&                domain,
+                       BCTags                   tags) {
       /**
-       * absorbing boundaries
+       * matching boundaries
        */
-      const auto ds = m_params.template get<real_t>(
-        "grid.boundaries.absorb.ds");
+      const auto ds = m_params.template get<real_t>("grid.boundaries.match.ds");
       const auto dim = direction.get_dim();
       real_t     xg_min, xg_max, xg_edge;
       auto       sign = direction.get_sign();
@@ -652,40 +651,49 @@ namespace ntt {
         range_min[d] = intersect_range[d].first;
         range_max[d] = intersect_range[d].second;
       }
-      if (dim == in::x1) {
-        Kokkos::parallel_for(
-          "AbsorbFields",
-          CreateRangePolicy<M::Dim>(range_min, range_max),
-          kernel::AbsorbBoundaries_kernel<M, 1>(domain.fields.em,
-                                                domain.mesh.metric,
-                                                xg_edge,
-                                                ds,
-                                                tags));
-      } else if (dim == in::x2) {
-        if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) {
-          Kokkos::parallel_for(
-            "AbsorbFields",
-            CreateRangePolicy<M::Dim>(range_min, range_max),
-            kernel::AbsorbBoundaries_kernel<M, 2>(domain.fields.em,
-                                                  domain.mesh.metric,
-                                                  xg_edge,
-                                                  ds,
-                                                  tags));
-        } else {
-          raise::Error("Invalid dimension", HERE);
-        }
-      } else if (dim == in::x3) {
-        if constexpr (M::Dim == Dim::_3D) {
+      if constexpr (traits::has_member<traits::pgen::match_fields_t, pgen_t>::value) {
+        auto match_fields = m_pgen.MatchFields(time);
+        if (dim == in::x1) {
           Kokkos::parallel_for(
-            "AbsorbFields",
+            "MatchFields",
             CreateRangePolicy<M::Dim>(range_min, range_max),
-            kernel::AbsorbBoundaries_kernel<M, 3>(domain.fields.em,
-                                                  domain.mesh.metric,
-                                                  xg_edge,
-                                                  ds,
-                                                  tags));
-        } else {
-          raise::Error("Invalid dimension", HERE);
+            kernel::bc::MatchBoundaries_kernel<SimEngine::SRPIC, decltype(match_fields), M, in::x1>(
+              domain.fields.em,
+              match_fields,
+              domain.mesh.metric,
+              xg_edge,
+              ds,
+              tags));
+        } else if (dim == in::x2) {
+          if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) {
+            Kokkos::parallel_for(
+              "MatchFields",
+              CreateRangePolicy<M::Dim>(range_min, range_max),
+              kernel::bc::MatchBoundaries_kernel<SimEngine::SRPIC, decltype(match_fields), M, in::x2>(
+                domain.fields.em,
+                match_fields,
+                domain.mesh.metric,
+                xg_edge,
+                ds,
+                tags));
+          } else {
+            raise::Error("Invalid dimension", HERE);
+          }
+        } else if (dim == in::x3) {
+          if constexpr (M::Dim == Dim::_3D) {
+            Kokkos::parallel_for(
+              "MatchFields",
+              CreateRangePolicy<M::Dim>(range_min, range_max),
+              kernel::bc::MatchBoundaries_kernel<SimEngine::SRPIC, decltype(match_fields), M, in::x3>(
+                domain.fields.em,
+                match_fields,
+                domain.mesh.metric,
+                xg_edge,
+                ds,
+                tags));
+          } else {
+            raise::Error("Invalid dimension", HERE);
+          }
         }
       }
     }
@@ -708,12 +716,16 @@ namespace ntt {
         Kokkos::parallel_for(
           "AxisBCFields",
           domain.mesh.n_all(in::x1),
-          kernel::AxisBoundaries_kernel<M::Dim, false>(domain.fields.em, i2_min, tags));
+          kernel::bc::AxisBoundaries_kernel<M::Dim, false>(domain.fields.em,
+                                                           i2_min,
+                                                           tags));
       } else {
         Kokkos::parallel_for(
           "AxisBCFields",
           domain.mesh.n_all(in::x1),
-          kernel::AxisBoundaries_kernel<M::Dim, true>(domain.fields.em, i2_max, tags));
+          kernel::bc::AxisBoundaries_kernel<M::Dim, true>(domain.fields.em,
+                                                          i2_max,
+                                                          tags));
       }
     }
 
@@ -774,40 +786,51 @@ namespace ntt {
       if (tags & BC::B) {
         comps.push_back(normal_b_comp);
       }
-      if constexpr (traits::has_member<traits::pgen::field_driver_t, pgen_t>::value) {
-        raise::Error("Field driver for fixed fields not implemented", HERE);
-      } else {
-        // if field driver not present, set fields to fixed values
+      if constexpr (traits::has_member<traits::pgen::fix_fields_t, pgen_t>::value) {
+        raise::Error("Non-const fixed fields not implemented", HERE);
+      } else if constexpr (
+        traits::has_member<traits::pgen::fix_fields_const_t, pgen_t>::value) {
         for (const auto& comp : comps) {
-          auto value = ZERO;
+          auto value     = ZERO;
+          bool shouldset = false;
           if constexpr (
-            traits::has_member<traits::pgen::fix_field_t, pgen_t>::value) {
+            traits::has_member<traits::pgen::fix_fields_const_t, pgen_t>::value) {
             // if fix field function present, read from it
-            value = m_pgen.FixField((em)comp);
+            const auto newset = m_pgen.FixFieldsConst(
+              (bc_in)(sign * ((short)dim + 1)),
+              (em)comp);
+            value     = newset.first;
+            shouldset = newset.second;
           }
-          if constexpr (M::Dim == Dim::_1D) {
-            Kokkos::deep_copy(Kokkos::subview(domain.fields.em,
-                                              std::make_pair(xi_min[0], xi_max[0]),
-                                              comp),
-                              value);
-          } else if constexpr (M::Dim == Dim::_2D) {
-            Kokkos::deep_copy(Kokkos::subview(domain.fields.em,
-                                              std::make_pair(xi_min[0], xi_max[0]),
-                                              std::make_pair(xi_min[1], xi_max[1]),
-                                              comp),
-                              value);
-          } else if constexpr (M::Dim == Dim::_3D) {
-            Kokkos::deep_copy(
-              Kokkos::subview(domain.fields.em,
-                              std::make_pair(xi_min[0], xi_max[0]),
-                              std::make_pair(xi_min[1], xi_max[1]),
-                              std::make_pair(xi_min[2], xi_max[2]),
-                              comp),
-              value);
-          } else {
-            raise::Error("Invalid dimension", HERE);
+          if (shouldset) {
+            if constexpr (M::Dim == Dim::_1D) {
+              Kokkos::deep_copy(
+                Kokkos::subview(domain.fields.em,
+                                std::make_pair(xi_min[0], xi_max[0]),
+                                comp),
+                value);
+            } else if constexpr (M::Dim == Dim::_2D) {
+              Kokkos::deep_copy(
+                Kokkos::subview(domain.fields.em,
+                                std::make_pair(xi_min[0], xi_max[0]),
+                                std::make_pair(xi_min[1], xi_max[1]),
+                                comp),
+                value);
+            } else if constexpr (M::Dim == Dim::_3D) {
+              Kokkos::deep_copy(
+                Kokkos::subview(domain.fields.em,
+                                std::make_pair(xi_min[0], xi_max[0]),
+                                std::make_pair(xi_min[1], xi_max[1]),
+                                std::make_pair(xi_min[2], xi_max[2]),
+                                comp),
+                value);
+            } else {
+              raise::Error("Invalid dimension", HERE);
+            }
           }
         }
+      } else {
+        raise::Error("Fixed fields not present (both const and non-const)", HERE);
       }
     }
 
@@ -817,7 +840,7 @@ namespace ntt {
       /**
        * atmosphere field boundaries
        */
-      if constexpr (traits::has_member<traits::pgen::field_driver_t, pgen_t>::value) {
+      if constexpr (traits::has_member<traits::pgen::atm_fields_t, pgen_t>::value) {
         const auto [sign, dim, xg_min, xg_max] = get_atm_extent(direction);
         const auto           dd = static_cast<unsigned short>(dim);
         boundaries_t<real_t> box;
@@ -846,7 +869,7 @@ namespace ntt {
           range_min[d] = intersect_range[d].first;
           range_max[d] = intersect_range[d].second;
         }
-        auto        field_driver = m_pgen.FieldDriver(time);
+        auto        atm_fields = m_pgen.AtmFields(time);
         std::size_t il_edge;
         if (sign > 0) {
           il_edge = range_min[dd] - N_GHOSTS;
@@ -859,9 +882,9 @@ namespace ntt {
             Kokkos::parallel_for(
               "AtmosphereBCFields",
               range,
-              kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, true, in::x1>(
+              kernel::bc::EnforcedBoundaries_kernel<decltype(atm_fields), M, true, in::x1>(
                 domain.fields.em,
-                field_driver,
+                atm_fields,
                 domain.mesh.metric,
                 il_edge,
                 tags));
@@ -869,9 +892,9 @@ namespace ntt {
             Kokkos::parallel_for(
               "AtmosphereBCFields",
               range,
-              kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, false, in::x1>(
+              kernel::bc::EnforcedBoundaries_kernel<decltype(atm_fields), M, false, in::x1>(
                 domain.fields.em,
-                field_driver,
+                atm_fields,
                 domain.mesh.metric,
                 il_edge,
                 tags));
@@ -882,9 +905,9 @@ namespace ntt {
               Kokkos::parallel_for(
                 "AtmosphereBCFields",
                 range,
-                kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, true, in::x2>(
+                kernel::bc::EnforcedBoundaries_kernel<decltype(atm_fields), M, true, in::x2>(
                   domain.fields.em,
-                  field_driver,
+                  atm_fields,
                   domain.mesh.metric,
                   il_edge,
                   tags));
@@ -892,9 +915,9 @@ namespace ntt {
               Kokkos::parallel_for(
                 "AtmosphereBCFields",
                 range,
-                kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, false, in::x2>(
+                kernel::bc::EnforcedBoundaries_kernel<decltype(atm_fields), M, false, in::x2>(
                   domain.fields.em,
-                  field_driver,
+                  atm_fields,
                   domain.mesh.metric,
                   il_edge,
                   tags));
@@ -908,9 +931,9 @@ namespace ntt {
               Kokkos::parallel_for(
                 "AtmosphereBCFields",
                 range,
-                kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, true, in::x3>(
+                kernel::bc::EnforcedBoundaries_kernel<decltype(atm_fields), M, true, in::x3>(
                   domain.fields.em,
-                  field_driver,
+                  atm_fields,
                   domain.mesh.metric,
                   il_edge,
                   tags));
@@ -918,9 +941,9 @@ namespace ntt {
               Kokkos::parallel_for(
                 "AtmosphereBCFields",
                 range,
-                kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, false, in::x3>(
+                kernel::bc::EnforcedBoundaries_kernel<decltype(atm_fields), M, false, in::x3>(
                   domain.fields.em,
-                  field_driver,
+                  atm_fields,
                   domain.mesh.metric,
                   il_edge,
                   tags));
@@ -932,8 +955,7 @@ namespace ntt {
           raise::Error("Invalid dimension", HERE);
         }
       } else {
-        raise::Error("Field driver not implemented in PGEN for atmosphere BCs",
-                     HERE);
+        raise::Error("Atm fields not implemented in PGEN for atmosphere BCs", HERE);
       }
     }
 
diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp
index 0a605651c..4a9b3056a 100644
--- a/src/framework/parameters.cpp
+++ b/src/framework/parameters.cpp
@@ -47,7 +47,7 @@ namespace ntt {
 
   /*
    * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-   * Parameters that must not be changed during after the checkpoint restart
+   * Parameters that must not be changed during the checkpoint restart
    * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
    */
   void SimulationParams::setImmutableParams(const toml::value& toml_data) {
@@ -322,7 +322,7 @@ namespace ntt {
 
   /*
    * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
-   * Parameters that may be changed during after the checkpoint restart
+   * Parameters that may be changed during the checkpoint restart
    * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
    */
   void SimulationParams::setMutableParams(const toml::value& toml_data) {
@@ -351,9 +351,9 @@ namespace ntt {
       auto atm_defined = false;
       for (const auto& bcs : flds_bc) {
         for (const auto& bc : bcs) {
-          if (fmt::toLower(bc) == "absorb") {
-            promiseToDefine("grid.boundaries.absorb.ds");
-            promiseToDefine("grid.boundaries.absorb.coeff");
+          if (fmt::toLower(bc) == "match") {
+            promiseToDefine("grid.boundaries.match.ds");
+            promiseToDefine("grid.boundaries.match.coeff");
           }
           if (fmt::toLower(bc) == "atmosphere") {
             raise::ErrorIf(atm_defined,
@@ -386,7 +386,6 @@ namespace ntt {
         for (const auto& bc : bcs) {
           if (fmt::toLower(bc) == "absorb") {
             promiseToDefine("grid.boundaries.absorb.ds");
-            promiseToDefine("grid.boundaries.absorb.coeff");
           }
           if (fmt::toLower(bc) == "atmosphere") {
             raise::ErrorIf(atm_defined,
@@ -731,6 +730,38 @@ namespace ntt {
     set("grid.boundaries.fields", flds_bc_pairwise);
     set("grid.boundaries.particles", prtl_bc_pairwise);
 
+    if (isPromised("grid.boundaries.match.ds")) {
+      if (coord_enum == Coord::Cart) {
+        auto min_extent = std::numeric_limits<real_t>::max();
+        for (const auto& e : extent_pairwise) {
+          min_extent = std::min(min_extent, e.second - e.first);
+        }
+        set("grid.boundaries.match.ds",
+            toml::find_or(toml_data,
+                          "grid",
+                          "boundaries",
+                          "match",
+                          "ds",
+                          min_extent * defaults::bc::match::ds_frac));
+      } else {
+        auto r_extent = extent_pairwise[0].second - extent_pairwise[0].first;
+        set("grid.boundaries.match.ds",
+            toml::find_or(toml_data,
+                          "grid",
+                          "boundaries",
+                          "match",
+                          "ds",
+                          r_extent * defaults::bc::match::ds_frac));
+      }
+      set("grid.boundaries.match.coeff",
+          toml::find_or(toml_data,
+                        "grid",
+                        "boundaries",
+                        "match",
+                        "coeff",
+                        defaults::bc::match::coeff));
+    }
+
     if (isPromised("grid.boundaries.absorb.ds")) {
       if (coord_enum == Coord::Cart) {
         auto min_extent = std::numeric_limits<real_t>::max();
@@ -754,13 +785,6 @@ namespace ntt {
                           "ds",
                           r_extent * defaults::bc::absorb::ds_frac));
       }
-      set("grid.boundaries.absorb.coeff",
-          toml::find_or(toml_data,
-                        "grid",
-                        "boundaries",
-                        "absorb",
-                        "coeff",
-                        defaults::bc::absorb::coeff));
     }
 
     if (isPromised("grid.boundaries.atmosphere.temperature")) {
diff --git a/src/framework/tests/parameters.cpp b/src/framework/tests/parameters.cpp
index 1a4228642..7cd5ce46a 100644
--- a/src/framework/tests/parameters.cpp
+++ b/src/framework/tests/parameters.cpp
@@ -32,7 +32,7 @@ const auto mink_1d = u8R"(
     fields = [["PERIODIC"]]
     particles = [["ABSORB", "ABSORB"]]
 
-    [grid.boundaries.absorb]
+    [grid.boundaries.match]
       coeff = 10.0
       ds = 0.025
 
@@ -101,10 +101,10 @@ const auto sph_2d = u8R"(
     metric = "spherical"
 
   [grid.boundaries]
-    fields = [["ATMOSPHERE", "ABSORB"]]
+    fields = [["ATMOSPHERE", "MATCH"]]
     particles = [["ATMOSPHERE", "ABSORB"]]
 
-    [grid.boundaries.absorb]
+    [grid.boundaries.match]
       coeff = 10.0
 
     [grid.boundaries.atmosphere]
@@ -180,7 +180,7 @@ const auto qks_2d = u8R"(
     ks_a = 0.99
 
   [grid.boundaries]
-    fields = [["ABSORB"]]
+    fields = [["MATCH"]]
     particles = [["ABSORB"]]
 
 [scales]
@@ -345,8 +345,8 @@ auto main(int argc, char* argv[]) -> int {
                               "simulation.engine");
 
       boundaries_t<FldsBC> fbc = {
-        { FldsBC::ATMOSPHERE, FldsBC::ABSORB },
-        {       FldsBC::AXIS,   FldsBC::AXIS }
+        { FldsBC::ATMOSPHERE, FldsBC::MATCH },
+        {       FldsBC::AXIS,  FldsBC::AXIS }
       };
 
       assert_equal<real_t>(params_sph_2d.get<real_t>("scales.B0"),
@@ -381,16 +381,16 @@ auto main(int argc, char* argv[]) -> int {
         fbc.size(),
         "grid.boundaries.fields.size()");
 
-      // absorb coeffs
+      // match coeffs
       assert_equal<real_t>(
-        params_sph_2d.get<real_t>("grid.boundaries.absorb.ds"),
-        (real_t)(defaults::bc::absorb::ds_frac * 19.0),
-        "grid.boundaries.absorb.ds");
+        params_sph_2d.get<real_t>("grid.boundaries.match.ds"),
+        (real_t)(defaults::bc::match::ds_frac * 19.0),
+        "grid.boundaries.match.ds");
 
       assert_equal<real_t>(
-        params_sph_2d.get<real_t>("grid.boundaries.absorb.coeff"),
+        params_sph_2d.get<real_t>("grid.boundaries.match.coeff"),
         (real_t)10.0,
-        "grid.boundaries.absorb.coeff");
+        "grid.boundaries.match.coeff");
 
       assert_equal(params_sph_2d.get<bool>("particles.use_weights"),
                    true,
@@ -537,16 +537,16 @@ auto main(int argc, char* argv[]) -> int {
         pbc.size(),
         "grid.boundaries.particles.size()");
 
-      // absorb coeffs
+      // match coeffs
       assert_equal<real_t>(
-        params_qks_2d.get<real_t>("grid.boundaries.absorb.ds"),
-        (real_t)(defaults::bc::absorb::ds_frac * (100.0 - 0.8)),
-        "grid.boundaries.absorb.ds");
+        params_qks_2d.get<real_t>("grid.boundaries.match.ds"),
+        (real_t)(defaults::bc::match::ds_frac * (100.0 - 0.8)),
+        "grid.boundaries.match.ds");
 
       assert_equal<real_t>(
-        params_qks_2d.get<real_t>("grid.boundaries.absorb.coeff"),
-        defaults::bc::absorb::coeff,
-        "grid.boundaries.absorb.coeff");
+        params_qks_2d.get<real_t>("grid.boundaries.match.coeff"),
+        defaults::bc::match::coeff,
+        "grid.boundaries.match.coeff");
 
       const auto species = params_qks_2d.get<std::vector<ParticleSpecies>>(
         "particles.species");
diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h
index 9fd40e201..4cde4fca5 100644
--- a/src/global/arch/traits.h
+++ b/src/global/arch/traits.h
@@ -10,7 +10,11 @@
  *   - traits::run_t, traits::to_string_t
  *   - traits::pgen::init_flds_t
  *   - traits::pgen::ext_force_t
- *   - traits::pgen::field_driver_t
+ *   - traits::pgen::atm_fields_t
+ *   - traits::pgen::match_fields_const_t
+ *   - traits::pgen::match_fields_t
+ *   - traits::pgen::fix_fields_const_t
+ *   - traits::pgen::fix_fields_t
  *   - traits::pgen::init_prtls_t
  *   - traits::pgen::custom_fields_t
  *   - traits::pgen::custom_field_output_t
@@ -94,10 +98,19 @@ namespace traits {
     using ext_force_t = decltype(&T::ext_force);
 
     template <typename T>
-    using field_driver_t = decltype(&T::FieldDriver);
+    using atm_fields_t = decltype(&T::AtmFields);
 
     template <typename T>
-    using fix_field_t = decltype(&T::FixField);
+    using match_fields_t = decltype(&T::MatchFields);
+
+    template <typename T>
+    using match_fields_const_t = decltype(&T::MatchFieldsConst);
+
+    template <typename T>
+    using fix_fields_t = decltype(&T::FixFields);
+
+    template <typename T>
+    using fix_fields_const_t = decltype(&T::FixFieldsConst);
 
     template <typename T>
     using custom_fields_t = decltype(&T::CustomFields);
diff --git a/src/global/defaults.h b/src/global/defaults.h
index b7b0107e7..f44fd1844 100644
--- a/src/global/defaults.h
+++ b/src/global/defaults.h
@@ -41,9 +41,13 @@ namespace ntt::defaults {
   } // namespace gr
 
   namespace bc {
-    namespace absorb {
+    namespace match {
       const real_t ds_frac = 0.01;
       const real_t coeff   = 1.0;
+    } // namespace match
+
+    namespace absorb {
+      const real_t ds_frac = 0.01;
     } // namespace absorb
   } // namespace bc
 
diff --git a/src/global/enums.h b/src/global/enums.h
index 283cb456d..8f2495c13 100644
--- a/src/global/enums.h
+++ b/src/global/enums.h
@@ -8,7 +8,7 @@
  *   - enum ntt::SimEngine         // SRPIC, GRPIC
  *   - enum ntt::PrtlBC            // periodic, absorb, atmosphere, custom,
  *                                    reflect, horizon, axis, sync
- *   - enum ntt::FldsBC            // periodic, absorb, fixed, atmosphere,
+ *   - enum ntt::FldsBC            // periodic, match, fixed, atmosphere,
  *                                    custom, horizon, axis, sync
  *   - enum ntt::PrtlPusher        // boris, vay, photon, none
  *   - enum ntt::Cooling           // synchrotron, none
@@ -215,7 +215,7 @@ namespace ntt {
     enum type : uint8_t {
       INVALID    = 0,
       PERIODIC   = 1,
-      ABSORB     = 2,
+      MATCH      = 2,
       FIXED      = 3,
       ATMOSPHERE = 4,
       CUSTOM     = 5,
@@ -226,9 +226,9 @@ namespace ntt {
 
     constexpr FldsBC(uint8_t c) : enums_hidden::BaseEnum<FldsBC> { c } {}
 
-    static constexpr type variants[] = { PERIODIC, ABSORB,  FIXED, ATMOSPHERE,
+    static constexpr type variants[] = { PERIODIC, MATCH,   FIXED, ATMOSPHERE,
                                          CUSTOM,   HORIZON, AXIS,  SYNC };
-    static constexpr const char* lookup[] = { "periodic",   "absorb", "fixed",
+    static constexpr const char* lookup[] = { "periodic",   "match",  "fixed",
                                               "atmosphere", "custom", "horizon",
                                               "axis",       "sync" };
     static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]);
diff --git a/src/global/global.h b/src/global/global.h
index 6669981b6..77fa8c51c 100644
--- a/src/global/global.h
+++ b/src/global/global.h
@@ -13,6 +13,10 @@
  *   - enum PrepareOutput
  *   - enum CellLayer           // allLayer, activeLayer, minGhostLayer,
  *                                 minActiveLayer, maxActiveLayer, maxGhostLayer
+ *   - enum Idx                // U, D, T, XYZ, Sph, PU, PD
+ *   - enum Crd                // Cd, Ph, XYZ, Sph
+ *   - enum in                 // x1, x2, x3
+ *   - enum bc_in                // Px1, Mx1, Px2, Mx2, Px3, Mx3
  *   - type box_region_t
  *   - files::LogFile, files::ErrFile, files::InfoFile
  *   - type prtldx_t
@@ -184,6 +188,15 @@ enum class in : unsigned short {
   x3 = 2,
 };
 
+enum class bc_in : short {
+  Mx1 = -1,
+  Px1 = 1,
+  Mx2 = -2,
+  Px2 = 2,
+  Mx3 = -3,
+  Px3 = 3,
+};
+
 template <Dimension D>
 using box_region_t = CellLayer[D];
 
diff --git a/src/global/tests/enums.cpp b/src/global/tests/enums.cpp
index 4d678e85e..7785ec1a3 100644
--- a/src/global/tests/enums.cpp
+++ b/src/global/tests/enums.cpp
@@ -61,7 +61,7 @@ auto main() -> int {
   enum_str_t all_simulation_engines = { "srpic", "grpic" };
   enum_str_t all_particle_bcs = { "periodic", "absorb",  "atmosphere", "custom",
                                   "reflect",  "horizon", "axis",       "sync" };
-  enum_str_t all_fields_bcs   = { "periodic", "absorb",  "fixed", "atmosphere",
+  enum_str_t all_fields_bcs   = { "periodic", "match",   "fixed", "atmosphere",
                                   "custom",   "horizon", "axis",  "sync" };
   enum_str_t all_particle_pushers = { "boris", "vay", "photon", "none" };
   enum_str_t all_coolings         = { "synchrotron", "none" };
diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp
index 2f2a458bb..363ff3ad2 100644
--- a/src/kernels/fields_bcs.hpp
+++ b/src/kernels/fields_bcs.hpp
@@ -1,5 +1,12 @@
 /**
- * @brief: kernels/fields_bcs.hpp
+ * @file kernels/fields_bcs.hpp
+ * @brief Kernels used for field boundary conditions
+ * @implements
+ *   - kernel::bc::MatchBoundaries_kernel<>
+ *   - kernel::bc::AxisBoundaries_kernel<>
+ *   - kernel::bc::EnforcedBoundaries_kernel<>
+ * @namespaces:
+ *   - kernel::bc::
  */
 
 #ifndef KERNELS_FIELDS_BCS_HPP
@@ -12,64 +19,162 @@
 #include "utils/error.h"
 #include "utils/numeric.h"
 
-namespace kernel {
+namespace kernel::bc {
   using namespace ntt;
 
-  template <class M, idx_t i>
-  struct AbsorbBoundaries_kernel {
+  /*
+   * @tparam S: Simulation Engine
+   * @tparam I: Field Setter class
+   * @tparam M: Metric
+   * @tparam o: Orientation
+   *
+   * @brief Applies matching boundary conditions (with a smooth profile) in a specific direction.
+   * @note If a component is not specified in the field setter, it is ignored.
+   * @note It is supposed to only be called on the active side of the absorbing edge (so sign is not needed).
+   */
+  template <SimEngine::type S, class I, class M, in o>
+  struct MatchBoundaries_kernel {
     static_assert(M::is_metric, "M must be a metric class");
-    static_assert(i <= static_cast<unsigned short>(M::Dim),
+    static_assert(static_cast<unsigned short>(o) <
+                    static_cast<unsigned short>(M::Dim),
                   "Invalid component index");
+    static constexpr idx_t i = static_cast<idx_t>(o) + 1u;
+    static constexpr bool defines_dx1 = traits::has_method<traits::dx1_t, I>::value;
+    static constexpr bool defines_dx2 = traits::has_method<traits::dx2_t, I>::value;
+    static constexpr bool defines_dx3 = traits::has_method<traits::dx3_t, I>::value;
+    static constexpr bool defines_ex1 = traits::has_method<traits::ex1_t, I>::value;
+    static constexpr bool defines_ex2 = traits::has_method<traits::ex2_t, I>::value;
+    static constexpr bool defines_ex3 = traits::has_method<traits::ex3_t, I>::value;
+    static constexpr bool defines_bx1 = traits::has_method<traits::bx1_t, I>::value;
+    static constexpr bool defines_bx2 = traits::has_method<traits::bx2_t, I>::value;
+    static constexpr bool defines_bx3 = traits::has_method<traits::bx3_t, I>::value;
+    static_assert(
+      (S == SimEngine::SRPIC and (defines_ex1 or defines_ex2 or defines_ex3 or
+                                  defines_bx1 or defines_bx2 or defines_bx3)) or
+        ((S == SimEngine::GRPIC) and (defines_dx1 or defines_dx2 or defines_dx3 or
+                                      defines_bx1 or defines_bx2 or defines_bx3)),
+      "none of the components of E/D or B are specified in PGEN");
 
     ndfield_t<M::Dim, 6> Fld;
+    const I              fset;
     const M              metric;
     const real_t         xg_edge;
     const real_t         dx_abs;
     const BCTags         tags;
 
-    AbsorbBoundaries_kernel(ndfield_t<M::Dim, 6> Fld,
-                            const M&             metric,
-                            real_t               xg_edge,
-                            real_t               dx_abs,
-                            BCTags               tags)
+    MatchBoundaries_kernel(ndfield_t<M::Dim, 6> Fld,
+                           const I&             fset,
+                           const M&             metric,
+                           real_t               xg_edge,
+                           real_t               dx_abs,
+                           BCTags               tags)
       : Fld { Fld }
+      , fset { fset }
       , metric { metric }
       , xg_edge { xg_edge }
       , dx_abs { dx_abs }
       , tags { tags } {}
 
+    Inline auto shape(const real_t& dx) const -> real_t {
+      return math::tanh(dx * FOUR / dx_abs);
+    }
+
     Inline void operator()(index_t i1) const {
       if constexpr (M::Dim == Dim::_1D) {
         const auto i1_ = COORD(i1);
-        for (const auto comp :
-             { em::ex1, em::ex2, em::ex3, em::bx1, em::bx2, em::bx3 }) {
-          if ((comp == em::ex1) and not(tags & BC::Ex1)) {
-            continue;
-          } else if ((comp == em::ex2) and not(tags & BC::Ex2)) {
-            continue;
-          } else if ((comp == em::ex3) and not(tags & BC::Ex3)) {
-            continue;
-          } else if ((comp == em::bx1) and not(tags & BC::Bx1)) {
-            continue;
-          } else if ((comp == em::bx2) and not(tags & BC::Bx2)) {
-            continue;
-          } else if ((comp == em::bx3) and not(tags & BC::Bx3)) {
-            continue;
-          }
-          coord_t<M::Dim> x_Cd { ZERO };
-          if (comp == em::ex1 or comp == em::bx2 or comp == em::bx3) {
-            x_Cd[0] = i1_ + HALF;
-          } else if (comp == em::ex2 or comp == em::bx1 or comp == em::ex3) {
-            x_Cd[0] = i1_;
-          }
-          const auto dx = math::abs(
-            metric.template convert<i, Crd::Cd, Crd::Ph>(x_Cd[i - 1]) - xg_edge);
-          Fld(i1, comp) *= math::tanh(dx / (INV_4 * dx_abs));
+
+        if constexpr (S == SimEngine::SRPIC) {
+          coord_t<Dim::_1D> x_Ph_0 { ZERO };
+          coord_t<Dim::_1D> x_Ph_H { ZERO };
+          metric.template convert<Crd::Cd, Crd::Ph>({ i1_ }, x_Ph_0);
+          metric.template convert<Crd::Cd, Crd::Ph>({ i1_ + HALF }, x_Ph_H);
+
+          // SRPIC
+          auto ex1_U { ZERO }, ex2_U { ZERO }, ex3_U { ZERO }, bx1_U { ZERO },
+            bx2_U { ZERO }, bx3_U { ZERO };
+          if (tags & BC::E) {
+            if constexpr (defines_ex1) {
+              ex1_U = metric.template transform<1, Idx::T, Idx::U>(
+                { i1_ + HALF },
+                fset.ex1(x_Ph_H));
+            }
+            if constexpr (defines_ex2) {
+              ex2_U = metric.template transform<2, Idx::T, Idx::U>(
+                { i1_ },
+                fset.ex2(x_Ph_0));
+            }
+            if constexpr (defines_ex3) {
+              ex3_U = metric.template transform<3, Idx::T, Idx::U>(
+                { i1_ },
+                fset.ex3(x_Ph_0));
+            }
+          }
+          if (tags & BC::B) {
+            if constexpr (defines_bx1) {
+              bx1_U = metric.template transform<1, Idx::T, Idx::U>(
+                { i1_ },
+                fset.bx1(x_Ph_0));
+            }
+            if constexpr (defines_bx2) {
+              bx2_U = metric.template transform<2, Idx::T, Idx::U>(
+                { i1_ + HALF },
+                fset.bx2(x_Ph_H));
+            }
+            if constexpr (defines_bx3) {
+              bx3_U = metric.template transform<3, Idx::T, Idx::U>(
+                { i1_ + HALF },
+                fset.bx3(x_Ph_H));
+            }
+          }
+
+          if constexpr (defines_ex1 or defines_bx2 or defines_bx3) {
+            const auto dx = math::abs(
+              metric.template convert<i, Crd::Cd, Crd::Ph>(i1_ + HALF) - xg_edge);
+            const auto s = shape(dx);
+            if constexpr (defines_ex1) {
+              if (tags & BC::E) {
+                Fld(i1, em::ex1) = s * Fld(i1, em::ex1) + (ONE - s) * ex1_U;
+              }
+            }
+            if constexpr (defines_bx2 or defines_bx3) {
+              if (tags & BC::B) {
+                if constexpr (defines_bx2) {
+                  Fld(i1, em::bx2) = s * Fld(i1, em::bx2) + (ONE - s) * bx2_U;
+                }
+                if constexpr (defines_bx3) {
+                  Fld(i1, em::bx3) = s * Fld(i1, em::bx3) + (ONE - s) * bx3_U;
+                }
+              }
+            }
+          }
+          if constexpr (defines_bx1 or defines_ex2 or defines_ex3) {
+            const auto dx = math::abs(
+              metric.template convert<i, Crd::Cd, Crd::Ph>(i1_) - xg_edge);
+            const auto s = shape(dx);
+            if constexpr (defines_bx1) {
+              if (tags & BC::B) {
+                Fld(i1, em::bx1) = s * Fld(i1, em::bx1) + (ONE - s) * bx1_U;
+              }
+            }
+            if constexpr (defines_ex2 or defines_ex3) {
+              if (tags & BC::E) {
+                if constexpr (defines_ex2) {
+                  Fld(i1, em::ex2) = s * Fld(i1, em::ex2) + (ONE - s) * ex2_U;
+                }
+                if constexpr (defines_ex3) {
+                  Fld(i1, em::ex3) = s * Fld(i1, em::ex3) + (ONE - s) * ex3_U;
+                }
+              }
+            }
+          }
+        } else {
+          // GRPIC
+          raise::KernelError(HERE, "1D GRPIC not implemented");
         }
       } else {
         raise::KernelError(
           HERE,
-          "AbsorbFields_kernel: 1D implementation called for D != 1");
+          "MatchBoundaries_kernel: 1D implementation called for D != 1");
       }
     }
 
@@ -77,43 +182,129 @@ namespace kernel {
       if constexpr (M::Dim == Dim::_2D) {
         const auto i1_ = COORD(i1);
         const auto i2_ = COORD(i2);
-        for (const auto comp :
-             { em::ex1, em::ex2, em::ex3, em::bx1, em::bx2, em::bx3 }) {
-          if ((comp == em::ex1) and not(tags & BC::Ex1)) {
-            continue;
-          } else if ((comp == em::ex2) and not(tags & BC::Ex2)) {
-            continue;
-          } else if ((comp == em::ex3) and not(tags & BC::Ex3)) {
-            continue;
-          } else if ((comp == em::bx1) and not(tags & BC::Bx1)) {
-            continue;
-          } else if ((comp == em::bx2) and not(tags & BC::Bx2)) {
-            continue;
-          } else if ((comp == em::bx3) and not(tags & BC::Bx3)) {
-            continue;
-          }
-          coord_t<M::Dim> x_Cd { ZERO };
-          if (comp == em::ex1 or comp == em::bx2) {
-            x_Cd[0] = i1_ + HALF;
-            x_Cd[1] = i2_;
-          } else if (comp == em::ex2 or comp == em::bx1) {
-            x_Cd[0] = i1_;
-            x_Cd[1] = i2_ + HALF;
-          } else if (comp == em::ex3) {
-            x_Cd[0] = i1_;
-            x_Cd[1] = i2_;
-          } else if (comp == em::bx3) {
-            x_Cd[0] = i1_ + HALF;
-            x_Cd[1] = i2_ + HALF;
-          }
-          const auto dx = math::abs(
-            metric.template convert<i, Crd::Cd, Crd::Ph>(x_Cd[i - 1]) - xg_edge);
-          Fld(i1, i2, comp) *= math::tanh(dx / (INV_4 * dx_abs));
+
+        if constexpr (S == SimEngine::SRPIC) {
+          // SRPIC
+          if constexpr (defines_ex1 or defines_bx2) {
+            coord_t<Dim::_2D> x_Ph_H0 { ZERO };
+            metric.template convert<Crd::Cd, Crd::Ph>({ i1_ + HALF, i2_ }, x_Ph_H0);
+            // i1 + 1/2, i2
+            real_t xi_Cd;
+            if constexpr (o == in::x1) {
+              xi_Cd = i1_ + HALF;
+            } else {
+              xi_Cd = i2_;
+            }
+
+            const auto dx = math::abs(
+              metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+            const auto s = shape(dx);
+
+            if constexpr (defines_ex1) {
+              if (tags & BC::E) {
+                const auto ex1_U = metric.template transform<1, Idx::T, Idx::U>(
+                  { i1_ + HALF, i2_ },
+                  fset.ex1(x_Ph_H0));
+                Fld(i1, i2, em::ex1) = s * Fld(i1, i2, em::ex1) + (ONE - s) * ex1_U;
+              }
+            }
+            if constexpr (defines_bx2) {
+              if (tags & BC::B) {
+                const auto bx2_U = metric.template transform<2, Idx::T, Idx::U>(
+                  { i1_ + HALF, i2_ },
+                  fset.bx2(x_Ph_H0));
+                Fld(i1, i2, em::bx2) = s * Fld(i1, i2, em::bx2) + (ONE - s) * bx2_U;
+              }
+            }
+          }
+
+          if constexpr (defines_ex2 or defines_bx1) {
+            coord_t<Dim::_2D> x_Ph_0H { ZERO };
+            metric.template convert<Crd::Cd, Crd::Ph>({ i1_, i2_ + HALF }, x_Ph_0H);
+            // i1, i2 + 1/2
+            real_t xi_Cd;
+            if constexpr (o == in::x1) {
+              xi_Cd = i1_;
+            } else {
+              xi_Cd = i2_ + HALF;
+            }
+
+            const auto dx = math::abs(
+              metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+            const auto s = shape(dx);
+            if constexpr (defines_ex2) {
+              if (tags & BC::E) {
+                auto ex2_U { ZERO };
+                ex2_U = metric.template transform<2, Idx::T, Idx::U>(
+                  { i1_, i2_ + HALF },
+                  fset.ex2(x_Ph_0H));
+                Fld(i1, i2, em::ex2) = s * Fld(i1, i2, em::ex2) + (ONE - s) * ex2_U;
+              }
+            }
+            if constexpr (defines_bx1) {
+              if (tags & BC::B) {
+                auto bx1_U { ZERO };
+                bx1_U = metric.template transform<1, Idx::T, Idx::U>(
+                  { i1_, i2_ + HALF },
+                  fset.bx1(x_Ph_0H));
+                Fld(i1, i2, em::bx1) = s * Fld(i1, i2, em::bx1) + (ONE - s) * bx1_U;
+              }
+            }
+          }
+
+          if constexpr (defines_ex3) {
+            if (tags & BC::E) {
+              auto              ex3_U { ZERO };
+              coord_t<Dim::_2D> x_Ph_00 { ZERO };
+              metric.template convert<Crd::Cd, Crd::Ph>({ i1_, i2_ }, x_Ph_00);
+              ex3_U = metric.template transform<3, Idx::T, Idx::U>(
+                { i1_, i2_ },
+                fset.ex3(x_Ph_00));
+              // i1, i2
+              real_t xi_Cd;
+              if constexpr (o == in::x1) {
+                xi_Cd = i1_;
+              } else {
+                xi_Cd = i2_;
+              }
+              const auto dx = math::abs(
+                metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+              const auto s = shape(dx);
+              Fld(i1, i2, em::ex3) = s * Fld(i1, i2, em::ex3) + (ONE - s) * ex3_U;
+            }
+          }
+
+          if constexpr (defines_bx3) {
+            if (tags & BC::B) {
+              auto              bx3_U { ZERO };
+              coord_t<Dim::_2D> x_Ph_HH { ZERO };
+              metric.template convert<Crd::Cd, Crd::Ph>({ i1_ + HALF, i2_ + HALF },
+                                                        x_Ph_HH);
+              bx3_U = metric.template transform<3, Idx::T, Idx::U>(
+                { i1_ + HALF, i2_ + HALF },
+                fset.bx3(x_Ph_HH));
+              // i1 + 1/2, i2 + 1/2
+              real_t xi_Cd;
+              if constexpr (o == in::x1) {
+                xi_Cd = i1_ + HALF;
+              } else {
+                xi_Cd = i2_ + HALF;
+              }
+              const auto dx = math::abs(
+                metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+              const auto s = shape(dx);
+              // bx3
+              Fld(i1, i2, em::bx3) = s * Fld(i1, i2, em::bx3) + (ONE - s) * bx3_U;
+            }
+          }
+        } else {
+          // GRPIC
+          raise::KernelError(HERE, "GRPIC not implemented");
         }
       } else {
         raise::KernelError(
           HERE,
-          "AbsorbFields_kernel: 2D implementation called for D != 2");
+          "MatchBoundaries_kernel: 2D implementation called for D != 2");
       }
     }
 
@@ -122,59 +313,184 @@ namespace kernel {
         const auto i1_ = COORD(i1);
         const auto i2_ = COORD(i2);
         const auto i3_ = COORD(i3);
-        for (const auto comp :
-             { em::ex1, em::ex2, em::ex3, em::bx1, em::bx2, em::bx3 }) {
-          if ((comp == em::ex1) and not(tags & BC::Ex1)) {
-            continue;
-          } else if ((comp == em::ex2) and not(tags & BC::Ex2)) {
-            continue;
-          } else if ((comp == em::ex3) and not(tags & BC::Ex3)) {
-            continue;
-          } else if ((comp == em::bx1) and not(tags & BC::Bx1)) {
-            continue;
-          } else if ((comp == em::bx2) and not(tags & BC::Bx2)) {
-            continue;
-          } else if ((comp == em::bx3) and not(tags & BC::Bx3)) {
-            continue;
-          }
-          coord_t<M::Dim> x_Cd { ZERO };
-          if (comp == em::ex1) {
-            x_Cd[0] = i1_ + HALF;
-            x_Cd[1] = i2_;
-            x_Cd[2] = i3_;
-          } else if (comp == em::ex2) {
-            x_Cd[0] = i1_;
-            x_Cd[1] = i2_ + HALF;
-            x_Cd[2] = i3_;
-          } else if (comp == em::ex3) {
-            x_Cd[0] = i1_;
-            x_Cd[1] = i2_;
-            x_Cd[2] = i3_ + HALF;
-          } else if (comp == em::bx1) {
-            x_Cd[0] = i1_;
-            x_Cd[1] = i2_ + HALF;
-            x_Cd[2] = i3_ + HALF;
-          } else if (comp == em::bx2) {
-            x_Cd[0] = i1_ + HALF;
-            x_Cd[1] = i2_;
-            x_Cd[2] = i3_ + HALF;
-          } else if (comp == em::bx3) {
-            x_Cd[0] = i1_ + HALF;
-            x_Cd[1] = i2_ + HALF;
-            x_Cd[2] = i3_;
-          }
-          const auto dx = math::abs(
-            metric.template convert<i, Crd::Cd, Crd::Ph>(x_Cd[i - 1]) - xg_edge);
-          Fld(i1, i2, i3, comp) *= math::tanh(dx / (INV_4 * dx_abs));
+
+        if constexpr (S == SimEngine::SRPIC) {
+          // SRPIC
+          if constexpr (defines_ex1 or defines_ex2 or defines_ex3) {
+            if (tags & BC::E) {
+              if constexpr (defines_ex1) {
+                // i1 + 1/2, i2, i3
+                real_t xi_Cd;
+                if constexpr (o == in::x1) {
+                  xi_Cd = i1_ + HALF;
+                } else if constexpr (o == in::x2) {
+                  xi_Cd = i2_;
+                } else {
+                  xi_Cd = i3_;
+                }
+                const auto dx = math::abs(
+                  metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+                const auto        s = shape(dx);
+                auto              ex1_U { ZERO };
+                coord_t<Dim::_3D> x_Ph_H00 { ZERO };
+                metric.template convert<Crd::Cd, Crd::Ph>({ i1_ + HALF, i2_, i3_ },
+                                                          x_Ph_H00);
+                ex1_U = metric.template transform<1, Idx::T, Idx::U>(
+                  { i1_ + HALF, i2_, i3_ },
+                  fset.ex1(x_Ph_H00));
+                Fld(i1, i2, i3, em::ex1) = s * Fld(i1, i2, i3, em::ex1) +
+                                           (ONE - s) * ex1_U;
+              }
+
+              if constexpr (defines_ex2) {
+                // i1, i2 + 1/2, i3
+                real_t xi_Cd;
+                if constexpr (o == in::x1) {
+                  xi_Cd = i1_;
+                } else if constexpr (o == in::x2) {
+                  xi_Cd = i2_ + HALF;
+                } else {
+                  xi_Cd = i3_;
+                }
+                const auto dx = math::abs(
+                  metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+                const auto        s = shape(dx);
+                auto              ex2_U { ZERO };
+                coord_t<Dim::_3D> x_Ph_0H0 { ZERO };
+                metric.template convert<Crd::Cd, Crd::Ph>({ i1_, i2_ + HALF, i3_ },
+                                                          x_Ph_0H0);
+                ex2_U = metric.template transform<2, Idx::T, Idx::U>(
+                  { i1_, i2_ + HALF, i3_ },
+                  fset.ex2(x_Ph_0H0));
+                Fld(i1, i2, i3, em::ex2) = s * Fld(i1, i2, i3, em::ex2) +
+                                           (ONE - s) * ex2_U;
+              }
+
+              if constexpr (defines_ex3) {
+                // i1, i2, i3 + 1/2
+                real_t xi_Cd;
+                if constexpr (o == in::x1) {
+                  xi_Cd = i1_;
+                } else if constexpr (o == in::x2) {
+                  xi_Cd = i2_;
+                } else {
+                  xi_Cd = i3_ + HALF;
+                }
+                const auto dx = math::abs(
+                  metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+                const auto        s = shape(dx);
+                auto              ex3_U { ZERO };
+                coord_t<Dim::_3D> x_Ph_00H { ZERO };
+                metric.template convert<Crd::Cd, Crd::Ph>({ i1_, i2_, i3_ + HALF },
+                                                          x_Ph_00H);
+                ex3_U = metric.template transform<3, Idx::T, Idx::U>(
+                  { i1_, i2_, i3_ + HALF },
+                  fset.ex3(x_Ph_00H));
+                Fld(i1, i2, i3, em::ex3) = s * Fld(i1, i2, i3, em::ex3) +
+                                           (ONE - s) * ex3_U;
+              }
+            }
+          }
+
+          if constexpr (defines_bx1 or defines_bx2 or defines_bx3) {
+            if (tags & BC::B) {
+              if constexpr (defines_bx1) {
+                // i1, i2 + 1/2, i3 + 1/2
+                real_t xi_Cd;
+                if constexpr (o == in::x1) {
+                  xi_Cd = i1_;
+                } else if constexpr (o == in::x2) {
+                  xi_Cd = i2_ + HALF;
+                } else {
+                  xi_Cd = i3_ + HALF;
+                }
+                const auto dx = math::abs(
+                  metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+                const auto s = shape(dx);
+                auto       bx1_U { ZERO };
+                if constexpr (defines_bx1) {
+                  coord_t<Dim::_3D> x_Ph_0HH { ZERO };
+                  metric.template convert<Crd::Cd, Crd::Ph>(
+                    { i1_, i2_ + HALF, i3_ + HALF },
+                    x_Ph_0HH);
+                  bx1_U = metric.template transform<1, Idx::T, Idx::U>(
+                    { i1_, i2_ + HALF, i3_ + HALF },
+                    fset.bx1(x_Ph_0HH));
+                }
+                // bx1
+                Fld(i1, i2, i3, em::bx1) = s * Fld(i1, i2, i3, em::bx1) +
+                                           (ONE - s) * bx1_U;
+              }
+
+              if constexpr (defines_bx2) {
+                // i1 + 1/2, i2, i3 + 1/2
+                real_t xi_Cd;
+                if constexpr (o == in::x1) {
+                  xi_Cd = i1_ + HALF;
+                } else if constexpr (o == in::x2) {
+                  xi_Cd = i2_;
+                } else {
+                  xi_Cd = i3_ + HALF;
+                }
+                const auto dx = math::abs(
+                  metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+                const auto        s = shape(dx);
+                auto              bx2_U { ZERO };
+                coord_t<Dim::_3D> x_Ph_H0H { ZERO };
+                metric.template convert<Crd::Cd, Crd::Ph>(
+                  { i1_ + HALF, i2_, i3_ + HALF },
+                  x_Ph_H0H);
+                bx2_U = metric.template transform<2, Idx::T, Idx::U>(
+                  { i1_ + HALF, i2_, i3_ + HALF },
+                  fset.bx2(x_Ph_H0H));
+                Fld(i1, i2, i3, em::bx2) = s * Fld(i1, i2, i3, em::bx2) +
+                                           (ONE - s) * bx2_U;
+              }
+
+              if constexpr (defines_bx3) {
+                // i1 + 1/2, i2 + 1/2, i3
+                real_t xi_Cd;
+                if constexpr (o == in::x1) {
+                  xi_Cd = i1_ + HALF;
+                } else if constexpr (o == in::x2) {
+                  xi_Cd = i2_ + HALF;
+                } else {
+                  xi_Cd = i3_;
+                }
+                const auto dx = math::abs(
+                  metric.template convert<i, Crd::Cd, Crd::Ph>(xi_Cd) - xg_edge);
+                const auto        s = shape(dx);
+                auto              bx3_U { ZERO };
+                coord_t<Dim::_3D> x_Ph_HH0 { ZERO };
+                metric.template convert<Crd::Cd, Crd::Ph>(
+                  { i1_ + HALF, i2_ + HALF, i3_ },
+                  x_Ph_HH0);
+                bx3_U = metric.template transform<3, Idx::T, Idx::U>(
+                  { i1_ + HALF, i2_ + HALF, i3_ },
+                  fset.bx3(x_Ph_HH0));
+                Fld(i1, i2, i3, em::bx3) = s * Fld(i1, i2, i3, em::bx3) +
+                                           (ONE - s) * bx3_U;
+              }
+            }
+          }
+        } else {
+          // GRPIC
+          raise::KernelError(HERE, "GRPIC not implemented");
         }
       } else {
         raise::KernelError(
           HERE,
-          "AbsorbFields_kernel: 3D implementation called for D != 3");
+          "MatchBoundaries_kernel: 3D implementation called for D != 3");
       }
     }
   };
 
+  /*
+   * @tparam D: Dimension
+   * @tparam P: Positive/Negative direction
+   *
+   * @brief Applies boundary conditions near the polar axis
+   */
   template <Dimension D, bool P>
   struct AxisBoundaries_kernel {
     ndfield_t<D, 6>   Fld;
@@ -216,6 +532,14 @@ namespace kernel {
     }
   };
 
+  /*
+   * @tparam I: Field Setter class
+   * @tparam M: Metric
+   * @tparam P: Positive/Negative direction
+   * @tparam O: Orientation
+   *
+   * @brief Applies enforced boundary conditions (fixed value)
+   */
   template <class I, class M, bool P, in O>
   struct EnforcedBoundaries_kernel {
     static constexpr Dimension D = M::Dim;
@@ -226,31 +550,28 @@ namespace kernel {
     static constexpr bool defines_bx2 = traits::has_method<traits::bx2_t, I>::value;
     static constexpr bool defines_bx3 = traits::has_method<traits::bx3_t, I>::value;
 
-    static_assert(defines_ex1 and defines_ex2 and defines_ex3 and
-                    defines_bx1 and defines_bx2 and defines_bx3,
-                  "not all components of E or B are specified in PGEN");
+    static_assert(defines_ex1 or defines_ex2 or defines_ex3 or defines_bx1 or
+                    defines_bx2 or defines_bx3,
+                  "none of the components of E or B are specified in PGEN");
     static_assert(M::is_metric, "M must be a metric class");
     static_assert(static_cast<unsigned short>(O) <
                     static_cast<unsigned short>(M::Dim),
                   "Invalid Orientation");
 
     ndfield_t<D, 6>   Fld;
-    const I           finit;
+    const I           fset;
     const M           metric;
     const std::size_t i_edge;
-    const bool        setE, setB;
 
     EnforcedBoundaries_kernel(ndfield_t<M::Dim, 6>& Fld,
-                              const I&              finit,
+                              const I&              fset,
                               const M&              metric,
                               std::size_t           i_edge,
                               BCTags                tags)
       : Fld { Fld }
-      , finit { finit }
+      , fset { fset }
       , metric { metric }
-      , i_edge { i_edge + N_GHOSTS }
-      , setE { tags & BC::Ex1 or tags & BC::Ex2 or tags & BC::Ex3 }
-      , setB { tags & BC::Bx1 or tags & BC::Bx2 or tags & BC::Bx3 } {}
+      , i_edge { i_edge + N_GHOSTS } {}
 
     Inline void operator()(index_t i1) const {
       if constexpr (D == Dim::_1D) {
@@ -259,8 +580,8 @@ namespace kernel {
         coord_t<Dim::_1D> x_Ph_H { ZERO };
         metric.template convert<Crd::Cd, Crd::Ph>({ i1_ }, x_Ph_0);
         metric.template convert<Crd::Cd, Crd::Ph>({ i1_ + HALF }, x_Ph_H);
-        bool setEx1 = setE, setEx2 = setE, setEx3 = setE, setBx1 = setB,
-             setBx2 = setB, setBx3 = setB;
+        bool setEx1 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3,
+             setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3;
         if constexpr (O == in::x1) {
           // x1 -- normal
           // x2,x3 -- tangential
@@ -276,35 +597,47 @@ namespace kernel {
         } else {
           raise::KernelError(HERE, "Invalid Orientation");
         }
-        if (setEx1) {
-          Fld(i1, em::ex1) = metric.template transform<1, Idx::T, Idx::U>(
-            { i1_ + HALF },
-            finit.ex1(x_Ph_H));
+        if constexpr (defines_ex1) {
+          if (setEx1) {
+            Fld(i1, em::ex1) = metric.template transform<1, Idx::T, Idx::U>(
+              { i1_ + HALF },
+              fset.ex1(x_Ph_H));
+          }
         }
-        if (setEx2) {
-          Fld(i1, em::ex2) = metric.template transform<2, Idx::T, Idx::U>(
-            { i1_ },
-            finit.ex2(x_Ph_0));
+        if constexpr (defines_ex2) {
+          if (setEx2) {
+            Fld(i1, em::ex2) = metric.template transform<2, Idx::T, Idx::U>(
+              { i1_ },
+              fset.ex2(x_Ph_0));
+          }
         }
-        if (setEx3) {
-          Fld(i1, em::ex3) = metric.template transform<3, Idx::T, Idx::U>(
-            { i1_ },
-            finit.ex3(x_Ph_0));
+        if constexpr (defines_ex3) {
+          if (setEx3) {
+            Fld(i1, em::ex3) = metric.template transform<3, Idx::T, Idx::U>(
+              { i1_ },
+              fset.ex3(x_Ph_0));
+          }
         }
-        if (setBx1) {
-          Fld(i1, em::bx1) = metric.template transform<1, Idx::T, Idx::U>(
-            { i1_ },
-            finit.bx1(x_Ph_0));
+        if constexpr (defines_bx1) {
+          if (setBx1) {
+            Fld(i1, em::bx1) = metric.template transform<1, Idx::T, Idx::U>(
+              { i1_ },
+              fset.bx1(x_Ph_0));
+          }
         }
-        if (setBx2) {
-          Fld(i1, em::bx2) = metric.template transform<2, Idx::T, Idx::U>(
-            { i1_ + HALF },
-            finit.bx2(x_Ph_H));
+        if constexpr (defines_bx2) {
+          if (setBx2) {
+            Fld(i1, em::bx2) = metric.template transform<2, Idx::T, Idx::U>(
+              { i1_ + HALF },
+              fset.bx2(x_Ph_H));
+          }
         }
-        if (setBx3) {
-          Fld(i1, em::bx3) = metric.template transform<3, Idx::T, Idx::U>(
-            { i1_ + HALF },
-            finit.bx3(x_Ph_H));
+        if constexpr (defines_bx3) {
+          if (setBx3) {
+            Fld(i1, em::bx3) = metric.template transform<3, Idx::T, Idx::U>(
+              { i1_ + HALF },
+              fset.bx3(x_Ph_H));
+          }
         }
       } else {
         raise::KernelError(HERE, "Invalid Dimension");
@@ -324,8 +657,8 @@ namespace kernel {
         metric.template convert<Crd::Cd, Crd::Ph>({ i1_ + HALF, i2_ }, x_Ph_H0);
         metric.template convert<Crd::Cd, Crd::Ph>({ i1_ + HALF, i2_ + HALF },
                                                   x_Ph_HH);
-        bool setEx1 = setE, setEx2 = setE, setEx3 = setE, setBx1 = setB,
-             setBx2 = setB, setBx3 = setB;
+        bool setEx1 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3,
+             setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3;
         if constexpr (O == in::x1) {
           // x1 -- normal
           // x2,x3 -- tangential
@@ -353,35 +686,47 @@ namespace kernel {
         } else {
           raise::KernelError(HERE, "Invalid Orientation");
         }
-        if (setEx1) {
-          Fld(i1, i2, em::ex1) = metric.template transform<1, Idx::T, Idx::U>(
-            { i1_ + HALF, i2_ },
-            finit.ex1(x_Ph_H0));
+        if constexpr (defines_ex1) {
+          if (setEx1) {
+            Fld(i1, i2, em::ex1) = metric.template transform<1, Idx::T, Idx::U>(
+              { i1_ + HALF, i2_ },
+              fset.ex1(x_Ph_H0));
+          }
         }
-        if (setEx2) {
-          Fld(i1, i2, em::ex2) = metric.template transform<2, Idx::T, Idx::U>(
-            { i1_, i2_ + HALF },
-            finit.ex2(x_Ph_0H));
+        if constexpr (defines_ex2) {
+          if (setEx2) {
+            Fld(i1, i2, em::ex2) = metric.template transform<2, Idx::T, Idx::U>(
+              { i1_, i2_ + HALF },
+              fset.ex2(x_Ph_0H));
+          }
         }
-        if (setEx3) {
-          Fld(i1, i2, em::ex3) = metric.template transform<3, Idx::T, Idx::U>(
-            { i1_, i2_ },
-            finit.ex3(x_Ph_00));
+        if constexpr (defines_ex3) {
+          if (setEx3) {
+            Fld(i1, i2, em::ex3) = metric.template transform<3, Idx::T, Idx::U>(
+              { i1_, i2_ },
+              fset.ex3(x_Ph_00));
+          }
         }
-        if (setBx1) {
-          Fld(i1, i2, em::bx1) = metric.template transform<1, Idx::T, Idx::U>(
-            { i1_, i2_ + HALF },
-            finit.bx1(x_Ph_0H));
+        if constexpr (defines_bx1) {
+          if (setBx1) {
+            Fld(i1, i2, em::bx1) = metric.template transform<1, Idx::T, Idx::U>(
+              { i1_, i2_ + HALF },
+              fset.bx1(x_Ph_0H));
+          }
         }
-        if (setBx2) {
-          Fld(i1, i2, em::bx2) = metric.template transform<2, Idx::T, Idx::U>(
-            { i1_ + HALF, i2_ },
-            finit.bx2(x_Ph_H0));
+        if constexpr (defines_bx2) {
+          if (setBx2) {
+            Fld(i1, i2, em::bx2) = metric.template transform<2, Idx::T, Idx::U>(
+              { i1_ + HALF, i2_ },
+              fset.bx2(x_Ph_H0));
+          }
         }
-        if (setBx3) {
-          Fld(i1, i2, em::bx3) = metric.template transform<3, Idx::T, Idx::U>(
-            { i1_ + HALF, i2_ + HALF },
-            finit.bx3(x_Ph_HH));
+        if constexpr (defines_bx3) {
+          if (setBx3) {
+            Fld(i1, i2, em::bx3) = metric.template transform<3, Idx::T, Idx::U>(
+              { i1_ + HALF, i2_ + HALF },
+              fset.bx3(x_Ph_HH));
+          }
         }
       } else {
         raise::KernelError(HERE, "Invalid Dimension");
@@ -412,8 +757,8 @@ namespace kernel {
                                                   x_Ph_H0H);
         metric.template convert<Crd::Cd, Crd::Ph>({ i1_, i2_ + HALF, i3_ + HALF },
                                                   x_Ph_0HH);
-        bool setEx1 = setE, setEx2 = setE, setEx3 = setE, setBx1 = setB,
-             setBx2 = setB, setBx3 = setB;
+        bool setEx1 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3,
+             setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3;
         if constexpr (O == in::x1) {
           // x1 -- normal
           // x2,x3 -- tangential
@@ -453,35 +798,47 @@ namespace kernel {
         } else {
           raise::KernelError(HERE, "Invalid Orientation");
         }
-        if (setEx1) {
-          Fld(i1, i2, i3, em::ex1) = metric.template transform<1, Idx::T, Idx::U>(
-            { i1_ + HALF, i2_, i3_ },
-            finit.ex1(x_Ph_H00));
+        if constexpr (defines_ex1) {
+          if (setEx1) {
+            Fld(i1, i2, i3, em::ex1) = metric.template transform<1, Idx::T, Idx::U>(
+              { i1_ + HALF, i2_, i3_ },
+              fset.ex1(x_Ph_H00));
+          }
         }
-        if (setEx2) {
-          Fld(i1, i2, i3, em::ex2) = metric.template transform<2, Idx::T, Idx::U>(
-            { i1_, i2_ + HALF, i3_ },
-            finit.ex2(x_Ph_0H0));
+        if constexpr (defines_ex2) {
+          if (setEx2) {
+            Fld(i1, i2, i3, em::ex2) = metric.template transform<2, Idx::T, Idx::U>(
+              { i1_, i2_ + HALF, i3_ },
+              fset.ex2(x_Ph_0H0));
+          }
         }
-        if (setEx3) {
-          Fld(i1, i2, i3, em::ex3) = metric.template transform<3, Idx::T, Idx::U>(
-            { i1_, i2_, i3_ + HALF },
-            finit.ex3(x_Ph_00H));
+        if constexpr (defines_ex3) {
+          if (setEx3) {
+            Fld(i1, i2, i3, em::ex3) = metric.template transform<3, Idx::T, Idx::U>(
+              { i1_, i2_, i3_ + HALF },
+              fset.ex3(x_Ph_00H));
+          }
         }
-        if (setBx1) {
-          Fld(i1, i2, i3, em::bx1) = metric.template transform<1, Idx::T, Idx::U>(
-            { i1_, i2_ + HALF, i3_ + HALF },
-            finit.bx1(x_Ph_0HH));
+        if constexpr (defines_bx1) {
+          if (setBx1) {
+            Fld(i1, i2, i3, em::bx1) = metric.template transform<1, Idx::T, Idx::U>(
+              { i1_, i2_ + HALF, i3_ + HALF },
+              fset.bx1(x_Ph_0HH));
+          }
         }
-        if (setBx2) {
-          Fld(i1, i2, i3, em::bx2) = metric.template transform<2, Idx::T, Idx::U>(
-            { i1_ + HALF, i2_, i3_ + HALF },
-            finit.bx2(x_Ph_H0H));
+        if constexpr (defines_bx2) {
+          if (setBx2) {
+            Fld(i1, i2, i3, em::bx2) = metric.template transform<2, Idx::T, Idx::U>(
+              { i1_ + HALF, i2_, i3_ + HALF },
+              fset.bx2(x_Ph_H0H));
+          }
         }
-        if (setBx3) {
-          Fld(i1, i2, i3, em::bx3) = metric.template transform<3, Idx::T, Idx::U>(
-            { i1_ + HALF, i2_ + HALF, i3_ },
-            finit.bx3(x_Ph_HH0));
+        if constexpr (defines_bx3) {
+          if (setBx3) {
+            Fld(i1, i2, i3, em::bx3) = metric.template transform<3, Idx::T, Idx::U>(
+              { i1_ + HALF, i2_ + HALF, i3_ },
+              fset.bx3(x_Ph_HH0));
+          }
         }
       } else {
         raise::KernelError(HERE, "Invalid Dimension");
@@ -489,6 +846,6 @@ namespace kernel {
     }
   };
 
-} // namespace kernel
+} // namespace kernel::bc
 
 #endif // KERNELS_FIELDS_BCS_HPP
diff --git a/src/kernels/tests/CMakeLists.txt b/src/kernels/tests/CMakeLists.txt
index 10e8bb944..a41ea43ef 100644
--- a/src/kernels/tests/CMakeLists.txt
+++ b/src/kernels/tests/CMakeLists.txt
@@ -31,3 +31,4 @@ gen_test(fields_to_phys)
 gen_test(prtls_to_phys)
 gen_test(gca_pusher)
 gen_test(prtl_bc)
+gen_test(flds_bc)
diff --git a/src/kernels/tests/flds_bc.cpp b/src/kernels/tests/flds_bc.cpp
new file mode 100644
index 000000000..aba829e8b
--- /dev/null
+++ b/src/kernels/tests/flds_bc.cpp
@@ -0,0 +1,210 @@
+#include "enums.h"
+#include "global.h"
+
+#include "arch/kokkos_aliases.h"
+#include "utils/comparators.h"
+#include "utils/error.h"
+
+#include "metrics/minkowski.h"
+
+#include "kernels/fields_bcs.hpp"
+
+#include <Kokkos_Core.hpp>
+
+#include <iostream>
+#include <stdexcept>
+#include <string>
+
+using namespace ntt;
+using namespace kernel::bc;
+using namespace metric;
+
+void errorIf(bool condition, const std::string& message) {
+  if (condition) {
+    throw std::runtime_error(message);
+  }
+}
+
+template <Dimension D>
+struct DummyFieldsBCs {
+  DummyFieldsBCs() {}
+
+  Inline auto ex1(const coord_t<D>&) const -> real_t {
+    return TWO;
+  }
+
+  Inline auto ex2(const coord_t<D>&) const -> real_t {
+    return THREE;
+  }
+
+  Inline auto bx2(const coord_t<D>&) const -> real_t {
+    return FOUR;
+  }
+
+  Inline auto bx3(const coord_t<D>&) const -> real_t {
+    return FIVE;
+  }
+};
+
+Inline auto equal(real_t a, real_t b, const char* msg, real_t acc) -> bool {
+  if (not(math::abs(a - b) < acc)) {
+    printf("%.12e != %.12e [%.12e] %s\n", a, b, math::abs(a - b), msg);
+    return false;
+  }
+  return true;
+}
+
+template <Dimension D>
+void testFldsBCs(const std::vector<std::size_t>& res) {
+  errorIf(res.size() != (unsigned short)D, "res.size() != D");
+  boundaries_t<real_t> sx;
+  for (const auto& r : res) {
+    sx.emplace_back(ZERO, r);
+  }
+  const auto      metric = Minkowski<D> { res, sx };
+  auto            fset   = DummyFieldsBCs<D> {};
+  ndfield_t<D, 6> flds;
+  if constexpr (D == Dim::_1D) {
+    flds = ndfield_t<D, 6> { "flds", res[0] + 2 * N_GHOSTS };
+  } else if constexpr (D == Dim::_2D) {
+    flds = ndfield_t<D, 6> { "flds", res[0] + 2 * N_GHOSTS, res[1] + 2 * N_GHOSTS };
+  } else if constexpr (D == Dim::_3D) {
+    flds = ndfield_t<D, 6> { "flds",
+                             res[0] + 2 * N_GHOSTS,
+                             res[1] + 2 * N_GHOSTS,
+                             res[2] + 2 * N_GHOSTS };
+  }
+
+  range_t<D> range;
+
+  if constexpr (D == Dim::_1D) {
+    range = CreateRangePolicy<D>({ res[0] / 2 + N_GHOSTS },
+                                 { res[0] + 2 * N_GHOSTS });
+  } else if constexpr (D == Dim::_2D) {
+    range = CreateRangePolicy<D>({ res[0] / 2 + N_GHOSTS, 0 },
+                                 { res[0] + 2 * N_GHOSTS, res[1] + N_GHOSTS });
+  } else if constexpr (D == Dim::_3D) {
+    range = CreateRangePolicy<D>(
+      { res[0] / 2 + N_GHOSTS, 0, 0 },
+      { res[0] + 2 * N_GHOSTS, res[1] + N_GHOSTS, res[2] + N_GHOSTS });
+  }
+
+  const auto xg_edge = (real_t)(sx[0].second);
+  const auto dx_abs  = (real_t)(res[0] / 10.0);
+
+  Kokkos::parallel_for(
+    "MatchBoundaries_kernel",
+    range,
+    MatchBoundaries_kernel<SimEngine::SRPIC, decltype(fset), decltype(metric), in::x1>(
+      flds,
+      fset,
+      metric,
+      xg_edge,
+      dx_abs,
+      BC::E | BC::B));
+
+  if constexpr (D == Dim::_1D) {
+    Kokkos::parallel_for(
+      "MatchBoundaries_kernel",
+      CreateRangePolicy<Dim::_1D>({ N_GHOSTS }, { res[0] + N_GHOSTS }),
+      Lambda(index_t i1) {
+        const auto x       = static_cast<real_t>(i1 - N_GHOSTS);
+        const auto factor1 = math::tanh(
+          FOUR * math::abs(x + HALF - xg_edge) / dx_abs);
+        const auto factor2 = math::tanh(FOUR * math::abs(x - xg_edge) / dx_abs);
+        if (not cmp::AlmostEqual(flds(i1, em::ex1), TWO * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, em::ex1), TWO * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect ex1");
+        }
+        if (not cmp::AlmostEqual(flds(i1, em::ex2), THREE * (ONE - factor2))) {
+          printf("%f != %f\n", flds(i1, em::ex2), THREE * (ONE - factor2));
+          raise::KernelError(HERE, "incorrect ex2");
+        }
+        if (not cmp::AlmostEqual(flds(i1, em::bx2), FOUR * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, em::bx2), FOUR * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect bx2");
+        }
+        if (not cmp::AlmostEqual(flds(i1, em::bx3), FIVE * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, em::bx3), FIVE * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect bx3");
+        }
+      });
+  } else if constexpr (D == Dim::_2D) {
+    Kokkos::parallel_for(
+      "MatchBoundaries_kernel",
+      CreateRangePolicy<Dim::_2D>({ N_GHOSTS, N_GHOSTS },
+                                  { res[0] + N_GHOSTS, res[1] + N_GHOSTS }),
+      Lambda(index_t i1, index_t i2) {
+        const auto x       = static_cast<real_t>(i1 - N_GHOSTS);
+        const auto factor1 = math::tanh(
+          FOUR * math::abs(x + HALF - xg_edge) / dx_abs);
+        const auto factor2 = math::tanh(FOUR * math::abs(x - xg_edge) / dx_abs);
+        if (not cmp::AlmostEqual(flds(i1, i2, em::ex1), TWO * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, i2, em::ex1), TWO * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect ex1");
+        }
+        if (not cmp::AlmostEqual(flds(i1, i2, em::ex2), THREE * (ONE - factor2))) {
+          printf("%f != %f\n", flds(i1, i2, em::ex2), THREE * (ONE - factor2));
+          raise::KernelError(HERE, "incorrect ex2");
+        }
+        if (not cmp::AlmostEqual(flds(i1, i2, em::bx2), FOUR * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, i2, em::bx2), FOUR * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect bx2");
+        }
+        if (not cmp::AlmostEqual(flds(i1, i2, em::bx3), FIVE * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, i2, em::bx3), FIVE * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect bx3");
+        }
+      });
+  } else if constexpr (D == Dim::_3D) {
+    Kokkos::parallel_for(
+      "MatchBoundaries_kernel",
+      CreateRangePolicy<Dim::_3D>(
+        { N_GHOSTS, N_GHOSTS, N_GHOSTS },
+        { res[0] + N_GHOSTS, res[1] + N_GHOSTS, res[2] + N_GHOSTS }),
+      Lambda(index_t i1, index_t i2, index_t i3) {
+        const auto x       = static_cast<real_t>(i1 - N_GHOSTS);
+        const auto factor1 = math::tanh(
+          FOUR * math::abs(x + HALF - xg_edge) / dx_abs);
+        const auto factor2 = math::tanh(FOUR * math::abs(x - xg_edge) / dx_abs);
+        if (not cmp::AlmostEqual(flds(i1, i2, i3, em::ex1), TWO * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, i2, i3, em::ex1), TWO * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect ex1");
+        }
+        if (not cmp::AlmostEqual(flds(i1, i2, i3, em::ex2),
+                                 THREE * (ONE - factor2))) {
+          printf("%f != %f\n", flds(i1, i2, i3, em::ex2), THREE * (ONE - factor2));
+          raise::KernelError(HERE, "incorrect ex2");
+        }
+        if (not cmp::AlmostEqual(flds(i1, i2, i3, em::bx2),
+                                 FOUR * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, i2, i3, em::bx2), FOUR * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect bx2");
+        }
+        if (not cmp::AlmostEqual(flds(i1, i2, i3, em::bx3),
+                                 FIVE * (ONE - factor1))) {
+          printf("%f != %f\n", flds(i1, i2, i3, em::bx3), FIVE * (ONE - factor1));
+          raise::KernelError(HERE, "incorrect bx3");
+        }
+      });
+  }
+}
+
+auto main(int argc, char* argv[]) -> int {
+  Kokkos::initialize(argc, argv);
+
+  try {
+    using namespace ntt;
+
+    testFldsBCs<Dim::_1D>({ 24 });
+    testFldsBCs<Dim::_2D>({ 64, 32 });
+    testFldsBCs<Dim::_3D>({ 14, 22, 15 });
+
+  } catch (std::exception& e) {
+    std::cerr << e.what() << std::endl;
+    Kokkos::finalize();
+    return 1;
+  }
+  Kokkos::finalize();
+  return 0;
+}
diff --git a/src/output/tests/writer-nompi.cpp b/src/output/tests/writer-nompi.cpp
index 08200d804..8fb2ac026 100644
--- a/src/output/tests/writer-nompi.cpp
+++ b/src/output/tests/writer-nompi.cpp
@@ -70,7 +70,7 @@ auto main(int argc, char* argv[]) -> int {
     {
       // write
       auto writer = out::Writer();
-      writer.init(&adios, "hdf5", "test");
+      writer.init(&adios, "hdf5", "test", false);
       writer.defineMeshLayout({ nx1, nx2, nx3 },
                               { 0, 0, 0 },
                               { nx1, nx2, nx3 },
@@ -84,13 +84,13 @@ auto main(int argc, char* argv[]) -> int {
         field_names.push_back(writer.fieldWriters()[0].name(i));
         addresses.push_back(i);
       }
-      writer.beginWriting(10, 123.0);
+      writer.beginWriting(WriteMode::Fields, 10, 123.0);
       writer.writeField<Dim::_3D, 3>(field_names, field, addresses);
-      writer.endWriting();
+      writer.endWriting(WriteMode::Fields);
 
-      writer.beginWriting(20, 123.4);
+      writer.beginWriting(WriteMode::Fields, 20, 123.4);
       writer.writeField<Dim::_3D, 3>(field_names, field, addresses);
-      writer.endWriting();
+      writer.endWriting(WriteMode::Fields);
     }
 
     adios.FlushAll();