diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml
index f60ee9061..a5c546328 100644
--- a/.github/workflows/actions.yml
+++ b/.github/workflows/actions.yml
@@ -27,9 +27,14 @@ jobs:
       matrix:
         device: [cpu, amd-gpu, nvidia-gpu]
         precision: [double, single]
+        mpi: [serial, parallel]
         exclude: # my AMD GPU doesn't support fp64 atomics : (
         - device: amd-gpu 
           precision: double
+        - device: amd-gpu
+          mpi: parallel
+        - device: nvidia-gpu
+          mpi: parallel
     runs-on: [self-hosted, "${{ matrix.device }}"]
     steps:
       - name: Checkout
@@ -47,7 +52,7 @@ jobs:
             fi
           elif [ "${{ matrix.device }}" = "amd-gpu" ]; then
             FLAGS="-D Kokkos_ENABLE_HIP=ON -D Kokkos_ARCH_AMD_GFX1100=ON"
-          elif [ "${{ matrix.device }}" = "cpu" ]; then
+          elif [ "${{ matrix.mpi }}" = "parallel" ]; then
             FLAGS="-D mpi=ON"
           fi
           cmake -B build -D TESTS=ON -D output=ON -D precision=${{ matrix.precision }} $FLAGS
diff --git a/.gitmodules b/.gitmodules
index e06c332fe..577d08ea6 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -4,7 +4,6 @@
 [submodule "extern/adios2"]
 	path = extern/adios2
 	url = https://github.com/ornladios/ADIOS2.git
-	branch = master
 [submodule "extern/Kokkos"]
 	path = extern/Kokkos
 	url = https://github.com/kokkos/kokkos.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index dd22b9308..f83e6637c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,3 +1,5 @@
+# cmake-lint: disable=C0103,C0111,E1120,R0913,R0915
+
 cmake_minimum_required(VERSION 3.16)
 cmake_policy(SET CMP0110 NEW)
 
@@ -8,10 +10,10 @@ project(
   VERSION 1.2.0
   LANGUAGES CXX C)
 add_compile_options("-D ENTITY_VERSION=\"${PROJECT_VERSION}\"")
+set(hash_cmd "git diff --quiet src/ && echo $(git rev-parse HEAD) ")
+string(APPEND hash_cmd "|| echo $(git rev-parse HEAD)-mod")
 execute_process(
-  COMMAND
-    bash -c
-    "git diff --quiet src/ && echo $(git rev-parse HEAD) || echo $(git rev-parse HEAD)-mod"
+  COMMAND bash -c ${hash_cmd}
   WORKING_DIRECTORY "${CMAKE_SOURCE_DIR}"
   OUTPUT_VARIABLE GIT_HASH
   ERROR_QUIET OUTPUT_STRIP_TRAILING_WHITESPACE)
@@ -55,9 +57,8 @@ if(${DEBUG} STREQUAL "OFF")
   set(CMAKE_BUILD_TYPE
       Release
       CACHE STRING "CMake build type")
-  set(CMAKE_CXX_FLAGS
-      "${CMAKE_CXX_FLAGS} -DNDEBUG -Wno-unused-local-typedefs -Wno-unknown-cuda-version"
-  )
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DNDEBUG "
+                      "-Wno-unused-local-typedefs -Wno-unknown-cuda-version")
 else()
   set(CMAKE_BUILD_TYPE
       Debug
@@ -79,11 +80,10 @@ set(BUILD_TESTING
     CACHE BOOL "Build tests")
 
 # ------------------------ Third-party dependencies ------------------------ #
-include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/kokkosConfig.cmake)
 include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/dependencies.cmake)
 
-find_or_fetch_dependency(Kokkos FALSE)
-find_or_fetch_dependency(plog TRUE)
+find_or_fetch_dependency(Kokkos FALSE QUIET)
+find_or_fetch_dependency(plog TRUE QUIET)
 set(DEPENDENCIES Kokkos::kokkos)
 include_directories(${plog_SRC}/include)
 
@@ -92,25 +92,16 @@ set_precision(${precision})
 
 # MPI
 if(${mpi})
-  include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/MPIConfig.cmake)
+  find_or_fetch_dependency(MPI FALSE REQUIRED)
+  include_directories(${MPI_CXX_INCLUDE_PATH})
+  add_compile_options("-D MPI_ENABLED")
   set(DEPENDENCIES ${DEPENDENCIES} MPI::MPI_CXX)
 endif()
 
 # Output
 if(${output})
-  include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/adios2Config.cmake)
-  find_or_fetch_dependency(adios2 FALSE)
-  if(NOT DEFINED ENV{HDF5_ROOT})
-    if(DEFINED ENV{CONDA_PREFIX})
-      execute_process(COMMAND bash -c "conda list | grep \"hdf5\" -q"
-                      RESULT_VARIABLE HDF5_INSTALLED)
-      if(HDF5_INSTALLED EQUAL 0)
-        set(HDF5_ROOT $ENV{CONDA_PREFIX})
-      endif()
-    endif()
-  endif()
-  find_package(HDF5 REQUIRED)
-
+  find_or_fetch_dependency(adios2 FALSE QUIET)
+  add_compile_options("-D OUTPUT_ENABLED")
   if(${mpi})
     set(DEPENDENCIES ${DEPENDENCIES} adios2::cxx11_mpi)
   else()
@@ -129,11 +120,12 @@ elseif(BENCHMARK)
 else()
   # ----------------------------------- GUI ---------------------------------- #
   if(${gui})
-    find_or_fetch_dependency(nttiny FALSE)
+    find_or_fetch_dependency(nttiny FALSE QUIET)
   endif()
 
   # ------------------------------- Main source ------------------------------ #
   set_problem_generator(${pgen})
   add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/src src)
-  include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/report.cmake)
 endif()
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/report.cmake)
diff --git a/cmake/MPIConfig.cmake b/cmake/MPIConfig.cmake
deleted file mode 100644
index d1bfeaab2..000000000
--- a/cmake/MPIConfig.cmake
+++ /dev/null
@@ -1,4 +0,0 @@
-find_package(MPI REQUIRED)
-include_directories(${MPI_CXX_INCLUDE_PATH})
-add_compile_options("-D MPI_ENABLED")
-
diff --git a/cmake/adios2Config.cmake b/cmake/adios2Config.cmake
index 5c480f3d8..a4ce46179 100644
--- a/cmake/adios2Config.cmake
+++ b/cmake/adios2Config.cmake
@@ -12,16 +12,18 @@ set(ADIOS2_USE_Fortran
     CACHE BOOL "Use Fortran for ADIOS2")
 
 # Format/compression support
-set(ADIOS2_USE_ZeroMQ
-    OFF
-    CACHE BOOL "Use ZeroMQ for ADIOS2")
+set(ADIOS2_USE_HDF5
+    ON
+    CACHE BOOL "Use HDF5 for ADIOS2")
 
 set(ADIOS2_USE_MPI
     ${mpi}
     CACHE BOOL "Use MPI for ADIOS2")
 
+set(ADIOS2_USE_ZeroMQ
+    OFF
+    CACHE BOOL "Use ZeroMQ for ADIOS2")
+
 set(ADIOS2_USE_CUDA
     OFF
     CACHE BOOL "Use CUDA for ADIOS2")
-
-add_compile_options("-D OUTPUT_ENABLED")
diff --git a/cmake/benchmark.cmake b/cmake/benchmark.cmake
index d2e8ca47c..39b075716 100644
--- a/cmake/benchmark.cmake
+++ b/cmake/benchmark.cmake
@@ -1,3 +1,5 @@
+# cmake-lint: disable=C0103
+
 set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
 
 add_subdirectory(${SRC_DIR}/global ${CMAKE_CURRENT_BINARY_DIR}/global)
diff --git a/cmake/config.cmake b/cmake/config.cmake
index 58dd467e9..7214812cd 100644
--- a/cmake/config.cmake
+++ b/cmake/config.cmake
@@ -1,3 +1,5 @@
+# cmake-lint: disable=C0103
+
 # -------------------------------- Precision ------------------------------- #
 function(set_precision precision_name)
   list(FIND precisions ${precision_name} PRECISION_FOUND)
@@ -28,10 +30,8 @@ function(set_problem_generator pgen_name)
   endforeach()
   list(FIND PGEN_NAMES ${pgen_name} PGEN_FOUND)
   if(NOT ${pgen_name} STREQUAL "." AND ${PGEN_FOUND} EQUAL -1)
-    message(
-      FATAL_ERROR
-        "Invalid problem generator: ${pgen_name}\nValid options are: ${PGEN_NAMES}"
-    )
+    message(FATAL_ERROR "Invalid problem generator: "
+                        "${pgen_name}\nValid options are: ${PGEN_NAMES}")
   endif()
   set(PGEN
       ${pgen_name}
diff --git a/cmake/defaults.cmake b/cmake/defaults.cmake
index 46b4609c5..30e605a5c 100644
--- a/cmake/defaults.cmake
+++ b/cmake/defaults.cmake
@@ -1,3 +1,5 @@
+# cmake-lint: disable=C0103
+
 # ----------------------------- Defaults ---------------------------------- #
 if(DEFINED ENV{Entity_ENABLE_DEBUG})
   set(default_debug
@@ -33,7 +35,7 @@ if(DEFINED ENV{Entity_ENABLE_OUTPUT})
       CACHE INTERNAL "Default flag for output")
 else()
   set(default_output
-      OFF
+      ON
       CACHE INTERNAL "Default flag for output")
 endif()
 
@@ -51,42 +53,6 @@ endif()
 
 set_property(CACHE default_gui PROPERTY TYPE BOOL)
 
-if(DEFINED ENV{Kokkos_ENABLE_CUDA})
-  set(default_KOKKOS_ENABLE_CUDA
-      $ENV{Kokkos_ENABLE_CUDA}
-      CACHE INTERNAL "Default flag for CUDA")
-else()
-  set(default_KOKKOS_ENABLE_CUDA
-      OFF
-      CACHE INTERNAL "Default flag for CUDA")
-endif()
-
-set_property(CACHE default_KOKKOS_ENABLE_CUDA PROPERTY TYPE BOOL)
-
-if(DEFINED ENV{Kokkos_ENABLE_HIP})
-  set(default_KOKKOS_ENABLE_HIP
-      $ENV{Kokkos_ENABLE_HIP}
-      CACHE INTERNAL "Default flag for HIP")
-else()
-  set(default_KOKKOS_ENABLE_HIP
-      OFF
-      CACHE INTERNAL "Default flag for HIP")
-endif()
-
-set_property(CACHE default_KOKKOS_ENABLE_HIP PROPERTY TYPE BOOL)
-
-if(DEFINED ENV{Kokkos_ENABLE_OPENMP})
-  set(default_KOKKOS_ENABLE_OPENMP
-      $ENV{Kokkos_ENABLE_OPENMP}
-      CACHE INTERNAL "Default flag for OpenMP")
-else()
-  set(default_KOKKOS_ENABLE_OPENMP
-      OFF
-      CACHE INTERNAL "Default flag for OpenMP")
-endif()
-
-set_property(CACHE default_KOKKOS_ENABLE_OPENMP PROPERTY TYPE BOOL)
-
 if(DEFINED ENV{Entity_ENABLE_MPI})
   set(default_mpi
       $ENV{Entity_ENABLE_MPI}
diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake
index 06a3e6a1f..a21ea00e8 100644
--- a/cmake/dependencies.cmake
+++ b/cmake/dependencies.cmake
@@ -1,12 +1,15 @@
+# cmake-lint: disable=C0103,C0111,R0915,R0912
+
 set(Kokkos_REPOSITORY
     https://github.com/kokkos/kokkos.git
     CACHE STRING "Kokkos repository")
 set(plog_REPOSITORY
     https://github.com/SergiusTheBest/plog.git
     CACHE STRING "plog repository")
+set(adios2_REPOSITORY
+    https://github.com/ornladios/ADIOS2.git
+    CACHE STRING "ADIOS2 repository")
 
-# set (adios2_REPOSITORY https://github.com/ornladios/ADIOS2.git CACHE STRING
-# "ADIOS2 repository")
 function(check_internet_connection)
   if(OFFLINE STREQUAL "ON")
     set(FETCHCONTENT_FULLY_DISCONNECTED
@@ -34,25 +37,28 @@ function(check_internet_connection)
   endif()
 endfunction()
 
-function(find_or_fetch_dependency package_name header_only)
+function(find_or_fetch_dependency package_name header_only mode)
   if(NOT header_only)
-    find_package(${package_name} QUIET)
+    find_package(${package_name} ${mode})
   endif()
 
   if(NOT ${package_name}_FOUND)
+    if(${package_name} STREQUAL "Kokkos")
+      include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/kokkosConfig.cmake)
+    elseif(${package_name} STREQUAL "adios2")
+      include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/adios2Config.cmake)
+    endif()
     if(DEFINED ${package_name}_REPOSITORY AND NOT
                                               FETCHCONTENT_FULLY_DISCONNECTED)
       # fetching package
-      message(
-        STATUS
-          "${Blue}${package_name} not found. Fetching from ${${package_name}_REPOSITORY}${ColorReset}"
-      )
+      message(STATUS "${Blue}${package_name} not found. "
+                     "Fetching from ${${package_name}_REPOSITORY}${ColorReset}")
       include(FetchContent)
       if(${package_name} STREQUAL "Kokkos")
         FetchContent_Declare(
           ${package_name}
           GIT_REPOSITORY ${${package_name}_REPOSITORY}
-          GIT_TAG 4.3.00)
+          GIT_TAG 4.5.01)
       else()
         FetchContent_Declare(${package_name}
                              GIT_REPOSITORY ${${package_name}_REPOSITORY})
@@ -65,6 +71,9 @@ function(find_or_fetch_dependency package_name header_only)
       set(${package_name}_SRC
           ${CMAKE_CURRENT_BINARY_DIR}/_deps/${lower_pckg_name}-src
           CACHE PATH "Path to ${package_name} src")
+      set(${package_name}_BUILD_DIR
+          ${CMAKE_CURRENT_BINARY_DIR}/_deps/${lower_pckg_name}-build
+          CACHE PATH "Path to ${package_name} build")
       set(${package_name}_FETCHED
           TRUE
           CACHE BOOL "Whether ${package_name} was fetched")
@@ -127,7 +136,35 @@ function(find_or_fetch_dependency package_name header_only)
           ${Kokkos_VERSION}
           CACHE INTERNAL "Kokkos version")
     endif()
+    if(NOT DEFINED Kokkos_ARCH
+       OR Kokkos_ARCH STREQUAL ""
+       OR NOT DEFINED Kokkos_DEVICES
+       OR Kokkos_DEVICES STREQUAL "")
+      if(${Kokkos_FOUND})
+        include(${Kokkos_DIR}/KokkosConfigCommon.cmake)
+      elseif(NOT ${Kokkos_BUILD_DIR} STREQUAL "")
+        include(${Kokkos_BUILD_DIR}/KokkosConfigCommon.cmake)
+      else()
+        message(
+          STATUS "${Red}Kokkos_DIR and Kokkos_BUILD_DIR not set.${ColorReset}")
+      endif()
+    endif()
+    set(Kokkos_ARCH
+        ${Kokkos_ARCH}
+        PARENT_SCOPE)
+    set(Kokkos_DEVICES
+        ${Kokkos_DEVICES}
+        PARENT_SCOPE)
   endif()
+  set(${package_name}_FOUND
+      ${${package_name}_FOUND}
+      PARENT_SCOPE)
+  set(${package_name}_FETCHED
+      ${${package_name}_FETCHED}
+      PARENT_SCOPE)
+  set(${package_name}_BUILD_DIR
+      ${${package_name}_BUILD_DIR}
+      PARENT_SCOPE)
 endfunction()
 
 check_internet_connection()
diff --git a/cmake/kokkosConfig.cmake b/cmake/kokkosConfig.cmake
index 63c32622d..f1a1cf207 100644
--- a/cmake/kokkosConfig.cmake
+++ b/cmake/kokkosConfig.cmake
@@ -1,3 +1,5 @@
+# cmake-lint: disable=C0103
+
 # ----------------------------- Kokkos settings ---------------------------- #
 if(${DEBUG} STREQUAL "OFF")
   set(Kokkos_ENABLE_AGGRESSIVE_VECTORIZATION
@@ -27,51 +29,6 @@ else()
       CACHE BOOL "Kokkos debug bounds check")
 endif()
 
