diff --git a/.gitignore b/.gitignore
index a1b05e75..9a167b9d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -51,6 +51,7 @@ venv/
 # CMake testing files
 Testing/
 
+tags
 .clangd
 .schema.json
 *_old/
diff --git a/CMakeLists.txt b/CMakeLists.txt
index efd24099..dd22b930 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -55,7 +55,9 @@ if(${DEBUG} STREQUAL "OFF")
   set(CMAKE_BUILD_TYPE
       Release
       CACHE STRING "CMake build type")
-  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNDEBUG")
+  set(CMAKE_CXX_FLAGS
+      "${CMAKE_CXX_FLAGS} -DNDEBUG -Wno-unused-local-typedefs -Wno-unknown-cuda-version"
+  )
 else()
   set(CMAKE_BUILD_TYPE
       Debug
@@ -64,8 +66,6 @@ else()
       "${CMAKE_CXX_FLAGS} -DDEBUG -Wall -Wextra -Wno-unknown-pragmas")
 endif()
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-unused-local-typedefs")
-
 # options
 set(precisions
     "single" "double"
diff --git a/input.example.toml b/input.example.toml
index 2f6d2b28..c5622d65 100644
--- a/input.example.toml
+++ b/input.example.toml
@@ -90,10 +90,10 @@
     # Boundary conditions for fields:
     #   @required
     #   @type: 1/2/3-size array of string tuples, each of size 1 or 2
-    #   @valid: "PERIODIC", "ABSORB", "ATMOSPHERE", "CUSTOM", "HORIZON"
+    #   @valid: "PERIODIC", "ABSORB", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON"
     #   @example: [["CUSTOM", "ABSORB"]] (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: 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"]]
     fields = ""
     # Boundary conditions for fields:
@@ -319,6 +319,10 @@
   #   @default: -1.0 (disabled)
   #   @note: When `interval_time` < 0, the output is controlled by `interval`, otherwise by `interval_time`
   interval_time = ""
+  # Whether to output each timestep into separate files:
+  #   @type: bool
+  #   @default: true
+  separate_files = ""
 
   [output.fields]
     # Toggle for the field output:
diff --git a/setups/srpic/shock/pgen.hpp b/setups/srpic/shock/pgen.hpp
index f07b9987..1eedb3a0 100644
--- a/setups/srpic/shock/pgen.hpp
+++ b/setups/srpic/shock/pgen.hpp
@@ -5,6 +5,8 @@
 #include "global.h"
 
 #include "arch/traits.h"
+#include "utils/error.h"
+#include "utils/numeric.h"
 
 #include "archetypes/energy_dist.h"
 #include "archetypes/particle_injector.h"
@@ -14,6 +16,53 @@
 namespace user {
   using namespace ntt;
 
+  template <Dimension D>
+  struct InitFields {
+    /*
+      Sets up magnetic and electric field components for the simulation.
+      Must satisfy E = -v x B for Lorentz Force to be zero.
+
+      @param bmag: magnetic field scaling
+      @param btheta: magnetic field polar angle
+      @param bphi: magnetic field azimuthal angle
+      @param drift_ux: drift velocity in the x direction
+    */
+    InitFields(real_t bmag, real_t btheta, real_t bphi, real_t drift_ux)
+      : Bmag { bmag * static_cast<real_t>(convert::deg2rad) }
+      , Btheta { btheta * static_cast<real_t>(convert::deg2rad) }
+      , Bphi { bphi * static_cast<real_t>(convert::deg2rad) }
+      , Vx { drift_ux } {}
+
+    // magnetic field components
+    Inline auto bx1(const coord_t<D>&) const -> real_t {
+      return Bmag * math::cos(Btheta);
+    }
+
+    Inline auto bx2(const coord_t<D>&) const -> real_t {
+      return Bmag * math::sin(Btheta) * math::sin(Bphi);
+    }
+
+    Inline auto bx3(const coord_t<D>&) const -> real_t {
+      return Bmag * math::sin(Btheta) * math::cos(Bphi);
+    }
+
+    // electric field components
+    Inline auto ex1(const coord_t<D>&) const -> real_t {
+      return ZERO;
+    }
+
+    Inline auto ex2(const coord_t<D>&) const -> real_t {
+      return -Vx * Bmag * math::sin(Btheta) * math::cos(Bphi);
+    }
+
+    Inline auto ex3(const coord_t<D>&) const -> real_t {
+      return Vx * Bmag * math::sin(Btheta) * math::sin(Bphi);
+    }
+
+  private:
+    const real_t Btheta, Bphi, Vx, Bmag;
+  };
+
   template <SimEngine::type S, class M>
   struct PGen : public arch::ProblemGenerator<S, M> {
     // compatibility traits for the problem generator
@@ -30,20 +79,42 @@ namespace user {
 
     const real_t drift_ux, temperature;
 
+    const real_t  Btheta, Bphi, Bmag;
+    InitFields<D> init_flds;
+
     inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
-      : arch::ProblemGenerator<S, M>(p)
+      : arch::ProblemGenerator<S, M> { p }
       , drift_ux { p.template get<real_t>("setup.drift_ux") }
-      , temperature { p.template get<real_t>("setup.temperature") } {}
+      , temperature { p.template get<real_t>("setup.temperature") }
+      , Bmag { p.template get<real_t>("setup.Bmag", ZERO) }
+      , Btheta { p.template get<real_t>("setup.Btheta", ZERO) }
+      , Bphi { p.template get<real_t>("setup.Bphi", ZERO) }
+      , init_flds { Bmag, Btheta, Bphi, drift_ux } {}
 
     inline PGen() {}
 
+    auto FixField(const em& comp) const -> real_t {
+      if (comp == em::ex2) {
+        return init_flds.ex2({ ZERO });
+      } else if (comp == em::ex3) {
+        return init_flds.ex3({ ZERO });
+      } else if (comp == em::bx1) {
+        return init_flds.bx1({ ZERO });
+      } else {
+        raise::Error("Other components should not be requested when BC is in X",
+                     HERE);
+        return ZERO;
+      }
+    }
+
     inline void InitPrtls(Domain<S, M>& local_domain) {
       const auto energy_dist = arch::Maxwellian<S, M>(local_domain.mesh.metric,
                                                       local_domain.random_pool,
                                                       temperature,
                                                       -drift_ux,
                                                       in::x1);
-      const auto injector    = arch::UniformInjector<S, M, arch::Maxwellian>(
+
+      const auto injector = arch::UniformInjector<S, M, arch::Maxwellian>(
         energy_dist,
         { 1, 2 });
       arch::InjectUniform<S, M, arch::UniformInjector<S, M, arch::Maxwellian>>(
diff --git a/setups/srpic/shock/shock.py b/setups/srpic/shock/shock.py
index 64224c72..dc156557 100644
--- a/setups/srpic/shock/shock.py
+++ b/setups/srpic/shock/shock.py
@@ -2,7 +2,7 @@
 import matplotlib.pyplot as plt
 import matplotlib as mpl
 
-data = nt2r.Data("shock-03.h5")
+data = nt2r.Data("shock.h5")
 
 
 def frame(ti, f):
@@ -55,7 +55,7 @@ def frame(ti, f):
     axs = [fig.add_subplot(gs[i]) for i in range(len(quantities))]
 
     for ax, q in zip(axs, quantities):
-        q["compute"](f).coarsen(x=2, y=2).mean().plot(
+        q["compute"](f.isel(t=ti)).plot(
             ax=ax,
             cmap=q["cmap"],
             norm=q["norm"],
diff --git a/setups/srpic/shock/shock.toml b/setups/srpic/shock/shock.toml
index 7b2cdde2..c8a6f7c9 100644
--- a/setups/srpic/shock/shock.toml
+++ b/setups/srpic/shock/shock.toml
@@ -11,7 +11,7 @@
     metric = "minkowski"
 
   [grid.boundaries]
-    fields    = [["CONDUCTOR", "ABSORB"], ["PERIODIC"]]
+    fields = [["FIXED", "ABSORB"], ["PERIODIC"]]
     particles = [["REFLECT", "ABSORB"], ["PERIODIC"]]
 
 [scales]
@@ -42,6 +42,9 @@
 [setup]
   drift_ux    = 0.1
   temperature = 1e-3
+  Bmag = 1.0
+  Btheta = 0.0
+  Bphi = 0.0
 
 [output]
   interval_time = 0.1
diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp
index b5429154..d8cb1e64 100644
--- a/src/engines/srpic.hpp
+++ b/src/engines/srpic.hpp
@@ -42,7 +42,6 @@
 #include <Kokkos_Core.hpp>
 #include <Kokkos_ScatterView.hpp>
 
-#include <string>
 #include <utility>
 
 namespace ntt {
@@ -593,9 +592,9 @@ namespace ntt {
           }
         } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::ATMOSPHERE) {
           AtmosphereFieldsIn(direction, domain, tags);
-        } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::CONDUCTOR) {
-          if (domain.mesh.flds_bc_in(direction) == FldsBC::CONDUCTOR) {
-            ConductorFieldsIn(direction, domain, tags);
+        } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::FIXED) {
+          if (domain.mesh.flds_bc_in(direction) == FldsBC::FIXED) {
+            FixedFieldsIn(direction, domain, tags);
           }
         } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::CUSTOM) {
           if (domain.mesh.flds_bc_in(direction) == FldsBC::CUSTOM) {
@@ -718,11 +717,105 @@ namespace ntt {
       }
     }
 
+    void FixedFieldsIn(dir::direction_t<M::Dim> direction,
+                       domain_t&                domain,
+                       BCTags                   tags) {
+      /**
+       * fixed field boundaries
+       */
+      const auto sign = direction.get_sign();
+      const auto dim  = direction.get_dim();
+      raise::ErrorIf(dim != in::x1 and M::CoordType != Coord::Cart,
+                     "Fixed BCs only implemented for x1 in "
+                     "non-cartesian coordinates",
+                     HERE);
+      em normal_b_comp, tang_e_comp1, tang_e_comp2;
+      if (dim == in::x1) {
+        normal_b_comp = em::bx1;
+        tang_e_comp1  = em::ex2;
+        tang_e_comp2  = em::ex3;
+      } else if (dim == in::x2) {
+        normal_b_comp = em::bx2;
+        tang_e_comp1  = em::ex1;
+        tang_e_comp2  = em::ex3;
+      } else if (dim == in::x3) {
+        normal_b_comp = em::bx3;
+        tang_e_comp1  = em::ex1;
+        tang_e_comp2  = em::ex2;
+      } else {
+        raise::Error("Invalid dimension", HERE);
+      }
+      std::vector<std::size_t> xi_min, xi_max;
+      const std::vector<in>    all_dirs { in::x1, in::x2, in::x3 };
+      for (unsigned short d { 0 }; d < static_cast<unsigned short>(M::Dim); ++d) {
+        const auto dd = all_dirs[d];
+        if (dim == dd) {
+          if (sign > 0) { // + direction
+            xi_min.push_back(domain.mesh.n_all(dd) - N_GHOSTS);
+            xi_max.push_back(domain.mesh.n_all(dd));
+          } else { // - direction
+            xi_min.push_back(0);
+            xi_max.push_back(N_GHOSTS);
+          }
+        } else {
+          xi_min.push_back(0);
+          xi_max.push_back(domain.mesh.n_all(dd));
+        }
+      }
+      raise::ErrorIf(xi_min.size() != xi_max.size() or
+                       xi_min.size() != static_cast<std::size_t>(M::Dim),
+                     "Invalid range size",
+                     HERE);
+      std::vector<unsigned short> comps;
+      if (tags & BC::E) {
+        comps.push_back(tang_e_comp1);
+        comps.push_back(tang_e_comp2);
+      }
+      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
+        for (const auto& comp : comps) {
+          auto value = ZERO;
+          if constexpr (
+            traits::has_member<traits::pgen::fix_field_t, pgen_t>::value) {
+            // if fix field function present, read from it
+            value = m_pgen.FixField((em)comp);
+          }
+          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);
+          }
+        }
+      }
+    }
+
     void AtmosphereFieldsIn(dir::direction_t<M::Dim> direction,
                             domain_t&                domain,
                             BCTags                   tags) {
       /**
-       * atmosphere boundaries
+       * atmosphere field boundaries
        */
       if constexpr (traits::has_member<traits::pgen::field_driver_t, pgen_t>::value) {
         const auto [sign, dim, xg_min, xg_max] = get_atm_extent(direction);
@@ -766,7 +859,7 @@ namespace ntt {
             Kokkos::parallel_for(
               "AtmosphereBCFields",
               range,
-              kernel::AtmosphereBoundaries_kernel<decltype(field_driver), M, true, in::x1>(
+              kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, true, in::x1>(
                 domain.fields.em,
                 field_driver,
                 domain.mesh.metric,
@@ -776,7 +869,7 @@ namespace ntt {
             Kokkos::parallel_for(
               "AtmosphereBCFields",
               range,
-              kernel::AtmosphereBoundaries_kernel<decltype(field_driver), M, false, in::x1>(
+              kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, false, in::x1>(
                 domain.fields.em,
                 field_driver,
                 domain.mesh.metric,
@@ -789,7 +882,7 @@ namespace ntt {
               Kokkos::parallel_for(
                 "AtmosphereBCFields",
                 range,
-                kernel::AtmosphereBoundaries_kernel<decltype(field_driver), M, true, in::x2>(
+                kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, true, in::x2>(
                   domain.fields.em,
                   field_driver,
                   domain.mesh.metric,
@@ -799,7 +892,7 @@ namespace ntt {
               Kokkos::parallel_for(
                 "AtmosphereBCFields",
                 range,
-                kernel::AtmosphereBoundaries_kernel<decltype(field_driver), M, false, in::x2>(
+                kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, false, in::x2>(
                   domain.fields.em,
                   field_driver,
                   domain.mesh.metric,
@@ -815,7 +908,7 @@ namespace ntt {
               Kokkos::parallel_for(
                 "AtmosphereBCFields",
                 range,
-                kernel::AtmosphereBoundaries_kernel<decltype(field_driver), M, true, in::x3>(
+                kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, true, in::x3>(
                   domain.fields.em,
                   field_driver,
                   domain.mesh.metric,
@@ -825,7 +918,7 @@ namespace ntt {
               Kokkos::parallel_for(
                 "AtmosphereBCFields",
                 range,
-                kernel::AtmosphereBoundaries_kernel<decltype(field_driver), M, false, in::x3>(
+                kernel::EnforcedBoundaries_kernel<decltype(field_driver), M, false, in::x3>(
                   domain.fields.em,
                   field_driver,
                   domain.mesh.metric,
@@ -844,85 +937,6 @@ namespace ntt {
       }
     }
 
-    void ConductorFieldsIn(dir::direction_t<M::Dim> direction,
-                           domain_t&                domain,
-                           BCTags                   tags) {
-      const auto sign = direction.get_sign();
-      const auto dim  = direction.get_dim();
-      raise::ErrorIf(
-        dim != in::x1 and M::CoordType != Coord::Cart,
-        "Conductor BCs only implemented for x1 in non-cartesian coordinates",
-        HERE);
-      em normal_b_comp, tang_e_comp1, tang_e_comp2;
-      if (dim == in::x1) {
-        normal_b_comp = em::bx1;
-        tang_e_comp1  = em::ex2;
-        tang_e_comp2  = em::ex3;
-      } else if (dim == in::x2) {
-        normal_b_comp = em::bx2;
-        tang_e_comp1  = em::ex1;
-        tang_e_comp2  = em::ex3;
-      } else if (dim == in::x3) {
-        normal_b_comp = em::bx3;
-        tang_e_comp1  = em::ex1;
-        tang_e_comp2  = em::ex2;
-      } else {
-        raise::Error("Invalid dimension", HERE);
-      }
-      std::vector<std::size_t> xi_min, xi_max;
-      const std::vector<in>    all_dirs { in::x1, in::x2, in::x3 };
-      for (unsigned short d { 0 }; d < static_cast<unsigned short>(M::Dim); ++d) {
-        const auto dd = all_dirs[d];
-        if (dim == dd) {
-          if (sign > 0) { // + direction
-            xi_min.push_back(domain.mesh.n_all(dd) - N_GHOSTS);
-            xi_max.push_back(domain.mesh.n_all(dd));
-          } else { // - direction
-            xi_min.push_back(0);
-            xi_max.push_back(N_GHOSTS);
-          }
-        } else {
-          xi_min.push_back(0);
-          xi_max.push_back(domain.mesh.n_all(dd));
-        }
-      }
-      raise::ErrorIf(xi_min.size() != xi_max.size() or
-                       xi_min.size() != static_cast<std::size_t>(M::Dim),
-                     "Invalid range size",
-                     HERE);
-      std::vector<unsigned short> comps;
-      if (tags & BC::E) {
-        comps.push_back(tang_e_comp1);
-        comps.push_back(tang_e_comp2);
-      }
-      if (tags & BC::B) {
-        comps.push_back(normal_b_comp);
-      }
-      for (const auto& comp : comps) {
-        if constexpr (M::Dim == Dim::_1D) {
-          Kokkos::deep_copy(Kokkos::subview(domain.fields.em,
-                                            std::make_pair(xi_min[0], xi_max[0]),
-                                            comp),
-                            ZERO);
-        } 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),
-                            ZERO);
-        } 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),
-                            ZERO);
-        } else {
-          raise::Error("Invalid dimension", HERE);
-        }
-      }
-    }
-
     void CustomFieldsIn(dir::direction_t<M::Dim> direction,
                         domain_t&                domain,
                         BCTags                   tags) {
diff --git a/src/framework/domain/metadomain.cpp b/src/framework/domain/metadomain.cpp
index ec8561a9..ed4373df 100644
--- a/src/framework/domain/metadomain.cpp
+++ b/src/framework/domain/metadomain.cpp
@@ -384,7 +384,7 @@ namespace ntt {
 #if defined(MPI_ENABLED)
     auto dx_mins        = std::vector<real_t>(g_ndomains);
     dx_mins[g_mpi_rank] = dx_min;
-    MPI_Allgather(&dx_mins[g_mpi_rank],
+    MPI_Allgather(&dx_min,
                   1,
                   mpi::get_type<real_t>(),
                   dx_mins.data(),
diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp
index c39f0c67..6961d282 100644
--- a/src/framework/domain/output.cpp
+++ b/src/framework/domain/output.cpp
@@ -65,7 +65,8 @@ namespace ntt {
 
     g_writer.init(ptr_adios,
                   params.template get<std::string>("output.format"),
-                  params.template get<std::string>("simulation.name"));
+                  params.template get<std::string>("simulation.name"),
+                  params.template get<bool>("output.separate_files"));
     g_writer.defineMeshLayout(glob_shape_with_ghosts,
                               off_ncells_with_ghosts,
                               loc_shape_with_ghosts,
@@ -215,8 +216,8 @@ namespace ntt {
                    "local_domain is a placeholder",
                    HERE);
     logger::Checkpoint("Writing output", HERE);
-    g_writer.beginWriting(current_step, current_time);
     if (write_fields) {
+      g_writer.beginWriting(WriteMode::Fields, current_step, current_time);
       const auto incl_ghosts = params.template get<bool>("output.debug.ghosts");
       const auto dwn         = params.template get<std::vector<unsigned int>>(
         "output.fields.downsampling");
@@ -466,9 +467,11 @@ namespace ntt {
         }
         g_writer.writeField<M::Dim, 6>(names, local_domain->fields.bckp, addresses);
       }
+      g_writer.endWriting(WriteMode::Fields);
     } // end shouldWrite("fields", step, time)
 
     if (write_particles) {
+      g_writer.beginWriting(WriteMode::Particles, current_step, current_time);
       const auto prtl_stride = params.template get<std::size_t>(
         "output.particles.stride");
       for (const auto& prtl : g_writer.speciesWriters()) {
@@ -544,9 +547,11 @@ namespace ntt {
           g_writer.writeParticleQuantity(buff_x3, glob_tot, offset, prtl.name("X", 3));
         }
       }
+      g_writer.endWriting(WriteMode::Particles);
     } // end shouldWrite("particles", step, time)
 
     if (write_spectra) {
+      g_writer.beginWriting(WriteMode::Spectra, current_step, current_time);
       const auto log_bins = params.template get<bool>(
         "output.spectra.log_bins");
       const auto n_bins = params.template get<std::size_t>(
@@ -610,9 +615,9 @@ namespace ntt {
         g_writer.writeSpectrum(dn, spec.name());
       }
       g_writer.writeSpectrumBins(energy, "sEbn");
+      g_writer.endWriting(WriteMode::Spectra);
     } // end shouldWrite("spectra", step, time)
 
-    g_writer.endWriting();
     return true;
   }
 
diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp
index af7a773e..0a605651 100644
--- a/src/framework/parameters.cpp
+++ b/src/framework/parameters.cpp
@@ -456,6 +456,9 @@ namespace ntt {
         toml::find_or(toml_data, "output", "interval", defaults::output::interval));
     set("output.interval_time",
         toml::find_or<long double>(toml_data, "output", "interval_time", -1.0));
+    set("output.separate_files",
+        toml::find_or<bool>(toml_data, "output", "separate_files", true));
+
     promiseToDefine("output.fields.interval");
     promiseToDefine("output.fields.interval_time");
     promiseToDefine("output.fields.enable");
@@ -502,11 +505,16 @@ namespace ntt {
     set("output.fields.downsampling", field_dwn);
 
     // particles
+    auto       all_specs = std::vector<unsigned short> {};
+    const auto nspec     = get<std::size_t>("particles.nspec");
+    for (auto i = 0u; i < nspec; ++i) {
+      all_specs.push_back(static_cast<unsigned short>(i + 1));
+    }
     const auto prtl_out = toml::find_or(toml_data,
                                         "output",
                                         "particles",
                                         "species",
-                                        std::vector<unsigned short> {});
+                                        all_specs);
     set("output.particles.species", prtl_out);
     set("output.particles.stride",
         toml::find_or(toml_data,
diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h
index e915bdf1..9fd40e20 100644
--- a/src/global/arch/traits.h
+++ b/src/global/arch/traits.h
@@ -96,6 +96,9 @@ namespace traits {
     template <typename T>
     using field_driver_t = decltype(&T::FieldDriver);
 
+    template <typename T>
+    using fix_field_t = decltype(&T::FixField);
+
     template <typename T>
     using custom_fields_t = decltype(&T::CustomFields);
 
diff --git a/src/global/enums.h b/src/global/enums.h
index 57822dec..283cb456 100644
--- a/src/global/enums.h
+++ b/src/global/enums.h
@@ -8,8 +8,8 @@
  *   - enum ntt::SimEngine         // SRPIC, GRPIC
  *   - enum ntt::PrtlBC            // periodic, absorb, atmosphere, custom,
  *                                    reflect, horizon, axis, sync
- *   - enum ntt::FldsBC            // periodic, absorb, atmosphere, custom,
- *                                    conductor, horizon, axis, sync
+ *   - enum ntt::FldsBC            // periodic, absorb, fixed, atmosphere,
+ *                                    custom, horizon, axis, sync
  *   - enum ntt::PrtlPusher        // boris, vay, photon, none
  *   - enum ntt::Cooling           // synchrotron, none
  *   - enum ntt::FldsID            // e, dive, d, divd, b, h, j,
@@ -216,9 +216,9 @@ namespace ntt {
       INVALID    = 0,
       PERIODIC   = 1,
       ABSORB     = 2,
-      ATMOSPHERE = 3,
-      CUSTOM     = 4,
-      CONDUCTOR  = 5,
+      FIXED      = 3,
+      ATMOSPHERE = 4,
+      CUSTOM     = 5,
       HORIZON    = 6,
       AXIS       = 7,
       SYNC       = 8, // <- SYNC means synchronization with other domains
@@ -226,11 +226,10 @@ namespace ntt {
 
     constexpr FldsBC(uint8_t c) : enums_hidden::BaseEnum<FldsBC> { c } {}
 
-    static constexpr type variants[] = { PERIODIC,  ABSORB,  ATMOSPHERE, CUSTOM,
-                                         CONDUCTOR, HORIZON, AXIS,       SYNC };
-    static constexpr const char* lookup[] = { "periodic",   "absorb",
-                                              "atmosphere", "custom",
-                                              "conductor",  "horizon",
+    static constexpr type variants[] = { PERIODIC, ABSORB,  FIXED, ATMOSPHERE,
+                                         CUSTOM,   HORIZON, AXIS,  SYNC };
+    static constexpr const char* lookup[] = { "periodic",   "absorb", "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 dad6afcc..6669981b 100644
--- a/src/global/global.h
+++ b/src/global/global.h
@@ -249,6 +249,17 @@ namespace Comm {
 
 typedef int CommTags;
 
+namespace WriteMode {
+  enum WriteModeTags_ {
+    None      = 0,
+    Fields    = 1 << 0,
+    Particles = 1 << 1,
+    Spectra   = 1 << 2,
+  };
+} // namespace WriteMode
+
+typedef int WriteModeTags;
+
 namespace BC {
   enum BCTags_ {
     None = 0,
diff --git a/src/global/tests/enums.cpp b/src/global/tests/enums.cpp
index 1fc57398..4d678e85 100644
--- a/src/global/tests/enums.cpp
+++ b/src/global/tests/enums.cpp
@@ -61,8 +61,8 @@ 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",    "atmosphere", "custom",
-                                "horizon",  "conductor", "axis",       "sync" };
+  enum_str_t all_fields_bcs   = { "periodic", "absorb",  "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/global/utils/numeric.h b/src/global/utils/numeric.h
index 0b09f6c1..719256d1 100644
--- a/src/global/utils/numeric.h
+++ b/src/global/utils/numeric.h
@@ -91,4 +91,8 @@ namespace constant {
   inline constexpr double        SQRT3      = 1.73205080756887729352;
 } // namespace constant
 
+namespace convert {
+  inline constexpr double deg2rad = constant::PI / 180.0;
+} // namespace convert
+
 #endif // GLOBAL_UTILS_NUMERIC_H
diff --git a/src/global/utils/timer.cpp b/src/global/utils/timer.cpp
index 7d5a9beb..a12d79b9 100644
--- a/src/global/utils/timer.cpp
+++ b/src/global/utils/timer.cpp
@@ -116,11 +116,12 @@ namespace timer {
         (timer.second / local_tot) * 100.0);
       timer_stats.insert(
         { name,
-          std::make_tuple(timer.second,
-                          timer.second / static_cast<duration_t>(npart),
-                          timer.second / static_cast<duration_t>(ncells),
-                          pcent,
-                          0u) });
+          std::make_tuple(
+            timer.second,
+            npart > 0 ? timer.second / static_cast<duration_t>(npart) : 0.0,
+            timer.second / static_cast<duration_t>(ncells),
+            pcent,
+            0u) });
     }
     timer_stats.insert({ "Total", std::make_tuple(local_tot, 0.0, 0.0, 100u, 0u) });
 #endif
diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp
index e617010b..2f2a458b 100644
--- a/src/kernels/fields_bcs.hpp
+++ b/src/kernels/fields_bcs.hpp
@@ -217,7 +217,7 @@ namespace kernel {
   };
 
   template <class I, class M, bool P, in O>
-  struct AtmosphereBoundaries_kernel {
+  struct EnforcedBoundaries_kernel {
     static constexpr Dimension D = M::Dim;
     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;
@@ -240,11 +240,11 @@ namespace kernel {
     const std::size_t i_edge;
     const bool        setE, setB;
 
-    AtmosphereBoundaries_kernel(ndfield_t<M::Dim, 6>& Fld,
-                                const I&              finit,
-                                const M&              metric,
-                                std::size_t           i_edge,
-                                BCTags                tags)
+    EnforcedBoundaries_kernel(ndfield_t<M::Dim, 6>& Fld,
+                              const I&              finit,
+                              const M&              metric,
+                              std::size_t           i_edge,
+                              BCTags                tags)
       : Fld { Fld }
       , finit { finit }
       , metric { metric }
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
index e6dbcc03..81333e9f 100644
--- a/src/output/CMakeLists.txt
+++ b/src/output/CMakeLists.txt
@@ -30,6 +30,7 @@ add_library(ntt_output ${SOURCES})
 set(libs ntt_global)
 add_dependencies(ntt_output ${libs})
 target_link_libraries(ntt_output PUBLIC ${libs})
+target_link_libraries(ntt_output PRIVATE stdc++fs)
 
 target_include_directories(
   ntt_output
diff --git a/src/output/writer.cpp b/src/output/writer.cpp
index 4ba0ea14..95965c86 100644
--- a/src/output/writer.cpp
+++ b/src/output/writer.cpp
@@ -18,6 +18,7 @@
   #include <mpi.h>
 #endif
 
+#include <filesystem>
 #include <string>
 #include <vector>
 
@@ -25,9 +26,11 @@ namespace out {
 
   void Writer::init(adios2::ADIOS*     ptr_adios,
                     const std::string& engine,
-                    const std::string& title) {
-    m_engine = engine;
-    p_adios  = ptr_adios;
+                    const std::string& title,
+                    bool               use_separate_files) {
+    m_separate_files = use_separate_files;
+    m_engine         = engine;
+    p_adios          = ptr_adios;
 
     raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE);
 
@@ -36,7 +39,7 @@ namespace out {
 
     m_io.DefineVariable<std::size_t>("Step");
     m_io.DefineVariable<long double>("Time");
-    m_fname = title + (m_engine == "hdf5" ? ".h5" : ".bp");
+    m_fname = title;
   }
 
   void Writer::addTracker(const std::string& type,
@@ -412,33 +415,75 @@ namespace out {
     m_writer.Put(vare, xe_h);
   }
 
-  void Writer::beginWriting(std::size_t tstep, long double time) {
+  void Writer::beginWriting(WriteModeTags write_mode,
+                            std::size_t   tstep,
+                            long double   time) {
+    raise::ErrorIf(write_mode == WriteMode::None, "None is not a valid mode", HERE);
     raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE);
-    p_adios->ExitComputationBlock();
-    if (m_writing_mode) {
+    if (m_active_mode != WriteMode::None) {
       raise::Fatal("Already writing", HERE);
     }
-    m_writing_mode = true;
+    m_active_mode = write_mode;
     try {
-      m_writer = m_io.Open(m_fname, m_mode);
+      std::string       filename;
+      const std::string ext = m_engine == "hdf5" ? "h5" : "bp";
+      if (m_separate_files) {
+        std::string mode_str;
+        if (m_active_mode == WriteMode::Fields) {
+          mode_str = "fields";
+        } else if (m_active_mode == WriteMode::Particles) {
+          mode_str = "particles";
+        } else if (m_active_mode == WriteMode::Spectra) {
+          mode_str = "spectra";
+        } else {
+          raise::Fatal("Unknown write mode", HERE);
+        }
+        CallOnce(
+          [](auto& main_path, auto& mode_path) {
+            const std::filesystem::path main { main_path };
+            const std::filesystem::path mode { mode_path };
+            if (!std::filesystem::exists(main_path)) {
+              std::filesystem::create_directory(main_path);
+            }
+            if (!std::filesystem::exists(main / mode)) {
+              std::filesystem::create_directory(main / mode);
+            }
+          },
+          m_fname,
+          mode_str);
+        filename = fmt::format("%s/%s/%s.%08lu.%s",
+                               m_fname.c_str(),
+                               mode_str.c_str(),
+                               mode_str.c_str(),
+                               tstep,
+                               ext.c_str());
+        m_mode   = adios2::Mode::Write;
+      } else {
+        filename = fmt::format("%s.%s", m_fname.c_str(), ext.c_str());
+        m_mode   = std::filesystem::exists(filename) ? adios2::Mode::Append
+                                                     : adios2::Mode::Write;
+      }
+      m_writer = m_io.Open(filename, m_mode);
     } catch (std::exception& e) {
       raise::Fatal(e.what(), HERE);
     }
-    m_mode = adios2::Mode::Append;
     m_writer.BeginStep();
     m_writer.Put(m_io.InquireVariable<std::size_t>("Step"), &tstep);
     m_writer.Put(m_io.InquireVariable<long double>("Time"), &time);
   }
 
-  void Writer::endWriting() {
+  void Writer::endWriting(WriteModeTags write_mode) {
+    raise::ErrorIf(write_mode == WriteMode::None, "None is not a valid mode", HERE);
     raise::ErrorIf(p_adios == nullptr, "ADIOS pointer is null", HERE);
-    if (!m_writing_mode) {
+    if (m_active_mode == WriteMode::None) {
       raise::Fatal("Not writing", HERE);
     }
-    m_writing_mode = false;
+    if (m_active_mode != write_mode) {
+      raise::Fatal("Writing mode mismatch", HERE);
+    }
+    m_active_mode = WriteMode::None;
     m_writer.EndStep();
     m_writer.Close();
-    p_adios->EnterComputationBlock();
   }
 
   template void Writer::writeField<Dim::_1D, 3>(const std::vector<std::string>&,
diff --git a/src/output/writer.h b/src/output/writer.h
index 566da44b..a8abf4b1 100644
--- a/src/output/writer.h
+++ b/src/output/writer.h
@@ -36,6 +36,8 @@ namespace out {
     adios2::Engine m_writer;
     adios2::Mode   m_mode { adios2::Mode::Write };
 
+    bool m_separate_files;
+
     // global shape of the fields array to output
     std::vector<std::size_t> m_flds_g_shape;
     // local corner of the fields array to output
@@ -63,7 +65,7 @@ namespace out {
     std::vector<OutputSpecies> m_prtl_writers;
     std::vector<OutputSpectra> m_spectra_writers;
 
-    bool m_writing_mode { false };
+    WriteModeTags m_active_mode { WriteMode::None };
 
   public:
     Writer() {}
@@ -72,7 +74,7 @@ namespace out {
 
     Writer(Writer&&) = default;
 
-    void init(adios2::ADIOS*, const std::string&, const std::string&);
+    void init(adios2::ADIOS*, const std::string&, const std::string&, bool);
 
     void setMode(adios2::Mode);
 
@@ -106,8 +108,8 @@ namespace out {
     void writeSpectrum(const array_t<real_t*>&, const std::string&);
     void writeSpectrumBins(const array_t<real_t*>&, const std::string&);
 
-    void beginWriting(std::size_t, long double);
-    void endWriting();
+    void beginWriting(WriteModeTags, std::size_t, long double);
+    void endWriting(WriteModeTags);
 
     /* getters -------------------------------------------------------------- */
     auto fname() const -> const std::string& {