-set(Kokkos_ENABLE_HIP
-    ${default_KOKKOS_ENABLE_HIP}
-    CACHE BOOL "Enable HIP")
-set(Kokkos_ENABLE_CUDA
-    ${default_KOKKOS_ENABLE_CUDA}
-    CACHE BOOL "Enable CUDA")
-set(Kokkos_ENABLE_OPENMP
-    ${default_KOKKOS_ENABLE_OPENMP}
-    CACHE BOOL "Enable OpenMP")
-
-# set memory space
-if(${Kokkos_ENABLE_CUDA})
-  add_compile_definitions(CUDA_ENABLED)
-  set(ACC_MEM_SPACE Kokkos::CudaSpace)
-elseif(${Kokkos_ENABLE_HIP})
-  add_compile_definitions(HIP_ENABLED)
-  set(ACC_MEM_SPACE Kokkos::HIPSpace)
-else()
-  set(ACC_MEM_SPACE Kokkos::HostSpace)
-endif()
-
-set(HOST_MEM_SPACE Kokkos::HostSpace)
-
-# set execution space
-if(${Kokkos_ENABLE_CUDA})
-  set(ACC_EXE_SPACE Kokkos::Cuda)
-elseif(${Kokkos_ENABLE_HIP})
-  set(ACC_EXE_SPACE Kokkos::HIP)
-elseif(${Kokkos_ENABLE_OPENMP})
-  set(ACC_EXE_SPACE Kokkos::OpenMP)
-else()
-  set(ACC_EXE_SPACE Kokkos::Serial)
-endif()
-
-if(${Kokkos_ENABLE_OPENMP})
-  set(HOST_EXE_SPACE Kokkos::OpenMP)
-else()
-  set(HOST_EXE_SPACE Kokkos::Serial)
-endif()
-
-add_compile_options("-D AccelExeSpace=${ACC_EXE_SPACE}")
-add_compile_options("-D AccelMemSpace=${ACC_MEM_SPACE}")
-add_compile_options("-D HostExeSpace=${HOST_EXE_SPACE}")
-add_compile_options("-D HostMemSpace=${HOST_MEM_SPACE}")
-
 if(${BUILD_TESTING} STREQUAL "OFF")
   set(Kokkos_ENABLE_TESTS
       OFF
diff --git a/cmake/report.cmake b/cmake/report.cmake
index 13dde63f7..5a38b0dd5 100644
--- a/cmake/report.cmake
+++ b/cmake/report.cmake
@@ -1,95 +1,3 @@
-function(PadTo Text Padding Target Result)
-  set(rt ${Text})
-  string(FIND ${rt} "${Magenta}" mg_fnd)
-
-  if(mg_fnd GREATER -1)
-    string(REGEX REPLACE "${Esc}\\[35m" "" rt ${rt})
-  endif()
-
-  string(LENGTH "${rt}" TextLength)
-  math(EXPR PaddingNeeded "${Target} - ${TextLength}")
-  set(rt ${Text})
-
-  if(PaddingNeeded GREATER 0)
-    foreach(i RANGE 0 ${PaddingNeeded})
-      set(rt "${rt}${Padding}")
-    endforeach()
-  else()
-    set(${rt} "${rt}")
-  endif()
-
-  set(${Result}
-      "${rt}"
-      PARENT_SCOPE)
-endfunction()
-
-function(
-  PrintChoices
-  Label
-  Flag
-  Choices
-  Value
-  Default
-  Color
-  OutputString
-  Multiline
-  Padding)
-  list(LENGTH "${Choices}" nchoices)
-  set(rstring "")
-  set(counter 0)
-
-  foreach(ch ${Choices})
-    if(${counter} EQUAL 0)
-      set(rstring_i "- ${Label}")
-
-      if(NOT "${Flag}" STREQUAL "")
-        set(rstring_i "${rstring_i} [${Magenta}${Flag}${ColorReset}]")
-      endif()
-
-      set(rstring_i "${rstring_i}:")
-      padto("${rstring_i}" " " ${Padding} rstring_i)
-    else()
-      set(rstring_i "")
-
-      if(NOT ${counter} EQUAL ${nchoices})
-        if(${Multiline} EQUAL 1)
-          set(rstring_i "${rstring_i}\n")
-          padto("${rstring_i}" " " ${Padding} rstring_i)
-        else()
-          set(rstring_i "${rstring_i}/")
-        endif()
-      endif()
-    endif()
-
-    if(${ch} STREQUAL ${Value})
-      if(${ch} STREQUAL "ON")
-        set(col ${Green})
-      elseif(${ch} STREQUAL "OFF")
-        set(col ${Red})
-      else()
-        set(col ${Color})
-      endif()
-    else()
-      set(col ${Dim})
-    endif()
-
-    if(${ch} STREQUAL ${Default})
-      set(col ${Underline}${col})
-    endif()
-
-    set(rstring_i "${rstring_i}${col}${ch}${ColorReset}")
-    math(EXPR counter "${counter} + 1")
-    set(rstring "${rstring}${rstring_i}")
-    set(rstring_i "")
-  endforeach()
-
-  set(${OutputString}
-      "${rstring}"
-      PARENT_SCOPE)
-endfunction()
-
-set(ON_OFF_VALUES "ON" "OFF")
-
 if(${PGEN_FOUND})
   printchoices(
     "Problem generator"
@@ -99,8 +7,25 @@ if(${PGEN_FOUND})
     ${default_pgen}
     "${Blue}"
     PGEN_REPORT
-    1
-    36)
+    0)
+elseif(${TESTS})
+  set(TEST_NAMES "")
+  foreach(test_dir IN LISTS TEST_DIRECTORIES)
+    get_property(
+      LOCAL_TEST_NAMES
+      DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/${test_dir}/tests
+      PROPERTY TESTS)
+    list(APPEND TEST_NAMES ${LOCAL_TEST_NAMES})
+  endforeach()
+  printchoices(
+    "Test cases"
+    ""
+    "${TEST_NAMES}"
+    ""
+    "${ColorReset}"
+    ""
+    TESTS_REPORT
+    0)
 endif()
 
 printchoices(
@@ -111,7 +36,6 @@ printchoices(
   ${default_precision}
   "${Blue}"
   PRECISION_REPORT
-  1
   36)
 printchoices(
   "Output"
@@ -121,17 +45,6 @@ printchoices(
   ${default_output}
   "${Green}"
   OUTPUT_REPORT
-  0
-  36)
-printchoices(
-  "GUI"
-  "gui"
-  "${ON_OFF_VALUES}"
-  ${gui}
-  ${default_gui}
-  "${Green}"
-  GUI_REPORT
-  0
   36)
 printchoices(
   "MPI"
@@ -141,8 +54,7 @@ printchoices(
   OFF
   "${Green}"
   MPI_REPORT
-  0
-  42)
+  36)
 printchoices(
   "Debug mode"
   "DEBUG"
@@ -151,200 +63,158 @@ printchoices(
   OFF
   "${Green}"
   DEBUG_REPORT
-  0
-  42)
+  36)
 
-printchoices(
-  "CUDA"
-  "Kokkos_ENABLE_CUDA"
-  "${ON_OFF_VALUES}"
-  ${Kokkos_ENABLE_CUDA}
-  "OFF"
-  "${Green}"
-  CUDA_REPORT
-  0
-  42)
-printchoices(
-  "HIP"
-  "Kokkos_ENABLE_HIP"
-  "${ON_OFF_VALUES}"
-  ${Kokkos_ENABLE_HIP}
-  "OFF"
-  "${Green}"
-  HIP_REPORT
-  0
-  42)
-printchoices(
-  "OpenMP"
-  "Kokkos_ENABLE_OPENMP"
-  "${ON_OFF_VALUES}"
-  ${Kokkos_ENABLE_OPENMP}
-  "OFF"
-  "${Green}"
-  OPENMP_REPORT
-  0
-  42)
+if(NOT ${PROJECT_VERSION_TWEAK} EQUAL 0)
+  set(VERSION_SYMBOL "v${PROJECT_VERSION_MAJOR}." "${PROJECT_VERSION_MINOR}.")
+  string(APPEND VERSION_SYMBOL
+         "${PROJECT_VERSION_PATCH}-rc${PROJECT_VERSION_TWEAK}")
+else()
+  set(VERSION_SYMBOL "v${PROJECT_VERSION_MAJOR}.")
+  string(APPEND VERSION_SYMBOL
+         "${PROJECT_VERSION_MINOR}.${PROJECT_VERSION_PATCH}    ")
+endif()
 
-printchoices(
-  "C++ compiler"
-  "CMAKE_CXX_COMPILER"
-  "${CMAKE_CXX_COMPILER} v${CMAKE_CXX_COMPILER_VERSION}"
-  "${CMAKE_CXX_COMPILER} v${CMAKE_CXX_COMPILER_VERSION}"
-  "N/A"
-  "${ColorReset}"
-  CXX_COMPILER_REPORT
-  0
-  42)
+set(REPORT_TEXT
+    "${Blue}              __        __
+             /\\ \\__  __/\\ \\__
+    __    ___\\ \\  _\\/\\_\\ \\  _\\  __  __
+  / __ \\ / __ \\ \\ \\/\\/\\ \\ \\ \\/ /\\ \\/\\ \\
+ /\\  __//\\ \\/\\ \\ \\ \\_\\ \\ \\ \\ \\_\\ \\ \\_\\ \\  __
+ \\ \\____\\ \\_\\ \\_\\ \\__\\\\ \\_\\ \\__\\\\ \\____ \\/\\_\\
+  \\/____/\\/_/\\/_/\\/__/ \\/_/\\/__/ \\/___/  \\/_/
+                                    /\\___/
+Entity ${VERSION_SYMBOL}\t\t    \\/__/")
+string(APPEND REPORT_TEXT ${ColorReset} "\n")
 
-printchoices(
-  "C compiler"
-  "CMAKE_C_COMPILER"
-  "${CMAKE_C_COMPILER} v${CMAKE_C_COMPILER_VERSION}"
-  "${CMAKE_C_COMPILER} v${CMAKE_C_COMPILER_VERSION}"
-  "N/A"
-  "${ColorReset}"
-  C_COMPILER_REPORT
-  0
-  42)
+string(APPEND REPORT_TEXT ${DASHED_LINE_SYMBOL} "\n" "Configurations" "\n")
 
-get_cmake_property(_variableNames VARIABLES)
-foreach(_variableName ${_variableNames})
-  string(REGEX MATCH "Kokkos_ARCH_*" _isMatched ${_variableName})
-  if(_isMatched)
-    get_property(
-      isSet
-      CACHE ${_variableName}
-      PROPERTY VALUE)
-    if(isSet STREQUAL "ON")
-      string(REGEX REPLACE "Kokkos_ARCH_" "" ARCH ${_variableName})
-      break()
-    endif()
-  endif()
-endforeach()
-printchoices(
-  "Architecture"
-  "Kokkos_ARCH_*"
-  "${ARCH}"
-  "${ARCH}"
-  "N/A"
-  "${ColorReset}"
-  ARCH_REPORT
-  0
-  42)
+if(${PGEN_FOUND})
+  string(APPEND REPORT_TEXT "  " ${PGEN_REPORT} "\n")
+elseif(${TESTS})
+  string(APPEND REPORT_TEXT "  " ${TESTS_REPORT} "\n")
+endif()
 
-if(${Kokkos_ENABLE_CUDA})
+string(
+  APPEND
+  REPORT_TEXT
+  "  "
+  ${PRECISION_REPORT}
+  "\n"
+  "  "
+  ${OUTPUT_REPORT}
+  "\n")
+
+string(REPLACE ";" "+" Kokkos_ARCH "${Kokkos_ARCH}")
+string(REPLACE ";" "+" Kokkos_DEVICES "${Kokkos_DEVICES}")
+
+string(
+  APPEND
+  REPORT_TEXT
+  "  - ARCH [${Magenta}Kokkos_ARCH_***${ColorReset}]:         ${Kokkos_ARCH}"
+  "\n"
+  "  - DEVICES [${Magenta}Kokkos_ENABLE_***${ColorReset}]:    ${Kokkos_DEVICES}"
+  "\n"
+  "  "
+  ${MPI_REPORT}
+  "\n"
+  "  "
+  ${DEBUG_REPORT}
+  "\n"
+  ${DASHED_LINE_SYMBOL}
+  "\n"
+  "Compilers & dependencies"
+  "\n")
+
+string(
+  APPEND
+  REPORT_TEXT
+  "  - C compiler [${Magenta}CMAKE_C_COMPILER${ColorReset}]: v"
+  ${CMAKE_C_COMPILER_VERSION}
+  "\n"
+  "    ${Dim}"
+  ${CMAKE_C_COMPILER}
+  "${ColorReset}\n"
+  "  - C++ compiler [${Magenta}CMAKE_CXX_COMPILER${ColorReset}]: v"
+  ${CMAKE_CXX_COMPILER_VERSION}
+  "\n"
+  "    ${Dim}"
+  ${CMAKE_CXX_COMPILER}
+  "${ColorReset}\n")
+
+if(${Kokkos_DEVICES} MATCHES "CUDA")
   if("${CMAKE_CUDA_COMPILER}" STREQUAL "")
     execute_process(COMMAND which nvcc OUTPUT_VARIABLE CUDACOMP)
   else()
     set(CUDACOMP ${CMAKE_CUDA_COMPILER})
   endif()
-
   string(STRIP ${CUDACOMP} CUDACOMP)
-
-  message(STATUS "CUDA compiler: ${CUDACOMP}")
+  set(cmd "${CUDACOMP} --version |")
+  string(APPEND cmd " grep release | sed -e 's/.*release //' -e 's/,.*//'")
   execute_process(
-    COMMAND
-      bash -c
-      "${CUDACOMP} --version | grep release | sed -e 's/.*release //' -e 's/,.*//'"
+    COMMAND bash -c ${cmd}
     OUTPUT_VARIABLE CUDACOMP_VERSION
     OUTPUT_STRIP_TRAILING_WHITESPACE)
-
-  printchoices(
-    "CUDA compiler"
-    "CMAKE_CUDA_COMPILER"
-    "${CUDACOMP}"
-    "${CUDACOMP}"
-    "N/A"
-    "${ColorReset}"
-    CUDA_COMPILER_REPORT
-    0
-    42)
-endif()
-
-if(${Kokkos_ENABLE_HIP})
+  message(STATUS "CUDACOMP: ${CUDACOMP_VERSION}")
+  string(
+    APPEND
+    REPORT_TEXT
+    "  - CUDA compiler: v"
+    ${CUDACOMP_VERSION}
+    "\n"
+    "    ${Dim}"
+    ${CUDACOMP}
+    "${ColorReset}\n")
+elseif(${Kokkos_DEVICES} MATCHES "HIP")
+  set(cmd "hipcc --version | grep HIP | cut -d ':' -f 2 | tr -d ' '")
   execute_process(
-    COMMAND bash -c "hipcc --version | grep HIP | cut -d ':' -f 2 | tr -d ' '"
+    COMMAND bash -c ${cmd}
     OUTPUT_VARIABLE ROCM_VERSION
     OUTPUT_STRIP_TRAILING_WHITESPACE)
+  string(APPEND REPORT_TEXT "  - ROCm: v" ${ROCM_VERSION} "\n")
 endif()
 
-set(DOT_SYMBOL "${ColorReset}.")
-set(DOTTED_LINE_SYMBOL
-    "${ColorReset}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . "
-)
-
-set(DASHED_LINE_SYMBOL
-    "${ColorReset}....................................................................... "
-)
-
-if(NOT ${PROJECT_VERSION_TWEAK} EQUAL 0)
-  set(VERSION_SYMBOL
-      "v${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}.${PROJECT_VERSION_PATCH}-rc${PROJECT_VERSION_TWEAK}"
-  )
+string(APPEND REPORT_TEXT "  - Kokkos: v" ${Kokkos_VERSION} "\n")
+if(${Kokkos_FOUND})
+  string(APPEND REPORT_TEXT "    " ${Dim}${Kokkos_DIR}${ColorReset} "\n")
 else()
-  set(VERSION_SYMBOL
-      "v${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR}.${PROJECT_VERSION_PATCH}    "
-  )
+  string(APPEND REPORT_TEXT "    " ${Dim}${Kokkos_BUILD_DIR}${ColorReset} "\n")
 endif()
 
-message(
-  "${Blue}              __        __
-             /\\ \\__  __/\\ \\__
-    __    ___\\ \\  _\\/\\_\\ \\  _\\  __  __
-  / __ \\ / __ \\ \\ \\/\\/\\ \\ \\ \\/ /\\ \\/\\ \\
- /\\  __//\\ \\/\\ \\ \\ \\_\\ \\ \\ \\ \\_\\ \\ \\_\\ \\  __
- \\ \\____\\ \\_\\ \\_\\ \\__\\\\ \\_\\ \\__\\\\ \\____ \\/\\_\\
-  \\/____/\\/_/\\/_/\\/__/ \\/_/\\/__/ \\/___/  \\/_/
-                                    /\\___/
-Entity ${VERSION_SYMBOL}\t\t    \\/__/")
-message("${DASHED_LINE_SYMBOL}
-Main configurations")
-
-if(${PGEN_FOUND})
-  message("  ${PGEN_REPORT}")
-endif()
-
-message("  ${PRECISION_REPORT}")
-message("  ${OUTPUT_REPORT}")
-message("${DASHED_LINE_SYMBOL}\nCompile configurations")
-
-if(NOT "${ARCH_REPORT}" STREQUAL "")
-  message("  ${ARCH_REPORT}")
+if(${output})
+  string(APPEND REPORT_TEXT "  - ADIOS2: v" ${adios2_VERSION} "\n")
+  if(${adios2_FOUND})
+    string(APPEND REPORT_TEXT "    " "${Dim}${adios2_DIR}${ColorReset}" "\n")
+  else()
+    string(APPEND REPORT_TEXT "    " "${Dim}${adios2_BUILD_DIR}${ColorReset}"
+           "\n")
+  endif()
 endif()
-message("  ${CUDA_REPORT}")
-message("  ${HIP_REPORT}")
-message("  ${OPENMP_REPORT}")
-
-message("  ${C_COMPILER_REPORT}")
 
-message("  ${CXX_COMPILER_REPORT}")
-
-if(NOT "${CUDA_COMPILER_REPORT}" STREQUAL "")
-  message("  ${CUDA_COMPILER_REPORT}")
+string(
+  APPEND
+  REPORT_TEXT
+  ${DASHED_LINE_SYMBOL}
+  "\n"
+  "Notes"
+  "\n"
+  "    ${Dim}: Set flags with `cmake ... -D "
+  "${Magenta}<FLAG>${ColorReset}${Dim}=<VALUE>`, "
+  "the ${Underline}default${ColorReset}${Dim} value"
+  "\n"
+  "    :   will be used unless the variable is explicitly set.${ColorReset}")
+
+if(${TESTS})
+  string(
+    APPEND
+    REPORT_TEXT
+    "\n"
+    "    ${Dim}: Run the tests with the following command:"
+    "\n"
+    "    :   ctest --test-dir ${CMAKE_CURRENT_BINARY_DIR}${ColorReset}"
+    "\n")
 endif()
 
-message("  ${MPI_REPORT}")
-
-message("  ${DEBUG_REPORT}")
-
-message("${DASHED_LINE_SYMBOL}\nDependencies")
-
-if(NOT "${CUDACOMP_VERSION}" STREQUAL "")
-  message("  - CUDA:\tv${CUDACOMP_VERSION}")
-elseif(NOT "${ROCM_VERSION}" STREQUAL "")
-  message("  - ROCm:\tv${ROCM_VERSION}")
-endif()
-message("  - Kokkos:\tv${Kokkos_VERSION}")
-if(${output})
-  message("  - ADIOS2:\tv${adios2_VERSION}")
-endif()
-if(${HDF5_FOUND})
-  message("  - HDF5:\tv${HDF5_VERSION}")
-endif()
+string(APPEND REPORT_TEXT "\n")
 
-message(
-  "${DASHED_LINE_SYMBOL}
-Notes
-   ${Dim}: Set flags with `cmake ... -D ${Magenta}<FLAG>${ColorReset}${Dim}=<VALUE>`, the ${Underline}default${ColorReset}${Dim} value
-   :   will be used unless the variable is explicitly set.${ColorReset}
-")
+message(${REPORT_TEXT})
diff --git a/cmake/styling.cmake b/cmake/styling.cmake
index 70c448fff..878cb44a4 100644
--- a/cmake/styling.cmake
+++ b/cmake/styling.cmake
@@ -1,3 +1,5 @@
+# cmake-lint: disable=C0103,C0301,C0111,E1120,R0913,R0915
+
 if(NOT WIN32)
   string(ASCII 27 Esc)
   set(ColorReset "${Esc}[m")
@@ -23,17 +25,142 @@ if(NOT WIN32)
   set(StrikeEnd "${Esc}[0m")
 endif()
 
-# message("This is normal") message("${Red}This is Red${ColorReset}")
-# message("${Green}This is Green${ColorReset}") message("${Yellow}This is
-# Yellow${ColorReset}") message("${Blue}This is Blue${ColorReset}")
-# message("${Magenta}This is Magenta${ColorReset}") message("${Cyan}This is
-# Cyan${ColorReset}") message("${White}This is White${ColorReset}")
-# message("${BoldRed}This is BoldRed${ColorReset}") message("${BoldGreen}This is
-# BoldGreen${ColorReset}") message("${BoldYellow}This is
-# BoldYellow${ColorReset}") message("${BoldBlue}This is BoldBlue${ColorReset}")
-# message("${BoldMagenta}This is BoldMagenta${ColorReset}")
-# message("${BoldCyan}This is BoldCyan${ColorReset}") message("${BoldWhite}This
-# is BoldWhite\n\n${ColorReset}")
+set(DOTTED_LINE_SYMBOL "${ColorReset}. . . . . . . . . . . . . . . .")
+string(APPEND DOTTED_LINE_SYMBOL " . . . . . . . . . . . . . . . . . . . . ")
+
+set(DASHED_LINE_SYMBOL "${ColorReset}.................................")
+string(APPEND DASHED_LINE_SYMBOL "...................................... ")
+
+set(ON_OFF_VALUES "ON" "OFF")
+
+function(PureLength Text Result)
+  set(rt ${Text})
+  string(FIND ${rt} "${Magenta}" mg_fnd)
+
+  if(mg_fnd GREATER -1)
+    string(REGEX REPLACE "${Esc}\\[35m" "" rt ${rt})
+  endif()
+
+  string(LENGTH "${rt}" TextLength)
+  set(${Result}
+      "${TextLength}"
+      PARENT_SCOPE)
+endfunction()
+
+function(PadTo Text Padding Target Result)
+  purelength("${Text}" TextLength)
+  math(EXPR PaddingNeeded "${Target} - ${TextLength}")
+  set(rt ${Text})
+
+  if(PaddingNeeded GREATER 0)
+    foreach(i RANGE 0 ${PaddingNeeded})
+      set(rt "${rt}${Padding}")
+    endforeach()
+  else()
+    set(rt "${rt}")
+  endif()
+
+  set(${Result}
+      "${rt}"
+      PARENT_SCOPE)
+endfunction()
+
+function(
+  PrintChoices
+  Label
+  Flag
+  Choices
+  Value
+  Default
+  Color
+  OutputString
+  Padding)
+  set(rstring "- ${Label}")
+
+  if(NOT "${Flag}" STREQUAL "")
+    string(APPEND rstring " [${Magenta}${Flag}${ColorReset}]")
+  endif()
+
+  string(APPEND rstring ":")
+
+  if(${Padding} EQUAL 0)
+    list(LENGTH "${Choices}" nchoices)
+    math(EXPR lastchoice "${nchoices} - 1")
+    set(ncols 4)
+    math(EXPR lastcol "${ncols} - 1")
+
+    set(longest 0)
+    foreach(ch IN LISTS Choices)
+      string(LENGTH ${ch} clen)
+      if(clen GREATER longest)
+        set(longest ${clen})
+      endif()
+    endforeach()
+
+    set(counter 0)
+    foreach(ch IN LISTS Choices)
+      if(NOT ${Value} STREQUAL "")
+        if(${ch} STREQUAL ${Value})
+          set(col ${Color})
+        else()
+          set(col ${Dim})
+        endif()
+      else()
+        set(col ${Dim})
+      endif()
+
+      if(NOT ${Default} STREQUAL "")
+        if(${ch} STREQUAL ${Default})
+          set(col ${Underline}${col})
+        endif()
+      endif()
+
+      string(LENGTH "${ch}" clen)
+      math(EXPR PaddingNeeded "${longest} - ${clen} + 4")
+
+      if(counter EQUAL ${lastcol} AND NOT ${counter} EQUAL ${lastchoice})
+        string(APPEND rstring "${col}~ ${ch}${ColorReset}")
+      else()
+        if(counter EQUAL 0)
+          string(APPEND rstring "\n    ")
+        endif()
+        string(APPEND rstring "${col}~ ${ch}${ColorReset}")
+        foreach(i RANGE 0 ${PaddingNeeded})
+          string(APPEND rstring " ")
+        endforeach()
+      endif()
+
+      math(EXPR counter "(${counter} + 1) % ${ncols}")
+    endforeach()
+  else()
+    padto("${rstring}" " " ${Padding} rstring)
 
-# message()
+    set(new_choices ${Choices})
+    foreach(ch IN LISTS new_choices)
+      string(REPLACE ${ch} "${Dim}${ch}${ColorReset}" new_choices
+                     "${new_choices}")
+    endforeach()
+    set(Choices ${new_choices})
+    if(${Value} STREQUAL "ON")
+      set(col ${Green})
+    elseif(${Value} STREQUAL "OFF")
+      set(col ${Red})
+    else()
+      set(col ${Color})
+    endif()
+    if(NOT "${Value}" STREQUAL "")
+      string(REPLACE ${Value} "${col}${Value}${ColorReset}" Choices
+                     "${Choices}")
+    endif()
+    if(NOT "${Default}" STREQUAL "")
+      string(REPLACE ${Default} "${Underline}${Default}${ColorReset}" Choices
+                     "${Choices}")
+    endif()
+    string(REPLACE ";" "/" Choices "${Choices}")
+    string(APPEND rstring "${Choices}")
+  endif()
 
+  set(${OutputString}
+      "${rstring}"
+      PARENT_SCOPE)
+endfunction()
diff --git a/cmake/tests.cmake b/cmake/tests.cmake
index ca8ee69c4..0e108d365 100644
--- a/cmake/tests.cmake
+++ b/cmake/tests.cmake
@@ -13,32 +13,24 @@ if(${output})
   add_subdirectory(${SRC_DIR}/checkpoint ${CMAKE_CURRENT_BINARY_DIR}/checkpoint)
 endif()
 
-if(${mpi})
-  # tests with mpi
-  if(${output})
-    add_subdirectory(${SRC_DIR}/output/tests
-                     ${CMAKE_CURRENT_BINARY_DIR}/output/tests)
-    add_subdirectory(${SRC_DIR}/checkpoint/tests
-                     ${CMAKE_CURRENT_BINARY_DIR}/checkpoint/tests)
-    add_subdirectory(${SRC_DIR}/framework/tests
-                     ${CMAKE_CURRENT_BINARY_DIR}/framework/tests)
-  endif()
-else()
-  # tests without mpi
-  add_subdirectory(${SRC_DIR}/global/tests
-                   ${CMAKE_CURRENT_BINARY_DIR}/global/tests)
-  add_subdirectory(${SRC_DIR}/metrics/tests
-                   ${CMAKE_CURRENT_BINARY_DIR}/metrics/tests)
-  add_subdirectory(${SRC_DIR}/kernels/tests
-                   ${CMAKE_CURRENT_BINARY_DIR}/kernels/tests)
-  add_subdirectory(${SRC_DIR}/archetypes/tests
-                   ${CMAKE_CURRENT_BINARY_DIR}/archetypes/tests)
-  add_subdirectory(${SRC_DIR}/framework/tests
-                   ${CMAKE_CURRENT_BINARY_DIR}/framework/tests)
-  if(${output})
-    add_subdirectory(${SRC_DIR}/output/tests
-                     ${CMAKE_CURRENT_BINARY_DIR}/output/tests)
-    add_subdirectory(${SRC_DIR}/checkpoint/tests
-                     ${CMAKE_CURRENT_BINARY_DIR}/checkpoint/tests)
-  endif()
+set(TEST_DIRECTORIES "")
+
+if(NOT ${mpi})
+  list(APPEND TEST_DIRECTORIES global)
+  list(APPEND TEST_DIRECTORIES metrics)
+  list(APPEND TEST_DIRECTORIES kernels)
+  list(APPEND TEST_DIRECTORIES archetypes)
+  list(APPEND TEST_DIRECTORIES framework)
+elseif(${mpi} AND ${output})
+  list(APPEND TEST_DIRECTORIES framework)
 endif()
+
+if(${output})
+  list(APPEND TEST_DIRECTORIES output)
+  list(APPEND TEST_DIRECTORIES checkpoint)
+endif()
+
+foreach(test_dir IN LISTS TEST_DIRECTORIES)
+  add_subdirectory(${SRC_DIR}/${test_dir}/tests
+                   ${CMAKE_CURRENT_BINARY_DIR}/${test_dir}/tests)
+endforeach()
diff --git a/input.example.toml b/input.example.toml
index a49067811..f519a2799 100644
--- a/input.example.toml
+++ b/input.example.toml
@@ -90,7 +90,7 @@
     # Boundary conditions for fields:
     #   @required
     #   @type: 1/2/3-size array of string tuples, each of size 1 or 2
-    #   @valid: "PERIODIC", "MATCH", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON"
+    #   @valid: "PERIODIC", "MATCH", "FIXED", "ATMOSPHERE", "CUSTOM", "HORIZON", "CONDUCTOR"
     #   @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", "MATCH"]]
@@ -478,3 +478,9 @@
   #   @type: bool
   #   @default: true
   colored_stdout = ""
+  # Specify the log level:
+  #   @type: string
+  #   @valid: "VERBOSE", "WARNING", "ERROR"
+  #   @default: "VERBOSE"
+  #   @note: VERBOSE prints all messages, WARNING prints only warnings and errors, ERROR prints only errors
+  log_level = ""
diff --git a/legacy/src/pic/particles/particle_pusher.hpp b/legacy/src/pic/particles/particle_pusher.hpp
index 4c0ec639a..7991a95a4 100644
--- a/legacy/src/pic/particles/particle_pusher.hpp
+++ b/legacy/src/pic/particles/particle_pusher.hpp
@@ -1,14 +1,13 @@
 #ifndef PIC_PARTICLE_PUSHER_H
 #define PIC_PARTICLE_PUSHER_H
 
-#include "wrapper.h"
-
-#include "pic.h"
+#include "utils/qmath.h"
 
 #include "io/output.h"
 #include "meshblock/meshblock.h"
 #include "meshblock/particles.h"
-#include "utils/qmath.h"
+#include "pic.h"
+#include "wrapper.h"
 #include METRIC_HEADER
 
 #include <typeindex>
@@ -73,35 +72,34 @@ namespace ntt {
                   real_t                          time,
                   real_t                          coeff,
                   real_t                          dt,
-                  ProblemGenerator<D, PICEngine>& pgen) :
-      EB { mblock.em },
-      i1 { particles.i1 },
-      i2 { particles.i2 },
-      i3 { particles.i3 },
-      i1_prev { particles.i1_prev },
-      i2_prev { particles.i2_prev },
-      i3_prev { particles.i3_prev },
-      dx1 { particles.dx1 },
-      dx2 { particles.dx2 },
-      dx3 { particles.dx3 },
-      dx1_prev { particles.dx1_prev },
-      dx2_prev { particles.dx2_prev },
-      dx3_prev { particles.dx3_prev },
-      ux1 { particles.ux1 },
-      ux2 { particles.ux2 },
-      ux3 { particles.ux3 },
-      phi { particles.phi },
-      tag { particles.tag },
-      metric { mblock.metric },
-      time { time },
-      coeff { coeff },
-      dt { dt },
-      ni1 { (int)mblock.Ni1() },
-      ni2 { (int)mblock.Ni2() },
-      ni3 { (int)mblock.Ni3() }
+                  ProblemGenerator<D, PICEngine>& pgen)
+      : EB { mblock.em }
+      , i1 { particles.i1 }
+      , i2 { particles.i2 }
+      , i3 { particles.i3 }
+      , i1_prev { particles.i1_prev }
+      , i2_prev { particles.i2_prev }
+      , i3_prev { particles.i3_prev }
+      , dx1 { particles.dx1 }
+      , dx2 { particles.dx2 }
+      , dx3 { particles.dx3 }
+      , dx1_prev { particles.dx1_prev }
+      , dx2_prev { particles.dx2_prev }
+      , dx3_prev { particles.dx3_prev }
+      , ux1 { particles.ux1 }
+      , ux2 { particles.ux2 }
+      , ux3 { particles.ux3 }
+      , phi { particles.phi }
+      , tag { particles.tag }
+      , metric { mblock.metric }
+      , time { time }
+      , coeff { coeff }
+      , dt { dt }
+      , ni1 { (int)mblock.Ni1() }
+      , ni2 { (int)mblock.Ni2() }
+      , ni3 { (int)mblock.Ni3() }
 #ifdef EXTERNAL_FORCE
-      ,
-      pgen { pgen }
+      , pgen { pgen }
 #endif
     {
       (void)pgen;
@@ -237,7 +235,7 @@ namespace ntt {
     const auto coeff           = charge_ovr_mass * HALF * dt * params.B0();
     Kokkos::parallel_for(
       "ParticlesPush",
-      Kokkos::RangePolicy<AccelExeSpace, P>(0, particles.npart()),
+      Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace, P>(0, particles.npart()),
       Pusher_kernel<D>(mblock, particles, time, coeff, dt, pgen));
   }
 
@@ -638,9 +636,9 @@ namespace ntt {
   template <Dimension D>
   template <typename T>
   Inline void Pusher_kernel<D>::get3VelCntrv(T,
-                                             index_t&      p,
+                                             index_t&           p,
                                              vec_t<PrtlCoordD>& xp,
-                                             vec_t<Dim3>&  v) const {
+                                             vec_t<Dim3>&       v) const {
     metric.v3_Cart2Cntrv(xp, { ux1(p), ux2(p), ux3(p) }, v);
     auto inv_energy { ONE / getEnergy(T {}, p) };
     v[0] *= inv_energy;
@@ -666,7 +664,8 @@ namespace ntt {
   }
 
   template <>
-  Inline void Pusher_kernel<Dim2>::getPrtlPos(index_t& p, coord_t<PrtlCoordD>& xp) const {
+  Inline void Pusher_kernel<Dim2>::getPrtlPos(index_t&             p,
+                                              coord_t<PrtlCoordD>& xp) const {
     xp[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
     xp[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
     xp[2] = phi(p);
@@ -1066,7 +1065,7 @@ namespace ntt {
   #else
   template <Dimension D>
   Inline void Pusher_kernel<D>::initForce(coord_t<PrtlCoordD>& xp,
-                                          vec_t<Dim3>&    force_Cart) const {
+                                          vec_t<Dim3>& force_Cart) const {
     coord_t<PrtlCoordD> xp_Ph { ZERO };
     coord_t<PrtlCoordD> xp_Code { ZERO };
     for (short d { 0 }; d < static_cast<short>(PrtlCoordD); ++d) {
diff --git a/legacy/tests/kernels-gr.cpp b/legacy/tests/kernels-gr.cpp
index 84a0c952b..6962f7c9f 100644
--- a/legacy/tests/kernels-gr.cpp
+++ b/legacy/tests/kernels-gr.cpp
@@ -1,16 +1,16 @@
-#include "wrapper.h"
-
 #include <Kokkos_Core.hpp>
 
 #include <iostream>
 #include <stdexcept>
 
-#include METRIC_HEADER
+#include "wrapper.h"
 
-#include "particle_macros.h"
+#include METRIC_HEADER
 
 #include "kernels/particle_pusher_gr.hpp"
 
+#include "particle_macros.h"
+
 template <typename T>
 void put_value(ntt::array_t<T*>& arr, T value, int i) {
   auto arr_h = Kokkos::create_mirror_view(arr);
@@ -154,9 +154,10 @@ auto main(int argc, char* argv[]) -> int {
         static_cast<real_t>(1.0e-5),
         10,
         boundaries);
-      Kokkos::parallel_for("ParticlesPush",
-                           Kokkos::RangePolicy<AccelExeSpace, ntt::Massive_t>(0, 1),
-                           kernel);
+      Kokkos::parallel_for(
+        "ParticlesPush",
+        Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace, ntt::Massive_t>(0, 1),
+        kernel);
       auto [ra, tha]   = get_physical_coord(0, i1, i2, dx1, dx2, metric);
       const real_t pha = get_value(phi, 0);
 
@@ -207,4 +208,4 @@ auto main(int argc, char* argv[]) -> int {
   ntt::GlobalFinalize();
 
   return 0;
-}
\ No newline at end of file
+}
diff --git a/legacy/tests/kernels-sr.cpp b/legacy/tests/kernels-sr.cpp
index 59ce0646b..3f64122cd 100644
--- a/legacy/tests/kernels-sr.cpp
+++ b/legacy/tests/kernels-sr.cpp
@@ -1,17 +1,17 @@
-#include "wrapper.h"
-
 #include <Kokkos_Core.hpp>
 
 #include <iostream>
 #include <stdexcept>
 
+#include "wrapper.h"
+
 #include METRIC_HEADER
 #include PGEN_HEADER
 
-#include "particle_macros.h"
-
 #include "kernels/particle_pusher_sr.hpp"
 
+#include "particle_macros.h"
+
 template <typename T>
 void put_value(ntt::array_t<T*>& arr, T value, int i) {
   auto arr_h = Kokkos::create_mirror_view(arr);
@@ -181,9 +181,10 @@ auto main(int argc, char* argv[]) -> int {
                                               ZERO,
                                               ZERO,
                                               ZERO);
-      Kokkos::parallel_for("ParticlesPush",
-                           Kokkos::RangePolicy<AccelExeSpace, ntt::Boris_t>(0, 1),
-                           kernel);
+      Kokkos::parallel_for(
+        "ParticlesPush",
+        Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace, ntt::Boris_t>(0, 1),
+        kernel);
       auto [xa, ya] = get_cartesian_coord(0, i1, i2, dx1, dx2, phi, metric);
 
       if (!ntt::AlmostEqual(xa,
@@ -221,4 +222,4 @@ auto main(int argc, char* argv[]) -> int {
   ntt::GlobalFinalize();
 
   return 0;
-}
\ No newline at end of file
+}
diff --git a/setups/srpic/shock/pgen.hpp b/setups/srpic/shock/pgen.hpp
index b8f169521..ad260bda0 100644
--- a/setups/srpic/shock/pgen.hpp
+++ b/setups/srpic/shock/pgen.hpp
@@ -9,6 +9,7 @@
 #include "utils/numeric.h"
 
 #include "archetypes/energy_dist.h"
+#include "archetypes/field_setter.h"
 #include "archetypes/particle_injector.h"
 #include "archetypes/problem_generator.h"
 #include "framework/domain/metadomain.h"
@@ -79,9 +80,13 @@ namespace user {
     using arch::ProblemGenerator<S, M>::C;
     using arch::ProblemGenerator<S, M>::params;
 
-    const real_t drift_ux, temperature;
-
-    const real_t  Btheta, Bphi, Bmag;
+    // gas properties
+    const real_t  drift_ux, temperature, filling_fraction;
+    // injector properties
+    const real_t  injector_velocity, injection_start, dt;
+    const int     injection_frequency;
+    // magnetic field properties
+    real_t        Btheta, Bphi, Bmag;
     InitFields<D> init_flds;
 
     inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
@@ -91,40 +96,236 @@ namespace user {
       , 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 } {}
+      , init_flds { Bmag, Btheta, Bphi, drift_ux }
+      , filling_fraction { p.template get<real_t>("setup.filling_fraction", 1.0) }
+      , injector_velocity { p.template get<real_t>("setup.injector_velocity", 1.0) }
+      , injection_start { p.template get<real_t>("setup.injection_start", 0.0) }
+      , injection_frequency { p.template get<int>("setup.injection_frequency", 100) }
+      , dt { p.template get<real_t>("algorithms.timestep.dt") } {}
 
     inline PGen() {}
 
-    auto FixFieldsConst(const bc_in&, const em& comp) const
-      -> std::pair<real_t, bool> {
-      if (comp == em::ex2) {
-        return { init_flds.ex2({ ZERO }), true };
+    auto MatchFields(real_t time) const -> InitFields<D> {
+      return init_flds;
+    }
+
+    auto ResetFields(const em& comp) const -> real_t {
+      if (comp == em::ex1) {
+        return init_flds.ex1({ ZERO });
+      } else if (comp == em::ex2) {
+        return init_flds.ex2({ ZERO });
       } else if (comp == em::ex3) {
-        return { init_flds.ex3({ ZERO }), true };
+        return init_flds.ex3({ ZERO });
+      } else if (comp == em::bx1) {
+        return init_flds.bx1({ ZERO });
+      } else if (comp == em::bx2) {
+        return init_flds.bx2({ ZERO });
+      } else if (comp == em::bx3) {
+        return init_flds.bx3({ ZERO });
       } else {
-        return { ZERO, false };
+        raise::Error("Invalid component", HERE);
+        return ZERO;
       }
     }
 
-    auto MatchFields(real_t time) const -> InitFields<D> {
-      return init_flds;
-    }
-
     inline void InitPrtls(Domain<S, M>& local_domain) {
+
+      // minimum and maximum position of particles
+      real_t xg_min = local_domain.mesh.extent(in::x1).first;
+      real_t xg_max = local_domain.mesh.extent(in::x1).first +
+                      filling_fraction * (local_domain.mesh.extent(in::x1).second -
+                                          local_domain.mesh.extent(in::x1).first);
+
+      // define box to inject into
+      boundaries_t<real_t> box;
+      // loop over all dimensions
+      for (unsigned short d { 0 }; d < static_cast<unsigned short>(M::Dim); ++d) {
+        // compute the range for the x-direction
+        if (d == static_cast<unsigned short>(in::x1)) {
+          box.push_back({ xg_min, xg_max });
+        } else {
+          // inject into full range in other directions
+          box.push_back(Range::All);
+        }
+      }
+
+      // spatial distribution of the particles
+      // -> hack to use the uniform distribution in NonUniformInjector
+      const auto spatial_dist = arch::Piston<S, M>(local_domain.mesh.metric,
+                                                   xg_min,
+                                                   xg_max,
+                                                   in::x1);
+
+      // energy distribution of the particles
       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::NonUniformInjector<S, M, arch::Maxwellian, arch::Piston>(
         energy_dist,
+        spatial_dist,
         { 1, 2 });
-      arch::InjectUniform<S, M, arch::UniformInjector<S, M, arch::Maxwellian>>(
+
+      arch::InjectNonUniform<S, M, arch::NonUniformInjector<S, M, arch::Maxwellian, arch::Piston>>(
         params,
         local_domain,
         injector,
-        1.0);
+        1.0,   // target density
+        false, // no weights
+        box);
+    }
+
+    void CustomPostStep(std::size_t step, long double time, Domain<S, M>& domain) {
+
+      // check if the injector should be active
+      if (step % injection_frequency != 0) {
+        return;
+      }
+
+      // initial position of injector
+      const auto x_init = domain.mesh.extent(in::x1).first +
+                          filling_fraction * (domain.mesh.extent(in::x1).second -
+                                              domain.mesh.extent(in::x1).first);
+
+      // check if injector is supposed to start moving already
+      const auto dt_inj = time - injection_start > ZERO ? time - injection_start
+                                                        : ZERO;
+
+      // compute the position of the injector
+      auto xmax = x_init + injector_velocity * (dt_inj + dt);
+      if (xmax >= domain.mesh.extent(in::x1).second) {
+        xmax = domain.mesh.extent(in::x1).second;
+      }
+
+      // define box to inject into
+      boundaries_t<real_t> box;
+      // loop over all dimension
+      for (auto d = 0u; d < M::Dim; ++d) {
+        if (d == 0) {
+          box.push_back({ xmax - drift_ux / math::sqrt(1 + SQR(drift_ux)) * dt -
+                            injection_frequency * dt,
+                          xmax });
+        } else {
+          box.push_back(Range::All);
+        }
+      }
+
+      // define indice range to reset fields
+      boundaries_t<bool> incl_ghosts;
+      for (auto d = 0; d < M::Dim; ++d) {
+        incl_ghosts.push_back({ true, true });
+      }
+      auto fields_box = box;
+      // check if the box is still inside the domain
+      if (xmax + injection_frequency * dt < domain.mesh.extent(in::x1).second) {
+        fields_box[0].second += injection_frequency * dt;
+      } else { 
+        // if right side of the box is outside of the domain -> truncate box
+        fields_box[0].second = domain.mesh.extent(in::x1).second;
+      }
+      const auto extent = domain.mesh.ExtentToRange(fields_box, incl_ghosts);
+      tuple_t<std::size_t, M::Dim> x_min { 0 }, x_max { 0 };
+      for (auto d = 0; d < M::Dim; ++d) {
+        x_min[d] = extent[d].first;
+        x_max[d] = extent[d].second;
+      }
+
+      Kokkos::parallel_for("ResetFields",
+                           CreateRangePolicy<M::Dim>(x_min, x_max),
+                           arch::SetEMFields_kernel<decltype(init_flds), S, M> {
+                             domain.fields.em,
+                             init_flds,
+                             domain.mesh.metric });
+
+
+      /*
+        tag particles inside the injection zone as dead
+      */
+
+      // loop over particle species
+      for (std::size_t s { 0 }; s < 2; ++s) {
+
+        // get particle properties
+        auto& species = domain.species[s];
+        auto  i1      = species.i1;
+        auto  tag     = species.tag;
+
+        // tag all particles with x > box[0].first as dead
+        Kokkos::parallel_for(
+          "RemoveParticles",
+          species.rangeActiveParticles(),
+          Lambda(index_t p) {
+            // check if the particle is already dead
+            if (tag(p) == ParticleTag::dead) {
+              return;
+            }
+            // select the x-coordinate index
+            auto x_i1 = i1(p);
+            // check if the particle is inside the box of new plasma
+            if (x_i1 >= x_min[0]) {
+              tag(p) = ParticleTag::dead;
+            }
+          });
+      }
+
+      /*
+          Inject piston of fresh plasma
+      */
+
+      // same maxwell distribution as above
+      const auto energy_dist  = arch::Maxwellian<S, M>(domain.mesh.metric,
+                                                      domain.random_pool,
+                                                      temperature,
+                                                      -drift_ux,
+                                                      in::x1);
+      // spatial distribution of the particles
+      // -> hack to use the uniform distribution in NonUniformInjector
+      const auto spatial_dist = arch::Piston<S, M>(domain.mesh.metric,
+                                                   box[0].first,
+                                                   box[0].second,
+                                                   in::x1);
+
+      // inject piston of fresh plasma
+      const auto injector = arch::NonUniformInjector<S, M, arch::Maxwellian, arch::Piston>(
+        energy_dist,
+        spatial_dist,
+        { 1, 2 });
+
+      // inject non-uniformly within the defined box
+      arch::InjectNonUniform<S, M, decltype(injector)>(params,
+                                                       domain,
+                                                       injector,
+                                                       ONE,
+                                                       false,
+                                                       box);
+
+      /*
+          I thought this option would be better, but I can't get it to work
+      */
+
+      //   const auto spatial_dist = arch::Replenish<S, M,
+      //   decltype(TargetDensityProfile)>(domain.mesh.metric,
+      //                                                   domain.fields.bckp,
+      //                                                   box,
+      //                                                   TargetDensity,
+      //                                                   1.0);
+
+      //   const auto injector = arch::NonUniformInjector<S, M, arch::Maxwellian, arch::Replenish>(
+      //       energy_dist,
+      //       spatial_dist,
+      //       {1, 2});
+
+      // const auto injector = arch::MovingInjector<S, M, in::x1> {
+      //   domain.mesh.metric,
+      //   domain.fields.bckp,
+      //   energy_dist,
+      //   box[0].first,
+      //   box[0].second,
+      //   1.0,
+      //   { 1, 2 }
+      // };
     }
   };
 
diff --git a/setups/srpic/shock/shock.toml b/setups/srpic/shock/shock.toml
index 4ed3a2b9e..ca19a4078 100644
--- a/setups/srpic/shock/shock.toml
+++ b/setups/srpic/shock/shock.toml
@@ -1,22 +1,22 @@
 [simulation]
   name    = "shock"
   engine  = "srpic"
-  runtime = 50.0
+  runtime = 30000.0
 
 [grid]
-  resolution = [2048, 128]
-  extent     = [[0.0, 10.0], [-0.3125, 0.3125]]
+  resolution = [4096, 128]
+  extent     = [[0.0, 2000.0], [-31.25, 31.25]]
 
   [grid.metric]
     metric = "minkowski"
 
   [grid.boundaries]
-    fields = [["ABSORB", "FIXED"], ["PERIODIC"]]
+    fields = [["CONDUCTOR", "MATCH"], ["PERIODIC"]]
     particles = [["REFLECT", "ABSORB"], ["PERIODIC"]]
-
+    
 [scales]
-  larmor0    = 1e-2
-  skindepth0 = 1e-2
+  larmor0    = 100.0
+  skindepth0 = 1.0
 
 [algorithms]
   current_filters = 8
@@ -25,7 +25,7 @@
     CFL = 0.5
 
 [particles]
-  ppc0 = 16.0
+  ppc0 = 8.0
 
   [[particles.species]]
     label    = "e-"
@@ -34,21 +34,32 @@
     maxnpart = 1e8
 
   [[particles.species]]
-    label    = "e+"
+    label    = "p+"
     mass     = 1.0
     charge   = 1.0
     maxnpart = 1e8
 
 [setup]
-  drift_ux    = 0.1
-  temperature = 1e-3
-  Bmag = 1.0
-  Btheta = 0.0
-  Bphi = 0.0
+  drift_ux    = 0.1       # speed towards the wall [c]
+  temperature = 1e-4      # temeperature of maxwell distribution [m_e c^2]
+  Bmag = 1.0              # magnetic field strength as fraction of magnetisation
+  Btheta = 0.0            # magnetic field angle in the plane
+  Bphi = 0.0              # magnetic field angle out of plane
+  filling_fraction = 0.1  # fraction of the shock piston filled with plasma
+  injector_velocity = 1.0 # speed of injector [c]
+  injection_start = 0.0   # start time of moving injector
 
 [output]
-  interval_time = 0.1
+  interval_time = 10.0
   format        = "hdf5"
-
+  
   [output.fields]
-    quantities = ["N_1", "N_2", "E", "B", "T0i_1", "T0i_2", "J"]
+    quantities = ["N_1", "N_2", "B", "E"]
+
+  [output.particles]
+    enable = true
+    stride = 10 
+
+  [output.spectra]
+    enable = false
+
diff --git a/setups/srpic/shocktest/pgen.hpp b/setups/srpic/shocktest/pgen.hpp
new file mode 100644
index 000000000..a7bcf5c3d
--- /dev/null
+++ b/setups/srpic/shocktest/pgen.hpp
@@ -0,0 +1,155 @@
+#ifndef PROBLEM_GENERATOR_H
+#define PROBLEM_GENERATOR_H
+
+#include "enums.h"
+#include "global.h"
+
+#include "arch/traits.h"
+#include "utils/error.h"
+#include "utils/numeric.h"
+
+#include "archetypes/energy_dist.h"
+#include "archetypes/particle_injector.h"
+#include "archetypes/problem_generator.h"
+#include "framework/domain/metadomain.h"
+
+#include <utility>
+#include <vector>
+
+namespace user {
+  using namespace ntt;
+
+  template <Dimension D>
+  struct InitFields {
+    /*
+      Sets up magnetic and electric field components for the simulation.
+      Must satisfy E = -v x B for Lorentz Force to be zero.
+
+      @param bmag: magnetic field scaling
+      @param btheta: magnetic field polar angle
+      @param bphi: magnetic field azimuthal angle
+      @param drift_ux: drift velocity in the x direction
+    */
+    InitFields(real_t bmag, real_t btheta, real_t bphi, real_t drift_ux)
+      : Bmag { bmag }
+      , Btheta { btheta * static_cast<real_t>(convert::deg2rad) }
+      , Bphi { bphi * static_cast<real_t>(convert::deg2rad) }
+      , Vx { drift_ux } {}
+
+    // magnetic field components
+    Inline auto bx1(const coord_t<D>& x_ph) const -> real_t {
+      return Bmag * math::cos(Btheta);
+    }
+
+    Inline auto bx2(const coord_t<D>& x_ph) const -> real_t {
+      return Bmag * math::sin(Btheta) * math::sin(Bphi);
+    }
+
+    Inline auto bx3(const coord_t<D>& x_ph) const -> real_t {
+      return Bmag * math::sin(Btheta) * math::cos(Bphi);
+    }
+
+    // electric field components
+    Inline auto ex1(const coord_t<D>& x_ph) const -> real_t {
+      return ZERO;
+    }
+
+    Inline auto ex2(const coord_t<D>& x_ph) const -> real_t {
+      return -Vx * Bmag * math::sin(Btheta) * math::cos(Bphi);
+    }
+
+    Inline auto ex3(const coord_t<D>& x_ph) const -> real_t {
+      return Vx * Bmag * math::sin(Btheta) * math::sin(Bphi);
+    }
+
+  private:
+    const real_t Btheta, Bphi, Vx, Bmag;
+  };
+
+  template <SimEngine::type S, class M>
+  struct PGen : public arch::ProblemGenerator<S, M> {
+    // compatibility traits for the problem generator
+    static constexpr auto engines { traits::compatible_with<SimEngine::SRPIC>::value };
+    static constexpr auto metrics { traits::compatible_with<Metric::Minkowski>::value };
+    static constexpr auto dimensions {
+      traits::compatible_with<Dim::_1D, Dim::_2D, Dim::_3D>::value
+    };
+
+    // for easy access to variables in the child class
+    using arch::ProblemGenerator<S, M>::D;
+    using arch::ProblemGenerator<S, M>::C;
+    using arch::ProblemGenerator<S, M>::params;
+
+    const real_t drift_ux, temperature;
+
+    const std::vector<real_t> x1arr_e, x2arr_e, ux1arr_e, ux2arr_e, ux3arr_e;
+    const std::vector<real_t> x1arr_i, x2arr_i, ux1arr_i, ux2arr_i, ux3arr_i;
+
+    const real_t            Btheta, Bphi, Bmag;
+    InitFields<D>           init_flds;
+    const Metadomain<S, M>* metadomain;
+
+    inline PGen(const SimulationParams& p, const Metadomain<S, M>& m)
+      : arch::ProblemGenerator<S, M> { p }
+      , drift_ux { p.template get<real_t>("setup.drift_ux") }
+      , temperature { p.template get<real_t>("setup.temperature") }
+      , x1arr_e { p.template get<std::vector<real_t>>("setup.x_e") }
+      , x2arr_e { p.template get<std::vector<real_t>>("setup.y_e") }
+      , ux1arr_e { p.template get<std::vector<real_t>>("setup.ux_e") }
+      , ux2arr_e { p.template get<std::vector<real_t>>("setup.uy_e") }
+      , ux3arr_e { p.template get<std::vector<real_t>>("setup.uz_e") }
+      , x1arr_i { p.template get<std::vector<real_t>>("setup.x_i") }
+      , x2arr_i { p.template get<std::vector<real_t>>("setup.y_i") }
+      , ux1arr_i { p.template get<std::vector<real_t>>("setup.ux_i") }
+      , ux2arr_i { p.template get<std::vector<real_t>>("setup.uy_i") }
+      , ux3arr_i { p.template get<std::vector<real_t>>("setup.uz_i") }
+      , Btheta { p.template get<real_t>("setup.Btheta", ZERO) }
+      , Bmag { p.template get<real_t>("setup.Bmag", ZERO) }
+      , Bphi { p.template get<real_t>("setup.Bphi", ZERO) }
+      , init_flds { Bmag, Btheta, Bphi, drift_ux }
+      , metadomain { &m } {}
+
+    inline PGen() {}
+
+    auto FixFieldsConst(const bc_in&, const em& comp) const
+      -> std::pair<real_t, bool> {
+      if (comp == em::ex2) {
+        return { init_flds.ex2({ ZERO }), true };
+      } else if (comp == em::ex3) {
+        return { init_flds.ex3({ ZERO }), true };
+      } else {
+        return { ZERO, false };
+      }
+    }
+
+    auto MatchFields(real_t time) const -> InitFields<D> {
+      return init_flds;
+    }
+
+    inline void InitPrtls(Domain<S, M>& domain) {
+      arch::InjectGlobally(*metadomain,
+                           domain,
+                           1,
+                           {
+                             {  "x1",  x1arr_e },
+                             {  "x2",  x2arr_e },
+                             { "ux1", ux1arr_e },
+                             { "ux2", ux1arr_e },
+                             { "ux3", ux3arr_e }
+      });
+      arch::InjectGlobally(*metadomain,
+                           domain,
+                           2,
+                           {
+                             {  "x1",  x1arr_i },
+                             {  "x2",  x2arr_i },
+                             { "ux1", ux1arr_i },
+                             { "ux2", ux1arr_i },
+                             { "ux3", ux3arr_i }
+      });
+    }
+  };
+
+} // namespace user
+
+#endif
diff --git a/setups/srpic/shocktest/shock.py b/setups/srpic/shocktest/shock.py
new file mode 100644
index 000000000..dc1565572
--- /dev/null
+++ b/setups/srpic/shocktest/shock.py
@@ -0,0 +1,75 @@
+import nt2.read as nt2r
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+data = nt2r.Data("shock.h5")
+
+
+def frame(ti, f):
+    quantities = [
+        {
+            "name": "density",
+            "compute": lambda f: f.N_2 + f.N_1,
+            "cmap": "inferno",
+            "norm": mpl.colors.Normalize(0, 5),
+        },
+        {
+            "name": r"$E_x$",
+            "compute": lambda f: f.Ex,
+            "cmap": "RdBu_r",
+            "norm": mpl.colors.Normalize(-0.05, 0.05),
+        },
+        {
+            "name": r"$E_y$",
+            "compute": lambda f: f.Ey,
+            "cmap": "RdBu_r",
+            "norm": mpl.colors.Normalize(-0.05, 0.05),
+        },
+        {
+            "name": r"$E_z$",
+            "compute": lambda f: f.Ez,
+            "cmap": "RdBu_r",
+            "norm": mpl.colors.Normalize(-0.05, 0.05),
+        },
+        {
+            "name": r"$B_x$",
+            "compute": lambda f: f.Bx,
+            "cmap": "BrBG",
+            "norm": mpl.colors.Normalize(-0.05, 0.05),
+        },
+        {
+            "name": r"$B_y$",
+            "compute": lambda f: f.By,
+            "cmap": "BrBG",
+            "norm": mpl.colors.Normalize(-0.05, 0.05),
+        },
+        {
+            "name": r"$B_z$",
+            "compute": lambda f: f.Bz,
+            "cmap": "BrBG",
+            "norm": mpl.colors.Normalize(-0.05, 0.05),
+        },
+    ]
+    fig = plt.figure(figsize=(12, 5.5), dpi=300)
+    gs = fig.add_gridspec(len(quantities), 1, hspace=0.02)
+    axs = [fig.add_subplot(gs[i]) for i in range(len(quantities))]
+
+    for ax, q in zip(axs, quantities):
+        q["compute"](f.isel(t=ti)).plot(
+            ax=ax,
+            cmap=q["cmap"],
+            norm=q["norm"],
+            cbar_kwargs={"label": q["name"], "shrink": 0.8, "aspect": 10, "pad": 0.005},
+        )
+    for i, ax in enumerate(axs):
+        ax.set(aspect=1)
+        if i != 0:
+            ax.set(title=None)
+        if i != len(axs) - 1:
+            ax.set(
+                xticks=[],
+                xticklabels=[],
+                xlabel=None,
+                title=ax.get_title().split(",")[0],
+            )
+    return fig
diff --git a/setups/srpic/shocktest/shock.toml b/setups/srpic/shocktest/shock.toml
new file mode 100644
index 000000000..6c0c9a3a0
--- /dev/null
+++ b/setups/srpic/shocktest/shock.toml
@@ -0,0 +1,69 @@
+[simulation]
+  name    = "shock"
+  engine  = "srpic"
+  runtime = 50.0
+
+[grid]
+  resolution = [2048, 128]
+  extent     = [[0.0, 10.0], [-0.3125, 0.3125]]
+
+  [grid.metric]
+    metric = "minkowski"
+
+  [grid.boundaries]
+    fields    = [["CONDUCTOR", "FIXED"], ["PERIODIC"]]
+    particles = [["REFLECT", "ABSORB"], ["PERIODIC"]]
+
+[scales]
+  larmor0    = 1e-2
+  skindepth0 = 1e-2
+
+[algorithms]
+  current_filters = 8
+  fieldsolver     = "false"
+  deposit         = "false"
+
+  [algorithms.timestep]
+    CFL = 0.5
+
+[particles]
+  ppc0 = 16.0
+
+  [[particles.species]]
+    label    = "e-"
+    mass     = 1.0
+    charge   = -1.0
+    maxnpart = 1e8
+
+  [[particles.species]]
+    label    = "e+"
+    mass     = 1.0
+    charge   = 1.0
+    maxnpart = 1e8
+
+[setup]
+  drift_ux    = 0.1
+  temperature = 1e-3
+  Bmag        = 1.0
+  Btheta      = 0.0
+  Bphi        = 0.0
+  x_e         = [0.05]
+  y_e         = [0.0]
+  ux_e        = [-0.01]
+  uy_e        = [0.01]
+  uz_e        = [0.001]
+  x_i         = [0.05]
+  y_i         = [0.0]
+  ux_i        = [0.01]
+  uy_i        = [-0.01]
+  uz_i        = [-0.001]
+
+[output]
+  interval = 1
+  format   = "hdf5"
+
+  [output.fields]
+    quantities = ["N_1", "N_2", "E", "B", "T0i_1", "T0i_2", "J"]
+
+  [output.debug]
+    ghosts = true
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a41b84900..f9d921df0 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: entity [STATIC/SHARED]
 #
@@ -17,7 +18,6 @@
 #
 # * kokkos [required]
 # * plog [required]
-# * toml11 [required]
 # * ADIOS2 [optional]
 # * mpi [optional]
 # ------------------------------
diff --git a/src/archetypes/CMakeLists.txt b/src/archetypes/CMakeLists.txt
index 8e2f325af..93f8baaaa 100644
--- a/src/archetypes/CMakeLists.txt
+++ b/src/archetypes/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: ntt_archetypes [INTERFACE]
 #
@@ -24,4 +25,3 @@ target_link_libraries(ntt_archetypes INTERFACE ${libs})
 
 target_include_directories(ntt_archetypes
                            INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/../)
-
diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h
index 3884e6ced..62b9249c3 100644
--- a/src/archetypes/particle_injector.h
+++ b/src/archetypes/particle_injector.h
@@ -170,6 +170,68 @@ namespace arch {
     ~AtmosphereInjector() = default;
   };
 
+  template <SimEngine::type S, class M, in O>
+  struct MovingInjector {
+    struct TargetDensityProfile {
+      const real_t nmax, xinj, xdrift;
+
+      TargetDensityProfile(real_t xinj, real_t xdrift, real_t nmax)
+        : xinj { xinj }
+        , xdrift { xdrift }
+        , nmax { nmax } {}
+
+      Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t {
+        if constexpr ((O == in::x1) or
+                      (O == in::x2 and (M::Dim == Dim::_2D or M::Dim == Dim::_3D)) or
+                      (O == in::x3 and M::Dim == Dim::_3D)) {
+          const auto xi = x_Ph[static_cast<unsigned short>(O)];
+          // + direction
+          if (xi < xdrift or xi >= xinj) {
+            return ZERO;
+          } else {
+            if constexpr (M::CoordType == Coord::Cart) {
+              return nmax;
+            } else {
+              raise::KernelError(
+                HERE,
+                "Moving injector in +x cannot be applied for non-cartesian");
+              return ZERO;
+            }
+          }
+        } else {
+          raise::KernelError(HERE, "Wrong direction");
+          return ZERO;
+        }
+      }
+    };
+
+    using energy_dist_t  = Maxwellian<S, M>;
+    using spatial_dist_t = Replenish<S, M, TargetDensityProfile>;
+    static_assert(M::is_metric, "M must be a metric class");
+    static constexpr bool      is_nonuniform_injector { true };
+    static constexpr Dimension D { M::Dim };
+    static constexpr Coord     C { M::CoordType };
+
+    const energy_dist_t               energy_dist;
+    const TargetDensityProfile        target_density;
+    const spatial_dist_t              spatial_dist;
+    const std::pair<spidx_t, spidx_t> species;
+
+    MovingInjector(const M&                           metric,
+                   const ndfield_t<M::Dim, 6>&        density,
+                   const energy_dist_t&               energy_dist,
+                   real_t                             xinj,
+                   real_t                             xdrift,
+                   real_t                             nmax,
+                   const std::pair<spidx_t, spidx_t>& species)
+      : energy_dist { energy_dist }
+      , target_density { xinj, xdrift, nmax }
+      , spatial_dist { metric, density, 0, target_density, nmax }
+      , species { species } {}
+
+    ~MovingInjector() = default;
+  };
+
   /**
    * @brief Injects uniform number density of particles everywhere in the domain
    * @param domain Domain object
diff --git a/src/archetypes/spatial_dist.h b/src/archetypes/spatial_dist.h
index d036c0166..ad9404ea3 100644
--- a/src/archetypes/spatial_dist.h
+++ b/src/archetypes/spatial_dist.h
@@ -4,6 +4,7 @@
  * @implements
  *   - arch::SpatialDistribution<>
  *   - arch::Uniform<> : arch::SpatialDistribution<>
+ *   - arch::Piston<> : arch::SpatialDistribution<>
  *   - arch::Replenish<> : arch::SpatialDistribution<>
  * @namespace
  *   - arch::
@@ -49,6 +50,30 @@ namespace arch {
     }
   };
 
+  template <SimEngine::type S, class M>
+  struct Piston : public arch::SpatialDistribution<S, M> {
+    Piston(const M& metric, real_t xmin, real_t xmax, in piston_direction = in::x1)
+      : arch::SpatialDistribution<S, M> { metric }
+      , xmin { xmin }
+      , xmax { xmax }
+      , piston_direction { piston_direction } {}
+
+    Inline auto operator()(const coord_t<M::Dim>& x_Ph) const -> real_t override {
+      // dimentsion to fill
+      const auto fill_dim = static_cast<int>(piston_direction);
+
+      if (x_Ph[fill_dim] < xmin || x_Ph[fill_dim] > xmax) {
+        return ZERO;
+      } else {
+        return ONE;
+      }
+    }
+
+  private:
+    real_t xmin, xmax;
+    in     piston_direction;
+  };
+
   template <SimEngine::type S, class M, class T>
   struct Replenish : public SpatialDistribution<S, M> {
     using SpatialDistribution<S, M>::metric;
diff --git a/src/archetypes/tests/CMakeLists.txt b/src/archetypes/tests/CMakeLists.txt
index 694a6b4f9..9419847c5 100644
--- a/src/archetypes/tests/CMakeLists.txt
+++ b/src/archetypes/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103,C0111
 # ------------------------------
 # @brief: Generates tests for the `ntt_archetypes` module
 #
diff --git a/src/checkpoint/CMakeLists.txt b/src/checkpoint/CMakeLists.txt
index fa641bfb5..096aad690 100644
--- a/src/checkpoint/CMakeLists.txt
+++ b/src/checkpoint/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: ntt_checkpoint [STATIC/SHARED]
 #
diff --git a/src/checkpoint/tests/CMakeLists.txt b/src/checkpoint/tests/CMakeLists.txt
index 54400652e..cbfd63aa9 100644
--- a/src/checkpoint/tests/CMakeLists.txt
+++ b/src/checkpoint/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103,C0111
 # ------------------------------
 # @brief: Generates tests for the `ntt_checkpoint` module
 #
diff --git a/src/checkpoint/tests/checkpoint-mpi.cpp b/src/checkpoint/tests/checkpoint-mpi.cpp
index f97202ab1..bc6d6038a 100644
--- a/src/checkpoint/tests/checkpoint-mpi.cpp
+++ b/src/checkpoint/tests/checkpoint-mpi.cpp
@@ -144,8 +144,8 @@ auto main(int argc, char* argv[]) -> int {
       writer.saveField<Dim::_2D, 6>("em", field1);
       writer.saveField<Dim::_2D, 6>("em0", field2);
 
-      writer.savePerDomainVariable<std::size_t>("s1_npart", 1, 0, npart1);
-      writer.savePerDomainVariable<std::size_t>("s2_npart", 1, 0, npart2);
+      writer.savePerDomainVariable<npart_t>("s1_npart", 1, 0, npart1);
+      writer.savePerDomainVariable<npart_t>("s2_npart", 1, 0, npart2);
 
       writer.saveParticleQuantity<int>("s1_i1",
                                        npart1_globtot,
diff --git a/src/checkpoint/tests/checkpoint-nompi.cpp b/src/checkpoint/tests/checkpoint-nompi.cpp
index 23dbd8871..be5e97e24 100644
--- a/src/checkpoint/tests/checkpoint-nompi.cpp
+++ b/src/checkpoint/tests/checkpoint-nompi.cpp
@@ -105,8 +105,8 @@ auto main(int argc, char* argv[]) -> int {
       writer.saveField<Dim::_3D, 6>("em", field1);
       writer.saveField<Dim::_3D, 6>("em0", field2);
 
-      writer.savePerDomainVariable<std::size_t>("s1_npart", 1, 0, npart1);
-      writer.savePerDomainVariable<std::size_t>("s2_npart", 1, 0, npart2);
+      writer.savePerDomainVariable<npart_t>("s1_npart", 1, 0, npart1);
+      writer.savePerDomainVariable<npart_t>("s2_npart", 1, 0, npart2);
 
       writer.saveParticleQuantity<int>("s1_i1", npart1, 0, npart1, i1);
       writer.saveParticleQuantity<real_t>("s1_ux1", npart1, 0, npart1, u1);
diff --git a/src/checkpoint/writer.cpp b/src/checkpoint/writer.cpp
index a70d8de09..fe77d3b56 100644
--- a/src/checkpoint/writer.cpp
+++ b/src/checkpoint/writer.cpp
@@ -274,9 +274,9 @@ namespace checkpoint {
                                                       std::size_t,
                                                       double);
   template void Writer::savePerDomainVariable<npart_t>(const std::string&,
-                                                           std::size_t,
-                                                           std::size_t,
-                                                           npart_t);
+                                                       std::size_t,
+                                                       std::size_t,
+                                                       npart_t);
 
   template void Writer::saveField<Dim::_1D, 3>(const std::string&,
                                                const ndfield_t<Dim::_1D, 3>&);
diff --git a/src/engines/CMakeLists.txt b/src/engines/CMakeLists.txt
index 6da2f4efd..4cef18630 100644
--- a/src/engines/CMakeLists.txt
+++ b/src/engines/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: ntt_engines [STATIC/SHARED]
 #
@@ -25,7 +26,6 @@
 #
 # * kokkos [required]
 # * plog [required]
-# * toml11 [required]
 # * adios2 [optional]
 # * hdf5 [optional]
 # * mpi [optional]
@@ -39,7 +39,7 @@ add_library(ntt_engines ${SOURCES})
 set(libs ntt_global ntt_framework ntt_metrics ntt_archetypes ntt_kernels
          ntt_pgen)
 if(${output})
-  list(APPEND libs ntt_output hdf5::hdf5)
+  list(APPEND libs ntt_output)
 endif()
 add_dependencies(ntt_engines ${libs})
 target_link_libraries(ntt_engines PUBLIC ${libs})
diff --git a/src/engines/engine_printer.cpp b/src/engines/engine_printer.cpp
index f94715d09..eb8ff402d 100644
--- a/src/engines/engine_printer.cpp
+++ b/src/engines/engine_printer.cpp
@@ -20,7 +20,6 @@
 #endif
 
 #if defined(OUTPUT_ENABLED)
-  #include <H5public.h>
   #include <adios2.h>
 #endif
 
@@ -188,18 +187,11 @@ namespace ntt {
                                                 KOKKOS_VERSION % 100);
 
 #if defined(OUTPUT_ENABLED)
-        unsigned h5_major, h5_minor, h5_release;
-        H5get_libversion(&h5_major, &h5_minor, &h5_release);
-        const std::string hdf5_version   = fmt::format("%d.%d.%d",
-                                                     h5_major,
-                                                     h5_minor,
-                                                     h5_release);
         const std::string adios2_version = fmt::format("%d.%d.%d",
                                                        ADIOS2_VERSION / 10000,
                                                        ADIOS2_VERSION / 100 % 100,
                                                        ADIOS2_VERSION % 100);
 #else // not OUTPUT_ENABLED
-        const std::string hdf5_version   = "OFF";
         const std::string adios2_version = "OFF";
 #endif
 
@@ -217,7 +209,6 @@ namespace ntt {
         add_param(report, 4, "CXX", "%s [%s]", ccx.c_str(), cpp_standard.c_str());
         add_param(report, 4, "CUDA", "%s", cuda_version.c_str());
         add_param(report, 4, "MPI", "%s", mpi_version.c_str());
-        add_param(report, 4, "HDF5", "%s", hdf5_version.c_str());
         add_param(report, 4, "Kokkos", "%s", kokkos_version.c_str());
         add_param(report, 4, "ADIOS2", "%s", adios2_version.c_str());
         add_param(report, 4, "Precision", "%s", precision);
diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp
index 0a9cc311b..62cd23ad9 100644
--- a/src/engines/srpic.hpp
+++ b/src/engines/srpic.hpp
@@ -596,6 +596,10 @@ namespace ntt {
           if (domain.mesh.flds_bc_in(direction) == FldsBC::FIXED) {
             FixedFieldsIn(direction, domain, tags);
           }
+        } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::CONDUCTOR) {
+          if (domain.mesh.flds_bc_in(direction) == FldsBC::CONDUCTOR) {
+            PerfectConductorFieldsIn(direction, domain, tags);
+          }
         } else if (m_metadomain.mesh().flds_bc_in(direction) == FldsBC::CUSTOM) {
           if (domain.mesh.flds_bc_in(direction) == FldsBC::CUSTOM) {
             CustomFieldsIn(direction, domain, tags);
@@ -647,7 +651,7 @@ namespace ntt {
       tuple_t<ncells_t, M::Dim> range_min { 0 };
       tuple_t<ncells_t, M::Dim> range_max { 0 };
 
-      for (unsigned short d { 0 }; d < M::Dim; ++d) {
+      for (auto d { 0u }; d < M::Dim; ++d) {
         range_min[d] = intersect_range[d].first;
         range_max[d] = intersect_range[d].second;
       }
@@ -704,28 +708,32 @@ namespace ntt {
       /**
        * axis boundaries
        */
-      raise::ErrorIf(M::CoordType == Coord::Cart,
-                     "Invalid coordinate type for axis BCs",
-                     HERE);
-      raise::ErrorIf(direction.get_dim() != in::x2,
-                     "Invalid axis direction, should be x2",
-                     HERE);
-      const auto i2_min = domain.mesh.i_min(in::x2);
-      const auto i2_max = domain.mesh.i_max(in::x2);
-      if (direction.get_sign() < 0) {
-        Kokkos::parallel_for(
-          "AxisBCFields",
-          domain.mesh.n_all(in::x1),
-          kernel::bc::AxisBoundaries_kernel<M::Dim, false>(domain.fields.em,
-                                                           i2_min,
-                                                           tags));
+      if constexpr (M::CoordType != Coord::Cart) {
+        raise::ErrorIf(direction.get_dim() != in::x2,
+                       "Invalid axis direction, should be x2",
+                       HERE);
+        const auto i2_min = domain.mesh.i_min(in::x2);
+        const auto i2_max = domain.mesh.i_max(in::x2);
+        if (direction.get_sign() < 0) {
+          Kokkos::parallel_for(
+            "AxisBCFields",
+            domain.mesh.n_all(in::x1),
+            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::bc::AxisBoundaries_kernel<M::Dim, true>(domain.fields.em,
+                                                            i2_max,
+                                                            tags));
+        }
       } else {
-        Kokkos::parallel_for(
-          "AxisBCFields",
-          domain.mesh.n_all(in::x1),
-          kernel::bc::AxisBoundaries_kernel<M::Dim, true>(domain.fields.em,
-                                                          i2_max,
-                                                          tags));
+        (void)direction;
+        (void)domain;
+        (void)tags;
+        raise::Error("Invalid coordinate type for axis BCs", HERE);
       }
     }
 
@@ -830,10 +838,114 @@ namespace ntt {
           }
         }
       } else {
+        (void)direction;
+        (void)domain;
+        (void)tags;
         raise::Error("Fixed fields not present (both const and non-const)", HERE);
       }
     }
 
+    void PerfectConductorFieldsIn(dir::direction_t<M::Dim> direction,
+                                  domain_t&                domain,
+                                  BCTags                   tags) {
+      /**
+       * perfect conductor field boundaries
+       */
+      if constexpr (M::CoordType != Coord::Cart) {
+        (void)direction;
+        (void)domain;
+        (void)tags;
+        raise::Error(
+          "Perfect conductor BCs only applicable to cartesian coordinates",
+          HERE);
+      } else {
+        const auto sign = direction.get_sign();
+        const auto dim  = direction.get_dim();
+
+        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) {
+            xi_min.push_back(0);
+            xi_max.push_back((sign < 0) ? (N_GHOSTS + 1) : 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);
+
+        range_t<M::Dim> range;
+        if constexpr (M::Dim == Dim::_1D) {
+          range = CreateRangePolicy<M::Dim>({ xi_min[0] }, { xi_max[0] });
+        } else if constexpr (M::Dim == Dim::_2D) {
+          range = CreateRangePolicy<M::Dim>({ xi_min[0], xi_min[1] },
+                                            { xi_max[0], xi_max[1] });
+        } else if constexpr (M::Dim == Dim::_3D) {
+          range = CreateRangePolicy<M::Dim>({ xi_min[0], xi_min[1], xi_min[2] },
+                                            { xi_max[0], xi_max[1], xi_max[2] });
+        } else {
+          raise::Error("Invalid dimension", HERE);
+        }
+
+        if (dim == in::x1) {
+          if (sign > 0) {
+            Kokkos::parallel_for(
+              "ConductorFields",
+              range,
+              kernel::bc::ConductorBoundaries_kernel<M::Dim, in::x1, true>(
+                domain.fields.em,
+                tags));
+          } else {
+            Kokkos::parallel_for(
+              "ConductorFields",
+              range,
+              kernel::bc::ConductorBoundaries_kernel<M::Dim, in::x1, false>(
+                domain.fields.em,
+                tags));
+          }
+        } else if (dim == in::x2) {
+          if (sign > 0) {
+            Kokkos::parallel_for(
+              "ConductorFields",
+              range,
+              kernel::bc::ConductorBoundaries_kernel<M::Dim, in::x2, true>(
+                domain.fields.em,
+                tags));
+          } else {
+            Kokkos::parallel_for(
+              "ConductorFields",
+              range,
+              kernel::bc::ConductorBoundaries_kernel<M::Dim, in::x2, false>(
+                domain.fields.em,
+                tags));
+          }
+        } else {
+          if (sign > 0) {
+            Kokkos::parallel_for(
+              "ConductorFields",
+              range,
+              kernel::bc::ConductorBoundaries_kernel<M::Dim, in::x3, true>(
+                domain.fields.em,
+                tags));
+          } else {
+            Kokkos::parallel_for(
+              "ConductorFields",
+              range,
+              kernel::bc::ConductorBoundaries_kernel<M::Dim, in::x3, false>(
+                domain.fields.em,
+                tags));
+          }
+        }
+      }
+    }
+
     void AtmosphereFieldsIn(dir::direction_t<M::Dim> direction,
                             domain_t&                domain,
                             BCTags                   tags) {
@@ -862,15 +974,15 @@ namespace ntt {
           return;
         }
         const auto intersect_range = domain.mesh.ExtentToRange(box, incl_ghosts);
-        tuple_t<ncells_t, M::Dim> range_min { 0 };
-        tuple_t<ncells_t, M::Dim> range_max { 0 };
+        tuple_t<std::size_t, M::Dim> range_min { 0 };
+        tuple_t<std::size_t, M::Dim> range_max { 0 };
 
         for (unsigned short d { 0 }; d < M::Dim; ++d) {
           range_min[d] = intersect_range[d].first;
           range_max[d] = intersect_range[d].second;
         }
-        auto     atm_fields = m_pgen.AtmFields(time);
-        ncells_t il_edge;
+        auto        atm_fields = m_pgen.AtmFields(time);
+        std::size_t il_edge;
         if (sign > 0) {
           il_edge = range_min[dd] - N_GHOSTS;
         } else {
@@ -955,6 +1067,9 @@ namespace ntt {
           raise::Error("Invalid dimension", HERE);
         }
       } else {
+        (void)direction;
+        (void)domain;
+        (void)tags;
         raise::Error("Atm fields not implemented in PGEN for atmosphere BCs", HERE);
       }
     }
diff --git a/src/framework/CMakeLists.txt b/src/framework/CMakeLists.txt
index 8802f696b..4c407fb0c 100644
--- a/src/framework/CMakeLists.txt
+++ b/src/framework/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: ntt_framework [STATIC/SHARED]
 #
@@ -28,7 +29,6 @@
 #
 # * kokkos [required]
 # * plog [required]
-# * toml11 [required]
 # * ADIOS2 [optional]
 # * mpi [optional]
 # ------------------------------
diff --git a/src/framework/containers/particles.cpp b/src/framework/containers/particles.cpp
index 048d57cde..2f59a004d 100644
--- a/src/framework/containers/particles.cpp
+++ b/src/framework/containers/particles.cpp
@@ -84,7 +84,7 @@ namespace ntt {
       rangeActiveParticles(),
       Lambda(index_t p) {
         auto npptag_acc = npptag_scat.access();
-        if (this_tag(p) < 0 || this_tag(p) >= num_tags) {
+        if (this_tag(p) < 0 || this_tag(p) >= static_cast<short>(num_tags)) {
           raise::KernelError(HERE, "Invalid tag value");
         }
         npptag_acc(this_tag(p)) += 1;
@@ -144,9 +144,8 @@ namespace ntt {
 
   template <Dimension D, Coord::type C>
   void Particles<D, C>::RemoveDead() {
-    const auto n_part  = npart();
-    npart_t    n_alive = 0, n_dead = 0;
-    auto&      this_tag = tag;
+    npart_t n_alive = 0, n_dead = 0;
+    auto&   this_tag = tag;
 
     Kokkos::parallel_reduce(
       "CountDeadAlive",
@@ -218,13 +217,13 @@ namespace ntt {
 
     Kokkos::Experimental::fill(
       "TagAliveParticles",
-      AccelExeSpace(),
+      Kokkos::DefaultExecutionSpace(),
       Kokkos::subview(this_tag, std::make_pair(static_cast<npart_t>(0), n_alive)),
       ParticleTag::alive);
 
     Kokkos::Experimental::fill(
       "TagDeadParticles",
-      AccelExeSpace(),
+      Kokkos::DefaultExecutionSpace(),
       Kokkos::subview(this_tag, std::make_pair(n_alive, n_alive + n_dead)),
       ParticleTag::dead);
 
diff --git a/src/framework/domain/comm_mpi.hpp b/src/framework/domain/comm_mpi.hpp
index 8c6e532de..d9783bdb6 100644
--- a/src/framework/domain/comm_mpi.hpp
+++ b/src/framework/domain/comm_mpi.hpp
@@ -83,7 +83,7 @@ namespace comm {
                                  (long int)(send_slice[0].first);
           Kokkos::parallel_for(
             "CommunicateField-extract",
-            Kokkos::MDRangePolicy<Kokkos::Rank<2>, AccelExeSpace>(
+            Kokkos::MDRangePolicy<Kokkos::Rank<2>, Kokkos::DefaultExecutionSpace>(
               { recv_slice[0].first, comps.first },
               { recv_slice[0].second, comps.second }),
             Lambda(index_t i1, index_t ci) {
@@ -96,7 +96,7 @@ namespace comm {
                                  (long int)(send_slice[1].first);
           Kokkos::parallel_for(
             "CommunicateField-extract",
-            Kokkos::MDRangePolicy<Kokkos::Rank<3>, AccelExeSpace>(
+            Kokkos::MDRangePolicy<Kokkos::Rank<3>, Kokkos::DefaultExecutionSpace>(
               { recv_slice[0].first, recv_slice[1].first, comps.first },
               { recv_slice[0].second, recv_slice[1].second, comps.second }),
             Lambda(index_t i1, index_t i2, index_t ci) {
@@ -111,7 +111,7 @@ namespace comm {
                                  (long int)(send_slice[2].first);
           Kokkos::parallel_for(
             "CommunicateField-extract",
-            Kokkos::MDRangePolicy<Kokkos::Rank<4>, AccelExeSpace>(
+            Kokkos::MDRangePolicy<Kokkos::Rank<4>, Kokkos::DefaultExecutionSpace>(
               { recv_slice[0].first,
                 recv_slice[1].first,
                 recv_slice[2].first,
@@ -220,7 +220,6 @@ namespace comm {
 
       if (recv_rank >= 0) {
 
-        // !TODO: perhaps directly recv to the fld?
         if (not additive) {
           if constexpr (D == Dim::_1D) {
             Kokkos::deep_copy(Kokkos::subview(fld, recv_slice[0], comps), recv_fld);
@@ -239,7 +238,7 @@ namespace comm {
             const auto offset_c  = comps.first;
             Kokkos::parallel_for(
               "CommunicateField-extract",
-              Kokkos::MDRangePolicy<Kokkos::Rank<2>, AccelExeSpace>(
+              Kokkos::MDRangePolicy<Kokkos::Rank<2>, Kokkos::DefaultExecutionSpace>(
                 { recv_slice[0].first, comps.first },
                 { recv_slice[0].second, comps.second }),
               Lambda(index_t i1, index_t ci) {
@@ -251,7 +250,7 @@ namespace comm {
             const auto offset_c  = comps.first;
             Kokkos::parallel_for(
               "CommunicateField-extract",
-              Kokkos::MDRangePolicy<Kokkos::Rank<3>, AccelExeSpace>(
+              Kokkos::MDRangePolicy<Kokkos::Rank<3>, Kokkos::DefaultExecutionSpace>(
                 { recv_slice[0].first, recv_slice[1].first, comps.first },
                 { recv_slice[0].second, recv_slice[1].second, comps.second }),
               Lambda(index_t i1, index_t i2, index_t ci) {
@@ -266,7 +265,7 @@ namespace comm {
             const auto offset_c  = comps.first;
             Kokkos::parallel_for(
               "CommunicateField-extract",
-              Kokkos::MDRangePolicy<Kokkos::Rank<4>, AccelExeSpace>(
+              Kokkos::MDRangePolicy<Kokkos::Rank<4>, Kokkos::DefaultExecutionSpace>(
                 { recv_slice[0].first,
                   recv_slice[1].first,
                   recv_slice[2].first,
@@ -336,10 +335,8 @@ namespace comm {
     const unsigned short NPLDS   = species.npld();
 
     // buffers to store recv data
-    const auto       npart_alive = npptag_vec[ParticleTag::alive];
-    const auto       npart_dead  = npptag_vec[ParticleTag::dead];
-    const auto       npart_send  = outgoing_indices.extent(0) - npart_dead;
-    const auto       npart_recv  = std::accumulate(npptag_recv_vec.begin(),
+    const auto       npart_dead = npptag_vec[ParticleTag::dead];
+    const auto       npart_recv = std::accumulate(npptag_recv_vec.begin(),
                                             npptag_recv_vec.end(),
                                             static_cast<npart_t>(0));
     array_t<int*>    recv_buff_int { "recv_buff_int", npart_recv * NINTS };
diff --git a/src/framework/domain/comm_nompi.hpp b/src/framework/domain/comm_nompi.hpp
index 197d336fa..b477ac176 100644
--- a/src/framework/domain/comm_nompi.hpp
+++ b/src/framework/domain/comm_nompi.hpp
@@ -70,7 +70,7 @@ namespace comm {
                                  (long int)(send_slice[0].first);
           Kokkos::parallel_for(
             "CommunicateField-extract",
-            Kokkos::MDRangePolicy<Kokkos::Rank<2>, AccelExeSpace>(
+            Kokkos::MDRangePolicy<Kokkos::Rank<2>, Kokkos::DefaultExecutionSpace>(
               { recv_slice[0].first, comps.first },
               { recv_slice[0].second, comps.second }),
             Lambda(index_t i1, index_t ci) {
@@ -83,7 +83,7 @@ namespace comm {
                                  (long int)(send_slice[1].first);
           Kokkos::parallel_for(
             "CommunicateField-extract",
-            Kokkos::MDRangePolicy<Kokkos::Rank<3>, AccelExeSpace>(
+            Kokkos::MDRangePolicy<Kokkos::Rank<3>, Kokkos::DefaultExecutionSpace>(
               { recv_slice[0].first, recv_slice[1].first, comps.first },
               { recv_slice[0].second, recv_slice[1].second, comps.second }),
             Lambda(index_t i1, index_t i2, index_t ci) {
@@ -98,7 +98,7 @@ namespace comm {
                                  (long int)(send_slice[2].first);
           Kokkos::parallel_for(
             "CommunicateField-extract",
-            Kokkos::MDRangePolicy<Kokkos::Rank<4>, AccelExeSpace>(
+            Kokkos::MDRangePolicy<Kokkos::Rank<4>, Kokkos::DefaultExecutionSpace>(
               { recv_slice[0].first,
                 recv_slice[1].first,
                 recv_slice[2].first,
diff --git a/src/framework/domain/communications.cpp b/src/framework/domain/communications.cpp
index 946756f49..cf7e04974 100644
--- a/src/framework/domain/communications.cpp
+++ b/src/framework/domain/communications.cpp
@@ -36,10 +36,10 @@ namespace ntt {
   using comm_params_t = std::pair<address_t, std::vector<range_tuple_t>>;
 
   template <SimEngine::type S, class M>
-  auto GetSendRecvRanks(
-    Metadomain<S, M>*        metadomain,
-    Domain<S, M>&            domain,
-    dir::direction_t<M::Dim> direction) -> std::pair<address_t, address_t> {
+  auto GetSendRecvRanks(Metadomain<S, M>*        metadomain,
+                        Domain<S, M>&            domain,
+                        dir::direction_t<M::Dim> direction)
+    -> std::pair<address_t, address_t> {
     Domain<S, M>* send_to_nghbr_ptr   = nullptr;
     Domain<S, M>* recv_from_nghbr_ptr = nullptr;
     // set pointers to the correct send/recv domains
@@ -119,11 +119,11 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
-  auto GetSendRecvParams(
-    Metadomain<S, M>*        metadomain,
-    Domain<S, M>&            domain,
-    dir::direction_t<M::Dim> direction,
-    bool synchronize) -> std::pair<comm_params_t, comm_params_t> {
+  auto GetSendRecvParams(Metadomain<S, M>*        metadomain,
+                         Domain<S, M>&            domain,
+                         dir::direction_t<M::Dim> direction,
+                         bool                     synchronize)
+    -> std::pair<comm_params_t, comm_params_t> {
     const auto [send_indrank,
                 recv_indrank] = GetSendRecvRanks(metadomain, domain, direction);
     const auto [send_ind, send_rank] = send_indrank;
@@ -506,8 +506,7 @@ namespace ntt {
       const auto npart_dead          = npptag_vec[ParticleTag::dead];
       const auto npart_alive         = npptag_vec[ParticleTag::alive];
 
-      const auto npart       = species.npart();
-      const auto npart_holes = npart - npart_alive;
+      const auto npart = species.npart();
 
       // # of particles to receive per each tag (direction)
       std::vector<npart_t> npptag_recv_vec(ntags - 2, 0);
diff --git a/src/framework/domain/grid.cpp b/src/framework/domain/grid.cpp
index baa23fb5c..c022184b1 100644
--- a/src/framework/domain/grid.cpp
+++ b/src/framework/domain/grid.cpp
@@ -85,7 +85,6 @@ namespace ntt {
     return CreateRangePolicy<D>(imin, imax);
   }
 
-  // !TODO: too ugly, implement a better solution (combine with device)
   template <Dimension D>
   auto Grid<D>::rangeCellsOnHost(
     const box_region_t<D>& region) const -> range_h_t<D> {
diff --git a/src/framework/domain/metadomain.cpp b/src/framework/domain/metadomain.cpp
index 10e3e4fb0..8952d417d 100644
--- a/src/framework/domain/metadomain.cpp
+++ b/src/framework/domain/metadomain.cpp
@@ -46,7 +46,7 @@ namespace ntt {
 #if defined(MPI_ENABLED)
     MPI_Comm_size(MPI_COMM_WORLD, &g_mpi_size);
     MPI_Comm_rank(MPI_COMM_WORLD, &g_mpi_rank);
-    raise::ErrorIf(global_ndomains != g_mpi_size,
+    raise::ErrorIf(global_ndomains != (unsigned int)g_mpi_size,
                    "Exactly 1 domain per MPI rank is allowed",
                    HERE);
 #endif
@@ -234,7 +234,6 @@ namespace ntt {
 
   template <SimEngine::type S, class M>
   void Metadomain<S, M>::redefineBoundaries() {
-    // !TODO: not setting CommBC for now
     for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) {
       // offset of the subdomain[idx]
       auto&      current_domain = g_subdomains[idx];
diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp
index e3bddfd70..556d9f547 100644
--- a/src/framework/parameters.cpp
+++ b/src/framework/parameters.cpp
@@ -31,10 +31,10 @@
 namespace ntt {
 
   template <typename M>
-  auto get_dx0_V0(
-    const std::vector<ncells_t>&         resolution,
-    const boundaries_t<real_t>&          extent,
-    const std::map<std::string, real_t>& params) -> std::pair<real_t, real_t> {
+  auto get_dx0_V0(const std::vector<ncells_t>&         resolution,
+                  const boundaries_t<real_t>&          extent,
+                  const std::map<std::string, real_t>& params)
+    -> std::pair<real_t, real_t> {
     const auto      metric = M(resolution, extent, params);
     const auto      dx0    = metric.dxMin();
     coord_t<M::Dim> x_corner { ZERO };
@@ -601,6 +601,8 @@ namespace ntt {
         toml::find_or(toml_data, "diagnostics", "blocking_timers", false));
     set("diagnostics.colored_stdout",
         toml::find_or(toml_data, "diagnostics", "colored_stdout", false));
+    set("diagnostics.log_level",
+        toml::find_or(toml_data, "diagnostics", "log_level", defaults::diag::log_level));
 
     /* inferred variables --------------------------------------------------- */
     // fields/particle boundaries
@@ -839,7 +841,7 @@ namespace ntt {
 
   void SimulationParams::setSetupParams(const toml::value& toml_data) {
     /* [setup] -------------------------------------------------------------- */
-    const auto& setup = toml::find_or(toml_data, "setup", toml::table {});
+    const auto setup = toml::find_or(toml_data, "setup", toml::table {});
     for (const auto& [key, val] : setup) {
       if (val.is_boolean()) {
         set("setup." + key, (bool)(val.as_boolean()));
diff --git a/src/framework/simulation.cpp b/src/framework/simulation.cpp
index bea50ff09..6865140c9 100644
--- a/src/framework/simulation.cpp
+++ b/src/framework/simulation.cpp
@@ -32,7 +32,12 @@ namespace ntt {
 
     const auto raw_params = toml::parse(inputfname);
     const auto sim_name = toml::find<std::string>(raw_params, "simulation", "name");
-    logger::initPlog<files::LogFile, files::InfoFile, files::ErrFile>(sim_name);
+    const auto log_level = toml::find_or<std::string>(raw_params,
+                                                      "diagnostics",
+                                                      "log_level",
+                                                      defaults::diag::log_level);
+    logger::initPlog<files::LogFile, files::InfoFile, files::ErrFile>(sim_name,
+                                                                      log_level);
 
     m_requested_engine = SimEngine::pick(
       fmt::toLower(toml::find<std::string>(raw_params, "simulation", "engine")).c_str());
@@ -48,8 +53,6 @@ namespace ntt {
                    HERE);
     m_requested_dimension = static_cast<Dimension>(res.size());
 
-    // !TODO: when mixing checkpoint metadata with input,
-    // ... need to properly take care of the diffs
     m_params.setRawData(raw_params);
     timestep_t checkpoint_step = 0;
     if (is_resuming) {
diff --git a/src/framework/tests/CMakeLists.txt b/src/framework/tests/CMakeLists.txt
index ce188e9f1..92e327d80 100644
--- a/src/framework/tests/CMakeLists.txt
+++ b/src/framework/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103,C0111
 # ------------------------------
 # @brief: Generates tests for the `ntt_framework` module
 #
@@ -5,7 +6,6 @@
 #
 # * kokkos [required]
 # * plog [required]
-# * toml11 [required]
 # * mpi [optional]
 # * adios2 [optional]
 #
diff --git a/src/global/CMakeLists.txt b/src/global/CMakeLists.txt
index 97946f059..7546b6a98 100644
--- a/src/global/CMakeLists.txt
+++ b/src/global/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: ntt_global [STATIC/SHARED]
 #
diff --git a/src/global/arch/kokkos_aliases.cpp b/src/global/arch/kokkos_aliases.cpp
index d6622b0e5..25397af8e 100644
--- a/src/global/arch/kokkos_aliases.cpp
+++ b/src/global/arch/kokkos_aliases.cpp
@@ -9,73 +9,75 @@ auto CreateParticleRangePolicy(npart_t p1, npart_t p2) -> range_t<Dim::_1D> {
 }
 
 template <>
-auto CreateRangePolicy<Dim::_1D>(
-  const tuple_t<ncells_t, Dim::_1D>& i1,
-  const tuple_t<ncells_t, Dim::_1D>& i2) -> range_t<Dim::_1D> {
+auto CreateRangePolicy<Dim::_1D>(const tuple_t<ncells_t, Dim::_1D>& i1,
+                                 const tuple_t<ncells_t, Dim::_1D>& i2)
+  -> range_t<Dim::_1D> {
   index_t i1min = i1[0];
   index_t i1max = i2[0];
-  return Kokkos::RangePolicy<AccelExeSpace>(i1min, i1max);
+  return Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(i1min, i1max);
 }
 
 template <>
-auto CreateRangePolicy<Dim::_2D>(
-  const tuple_t<ncells_t, Dim::_2D>& i1,
-  const tuple_t<ncells_t, Dim::_2D>& i2) -> range_t<Dim::_2D> {
+auto CreateRangePolicy<Dim::_2D>(const tuple_t<ncells_t, Dim::_2D>& i1,
+                                 const tuple_t<ncells_t, Dim::_2D>& i2)
+  -> range_t<Dim::_2D> {
   index_t i1min = i1[0];
   index_t i1max = i2[0];
   index_t i2min = i1[1];
   index_t i2max = i2[1];
-  return Kokkos::MDRangePolicy<Kokkos::Rank<2>, AccelExeSpace>({ i1min, i2min },
-                                                               { i1max, i2max });
+  return Kokkos::MDRangePolicy<Kokkos::Rank<2>, Kokkos::DefaultExecutionSpace>(
+    { i1min, i2min },
+    { i1max, i2max });
 }
 
 template <>
-auto CreateRangePolicy<Dim::_3D>(
-  const tuple_t<ncells_t, Dim::_3D>& i1,
-  const tuple_t<ncells_t, Dim::_3D>& i2) -> range_t<Dim::_3D> {
+auto CreateRangePolicy<Dim::_3D>(const tuple_t<ncells_t, Dim::_3D>& i1,
+                                 const tuple_t<ncells_t, Dim::_3D>& i2)
+  -> range_t<Dim::_3D> {
   index_t i1min = i1[0];
   index_t i1max = i2[0];
   index_t i2min = i1[1];
   index_t i2max = i2[1];
   index_t i3min = i1[2];
   index_t i3max = i2[2];
-  return Kokkos::MDRangePolicy<Kokkos::Rank<3>, AccelExeSpace>(
+  return Kokkos::MDRangePolicy<Kokkos::Rank<3>, Kokkos::DefaultExecutionSpace>(
     { i1min, i2min, i3min },
     { i1max, i2max, i3max });
 }
 
 template <>
-auto CreateRangePolicyOnHost<Dim::_1D>(
-  const tuple_t<ncells_t, Dim::_1D>& i1,
-  const tuple_t<ncells_t, Dim::_1D>& i2) -> range_h_t<Dim::_1D> {
+auto CreateRangePolicyOnHost<Dim::_1D>(const tuple_t<ncells_t, Dim::_1D>& i1,
+                                       const tuple_t<ncells_t, Dim::_1D>& i2)
+  -> range_h_t<Dim::_1D> {
   index_t i1min = i1[0];
   index_t i1max = i2[0];
-  return Kokkos::RangePolicy<HostExeSpace>(i1min, i1max);
+  return Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(i1min, i1max);
 }
 
 template <>
-auto CreateRangePolicyOnHost<Dim::_2D>(
-  const tuple_t<ncells_t, Dim::_2D>& i1,
-  const tuple_t<ncells_t, Dim::_2D>& i2) -> range_h_t<Dim::_2D> {
+auto CreateRangePolicyOnHost<Dim::_2D>(const tuple_t<ncells_t, Dim::_2D>& i1,
+                                       const tuple_t<ncells_t, Dim::_2D>& i2)
+  -> range_h_t<Dim::_2D> {
   index_t i1min = i1[0];
   index_t i1max = i2[0];
   index_t i2min = i1[1];
   index_t i2max = i2[1];
-  return Kokkos::MDRangePolicy<Kokkos::Rank<2>, HostExeSpace>({ i1min, i2min },
-                                                              { i1max, i2max });
+  return Kokkos::MDRangePolicy<Kokkos::Rank<2>, Kokkos::DefaultHostExecutionSpace>(
+    { i1min, i2min },
+    { i1max, i2max });
 }
 
 template <>
-auto CreateRangePolicyOnHost<Dim::_3D>(
-  const tuple_t<ncells_t, Dim::_3D>& i1,
-  const tuple_t<ncells_t, Dim::_3D>& i2) -> range_h_t<Dim::_3D> {
+auto CreateRangePolicyOnHost<Dim::_3D>(const tuple_t<ncells_t, Dim::_3D>& i1,
+                                       const tuple_t<ncells_t, Dim::_3D>& i2)
+  -> range_h_t<Dim::_3D> {
   index_t i1min = i1[0];
   index_t i1max = i2[0];
   index_t i2min = i1[1];
   index_t i2max = i2[1];
   index_t i3min = i1[2];
   index_t i3max = i2[2];
-  return Kokkos::MDRangePolicy<Kokkos::Rank<3>, HostExeSpace>(
+  return Kokkos::MDRangePolicy<Kokkos::Rank<3>, Kokkos::DefaultHostExecutionSpace>(
     { i1min, i2min, i3min },
     { i1max, i2max, i3max });
 }
diff --git a/src/global/arch/kokkos_aliases.h b/src/global/arch/kokkos_aliases.h
index e1a759bed..712fc6eff 100644
--- a/src/global/arch/kokkos_aliases.h
+++ b/src/global/arch/kokkos_aliases.h
@@ -34,7 +34,7 @@
 namespace math = Kokkos;
 
 template <typename T>
-using array_t = Kokkos::View<T, AccelMemSpace>;
+using array_t = Kokkos::View<T>;
 
 // Array mirror alias of arbitrary type
 template <typename T>
@@ -174,17 +174,17 @@ namespace kokkos_aliases_hidden {
 
   template <>
   struct range_impl<Dim::_1D> {
-    using type = Kokkos::RangePolicy<AccelExeSpace>;
+    using type = Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>;
   };
 
   template <>
   struct range_impl<Dim::_2D> {
-    using type = Kokkos::MDRangePolicy<Kokkos::Rank<2>, AccelExeSpace>;
+    using type = Kokkos::MDRangePolicy<Kokkos::Rank<2>, Kokkos::DefaultExecutionSpace>;
   };
 
   template <>
   struct range_impl<Dim::_3D> {
-    using type = Kokkos::MDRangePolicy<Kokkos::Rank<3>, AccelExeSpace>;
+    using type = Kokkos::MDRangePolicy<Kokkos::Rank<3>, Kokkos::DefaultExecutionSpace>;
   };
 } // namespace kokkos_aliases_hidden
 
@@ -201,17 +201,17 @@ namespace kokkos_aliases_hidden {
 
   template <>
   struct range_h_impl<Dim::_1D> {
-    using type = Kokkos::RangePolicy<HostExeSpace>;
+    using type = Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>;
   };
 
   template <>
   struct range_h_impl<Dim::_2D> {
-    using type = Kokkos::MDRangePolicy<Kokkos::Rank<2>, HostExeSpace>;
+    using type = Kokkos::MDRangePolicy<Kokkos::Rank<2>, Kokkos::DefaultHostExecutionSpace>;
   };
 
   template <>
   struct range_h_impl<Dim::_3D> {
-    using type = Kokkos::MDRangePolicy<Kokkos::Rank<3>, HostExeSpace>;
+    using type = Kokkos::MDRangePolicy<Kokkos::Rank<3>, Kokkos::DefaultHostExecutionSpace>;
   };
 
 } // namespace kokkos_aliases_hidden
@@ -234,8 +234,8 @@ auto CreateParticleRangePolicy(npart_t, npart_t) -> range_t<Dim::_1D>;
  * @returns Kokkos::RangePolicy or Kokkos::MDRangePolicy in the accelerator execution space.
  */
 template <Dimension D>
-auto CreateRangePolicy(const tuple_t<ncells_t, D>&,
-                       const tuple_t<ncells_t, D>&) -> range_t<D>;
+auto CreateRangePolicy(const tuple_t<ncells_t, D>&, const tuple_t<ncells_t, D>&)
+  -> range_t<D>;
 
 /**
  * @brief Function template for generating ND Kokkos range policy on the host.
@@ -249,8 +249,8 @@ auto CreateRangePolicyOnHost(const tuple_t<ncells_t, D>&,
                              const tuple_t<ncells_t, D>&) -> range_h_t<D>;
 
 // Random number pool/generator type alias
-using random_number_pool_t = Kokkos::Random_XorShift1024_Pool<AccelExeSpace>;
-using random_generator_t   = typename random_number_pool_t::generator_type;
+using random_number_pool_t = Kokkos::Random_XorShift1024_Pool<Kokkos::DefaultExecutionSpace>;
+using random_generator_t = typename random_number_pool_t::generator_type;
 
 // Random number generator functions
 template <typename T>
diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h
index 4cde4fca5..65cc63cf8 100644
--- a/src/global/arch/traits.h
+++ b/src/global/arch/traits.h
@@ -112,6 +112,18 @@ namespace traits {
     template <typename T>
     using fix_fields_const_t = decltype(&T::FixFieldsConst);
 
+    template <typename T>
+    using perfect_conductor_fields_t = decltype(&T::PerfectConductorFields);
+
+    template <typename T>
+    using perfect_conductor_fields_const_t = decltype(&T::PerfectConductorFieldsConst);
+
+    template <typename T>
+    using perfect_conductor_currents_t = decltype(&T::PerfectConductorCurrents);
+
+    template <typename T>
+    using perfect_conductor_currents_const_t = decltype(&T::PerfectConductorCurrentsConst);
+
     template <typename T>
     using custom_fields_t = decltype(&T::CustomFields);
 
diff --git a/src/global/defaults.h b/src/global/defaults.h
index b81369022..d17647ccd 100644
--- a/src/global/defaults.h
+++ b/src/global/defaults.h
@@ -68,7 +68,8 @@ namespace ntt::defaults {
   } // namespace checkpoint
 
   namespace diag {
-    const timestep_t interval = 1;
+    const timestep_t  interval  = 1;
+    const std::string log_level = "VERBOSE";
   } // namespace diag
 
   namespace gca {
diff --git a/src/global/enums.h b/src/global/enums.h
index 3afa0497a..d80297d8d 100644
--- a/src/global/enums.h
+++ b/src/global/enums.h
@@ -9,7 +9,7 @@
  *   - enum ntt::PrtlBC            // periodic, absorb, atmosphere, custom,
  *                                    reflect, horizon, axis, sync
  *   - enum ntt::FldsBC            // periodic, match, fixed, atmosphere,
- *                                    custom, horizon, axis, sync
+ *                                    custom, horizon, axis, conductor, sync
  *   - enum ntt::PrtlPusher        // boris, vay, photon, none
  *   - enum ntt::Cooling           // synchrotron, none
  *   - enum ntt::FldsID            // e, dive, d, divd, b, h, j,
@@ -221,16 +221,20 @@ namespace ntt {
       CUSTOM     = 5,
       HORIZON    = 6,
       AXIS       = 7,
-      SYNC       = 8, // <- SYNC means synchronization with other domains
+      CONDUCTOR  = 8,
+      SYNC       = 9 // <- SYNC means synchronization with other domains
     };
 
     constexpr FldsBC(uint8_t c) : enums_hidden::BaseEnum<FldsBC> { c } {}
 
-    static constexpr type variants[] = { PERIODIC, MATCH,   FIXED, ATMOSPHERE,
-                                         CUSTOM,   HORIZON, AXIS,  SYNC };
-    static constexpr const char* lookup[] = { "periodic",   "match",  "fixed",
-                                              "atmosphere", "custom", "horizon",
-                                              "axis",       "sync" };
+    static constexpr type variants[] = {
+      PERIODIC, MATCH, FIXED,     ATMOSPHERE, CUSTOM,
+      HORIZON,  AXIS,  CONDUCTOR, SYNC,
+    };
+    static constexpr const char* lookup[] = {
+      "periodic", "match", "fixed",     "atmosphere", "custom",
+      "horizon",  "axis",  "conductor", "sync"
+    };
     static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]);
   };
 
diff --git a/src/global/tests/CMakeLists.txt b/src/global/tests/CMakeLists.txt
index e30da20a0..c206f85b0 100644
--- a/src/global/tests/CMakeLists.txt
+++ b/src/global/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103,C0111
 # ------------------------------
 # @brief: Generates tests for the `ntt_global` module
 #
diff --git a/src/global/tests/enums.cpp b/src/global/tests/enums.cpp
index 673efaf34..ebc074b72 100644
--- a/src/global/tests/enums.cpp
+++ b/src/global/tests/enums.cpp
@@ -61,8 +61,9 @@ 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", "match",   "fixed", "atmosphere",
-                                  "custom",   "horizon", "axis",  "sync" };
+  enum_str_t all_fields_bcs   = { "periodic",   "match",     "fixed",
+                                  "atmosphere", "custom",    "horizon",
+                                  "axis",       "conductor", "sync" };
   enum_str_t all_particle_pushers = { "boris", "vay", "photon", "none" };
   enum_str_t all_coolings         = { "synchrotron", "none" };
 
diff --git a/src/global/utils/formatting.h b/src/global/utils/formatting.h
index 8dc7b6ba8..c8c236cab 100644
--- a/src/global/utils/formatting.h
+++ b/src/global/utils/formatting.h
@@ -171,7 +171,7 @@ namespace fmt {
     std::size_t cntr { 0 };
     for (auto i { 0u }; i < columns.size(); ++i) {
       const auto  anch { static_cast<std::size_t>(anchors[i] < 0 ? -anchors[i]
-                                                                 : anchors[i]) };
+                                                                : anchors[i]) };
       const auto  leftalign { anchors[i] <= 0 };
       const auto  cmn { columns[i] };
       const auto  cmn_len { strlenUTF8(cmn) };
diff --git a/src/global/utils/plog.h b/src/global/utils/plog.h
index 03dc19319..2422c7cf9 100644
--- a/src/global/utils/plog.h
+++ b/src/global/utils/plog.h
@@ -13,6 +13,8 @@
 #ifndef GLOBAL_UTILS_PLOG_H
 #define GLOBAL_UTILS_PLOG_H
 
+#include "utils/formatting.h"
+
 #include <plog/Appenders/ColorConsoleAppender.h>
 #include <plog/Appenders/RollingFileAppender.h>
 #include <plog/Formatters/TxtFormatter.h>
@@ -57,7 +59,7 @@ namespace plog {
 namespace logger {
 
   template <int log_tag, int info_tag, int err_tag>
-  inline void initPlog(const std::string& fname) {
+  inline void initPlog(const std::string& fname, const std::string& log_level) {
     // setup logging
     const auto logfile_name  = fname + ".log";
     const auto infofile_name = fname + ".info";
@@ -77,7 +79,13 @@ namespace logger {
       infofile_name.c_str());
     static plog::RollingFileAppender<plog::NttInfoFormatter> errfileAppender(
       errfile_name.c_str());
-    plog::init<log_tag>(plog::verbose, &logfileAppender);
+    auto log_severity = plog::verbose;
+    if (fmt::toLower(log_level) == "WARNING") {
+      log_severity = plog::warning;
+    } else if (fmt::toLower(log_level) == "ERROR") {
+      log_severity = plog::error;
+    }
+    plog::init<log_tag>(log_severity, &logfileAppender);
     plog::init<info_tag>(plog::verbose, &infofileAppender);
     plog::init<err_tag>(plog::verbose, &errfileAppender);
 
@@ -93,4 +101,4 @@ namespace logger {
 
 } // namespace logger
 
-#endif // GLOBAL_UTILS_PLOG_H
\ No newline at end of file
+#endif // GLOBAL_UTILS_PLOG_H
diff --git a/src/global/utils/timer.cpp b/src/global/utils/timer.cpp
index bc133b4b7..3cca1b365 100644
--- a/src/global/utils/timer.cpp
+++ b/src/global/utils/timer.cpp
@@ -68,7 +68,7 @@ namespace timer {
       return {};
     }
     std::vector<duration_t> all_totals(size, 0.0);
-    for (auto i { 0u }; i < size; ++i) {
+    for (auto i { 0 }; i < size; ++i) {
       for (auto& [name, timer] : m_timers) {
         if (std::find(ignore_in_tot.begin(), ignore_in_tot.end(), name) ==
             ignore_in_tot.end()) {
diff --git a/src/kernels/CMakeLists.txt b/src/kernels/CMakeLists.txt
index c8a1f409f..60eda24cb 100644
--- a/src/kernels/CMakeLists.txt
+++ b/src/kernels/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: ntt_kernels [INTERFACE]
 #
@@ -24,4 +25,3 @@ target_link_libraries(ntt_kernels INTERFACE ${libs})
 
 target_include_directories(ntt_kernels
                            INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/../)
-
diff --git a/src/kernels/digital_filter.hpp b/src/kernels/digital_filter.hpp
index f01eb0fb8..5ac60327d 100644
--- a/src/kernels/digital_filter.hpp
+++ b/src/kernels/digital_filter.hpp
@@ -17,10 +17,41 @@
 #include "utils/error.h"
 #include "utils/numeric.h"
 
-#define FILTER_IN_I1(ARR, COMP, I, J)                                          \
+#define FILTER2D_IN_I1(ARR, COMP, I, J)                                        \
   INV_2*(ARR)((I), (J), (COMP)) +                                              \
     INV_4*((ARR)((I) - 1, (J), (COMP)) + (ARR)((I) + 1, (J), (COMP)))
 
+#define FILTER2D_IN_I2(ARR, COMP, I, J)                                        \
+  INV_2*(ARR)((I), (J), (COMP)) +                                              \
+    INV_4*((ARR)((I), (J) - 1, (COMP)) + (ARR)((I), (J) + 1, (COMP)))
+
+#define FILTER3D_IN_I1_I2(ARR, COMP, I, J, K)                                     \
+  INV_4*(ARR)(I, J, K, (COMP)) +                                                  \
+    INV_8*((ARR)((I) - 1, (J), (K), (COMP)) + (ARR)((I) + 1, (J), (K), (COMP)) +  \
+           (ARR)((I), (J) - 1, (K), (COMP)) + (ARR)((I), (J) + 1, (K), (COMP))) + \
+    INV_16*((ARR)((I) - 1, (J) - 1, (K), (COMP)) +                                \
+            (ARR)((I) + 1, (J) + 1, (K), (COMP)) +                                \
+            (ARR)((I) - 1, (J) + 1, (K), (COMP)) +                                \
+            (ARR)((I) + 1, (J) - 1, (K), (COMP)))
+
+#define FILTER3D_IN_I2_I3(ARR, COMP, I, J, K)                                     \
+  INV_4*(ARR)(I, J, K, (COMP)) +                                                  \
+    INV_8*((ARR)((I), (J) - 1, (K), (COMP)) + (ARR)((I), (J) + 1, (K), (COMP)) +  \
+           (ARR)((I), (J), (K) - 1, (COMP)) + (ARR)((I), (J), (K) + 1, (COMP))) + \
+    INV_16*((ARR)((I), (J) - 1, (K) - 1, (COMP)) +                                \
+            (ARR)((I), (J) + 1, (K) + 1, (COMP)) +                                \
+            (ARR)((I), (J) - 1, (K) + 1, (COMP)) +                                \
+            (ARR)((I), (J) + 1, (K) - 1, (COMP)))
+
+#define FILTER3D_IN_I1_I3(ARR, COMP, I, J, K)                                     \
+  INV_4*(ARR)(I, J, K, (COMP)) +                                                  \
+    INV_8*((ARR)((I) - 1, (J), (K), (COMP)) + (ARR)((I) + 1, (J), (K), (COMP)) +  \
+           (ARR)((I), (J), (K) - 1, (COMP)) + (ARR)((I), (J), (K) + 1, (COMP))) + \
+    INV_16*((ARR)((I) - 1, (J), (K) - 1, (COMP)) +                                \
+            (ARR)((I) + 1, (J), (K) + 1, (COMP)) +                                \
+            (ARR)((I) - 1, (J), (K) + 1, (COMP)) +                                \
+            (ARR)((I) + 1, (J), (K) - 1, (COMP)))
+
 namespace kernel {
   using namespace ntt;
 
@@ -28,9 +59,15 @@ namespace kernel {
   class DigitalFilter_kernel {
     ndfield_t<D, 3>       array;
     const ndfield_t<D, 3> buffer;
-    bool                  is_axis_i2min { false }, is_axis_i2max { false };
-    static constexpr auto i2_min = N_GHOSTS;
-    const ncells_t        i2_max;
+    const bool            is_axis_i2min, is_axis_i2max;
+    const bool            is_conductor_i1min, is_conductor_i1max;
+    const bool            is_conductor_i2min, is_conductor_i2max;
+    const bool            is_conductor_i3min, is_conductor_i3max;
+    static constexpr auto i1_min = N_GHOSTS, i2_min = N_GHOSTS, i3_min = N_GHOSTS;
+    const ncells_t i1_max, i2_max, i3_max;
+
+    // @TODO: Current implementation might have issues
+    // ... at the corners between two conductors
 
   public:
     DigitalFilter_kernel(ndfield_t<D, 3>&       array,
@@ -39,20 +76,50 @@ namespace kernel {
                          const boundaries_t<FldsBC>& boundaries)
       : array { array }
       , buffer { buffer }
-      , i2_max { (short)D > 1 ? size_[1] + N_GHOSTS : 0 } {
-      if constexpr ((C != Coord::Cart) && (D != Dim::_1D)) {
-        raise::ErrorIf(boundaries.size() < 2, "boundaries defined incorrectly", HERE);
-        is_axis_i2min = (boundaries[1].first == FldsBC::AXIS);
-        is_axis_i2max = (boundaries[1].second == FldsBC::AXIS);
-      }
-    }
+      , is_axis_i2min { (D == Dim::_2D) and (boundaries[1].first == FldsBC::AXIS) }
+      , is_axis_i2max { (D == Dim::_2D) and (boundaries[1].second == FldsBC::AXIS) }
+      , is_conductor_i1min { boundaries[0].first == FldsBC::CONDUCTOR }
+      , is_conductor_i1max { boundaries[0].second == FldsBC::CONDUCTOR }
+      , is_conductor_i2min { (short)D > 1
+                               ? (boundaries[1].first == FldsBC::CONDUCTOR)
+                               : false }
+      , is_conductor_i2max { (short)D > 1
+                               ? (boundaries[1].second == FldsBC::CONDUCTOR)
+                               : false }
+      , is_conductor_i3min { (short)D > 2
+                               ? (boundaries[2].first == FldsBC::CONDUCTOR)
+                               : false }
+      , is_conductor_i3max { (short)D > 2
+                               ? (boundaries[2].second == FldsBC::CONDUCTOR)
+                               : false }
+      , i1_max { size_[0] + N_GHOSTS }
+      , i2_max { (short)D > 1 ? (size_[1] + N_GHOSTS) : 0 }
+      , i3_max { (short)D > 2 ? (size_[2] + N_GHOSTS) : 0 } {}
 
     Inline void operator()(index_t i1) const {
       if constexpr ((D == Dim::_1D) && (C == Coord::Cart)) {
+        if ((is_conductor_i1min and i1 == i1_min) or
+            (is_conductor_i1max and i1 == i1_max - 1)) {
+          const auto i1side   = is_conductor_i1min ? (i1 + 1) : (i1 - 1);
+          array(i1, cur::jx1) = (THREE * INV_4) * buffer(i1, cur::jx1) +
+                                (INV_4)*buffer(i1side, cur::jx1);
+        } else if ((is_conductor_i1min and i1 == i1_min + 1) or
+                   (is_conductor_i1max and i1 == i1_max - 2)) {
+          const auto i1side   = is_conductor_i1min ? (i1 + 1) : (i1 - 1);
+          array(i1, cur::jx1) = INV_2 * buffer(i1, cur::jx1) +
+                                INV_4 * (buffer(i1 - 1, cur::jx1) +
+                                         buffer(i1 + 1, cur::jx1));
+          array(i1, cur::jx2) = (INV_2)*buffer(i1, cur::jx2) +
+                                (INV_4)*buffer(i1side, cur::jx2);
+          array(i1, cur::jx3) = (INV_2)*buffer(i1, cur::jx3) +
+                                (INV_4)*buffer(i1side, cur::jx3);
+        } else {
 #pragma unroll
-        for (const auto& comp : { cur::jx1, cur::jx2, cur::jx3 }) {
-          array(i1, comp) = INV_2 * buffer(i1, comp) +
-                            INV_4 * (buffer(i1 - 1, comp) + buffer(i1 + 1, comp));
+          for (const auto& comp : { cur::jx1, cur::jx2, cur::jx3 }) {
+            array(i1, comp) = INV_2 * buffer(i1, comp) +
+                              INV_4 *
+                                (buffer(i1 - 1, comp) + buffer(i1 + 1, comp));
+          }
         }
       } else {
         raise::KernelError(HERE, "DigitalFilter_kernel: 1D implementation called for D != 1 or for non-Cartesian metric");
@@ -62,25 +129,78 @@ namespace kernel {
     Inline void operator()(index_t i1, index_t i2) const {
       if constexpr (D == Dim::_2D) {
         if constexpr (C == Coord::Cart) {
+          if ((is_conductor_i1min and i1 == i1_min) or
+              (is_conductor_i1max and i1 == i1_max - 1)) {
+            const auto i1side = is_conductor_i1min ? (i1 + 1) : (i1 - 1);
+            array(i1, i2, cur::jx1) =
+              (THREE * INV_4) * (FILTER2D_IN_I2(buffer, cur::jx1, i1, i2)) +
+              (INV_4) * (FILTER2D_IN_I2(buffer, cur::jx1, i1side, i2));
+          } else if ((is_conductor_i1min and i1 == i1_min + 1) or
+                     (is_conductor_i1max and i1 == i1_max - 2)) {
+            const auto i1side = is_conductor_i1min ? (i1 + 1) : (i1 - 1);
+            array(i1,
+                  i2,
+                  cur::jx1) = INV_2 * (FILTER2D_IN_I2(buffer, cur::jx1, i1, i2)) +
+                              INV_4 *
+                                ((FILTER2D_IN_I2(buffer, cur::jx1, i1 - 1, i2)) +
+                                 (FILTER2D_IN_I2(buffer, cur::jx1, i1 + 1, i2)));
+            array(i1,
+                  i2,
+                  cur::jx2) = INV_2 * (FILTER2D_IN_I2(buffer, cur::jx2, i1, i2)) +
+                              INV_4 *
+                                (FILTER2D_IN_I2(buffer, cur::jx2, i1side, i2));
+            array(i1,
+                  i2,
+                  cur::jx3) = INV_2 * (FILTER2D_IN_I2(buffer, cur::jx3, i1, i2)) +
+                              INV_4 *
+                                (FILTER2D_IN_I2(buffer, cur::jx3, i1side, i2));
+          } else if ((is_conductor_i2min and i2 == i2_min) or
+                     (is_conductor_i2max and i2 == i2_max - 1)) {
+            const auto i2side = is_conductor_i2min ? (i2 + 1) : (i2 - 1);
+            array(i1, i2, cur::jx2) =
+              (THREE * INV_4) * (FILTER2D_IN_I1(buffer, cur::jx2, i1, i2)) +
+              (INV_4) * (FILTER2D_IN_I1(buffer, cur::jx2, i1, i2side));
+          } else if ((is_conductor_i2min and i2 == i2_min + 1) or
+                     (is_conductor_i2max and i2 == i2_max - 2)) {
+            const auto i2side = is_conductor_i2min ? (i2 + 1) : (i2 - 1);
+            array(i1,
+                  i2,
+                  cur::jx1) = INV_2 * (FILTER2D_IN_I1(buffer, cur::jx1, i1, i2)) +
+                              INV_4 *
+                                (FILTER2D_IN_I1(buffer, cur::jx1, i1, i2side));
+            array(i1,
+                  i2,
+                  cur::jx2) = INV_2 * (FILTER2D_IN_I1(buffer, cur::jx2, i1, i2)) +
+                              INV_4 *
+                                ((FILTER2D_IN_I1(buffer, cur::jx2, i1, i2 - 1)) +
+                                 (FILTER2D_IN_I1(buffer, cur::jx2, i1, i2 + 1)));
+            array(i1,
+                  i2,
+                  cur::jx3) = INV_2 * (FILTER2D_IN_I1(buffer, cur::jx3, i1, i2)) +
+                              INV_4 *
+                                (FILTER2D_IN_I1(buffer, cur::jx3, i1, i2side));
+          } else {
 #pragma unroll
-          for (const auto comp : { cur::jx1, cur::jx2, cur::jx3 }) {
-            array(i1, i2, comp) = INV_4 * buffer(i1, i2, comp) +
-                                  INV_8 * (buffer(i1 - 1, i2, comp) +
-                                           buffer(i1 + 1, i2, comp) +
-                                           buffer(i1, i2 - 1, comp) +
-                                           buffer(i1, i2 + 1, comp)) +
-                                  INV_16 * (buffer(i1 - 1, i2 - 1, comp) +
-                                            buffer(i1 + 1, i2 + 1, comp) +
-                                            buffer(i1 - 1, i2 + 1, comp) +
-                                            buffer(i1 + 1, i2 - 1, comp));
+            for (const auto comp : { cur::jx1, cur::jx2, cur::jx3 }) {
+              array(i1, i2, comp) = INV_4 * buffer(i1, i2, comp) +
+                                    INV_8 * (buffer(i1 - 1, i2, comp) +
+                                             buffer(i1 + 1, i2, comp) +
+                                             buffer(i1, i2 - 1, comp) +
+                                             buffer(i1, i2 + 1, comp)) +
+                                    INV_16 * (buffer(i1 - 1, i2 - 1, comp) +
+                                              buffer(i1 + 1, i2 + 1, comp) +
+                                              buffer(i1 - 1, i2 + 1, comp) +
+                                              buffer(i1 + 1, i2 - 1, comp));
+            }
           }
         } else { // spherical
+          // @TODO: get rid of temporary variables
           real_t cur_00, cur_0p1, cur_0m1;
           if (is_axis_i2min && (i2 == i2_min)) {
             /* --------------------------------- r, phi --------------------------------- */
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx1, i1, i2);
-            cur_0p1 = FILTER_IN_I1(buffer, cur::jx1, i1, i2 + 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2);
+            cur_0p1 = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2 + 1);
             // ... filter in theta
             array(i1, i2, cur::jx1) = INV_2 * cur_00 + INV_2 * cur_0p1;
 
@@ -88,58 +208,58 @@ namespace kernel {
 
             /* ---------------------------------- theta --------------------------------- */
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx2, i1, i2);
-            cur_0p1 = FILTER_IN_I1(buffer, cur::jx2, i1, i2 + 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx2, i1, i2);
+            cur_0p1 = FILTER2D_IN_I1(buffer, cur::jx2, i1, i2 + 1);
             // ... filter in theta
             array(i1, i2, cur::jx2) = INV_4 * (cur_00 + cur_0p1);
           } else if (is_axis_i2min && (i2 == i2_min + 1)) {
             /* --------------------------------- r, phi --------------------------------- */
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx1, i1, i2);
-            cur_0p1 = FILTER_IN_I1(buffer, cur::jx1, i1, i2 + 1);
-            cur_0m1 = FILTER_IN_I1(buffer, cur::jx1, i1, i2 - 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2);
+            cur_0p1 = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2 + 1);
+            cur_0m1 = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2 - 1);
             // ... filter in theta
             array(i1, i2, cur::jx1) = INV_2 * cur_00 + INV_4 * (cur_0p1 + cur_0m1);
 
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx3, i1, i2);
-            cur_0p1 = FILTER_IN_I1(buffer, cur::jx3, i1, i2 + 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx3, i1, i2);
+            cur_0p1 = FILTER2D_IN_I1(buffer, cur::jx3, i1, i2 + 1);
             // ... filter in theta
             array(i1, i2, cur::jx3) = INV_2 * cur_00 + INV_4 * cur_0p1;
 
             /* ---------------------------------- theta --------------------------------- */
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx2, i1, i2);
-            cur_0p1 = FILTER_IN_I1(buffer, cur::jx2, i1, i2 + 1);
-            cur_0m1 = FILTER_IN_I1(buffer, cur::jx2, i1, i2 - 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx2, i1, i2);
+            cur_0p1 = FILTER2D_IN_I1(buffer, cur::jx2, i1, i2 + 1);
+            cur_0m1 = FILTER2D_IN_I1(buffer, cur::jx2, i1, i2 - 1);
             // ... filter in theta
             array(i1, i2, cur::jx2) = INV_2 * cur_00 + INV_4 * (cur_0m1 + cur_0p1);
           } else if (is_axis_i2max && (i2 == i2_max - 1)) {
             /* --------------------------------- r, phi --------------------------------- */
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx1, i1, i2);
-            cur_0p1 = FILTER_IN_I1(buffer, cur::jx1, i1, i2 + 1);
-            cur_0m1 = FILTER_IN_I1(buffer, cur::jx1, i1, i2 - 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2);
+            cur_0p1 = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2 + 1);
+            cur_0m1 = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2 - 1);
             // ... filter in theta
             array(i1, i2, cur::jx1) = INV_2 * cur_00 + INV_4 * (cur_0m1 + cur_0p1);
 
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx3, i1, i2);
-            cur_0m1 = FILTER_IN_I1(buffer, cur::jx3, i1, i2 - 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx3, i1, i2);
+            cur_0m1 = FILTER2D_IN_I1(buffer, cur::jx3, i1, i2 - 1);
             // ... filter in theta
             array(i1, i2, cur::jx3) = INV_2 * cur_00 + INV_4 * cur_0m1;
 
             /* ---------------------------------- theta --------------------------------- */
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx2, i1, i2);
-            cur_0m1 = FILTER_IN_I1(buffer, cur::jx2, i1, i2 - 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx2, i1, i2);
+            cur_0m1 = FILTER2D_IN_I1(buffer, cur::jx2, i1, i2 - 1);
             // ... filter in theta
             array(i1, i2, cur::jx2) = INV_4 * (cur_00 + cur_0m1);
           } else if (is_axis_i2max && (i2 == i2_max)) {
             /* --------------------------------- r, phi --------------------------------- */
             // ... filter in r
-            cur_00  = FILTER_IN_I1(buffer, cur::jx1, i1, i2);
-            cur_0m1 = FILTER_IN_I1(buffer, cur::jx1, i1, i2 - 1);
+            cur_00  = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2);
+            cur_0m1 = FILTER2D_IN_I1(buffer, cur::jx1, i1, i2 - 1);
             // ... filter in theta
             array(i1, i2, cur::jx1) = INV_2 * cur_00 + INV_2 * cur_0m1;
 
@@ -169,33 +289,93 @@ namespace kernel {
     Inline void operator()(index_t i1, index_t i2, index_t i3) const {
       if constexpr (D == Dim::_3D) {
         if constexpr (C == Coord::Cart) {
+          if ((is_conductor_i1min and i1 == i1_min) or
+              (is_conductor_i1max and i1 == i1_max - 1)) {
+            const auto i1side = is_conductor_i1min ? (i1 + 1) : (i1 - 1);
+            array(i1, i2, i3, cur::jx1) =
+              (THREE * INV_4) * (FILTER3D_IN_I2_I3(buffer, cur::jx1, i1, i2, i3)) +
+              (INV_4) * (FILTER3D_IN_I2_I3(buffer, cur::jx1, i1side, i2, i3));
+          } else if ((is_conductor_i1min and i1 == i1_min + 1) or
+                     (is_conductor_i1max and i1 == i1_max - 2)) {
+            const auto i1side = is_conductor_i1min ? (i1 + 1) : (i1 - 1);
+            array(i1, i2, i3, cur::jx1) =
+              INV_2 * (FILTER3D_IN_I2_I3(buffer, cur::jx1, i1, i2, i3)) +
+              INV_4 * ((FILTER3D_IN_I2_I3(buffer, cur::jx1, i1 - 1, i2, i3)) +
+                       (FILTER3D_IN_I2_I3(buffer, cur::jx1, i1 + 1, i2, i3)));
+            array(i1, i2, i3, cur::jx2) =
+              INV_2 * (FILTER3D_IN_I2_I3(buffer, cur::jx2, i1, i2, i3)) +
+              INV_4 * (FILTER3D_IN_I2_I3(buffer, cur::jx2, i1side, i2, i3));
+            array(i1, i2, i3, cur::jx3) =
+              INV_2 * (FILTER3D_IN_I2_I3(buffer, cur::jx3, i1, i2, i3)) +
+              INV_4 * (FILTER3D_IN_I2_I3(buffer, cur::jx3, i1side, i2, i3));
+          } else if ((is_conductor_i2min and i2 == i2_min) or
+                     (is_conductor_i2max and i2 == i2_max - 1)) {
+            const auto i2side = is_conductor_i2min ? (i2 + 1) : (i2 - 1);
+            array(i1, i2, i3, cur::jx2) =
+              (THREE * INV_4) * (FILTER3D_IN_I1_I3(buffer, cur::jx2, i1, i2, i3)) +
+              (INV_4) * (FILTER3D_IN_I1_I3(buffer, cur::jx2, i1, i2side, i3));
+          } else if ((is_conductor_i2min and i2 == i2_min + 1) or
+                     (is_conductor_i2max and i2 == i2_max - 2)) {
+            const auto i2side = is_conductor_i2min ? (i2 + 1) : (i2 - 1);
+            array(i1, i2, i3, cur::jx1) =
+              INV_2 * (FILTER3D_IN_I1_I3(buffer, cur::jx1, i1, i2, i3)) +
+              INV_4 * (FILTER3D_IN_I1_I3(buffer, cur::jx1, i1, i2side, i3));
+            array(i1, i2, i3, cur::jx2) =
+              INV_2 * (FILTER3D_IN_I1_I3(buffer, cur::jx2, i1, i2, i3)) +
+              INV_4 * ((FILTER3D_IN_I1_I3(buffer, cur::jx2, i1, i2 - 1, i3)) +
+                       (FILTER3D_IN_I1_I3(buffer, cur::jx2, i1, i2 + 1, i3)));
+            array(i1, i2, i3, cur::jx3) =
+              INV_2 * (FILTER3D_IN_I1_I3(buffer, cur::jx3, i1, i2, i3)) +
+              INV_4 * (FILTER3D_IN_I1_I3(buffer, cur::jx3, i1, i2side, i3));
+          } else if ((is_conductor_i3min and i3 == i3_min) or
+                     (is_conductor_i3max and i3 == i3_max - 1)) {
+            const auto i3side = is_conductor_i3min ? (i3 + 1) : (i3 - 1);
+            array(i1, i2, i3, cur::jx3) =
+              (THREE * INV_4) * (FILTER3D_IN_I1_I2(buffer, cur::jx3, i1, i2, i3)) +
+              (INV_4) * (FILTER3D_IN_I1_I2(buffer, cur::jx3, i1, i2, i3side));
+          } else if ((is_conductor_i3min and i3 == i3_min + 1) or
+                     (is_conductor_i3max and i3 == i3_max - 2)) {
+            const auto i3side = is_conductor_i3min ? (i3 + 1) : (i3 - 1);
+            array(i1, i2, i3, cur::jx1) =
+              INV_2 * (FILTER3D_IN_I1_I2(buffer, cur::jx1, i1, i2, i3)) +
+              INV_4 * (FILTER3D_IN_I1_I2(buffer, cur::jx1, i1, i2, i3side));
+            array(i1, i2, i3, cur::jx2) =
+              INV_2 * (FILTER3D_IN_I1_I2(buffer, cur::jx2, i1, i2, i3)) +
+              INV_4 * (FILTER3D_IN_I1_I2(buffer, cur::jx2, i1, i2, i3side));
+            array(i1, i2, i3, cur::jx3) =
+              INV_2 * (FILTER3D_IN_I1_I2(buffer, cur::jx3, i1, i2, i3)) +
+              INV_4 * ((FILTER3D_IN_I1_I2(buffer, cur::jx3, i1, i2, i3 - 1)) +
+                       (FILTER3D_IN_I1_I2(buffer, cur::jx3, i1, i2, i3 + 1)));
+          } else {
 #pragma unroll
-          for (auto& comp : { cur::jx1, cur::jx2, cur::jx3 }) {
-            array(i1, i2, i3, comp) =
-              INV_8 * buffer(i1, i2, i3, comp) +
-              INV_16 *
-                (buffer(i1 - 1, i2, i3, comp) + buffer(i1 + 1, i2, i3, comp) +
-                 buffer(i1, i2 - 1, i3, comp) + buffer(i1, i2 + 1, i3, comp) +
-                 buffer(i1, i2, i3 - 1, comp) + buffer(i1, i2, i3 + 1, comp)) +
-              INV_32 *
-                (buffer(i1 - 1, i2 - 1, i3, comp) +
-                 buffer(i1 + 1, i2 + 1, i3, comp) +
-                 buffer(i1 - 1, i2 + 1, i3, comp) +
-                 buffer(i1 + 1, i2 - 1, i3, comp) +
-                 buffer(i1, i2 - 1, i3 - 1, comp) +
-                 buffer(i1, i2 + 1, i3 + 1, comp) + buffer(i1, i2, i3 - 1, comp) +
-                 buffer(i1, i2, i3 + 1, comp) + buffer(i1 - 1, i2, i3 - 1, comp) +
-                 buffer(i1 + 1, i2, i3 + 1, comp) +
-                 buffer(i1 - 1, i2, i3 + 1, comp) +
-                 buffer(i1 + 1, i2, i3 - 1, comp)) +
-              INV_64 * (buffer(i1 - 1, i2 - 1, i3 - 1, comp) +
-                        buffer(i1 + 1, i2 + 1, i3 + 1, comp) +
-                        buffer(i1 - 1, i2 + 1, i3 + 1, comp) +
-                        buffer(i1 + 1, i2 - 1, i3 - 1, comp) +
-                        buffer(i1 - 1, i2 - 1, i3 + 1, comp) +
-                        buffer(i1 + 1, i2 + 1, i3 - 1, comp) +
-                        buffer(i1 - 1, i2 + 1, i3 - 1, comp) +
-                        buffer(i1 + 1, i2 - 1, i3 + 1, comp));
+            for (auto& comp : { cur::jx1, cur::jx2, cur::jx3 }) {
+              array(i1, i2, i3, comp) =
+                INV_8 * buffer(i1, i2, i3, comp) +
+                INV_16 *
+                  (buffer(i1 - 1, i2, i3, comp) + buffer(i1 + 1, i2, i3, comp) +
+                   buffer(i1, i2 - 1, i3, comp) + buffer(i1, i2 + 1, i3, comp) +
+                   buffer(i1, i2, i3 - 1, comp) + buffer(i1, i2, i3 + 1, comp)) +
+                INV_32 *
+                  (buffer(i1 - 1, i2 - 1, i3, comp) +
+                   buffer(i1 + 1, i2 + 1, i3, comp) +
+                   buffer(i1 - 1, i2 + 1, i3, comp) +
+                   buffer(i1 + 1, i2 - 1, i3, comp) +
+                   buffer(i1, i2 - 1, i3 - 1, comp) +
+                   buffer(i1, i2 + 1, i3 + 1, comp) +
+                   buffer(i1, i2, i3 - 1, comp) + buffer(i1, i2, i3 + 1, comp) +
+                   buffer(i1 - 1, i2, i3 - 1, comp) +
+                   buffer(i1 + 1, i2, i3 + 1, comp) +
+                   buffer(i1 - 1, i2, i3 + 1, comp) +
+                   buffer(i1 + 1, i2, i3 - 1, comp)) +
+                INV_64 * (buffer(i1 - 1, i2 - 1, i3 - 1, comp) +
+                          buffer(i1 + 1, i2 + 1, i3 + 1, comp) +
+                          buffer(i1 - 1, i2 + 1, i3 + 1, comp) +
+                          buffer(i1 + 1, i2 - 1, i3 - 1, comp) +
+                          buffer(i1 - 1, i2 - 1, i3 + 1, comp) +
+                          buffer(i1 + 1, i2 + 1, i3 - 1, comp) +
+                          buffer(i1 - 1, i2 + 1, i3 - 1, comp) +
+                          buffer(i1 + 1, i2 - 1, i3 + 1, comp));
+            }
           }
         } else {
           raise::KernelNotImplementedError(HERE);
@@ -210,6 +390,11 @@ namespace kernel {
 
 } // namespace kernel
 
-#undef FILTER_IN_I1
+#undef FILTER3D_IN_I1_I3
+#undef FILTER3D_IN_I2_I3
+#undef FILTER3D_IN_I1_I2
+
+#undef FILTER2D_IN_I2
+#undef FILTER2D_IN_I1
 
 #endif // DIGITAL_FILTER_HPP
diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp
index dbb47f42c..782200e29 100644
--- a/src/kernels/fields_bcs.hpp
+++ b/src/kernels/fields_bcs.hpp
@@ -5,6 +5,7 @@
  *   - kernel::bc::MatchBoundaries_kernel<>
  *   - kernel::bc::AxisBoundaries_kernel<>
  *   - kernel::bc::EnforcedBoundaries_kernel<>
+ *   - kernel::bc::ConductorBoundaries_kernel<>
  * @namespaces:
  *   - kernel::bc::
  */
@@ -485,6 +486,301 @@ namespace kernel::bc {
     }
   };
 
+  template <Dimension D, in o, bool P>
+  struct ConductorBoundaries_kernel {
+    static_assert(static_cast<unsigned short>(o) < static_cast<unsigned short>(D),
+                  "Invalid component index");
+
+    ndfield_t<D, 6>   Fld;
+    const BCTags      tags;
+    const std::size_t i_edge;
+
+    ConductorBoundaries_kernel(ndfield_t<D, 6> Fld, std::size_t i_edge, BCTags tags)
+      : Fld { Fld }
+      , i_edge { i_edge }
+      , tags { tags } {}
+
+    Inline void operator()(index_t i1) const {
+      if constexpr (D == Dim::_1D) {
+        if (tags & BC::E) {
+          if (i1 == 0) {
+            Fld(i_edge, em::ex2) = ZERO;
+            Fld(i_edge, em::ex3) = ZERO;
+          } else {
+            if constexpr (not P) {
+              Fld(i_edge - i1, em::ex1) = Fld(i_edge + i1 - 1, em::ex1);
+              Fld(i_edge - i1, em::ex2) = -Fld(i_edge + i1, em::ex2);
+              Fld(i_edge - i1, em::ex3) = -Fld(i_edge + i1, em::ex3);
+            } else {
+              Fld(i_edge + i1 - 1, em::ex1) = Fld(i_edge - i1, em::ex1);
+              Fld(i_edge + i1, em::ex2)     = -Fld(i_edge - i1, em::ex2);
+              Fld(i_edge + i1, em::ex3)     = -Fld(i_edge - i1, em::ex3);
+            }
+          }
+        }
+
+        if (tags & BC::B) {
+          if (i1 == 0) {
+            Fld(i_edge, em::bx1) = ZERO;
+          } else {
+            if constexpr (not P) {
+              Fld(i_edge - i1, em::bx1) = -Fld(i_edge + i1, em::bx1);
+              Fld(i_edge - i1, em::bx2) = Fld(i_edge + i1 - 1, em::bx2);
+              Fld(i_edge - i1, em::bx3) = Fld(i_edge + i1 - 1, em::bx3);
+            } else {
+              Fld(i_edge + i1, em::bx1)     = -Fld(i_edge - i1, em::bx1);
+              Fld(i_edge + i1 - 1, em::bx2) = Fld(i_edge - i1, em::bx2);
+              Fld(i_edge + i1 - 1, em::bx3) = Fld(i_edge - i1, em::bx3);
+            }
+          }
+        }
+      } else {
+        raise::KernelError(
+          HERE,
+          "ConductorBoundaries_kernel: 1D implementation called for D != 1");
+      }
+    }
+
+    Inline void operator()(index_t i1, index_t i2) const {
+      if constexpr (D == Dim::_2D) {
+        if constexpr (o == in::x1) {
+          if (tags & BC::E) {
+            if (i1 == 0) {
+              Fld(i_edge, i2, em::ex2) = ZERO;
+              Fld(i_edge, i2, em::ex3) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i_edge - i1, i2, em::ex1) = Fld(i_edge + i1 - 1, i2, em::ex1);
+                Fld(i_edge - i1, i2, em::ex2) = -Fld(i_edge + i1, i2, em::ex2);
+                Fld(i_edge - i1, i2, em::ex3) = -Fld(i_edge + i1, i2, em::ex3);
+              } else {
+                Fld(i_edge + i1 - 1, i2, em::ex1) = Fld(i_edge - i1, i2, em::ex1);
+                Fld(i_edge + i1, i2, em::ex2) = -Fld(i_edge - i1, i2, em::ex2);
+                Fld(i_edge + i1, i2, em::ex3) = -Fld(i_edge - i1, i2, em::ex3);
+              }
+            }
+          }
+
+          if (tags & BC::B) {
+            if (i1 == 0) {
+              Fld(i_edge, i2, em::bx1) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i_edge - i1, i2, em::bx1) = -Fld(i_edge + i1, i2, em::bx1);
+                Fld(i_edge - i1, i2, em::bx2) = Fld(i_edge + i1 - 1, i2, em::bx2);
+                Fld(i_edge - i1, i2, em::bx3) = Fld(i_edge + i1 - 1, i2, em::bx3);
+              } else {
+                Fld(i_edge + i1, i2, em::bx1) = -Fld(i_edge - i1, i2, em::bx1);
+                Fld(i_edge + i1 - 1, i2, em::bx2) = Fld(i_edge - i1, i2, em::bx2);
+                Fld(i_edge + i1 - 1, i2, em::bx3) = Fld(i_edge - i1, i2, em::bx3);
+              }
+            }
+          }
+        } else {
+          if (tags & BC::E) {
+            if (i2 == 0) {
+              Fld(i1, i_edge, em::ex1) = ZERO;
+              Fld(i1, i_edge, em::ex3) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i1, i_edge - i2, em::ex1) = -Fld(i1, i_edge + i2, em::ex1);
+                Fld(i1, i_edge - i2, em::ex2) = Fld(i1, i_edge + i2 - 1, em::ex2);
+                Fld(i1, i_edge - i2, em::ex3) = -Fld(i1, i_edge + i2, em::ex3);
+              } else {
+                Fld(i1, i_edge + i2, em::ex1) = -Fld(i1, i_edge - i2, em::ex1);
+                Fld(i1, i_edge + i2 - 1, em::ex2) = Fld(i1, i_edge - i2, em::ex2);
+                Fld(i1, i_edge + i2, em::ex3) = -Fld(i1, i_edge - i2, em::ex3);
+              }
+            }
+          }
+
+          if (tags & BC::B) {
+            if (i2 == 0) {
+              Fld(i1, i_edge, em::bx2) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i1, i_edge - i2, em::bx1) = Fld(i1, i_edge + i2 - 1, em::bx1);
+                Fld(i1, i_edge - i2, em::bx2) = -Fld(i1, i_edge + i2, em::bx2);
+                Fld(i1, i_edge - i2, em::bx3) = Fld(i1, i_edge + i2 - 1, em::bx3);
+              } else {
+                Fld(i1, i_edge + i2 - 1, em::bx1) = Fld(i1, i_edge - i2, em::bx1);
+                Fld(i1, i_edge + i2, em::bx2) = -Fld(i1, i_edge - i2, em::bx2);
+                Fld(i1, i_edge + i2 - 1, em::bx3) = Fld(i1, i_edge - i2, em::bx3);
+              }
+            }
+          }
+        }
+      } else {
+        raise::KernelError(
+          HERE,
+          "ConductorBoundaries_kernel: 2D implementation called for D != 2");
+      }
+    }
+
+    Inline void operator()(index_t i1, index_t i2, index_t i3) const {
+      if constexpr (D == Dim::_3D) {
+        if constexpr (o == in::x1) {
+          if (tags & BC::E) {
+            if (i1 == 0) {
+              Fld(i_edge, i2, i3, em::ex2) = ZERO;
+              Fld(i_edge, i2, i3, em::ex3) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i_edge - i1, i2, i3, em::ex1) = Fld(i_edge + i1 - 1,
+                                                        i2,
+                                                        i3,
+                                                        em::ex1);
+                Fld(i_edge - i1, i2, i3, em::ex2) = -Fld(i_edge + i1, i2, i3, em::ex2);
+                Fld(i_edge - i1, i2, i3, em::ex3) = -Fld(i_edge + i1, i2, i3, em::ex3);
+              } else {
+                Fld(i_edge + i1 - 1, i2, i3, em::ex1) = Fld(i_edge - i1,
+                                                            i2,
+                                                            i3,
+                                                            em::ex1);
+                Fld(i_edge + i1, i2, i3, em::ex2) = -Fld(i_edge - i1, i2, i3, em::ex2);
+                Fld(i_edge + i1, i2, i3, em::ex3) = -Fld(i_edge - i1, i2, i3, em::ex3);
+              }
+            }
+          }
+
+          if (tags & BC::B) {
+            if (i1 == 0) {
+              Fld(i_edge, i2, i3, em::bx1) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i_edge - i1, i2, i3, em::bx1) = -Fld(i_edge + i1, i2, i3, em::bx1);
+                Fld(i_edge - i1, i2, i3, em::bx2) = Fld(i_edge + i1 - 1,
+                                                        i2,
+                                                        i3,
+                                                        em::bx2);
+                Fld(i_edge - i1, i2, i3, em::bx3) = Fld(i_edge + i1 - 1,
+                                                        i2,
+                                                        i3,
+                                                        em::bx3);
+              } else {
+                Fld(i_edge + i1, i2, i3, em::bx1) = -Fld(i_edge - i1, i2, i3, em::bx1);
+                Fld(i_edge + i1 - 1, i2, i3, em::bx2) = Fld(i_edge - i1,
+                                                            i2,
+                                                            i3,
+                                                            em::bx2);
+                Fld(i_edge + i1 - 1, i2, i3, em::bx3) = Fld(i_edge - i1,
+                                                            i2,
+                                                            i3,
+                                                            em::bx3);
+              }
+            }
+          }
+        } else if (o == in::x2) {
+          if (tags & BC::E) {
+            if (i2 == 0) {
+              Fld(i1, i_edge, i3, em::ex1) = ZERO;
+              Fld(i1, i_edge, i3, em::ex3) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i1, i_edge - i2, i3, em::ex1) = -Fld(i1, i_edge + i2, i3, em::ex1);
+                Fld(i1, i_edge - i2, i3, em::ex2) = Fld(i1,
+                                                        i_edge + i2 - 1,
+                                                        i3,
+                                                        em::ex2);
+                Fld(i1, i_edge - i2, i3, em::ex3) = -Fld(i1, i_edge + i2, i3, em::ex3);
+              } else {
+                Fld(i1, i_edge + i2, i3, em::ex1) = -Fld(i1, i_edge - i2, i3, em::ex1);
+                Fld(i1, i_edge + i2 - 1, i3, em::ex2) = Fld(i1,
+                                                            i_edge - i2,
+                                                            i3,
+                                                            em::ex2);
+                Fld(i1, i_edge + i2, i3, em::ex3) = -Fld(i1, i_edge - i2, i3, em::ex3);
+              }
+            }
+          }
+
+          if (tags & BC::B) {
+            if (i2 == 0) {
+              Fld(i1, i_edge, i3, em::bx2) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i1, i_edge - i2, i3, em::bx1) = Fld(i1,
+                                                        i_edge + i2 - 1,
+                                                        i3,
+                                                        em::bx1);
+                Fld(i1, i_edge - i2, i3, em::bx2) = -Fld(i1, i_edge + i2, i3, em::bx2);
+                Fld(i1, i_edge - i2, i3, em::bx3) = Fld(i1,
+                                                        i_edge + i2 - 1,
+                                                        i3,
+                                                        em::bx3);
+              } else {
+                Fld(i1, i_edge + i2 - 1, i3, em::bx1) = Fld(i1,
+                                                            i_edge - i2,
+                                                            i3,
+                                                            em::bx1);
+                Fld(i1, i_edge + i2, i3, em::bx2) = -Fld(i1, i_edge - i2, i3, em::bx2);
+                Fld(i1, i_edge + i2 - 1, i3, em::bx3) = Fld(i1,
+                                                            i_edge - i2,
+                                                            i3,
+                                                            em::bx3);
+              }
+            }
+          }
+        } else {
+          if (tags & BC::E) {
+            if (i3 == 0) {
+              Fld(i1, i2, i_edge, em::ex1) = ZERO;
+              Fld(i1, i2, i_edge, em::ex2) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i1, i2, i_edge - i3, em::ex1) = -Fld(i1, i2, i_edge + i3, em::ex1);
+                Fld(i1, i2, i_edge - i3, em::ex2) = -Fld(i1, i2, i_edge + i3, em::ex2);
+                Fld(i1, i2, i_edge - i3, em::ex3) = Fld(i1,
+                                                        i2,
+                                                        i_edge + i3 - 1,
+                                                        em::ex3);
+              } else {
+                Fld(i1, i2, i_edge + i3, em::ex1) = -Fld(i1, i2, i_edge - i3, em::ex1);
+                Fld(i1, i2, i_edge + i3, em::ex2) = -Fld(i1, i2, i_edge - i3, em::ex2);
+                Fld(i1, i2, i_edge + i3 - 1, em::ex3) = Fld(i1,
+                                                            i2,
+                                                            i_edge - i3,
+                                                            em::ex3);
+              }
+            }
+          }
+
+          if (tags & BC::B) {
+            if (i3 == 0) {
+              Fld(i1, i2, i_edge, em::bx3) = ZERO;
+            } else {
+              if constexpr (not P) {
+                Fld(i1, i2, i_edge - i3, em::bx1) = Fld(i1,
+                                                        i2,
+                                                        i_edge + i3 - 1,
+                                                        em::bx1);
+                Fld(i1, i2, i_edge - i3, em::bx2) = Fld(i1,
+                                                        i2,
+                                                        i_edge + i3 - 1,
+                                                        em::bx2);
+                Fld(i1, i2, i_edge - i3, em::bx3) = -Fld(i1, i2, i_edge + i3, em::bx3);
+              } else {
+                Fld(i1, i2, i_edge + i3 - 1, em::bx1) = Fld(i1,
+                                                            i2,
+                                                            i_edge - i3,
+                                                            em::bx1);
+                Fld(i1, i2, i_edge + i3 - 1, em::bx2) = Fld(i1,
+                                                            i2,
+                                                            i_edge - i3,
+                                                            em::bx2);
+                Fld(i1, i2, i_edge + i3, em::bx3) = -Fld(i1, i2, i_edge - i3, em::bx3);
+              }
+            }
+          }
+        }
+      } else {
+        raise::KernelError(
+          HERE,
+          "ConductorBoundaries_kernel: 3D implementation called for D != 3");
+      }
+    }
+  };
+
   /*
    * @tparam D: Dimension
    * @tparam P: Positive/Negative direction
@@ -562,6 +858,7 @@ namespace kernel::bc {
     const I         fset;
     const M         metric;
     const ncells_t  i_edge;
+    const BCTags    tags;
 
     EnforcedBoundaries_kernel(ndfield_t<M::Dim, 6>& Fld,
                               const I&              fset,
@@ -571,7 +868,8 @@ namespace kernel::bc {
       : Fld { Fld }
       , fset { fset }
       , metric { metric }
-      , i_edge { i_edge + N_GHOSTS } {}
+      , i_edge { i_edge + N_GHOSTS }
+      , tags { tags } {}
 
     Inline void operator()(index_t i1) const {
       if constexpr (D == Dim::_1D) {
@@ -580,8 +878,12 @@ namespace kernel::bc {
         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 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3,
-             setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3;
+        bool setEx1 = defines_ex1 and (tags & BC::E),
+             setEx2 = defines_ex2 and (tags & BC::E),
+             setEx3 = defines_ex3 and (tags & BC::E),
+             setBx1 = defines_bx1 and (tags & BC::B),
+             setBx2 = defines_bx2 and (tags & BC::B),
+             setBx3 = defines_bx3 and (tags & BC::B);
         if constexpr (O == in::x1) {
           // x1 -- normal
           // x2,x3 -- tangential
@@ -657,8 +959,12 @@ namespace kernel::bc {
         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 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3,
-             setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3;
+        bool setEx1 = defines_ex1 and (tags & BC::E),
+             setEx2 = defines_ex2 and (tags & BC::E),
+             setEx3 = defines_ex3 and (tags & BC::E),
+             setBx1 = defines_bx1 and (tags & BC::B),
+             setBx2 = defines_bx2 and (tags & BC::B),
+             setBx3 = defines_bx3 and (tags & BC::B);
         if constexpr (O == in::x1) {
           // x1 -- normal
           // x2,x3 -- tangential
@@ -757,8 +1063,12 @@ namespace kernel::bc {
                                                   x_Ph_H0H);
         metric.template convert<Crd::Cd, Crd::Ph>({ i1_, i2_ + HALF, i3_ + HALF },
                                                   x_Ph_0HH);
-        bool setEx1 = defines_ex1, setEx2 = defines_ex2, setEx3 = defines_ex3,
-             setBx1 = defines_bx1, setBx2 = defines_bx2, setBx3 = defines_bx3;
+        bool setEx1 = defines_ex1 and (tags & BC::E),
+             setEx2 = defines_ex2 and (tags & BC::E),
+             setEx3 = defines_ex3 and (tags & BC::E),
+             setBx1 = defines_bx1 and (tags & BC::B),
+             setBx2 = defines_bx2 and (tags & BC::B),
+             setBx3 = defines_bx3 and (tags & BC::B);
         if constexpr (O == in::x1) {
           // x1 -- normal
           // x2,x3 -- tangential
diff --git a/src/kernels/tests/CMakeLists.txt b/src/kernels/tests/CMakeLists.txt
index a41ea43ef..7579eb6d3 100644
--- a/src/kernels/tests/CMakeLists.txt
+++ b/src/kernels/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103,C0111
 # ------------------------------
 # @brief: Generates tests for the `ntt_kernels` module
 #
diff --git a/src/kernels/tests/digital_filter.cpp b/src/kernels/tests/digital_filter.cpp
index e0cc352f5..03a63753b 100644
--- a/src/kernels/tests/digital_filter.cpp
+++ b/src/kernels/tests/digital_filter.cpp
@@ -40,8 +40,13 @@ void testFilter(const std::vector<std::size_t>&      res,
   auto boundaries = boundaries_t<FldsBC> {};
   if constexpr (M::CoordType != Coord::Cart) {
     boundaries = {
-      {FldsBC::CUSTOM, FldsBC::CUSTOM},
-      {  FldsBC::AXIS,   FldsBC::AXIS}
+      { FldsBC::CUSTOM, FldsBC::CUSTOM },
+      {   FldsBC::AXIS,   FldsBC::AXIS }
+    };
+  } else {
+    boundaries = {
+      { FldsBC::PERIODIC, FldsBC::PERIODIC },
+      { FldsBC::PERIODIC, FldsBC::PERIODIC }
     };
   }
 
@@ -183,4 +188,4 @@ auto main(int argc, char* argv[]) -> int {
   }
   Kokkos::finalize();
   return 0;
-}
\ No newline at end of file
+}
diff --git a/src/kernels/tests/prtls_to_phys.cpp b/src/kernels/tests/prtls_to_phys.cpp
index 4719fe6a1..962c21b5c 100644
--- a/src/kernels/tests/prtls_to_phys.cpp
+++ b/src/kernels/tests/prtls_to_phys.cpp
@@ -132,7 +132,7 @@ void testPrtl2PhysSR(const std::vector<std::size_t>&      res,
 
   extent = {
     ext[0],
-    {ZERO, constant::PI}
+    { ZERO, constant::PI }
   };
 
   const M metric { res, extent, params };
@@ -177,28 +177,29 @@ void testPrtl2PhysSR(const std::vector<std::size_t>&      res,
   array_t<real_t*>  buff_ux3 { "buff_ux3", nprtl / stride };
   array_t<real_t*>  buff_wei { "buff_wei", nprtl / stride };
 
-  Kokkos::parallel_for("Init",
-                       Kokkos::RangePolicy<AccelExeSpace>(0, nprtl / stride),
-                       kernel::PrtlToPhys_kernel<SimEngine::SRPIC, M>(stride,
-                                                                      buff_x1,
-                                                                      buff_x2,
-                                                                      buff_x3,
-                                                                      buff_ux1,
-                                                                      buff_ux2,
-                                                                      buff_ux3,
-                                                                      buff_wei,
-                                                                      i1,
-                                                                      i2,
-                                                                      i3,
-                                                                      dx1,
-                                                                      dx2,
-                                                                      dx3,
-                                                                      ux1,
-                                                                      ux2,
-                                                                      ux3,
-                                                                      phi,
-                                                                      weight,
-                                                                      metric));
+  Kokkos::parallel_for(
+    "Init",
+    Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, nprtl / stride),
+    kernel::PrtlToPhys_kernel<SimEngine::SRPIC, M>(stride,
+                                                   buff_x1,
+                                                   buff_x2,
+                                                   buff_x3,
+                                                   buff_ux1,
+                                                   buff_ux2,
+                                                   buff_ux3,
+                                                   buff_wei,
+                                                   i1,
+                                                   i2,
+                                                   i3,
+                                                   dx1,
+                                                   dx2,
+                                                   dx3,
+                                                   ux1,
+                                                   ux2,
+                                                   ux3,
+                                                   phi,
+                                                   weight,
+                                                   metric));
   Kokkos::parallel_for("Check",
                        nprtl / stride,
                        Checker<M>(metric,
diff --git a/src/metrics/CMakeLists.txt b/src/metrics/CMakeLists.txt
index e053bb61c..0bb5b977c 100644
--- a/src/metrics/CMakeLists.txt
+++ b/src/metrics/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: ntt_metrics [INTERFACE]
 #
@@ -22,4 +23,3 @@ target_link_libraries(ntt_metrics INTERFACE ${libs})
 
 target_include_directories(ntt_metrics
                            INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/../)
-
diff --git a/src/metrics/tests/CMakeLists.txt b/src/metrics/tests/CMakeLists.txt
index c997ab079..0d661c318 100644
--- a/src/metrics/tests/CMakeLists.txt
+++ b/src/metrics/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103,C0111
 # ------------------------------
 # @brief: Generates tests for the `ntt_metrics` module
 #
@@ -28,4 +29,3 @@ gen_test(coord_trans)
 gen_test(sph-qsph)
 gen_test(ks-qks)
 gen_test(sr-cart-sph)
-
diff --git a/src/metrics/tests/ks-qks.cpp b/src/metrics/tests/ks-qks.cpp
index bed051e16..167f564ee 100644
--- a/src/metrics/tests/ks-qks.cpp
+++ b/src/metrics/tests/ks-qks.cpp
@@ -27,6 +27,7 @@ Inline auto equal(const vec_t<D>& a,
   const auto eps = epsilon * acc;
   for (unsigned short d = 0; d < D; ++d) {
     if (not cmp::AlmostEqual(a[d], b[d], eps)) {
+      printf("%s: %.12e : %.12e\n", msg, a[d], b[d]);
       return false;
     }
   }
diff --git a/src/output/CMakeLists.txt b/src/output/CMakeLists.txt
index 81333e9ff..8a2ea0f16 100644
--- a/src/output/CMakeLists.txt
+++ b/src/output/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103
 # ------------------------------
 # @defines: ntt_output [STATIC/SHARED]
 #
diff --git a/src/output/tests/CMakeLists.txt b/src/output/tests/CMakeLists.txt
index afc7950c4..835bb532f 100644
--- a/src/output/tests/CMakeLists.txt
+++ b/src/output/tests/CMakeLists.txt
@@ -1,3 +1,4 @@
+# cmake-lint: disable=C0103,C0111
 # ------------------------------
 # @brief: Generates tests for the `ntt_output` module
 #
diff --git a/src/output/tests/writer-mpi.cpp b/src/output/tests/writer-mpi.cpp
index f6d3ee88a..af8a38ef1 100644
--- a/src/output/tests/writer-mpi.cpp
+++ b/src/output/tests/writer-mpi.cpp
@@ -22,7 +22,8 @@ void cleanup() {
 }
 
 #define CEILDIV(a, b)                                                          \
-  (static_cast<int>(math::ceil(static_cast<real_t>(a) / static_cast<real_t>(b))))
+  (static_cast<std::size_t>(                                                   \
+    math::ceil(static_cast<real_t>(a) / static_cast<real_t>(b))))
 
 auto main(int argc, char* argv[]) -> int {
   Kokkos::initialize(argc, argv);
@@ -61,8 +62,8 @@ auto main(int argc, char* argv[]) -> int {
       // write
       auto writer = out::Writer();
       writer.init(&adios, "hdf5", "test", false);
-      writer.defineMeshLayout({ static_cast<unsigned long>(mpi_size) * nx1 },
-                              { static_cast<unsigned long>(mpi_rank) * nx1 },
+      writer.defineMeshLayout({ static_cast<ncells_t>(mpi_size) * nx1 },
+                              { static_cast<ncells_t>(mpi_rank) * nx1 },
                               { nx1 },
                               { dwn1 },
                               false,
@@ -89,9 +90,9 @@ auto main(int argc, char* argv[]) -> int {
     {
       // read
       adios2::IO io = adios.DeclareIO("read-test");
-      io.SetEngine("hdf5");
-      adios2::Engine reader = io.Open("test.h5", adios2::Mode::Read, MPI_COMM_SELF);
-      raise::ErrorIf(io.InquireAttribute<unsigned int>("NGhosts").Data()[0] != 0,
+      io.SetEngine("HDF5");
+      adios2::Engine reader = io.Open("test.h5", adios2::Mode::Read);
+      raise::ErrorIf(io.InquireAttribute<std::size_t>("NGhosts").Data()[0] != 0,
                      "NGhosts is not correct",
                      HERE);
       raise::ErrorIf(io.InquireAttribute<std::size_t>("Dimension").Data()[0] != 1,
@@ -99,13 +100,13 @@ auto main(int argc, char* argv[]) -> int {
                      HERE);
       for (std::size_t step = 0; reader.BeginStep() == adios2::StepStatus::OK;
            ++step) {
-        std::size_t step_read;
-        long double time_read;
+        timestep_t step_read;
+        simtime_t  time_read;
 
-        reader.Get(io.InquireVariable<std::size_t>("Step"),
+        reader.Get(io.InquireVariable<timestep_t>("Step"),
                    &step_read,
                    adios2::Mode::Sync);
-        reader.Get(io.InquireVariable<long double>("Time"),
+        reader.Get(io.InquireVariable<simtime_t>("Time"),
                    &time_read,
                    adios2::Mode::Sync);
         raise::ErrorIf(step_read != step, "Step is not correct", HERE);
@@ -115,16 +116,15 @@ auto main(int argc, char* argv[]) -> int {
 
         const auto l_size   = nx1;
         const auto l_offset = nx1 * mpi_rank;
-        const auto g_size   = nx1 * mpi_size;
 
         const double n = l_size;
         const double d = dwn1;
         const double l = l_offset;
         const double f = math::ceil(l / d) * d - l;
 
-        const auto first_cell = static_cast<std::size_t>(f);
-        const auto l_size_dwn = static_cast<std::size_t>(math::ceil((n - f) / d));
-        const auto l_corner_dwn = static_cast<std::size_t>(math::ceil(l / d));
+        const auto first_cell = static_cast<ncells_t>(f);
+        const auto l_size_dwn = static_cast<ncells_t>(math::ceil((n - f) / d));
+        const auto l_corner_dwn = static_cast<ncells_t>(math::ceil(l / d));
 
         array_t<real_t*> field_read {};
         int              cntr = 0;
diff --git a/src/output/tests/writer-nompi.cpp b/src/output/tests/writer-nompi.cpp
index b063539ca..593f37f92 100644
--- a/src/output/tests/writer-nompi.cpp
+++ b/src/output/tests/writer-nompi.cpp
@@ -23,7 +23,8 @@ void cleanup() {
 }
 
 #define CEILDIV(a, b)                                                          \
-  (static_cast<int>(math::ceil(static_cast<real_t>(a) / static_cast<real_t>(b))))
+  (static_cast<ncells_t>(                                                      \
+    math::ceil(static_cast<real_t>(a) / static_cast<real_t>(b))))
 
 auto main(int argc, char* argv[]) -> int {
   Kokkos::initialize(argc, argv);