diff --git a/src/core/grid_based_algorithms/lb_interface.cpp b/src/core/grid_based_algorithms/lb_interface.cpp
index 4922376883..076ede6468 100644
--- a/src/core/grid_based_algorithms/lb_interface.cpp
+++ b/src/core/grid_based_algorithms/lb_interface.cpp
@@ -26,15 +26,14 @@
 #include "grid.hpp"
 
 #include <utils/Vector.hpp>
-#include <utils/mpi/reduce_optional.hpp>
 
-#include <boost/optional.hpp>
-#include <boost/serialization/access.hpp>
-#include <boost/serialization/vector.hpp>
+#include <boost/mpi/collectives/reduce.hpp>
+#include <boost/mpi/communicator.hpp>
 
 #include <cmath>
 #include <functional>
 #include <limits>
+#include <optional>
 #include <stdexcept>
 #include <string>
 #include <vector>
@@ -147,7 +146,7 @@ Utils::Vector3d calc_fluid_momentum() {
   throw NoLBActive();
 }
 
-boost::optional<Utils::Vector3d>
+std::optional<Utils::Vector3d>
 get_interpolated_velocity(Utils::Vector3d const &pos) {
   if (lattice_switch == ActiveLB::WALBERLA_LB) {
 #ifdef WALBERLA
@@ -158,7 +157,7 @@ get_interpolated_velocity(Utils::Vector3d const &pos) {
   throw NoLBActive();
 }
 
-boost::optional<double> get_interpolated_density(Utils::Vector3d const &pos) {
+std::optional<double> get_interpolated_density(Utils::Vector3d const &pos) {
   if (lattice_switch == ActiveLB::WALBERLA_LB) {
 #ifdef WALBERLA
     auto const folded_pos = folded_position(pos, box_geo) / get_agrid();
diff --git a/src/core/grid_based_algorithms/lb_interface.hpp b/src/core/grid_based_algorithms/lb_interface.hpp
index 9759cb0272..892ded2ec3 100644
--- a/src/core/grid_based_algorithms/lb_interface.hpp
+++ b/src/core/grid_based_algorithms/lb_interface.hpp
@@ -21,8 +21,7 @@
 
 #include <utils/Vector.hpp>
 
-#include <boost/optional.hpp>
-
+#include <optional>
 #include <stdexcept>
 
 // Forward Declarations
@@ -116,7 +115,7 @@ Utils::Vector3d calc_fluid_momentum();
  * @param pos Position at which the velocity is to be calculated.
  * @retval optional interpolated fluid velocity.
  */
-boost::optional<Utils::Vector3d>
+std::optional<Utils::Vector3d>
 get_interpolated_velocity(Utils::Vector3d const &pos);
 
 /**
@@ -125,7 +124,7 @@ get_interpolated_velocity(Utils::Vector3d const &pos);
  * @param pos Position at which the density is to be calculated.
  * @retval optional interpolated fluid density.
  */
-boost::optional<double> get_interpolated_density(Utils::Vector3d const &pos);
+std::optional<double> get_interpolated_density(Utils::Vector3d const &pos);
 
 } // namespace LB
 
diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp
index 2a89b544c2..e0c6e6f99d 100644
--- a/src/script_interface/particle_data/ParticleHandle.cpp
+++ b/src/script_interface/particle_data/ParticleHandle.cpp
@@ -47,6 +47,7 @@
 #include <cmath>
 #include <cstddef>
 #include <memory>
+#include <optional>
 #include <set>
 #include <sstream>
 #include <stdexcept>
@@ -134,7 +135,7 @@ template <typename T, class F>
 T ParticleHandle::get_particle_property(F const &fun) const {
   auto const &comm = context()->get_comm();
   auto const ptr = const_cast<Particle const *>(get_real_particle(comm, m_pid));
-  boost::optional<T> ret;
+  std::optional<T> ret;
   if (ptr == nullptr) {
     ret = {};
   } else {
diff --git a/src/script_interface/walberla/EKContainer.hpp b/src/script_interface/walberla/EKContainer.hpp
index fbae513e74..e5f41554d2 100644
--- a/src/script_interface/walberla/EKContainer.hpp
+++ b/src/script_interface/walberla/EKContainer.hpp
@@ -32,7 +32,9 @@
 #include <script_interface/ObjectList.hpp>
 #include <script_interface/ScriptInterface.hpp>
 
+#include <cassert>
 #include <memory>
+#include <optional>
 #include <string>
 
 namespace ScriptInterface::walberla {
@@ -132,7 +134,7 @@ class EKContainer : public ObjectList<EKSpecies> {
   }
 
   void set_solver(Variant const &solver_variant) {
-    boost::optional<decltype(m_poisson_solver)> solver;
+    std::optional<decltype(m_poisson_solver)> solver;
     if (is_none(solver_variant)) {
       solver = std::shared_ptr<EKNone>();
     } else {
@@ -146,6 +148,7 @@ class EKContainer : public ObjectList<EKSpecies> {
         solver = get_value<std::shared_ptr<EKNone>>(solver_variant);
       }
     }
+    assert(solver.has_value());
     m_poisson_solver = *solver;
     auto const visitor = GetPoissonSolverInstance();
     auto const instance = boost::apply_visitor(visitor, m_poisson_solver);
diff --git a/src/script_interface/walberla/EKSpecies.cpp b/src/script_interface/walberla/EKSpecies.cpp
index 172f3c3bce..9f908ad15b 100644
--- a/src/script_interface/walberla/EKSpecies.cpp
+++ b/src/script_interface/walberla/EKSpecies.cpp
@@ -21,21 +21,21 @@
 #ifdef WALBERLA
 
 #include "EKSpecies.hpp"
+#include "EKWalberlaNodeState.hpp"
 #include "WalberlaCheckpoint.hpp"
 
 #include <walberla_bridge/LatticeWalberla.hpp>
-#include <walberla_bridge/electrokinetics/EKWalberlaNodeState.hpp>
 #include <walberla_bridge/electrokinetics/ek_walberla_init.hpp>
 
 #include <boost/mpi.hpp>
 #include <boost/mpi/collectives/all_reduce.hpp>
 #include <boost/mpi/collectives/broadcast.hpp>
-#include <boost/optional.hpp>
 
 #include <algorithm>
 #include <cassert>
 #include <functional>
 #include <memory>
+#include <optional>
 #include <sstream>
 #include <stdexcept>
 #include <string>
@@ -207,8 +207,8 @@ void EKSpecies::save_checkpoint(std::string const &filename, int mode) {
   auto const write_data = [&ek_obj,
                            mode](std::shared_ptr<CheckpointFile> cpfile_ptr,
                                  Context const &context) {
-    auto const get_node_checkpoint = [&](Utils::Vector3i const &ind)
-        -> boost::optional<EKWalberlaNodeState> {
+    auto const get_node_checkpoint =
+        [&](Utils::Vector3i const &ind) -> std::optional<EKWalberlaNodeState> {
       auto const density = ek_obj.get_node_density(ind);
       auto const is_b_d = ek_obj.get_node_is_density_boundary(ind);
       auto const dens_b = ek_obj.get_node_density_at_boundary(ind);
@@ -229,7 +229,7 @@ void EKSpecies::save_checkpoint(std::string const &filename, int mode) {
         }
         return cpnode;
       }
-      return {boost::none};
+      return std::nullopt;
     };
 
     auto failure = false;
diff --git a/src/script_interface/walberla/EKSpeciesNode.cpp b/src/script_interface/walberla/EKSpeciesNode.cpp
index b367deda65..61433ffafd 100644
--- a/src/script_interface/walberla/EKSpeciesNode.cpp
+++ b/src/script_interface/walberla/EKSpeciesNode.cpp
@@ -32,18 +32,18 @@
 #include <utils/mpi/reduce_optional.hpp>
 
 #include <boost/mpi/collectives/all_reduce.hpp>
-#include <boost/optional.hpp>
 #include <boost/serialization/vector.hpp>
 
 #include <cassert>
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <string>
 
 namespace ScriptInterface::walberla {
 
 static bool is_boundary_all_reduce(boost::mpi::communicator const &comm,
-                                   boost::optional<bool> const &is_boundary) {
+                                   std::optional<bool> const &is_boundary) {
   return boost::mpi::all_reduce(comm, is_boundary ? *is_boundary : false,
                                 std::logical_or<>());
 }
diff --git a/src/script_interface/walberla/EKSpeciesSlice.hpp b/src/script_interface/walberla/EKSpeciesSlice.hpp
index ce67ca31ee..d16688729f 100644
--- a/src/script_interface/walberla/EKSpeciesSlice.hpp
+++ b/src/script_interface/walberla/EKSpeciesSlice.hpp
@@ -38,6 +38,7 @@
 
 #include <cassert>
 #include <memory>
+#include <optional>
 #include <string>
 #include <type_traits>
 #include <unordered_map>
@@ -45,8 +46,8 @@
 
 namespace ScriptInterface::walberla {
 
-using DensityBoundaryType = boost::optional<double>;
-using FluxBoundaryType = boost::optional<Utils::Vector3d>;
+using DensityBoundaryType = std::optional<double>;
+using FluxBoundaryType = std::optional<Utils::Vector3d>;
 
 struct EKFieldSerializer {
 
@@ -77,7 +78,7 @@ struct EKFieldSerializer {
       auto const vector_variants = get_value<std::vector<Variant>>(variant);
       for (auto const &value : vector_variants) {
         if (is_none(value)) {
-          values.emplace_back(boost::none);
+          values.emplace_back(std::nullopt);
         } else {
           values.emplace_back(get_value<Utils::Vector3d>(value));
         }
@@ -86,7 +87,7 @@ struct EKFieldSerializer {
       auto const vector_variants = get_value<std::vector<Variant>>(variant);
       for (auto const &value : vector_variants) {
         if (is_none(value)) {
-          values.emplace_back(boost::none);
+          values.emplace_back(std::nullopt);
         } else {
           values.emplace_back(get_value<double>(value));
         }
diff --git a/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKWalberlaNodeState.hpp b/src/script_interface/walberla/EKWalberlaNodeState.hpp
similarity index 91%
rename from src/walberla_bridge/include/walberla_bridge/electrokinetics/EKWalberlaNodeState.hpp
rename to src/script_interface/walberla/EKWalberlaNodeState.hpp
index 6c4987c541..2dc0f8db04 100644
--- a/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKWalberlaNodeState.hpp
+++ b/src/script_interface/walberla/EKWalberlaNodeState.hpp
@@ -21,6 +21,8 @@
 
 #include <utils/Vector.hpp>
 
+#include <boost/serialization/access.hpp>
+
 /** Checkpoint data for a EK node. */
 struct EKWalberlaNodeState {
   double density;
@@ -32,7 +34,7 @@ struct EKWalberlaNodeState {
 private:
   friend boost::serialization::access;
   template <typename Archive>
-  void serialize(Archive &ar, long int /* version */) {
+  void serialize(Archive &ar, unsigned int /* version */) {
     ar &density &is_boundary_density &density_boundary &is_boundary_flux
         &flux_boundary;
   }
diff --git a/src/script_interface/walberla/LBFluid.cpp b/src/script_interface/walberla/LBFluid.cpp
index c70c3f82a3..39bd39f038 100644
--- a/src/script_interface/walberla/LBFluid.cpp
+++ b/src/script_interface/walberla/LBFluid.cpp
@@ -21,6 +21,7 @@
 #ifdef WALBERLA
 
 #include "LBFluid.hpp"
+#include "LBWalberlaNodeState.hpp"
 #include "WalberlaCheckpoint.hpp"
 
 #include "core/BoxGeometry.hpp"
@@ -34,7 +35,6 @@
 #include <script_interface/communication.hpp>
 
 #include <walberla_bridge/LatticeWalberla.hpp>
-#include <walberla_bridge/lattice_boltzmann/LBWalberlaNodeState.hpp>
 #include <walberla_bridge/lattice_boltzmann/LeesEdwardsPack.hpp>
 #include <walberla_bridge/lattice_boltzmann/lb_walberla_init.hpp>
 
@@ -45,13 +45,13 @@
 #include <boost/mpi.hpp>
 #include <boost/mpi/collectives/all_reduce.hpp>
 #include <boost/mpi/collectives/broadcast.hpp>
-#include <boost/optional.hpp>
 
 #include <algorithm>
 #include <cassert>
 #include <cstddef>
 #include <functional>
 #include <memory>
+#include <optional>
 #include <sstream>
 #include <stdexcept>
 #include <string>
@@ -299,8 +299,8 @@ void LBFluid::save_checkpoint(std::string const &filename, int mode) {
   auto const write_data = [&lb_obj,
                            mode](std::shared_ptr<CheckpointFile> cpfile_ptr,
                                  Context const &context) {
-    auto const get_node_checkpoint = [&](Utils::Vector3i const &ind)
-        -> boost::optional<LBWalberlaNodeState> {
+    auto const get_node_checkpoint =
+        [&](Utils::Vector3i const &ind) -> std::optional<LBWalberlaNodeState> {
       auto const pop = lb_obj.get_node_population(ind);
       auto const laf = lb_obj.get_node_last_applied_force(ind);
       auto const lbb = lb_obj.get_node_is_boundary(ind);
@@ -313,9 +313,9 @@ void LBFluid::save_checkpoint(std::string const &filename, int mode) {
         if (*lbb) {
           cpnode.slip_velocity = *vbb;
         }
-        return cpnode;
+        return {cpnode};
       }
-      return {boost::none};
+      return std::nullopt;
     };
 
     auto failure = false;
diff --git a/src/script_interface/walberla/LBFluid.hpp b/src/script_interface/walberla/LBFluid.hpp
index 43477e187e..a726108624 100644
--- a/src/script_interface/walberla/LBFluid.hpp
+++ b/src/script_interface/walberla/LBFluid.hpp
@@ -33,7 +33,6 @@
 
 #include <walberla_bridge/LatticeModel.hpp>
 #include <walberla_bridge/lattice_boltzmann/LBWalberlaBase.hpp>
-#include <walberla_bridge/lattice_boltzmann/LBWalberlaNodeState.hpp>
 
 #include <utils/Vector.hpp>
 #include <utils/math/int_pow.hpp>
diff --git a/src/script_interface/walberla/LBFluidNode.cpp b/src/script_interface/walberla/LBFluidNode.cpp
index aeb41837e6..4cab29bb23 100644
--- a/src/script_interface/walberla/LBFluidNode.cpp
+++ b/src/script_interface/walberla/LBFluidNode.cpp
@@ -29,15 +29,15 @@
 
 #include <boost/mpi/collectives/all_reduce.hpp>
 #include <boost/mpi/communicator.hpp>
-#include <boost/optional.hpp>
 
+#include <optional>
 #include <string>
 #include <vector>
 
 namespace ScriptInterface::walberla {
 
 static bool is_boundary_all_reduce(boost::mpi::communicator const &comm,
-                                   boost::optional<bool> const &is_boundary) {
+                                   std::optional<bool> const &is_boundary) {
   return boost::mpi::all_reduce(comm, is_boundary ? *is_boundary : false,
                                 std::logical_or<>());
 }
@@ -113,7 +113,7 @@ Variant LBFluidNode::do_call_method(std::string const &name,
   }
   if (name == "get_pressure_tensor" or name == "get_pressure_tensor_neq") {
     auto const result = m_lb_fluid->get_node_pressure_tensor(m_index);
-    auto value = boost::optional<std::vector<double>>{};
+    auto value = std::optional<std::vector<double>>{};
     if (result) {
       value = (*result / m_conv_press).as_vector();
     }
diff --git a/src/script_interface/walberla/LBFluidSlice.hpp b/src/script_interface/walberla/LBFluidSlice.hpp
index c73f28bfae..4701ee64bc 100644
--- a/src/script_interface/walberla/LBFluidSlice.hpp
+++ b/src/script_interface/walberla/LBFluidSlice.hpp
@@ -38,6 +38,7 @@
 
 #include <cassert>
 #include <memory>
+#include <optional>
 #include <string>
 #include <type_traits>
 #include <unordered_map>
@@ -45,7 +46,7 @@
 
 namespace ScriptInterface::walberla {
 
-using VelocityBounceBackType = boost::optional<Utils::Vector3d>;
+using VelocityBounceBackType = std::optional<Utils::Vector3d>;
 
 struct LBFieldSerializer {
   template <typename T> static Variant serialize(std::vector<T> const &values) {
@@ -74,7 +75,7 @@ struct LBFieldSerializer {
       auto const vector_variants = get_value<std::vector<Variant>>(variant);
       for (auto const &value : vector_variants) {
         if (is_none(value)) {
-          values.emplace_back(boost::none);
+          values.emplace_back(std::nullopt);
         } else {
           values.emplace_back(get_value<Utils::Vector3d>(value));
         }
diff --git a/src/walberla_bridge/include/walberla_bridge/lattice_boltzmann/LBWalberlaNodeState.hpp b/src/script_interface/walberla/LBWalberlaNodeState.hpp
similarity index 91%
rename from src/walberla_bridge/include/walberla_bridge/lattice_boltzmann/LBWalberlaNodeState.hpp
rename to src/script_interface/walberla/LBWalberlaNodeState.hpp
index 020ffa2cd3..a47099b3cc 100644
--- a/src/walberla_bridge/include/walberla_bridge/lattice_boltzmann/LBWalberlaNodeState.hpp
+++ b/src/script_interface/walberla/LBWalberlaNodeState.hpp
@@ -21,6 +21,8 @@
 
 #include <utils/Vector.hpp>
 
+#include <boost/serialization/access.hpp>
+
 #include <vector>
 
 /** Checkpoint data for a LB node. */
@@ -33,7 +35,7 @@ struct LBWalberlaNodeState {
 private:
   friend boost::serialization::access;
   template <typename Archive>
-  void serialize(Archive &ar, long int /* version */) {
+  void serialize(Archive &ar, unsigned int /* version */) {
     ar &populations &last_applied_force &slip_velocity &is_boundary;
   }
 };
diff --git a/src/script_interface/walberla/LatticeSlice.impl.hpp b/src/script_interface/walberla/LatticeSlice.impl.hpp
index 5d81b24773..0f0ac558e4 100644
--- a/src/script_interface/walberla/LatticeSlice.impl.hpp
+++ b/src/script_interface/walberla/LatticeSlice.impl.hpp
@@ -20,17 +20,17 @@
 #pragma once
 
 #include <utils/Vector.hpp>
+#include <utils/serialization/optional.hpp>
 
 #include <boost/mpi/collectives/gather.hpp>
 #include <boost/mpi/collectives/scatter.hpp>
 #include <boost/mpi/communicator.hpp>
 #include <boost/multi_array.hpp>
-#include <boost/optional.hpp>
-#include <boost/serialization/optional.hpp>
 #include <boost/serialization/vector.hpp>
 
 #include <algorithm>
 #include <functional>
+#include <optional>
 #include <stdexcept>
 #include <tuple>
 #include <type_traits>
@@ -40,10 +40,10 @@ namespace ScriptInterface::walberla {
 
 namespace detail {
 
-// boundary types are always boost::optional types
+// boundary types are always std::optional types
 template <class> struct is_optional_type : public std::false_type {};
 template <class T>
-struct is_optional_type<boost::optional<T>> : public std::true_type {};
+struct is_optional_type<std::optional<T>> : public std::true_type {};
 
 template <class ArrayView, typename T>
 void unflatten_grid(ArrayView &view, std::vector<T> const &values) {
@@ -84,7 +84,7 @@ void flatten_grid(ArrayView const &view, std::vector<T> &out,
             if (view[i][j][k][t]) {
               out.emplace_back(*(view[i][j][k][t]) * units_conversion);
             } else {
-              out.emplace_back(boost::none);
+              out.emplace_back(std::nullopt);
             }
           } else {
             out.emplace_back(view[i][j][k][t]);
diff --git a/src/utils/include/utils/mpi/reduce_optional.hpp b/src/utils/include/utils/mpi/reduce_optional.hpp
index 216a234c36..6ce76b4e4e 100644
--- a/src/utils/include/utils/mpi/reduce_optional.hpp
+++ b/src/utils/include/utils/mpi/reduce_optional.hpp
@@ -21,10 +21,10 @@
 
 #include <boost/mpi/collectives/all_reduce.hpp>
 #include <boost/mpi/communicator.hpp>
-#include <boost/optional.hpp>
 
 #include <cassert>
 #include <functional>
+#include <optional>
 
 namespace Utils::Mpi {
 
@@ -34,7 +34,7 @@ namespace Utils::Mpi {
  */
 template <typename T>
 T reduce_optional(boost::mpi::communicator const &comm,
-                  boost::optional<T> const &result) {
+                  std::optional<T> const &result) {
   assert(1 == boost::mpi::all_reduce(comm, static_cast<int>(!!result),
                                      std::plus<>()) &&
          "Incorrect number of return values");
diff --git a/src/utils/tests/reduce_optional_test.cpp b/src/utils/tests/reduce_optional_test.cpp
index 1a95151d8d..04ce72feda 100644
--- a/src/utils/tests/reduce_optional_test.cpp
+++ b/src/utils/tests/reduce_optional_test.cpp
@@ -25,13 +25,14 @@
 #include <utils/mpi/reduce_optional.hpp>
 
 #include <boost/mpi.hpp>
-#include <boost/optional.hpp>
+
+#include <optional>
 
 BOOST_AUTO_TEST_CASE(reduce_optional) {
   boost::mpi::communicator comm;
 
   for (int rank = 0; rank < comm.size(); ++rank) {
-    boost::optional<int> maybe;
+    std::optional<int> maybe;
     if (comm.rank() == rank) {
       maybe = 42;
     }
diff --git a/src/walberla_bridge/CMakeLists.txt b/src/walberla_bridge/CMakeLists.txt
index 78040cab5f..f7912560ff 100644
--- a/src/walberla_bridge/CMakeLists.txt
+++ b/src/walberla_bridge/CMakeLists.txt
@@ -51,8 +51,7 @@ endif()
 
 target_link_libraries(
   espresso_walberla PUBLIC MPI::MPI_CXX espresso::utils
-  PRIVATE Boost::boost espresso::cpp_flags espresso::walberla::cpp_flags
-          ${WALBERLA_LIBS})
+  PRIVATE espresso::cpp_flags espresso::walberla::cpp_flags ${WALBERLA_LIBS})
 target_include_directories(
   espresso_walberla PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
   PRIVATE ${WALBERLA_INCLUDE_DIRS} ${walberla_BINARY_DIR}/src)
diff --git a/src/walberla_bridge/include/walberla_bridge/BlockAndCell.hpp b/src/walberla_bridge/include/walberla_bridge/BlockAndCell.hpp
index 1ad71b75b5..908e5fba4b 100644
--- a/src/walberla_bridge/include/walberla_bridge/BlockAndCell.hpp
+++ b/src/walberla_bridge/include/walberla_bridge/BlockAndCell.hpp
@@ -26,9 +26,8 @@
 
 #include "LatticeWalberla.hpp"
 
-#include <boost/optional.hpp>
-
 #include <memory>
+#include <optional>
 
 namespace walberla {
 // Helpers to retrieve blocks and cells
@@ -52,7 +51,7 @@ IBlock *get_block_extended(std::shared_ptr<StructuredBlockForest> const &blocks,
   return nullptr;
 }
 
-inline boost::optional<BlockAndCell>
+inline std::optional<BlockAndCell>
 get_block_and_cell(::LatticeWalberla const &lattice,
                    Utils::Vector3i const &node, bool consider_ghost_layers) {
   // Get block and local cell
@@ -65,7 +64,7 @@ get_block_and_cell(::LatticeWalberla const &lattice,
     block = get_block_extended(blocks, node, lattice.get_ghost_layers());
   }
   if (!block)
-    return {boost::none};
+    return std::nullopt;
 
   // Transform coords to block local
   Cell local_cell;
diff --git a/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKinWalberlaBase.hpp b/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKinWalberlaBase.hpp
index baa6a65cba..56932530fe 100644
--- a/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKinWalberlaBase.hpp
+++ b/src/walberla_bridge/include/walberla_bridge/electrokinetics/EKinWalberlaBase.hpp
@@ -23,9 +23,8 @@
 
 #include <utils/Vector.hpp>
 
-#include <boost/optional.hpp>
-
 #include <cstddef>
+#include <optional>
 #include <vector>
 
 /** @brief Interface of a lattice-based electrokinetic model. */
@@ -46,7 +45,7 @@ class EKinWalberlaBase : public LatticeModel {
                                 double density) = 0;
 
   /** @brief Get node density. */
-  [[nodiscard]] virtual boost::optional<double>
+  [[nodiscard]] virtual std::optional<double>
   get_node_density(Utils::Vector3i const &node,
                    bool consider_ghosts = false) const = 0;
 
@@ -65,17 +64,17 @@ class EKinWalberlaBase : public LatticeModel {
                                       Utils::Vector3d const &flux) = 0;
 
   /** @brief Get node flux boundary conditions. */
-  [[nodiscard]] virtual boost::optional<Utils::Vector3d>
+  [[nodiscard]] virtual std::optional<Utils::Vector3d>
   get_node_flux_at_boundary(Utils::Vector3i const &node,
                             bool consider_ghosts = false) const = 0;
 
   /** @brief Set slice flux boundary conditions. */
   virtual void set_slice_flux_boundary(
       Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner,
-      std::vector<boost::optional<Utils::Vector3d>> const &flux) = 0;
+      std::vector<std::optional<Utils::Vector3d>> const &flux) = 0;
 
   /** @brief Get slice flux boundary conditions. */
-  [[nodiscard]] virtual std::vector<boost::optional<Utils::Vector3d>>
+  [[nodiscard]] virtual std::vector<std::optional<Utils::Vector3d>>
   get_slice_flux_at_boundary(Utils::Vector3i const &lower_corner,
                              Utils::Vector3i const &upper_corner) const = 0;
 
@@ -86,17 +85,17 @@ class EKinWalberlaBase : public LatticeModel {
                                          double density) = 0;
 
   /** @brief Get node density boundary conditions. */
-  [[nodiscard]] virtual boost::optional<double>
+  [[nodiscard]] virtual std::optional<double>
   get_node_density_at_boundary(Utils::Vector3i const &node,
                                bool consider_ghosts = false) const = 0;
 
   /** @brief Set slice density boundary conditions. */
   virtual void set_slice_density_boundary(
       Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner,
-      std::vector<boost::optional<double>> const &density) = 0;
+      std::vector<std::optional<double>> const &density) = 0;
 
   /** @brief Get slice density boundary conditions. */
-  [[nodiscard]] virtual std::vector<boost::optional<double>>
+  [[nodiscard]] virtual std::vector<std::optional<double>>
   get_slice_density_at_boundary(Utils::Vector3i const &lower_corner,
                                 Utils::Vector3i const &upper_corner) const = 0;
 
@@ -104,17 +103,17 @@ class EKinWalberlaBase : public LatticeModel {
   remove_node_from_density_boundary(Utils::Vector3i const &node) = 0;
 
   /** @brief Check if node has flux boundary conditions. */
-  [[nodiscard]] virtual boost::optional<bool>
+  [[nodiscard]] virtual std::optional<bool>
   get_node_is_flux_boundary(Utils::Vector3i const &node,
                             bool consider_ghosts = false) const = 0;
 
   /** @brief Check if node has density boundary conditions. */
-  [[nodiscard]] virtual boost::optional<bool>
+  [[nodiscard]] virtual std::optional<bool>
   get_node_is_density_boundary(Utils::Vector3i const &node,
                                bool consider_ghosts = false) const = 0;
 
   /** @brief Check if node has any boundary conditions. */
-  [[nodiscard]] virtual boost::optional<bool>
+  [[nodiscard]] virtual std::optional<bool>
   get_node_is_boundary(Utils::Vector3i const &node,
                        bool consider_ghosts = false) const = 0;
 
diff --git a/src/walberla_bridge/include/walberla_bridge/lattice_boltzmann/LBWalberlaBase.hpp b/src/walberla_bridge/include/walberla_bridge/lattice_boltzmann/LBWalberlaBase.hpp
index 1b74ea3863..622a256103 100644
--- a/src/walberla_bridge/include/walberla_bridge/lattice_boltzmann/LBWalberlaBase.hpp
+++ b/src/walberla_bridge/include/walberla_bridge/lattice_boltzmann/LBWalberlaBase.hpp
@@ -32,10 +32,9 @@
 
 #include <utils/Vector.hpp>
 
-#include <boost/optional.hpp>
-
 #include <cstddef>
 #include <memory>
+#include <optional>
 #include <vector>
 
 /** @brief Interface of a lattice-based fluid model. */
@@ -56,12 +55,12 @@ class LBWalberlaBase : public LatticeModel {
   [[nodiscard]] virtual bool is_double_precision() const noexcept = 0;
 
   /** @brief Get interpolated velocities at a position. */
-  virtual boost::optional<Utils::Vector3d>
+  virtual std::optional<Utils::Vector3d>
   get_velocity_at_pos(Utils::Vector3d const &position,
                       bool consider_points_in_halo = false) const = 0;
 
   /** @brief Get interpolated densities at a position. */
-  virtual boost::optional<double>
+  virtual std::optional<double>
   get_density_at_pos(Utils::Vector3d const &position,
                      bool consider_points_in_halo = false) const = 0;
 
@@ -73,11 +72,11 @@ class LBWalberlaBase : public LatticeModel {
                                 Utils::Vector3d const &force) = 0;
 
   /** @brief Get stored force to be applied on node in the next time step. */
-  virtual boost::optional<Utils::Vector3d>
+  virtual std::optional<Utils::Vector3d>
   get_node_force_to_be_applied(Utils::Vector3i const &node) const = 0;
 
   /** @brief Get stored force that was applied on node in the last time step. */
-  virtual boost::optional<Utils::Vector3d>
+  virtual std::optional<Utils::Vector3d>
   get_node_last_applied_force(Utils::Vector3i const &node,
                               bool consider_ghosts = false) const = 0;
 
@@ -99,7 +98,7 @@ class LBWalberlaBase : public LatticeModel {
                                std::vector<double> const &force) = 0;
 
   /** @brief Get node population. */
-  virtual boost::optional<std::vector<double>>
+  virtual std::optional<std::vector<double>>
   get_node_population(Utils::Vector3i const &node,
                       bool consider_ghosts = false) const = 0;
 
@@ -118,7 +117,7 @@ class LBWalberlaBase : public LatticeModel {
                                     std::vector<double> const &population) = 0;
 
   /** @brief Get node velocity. */
-  virtual boost::optional<Utils::Vector3d>
+  virtual std::optional<Utils::Vector3d>
   get_node_velocity(Utils::Vector3i const &node,
                     bool consider_ghosts = false) const = 0;
 
@@ -137,7 +136,7 @@ class LBWalberlaBase : public LatticeModel {
                                   std::vector<double> const &velocity) = 0;
 
   /** @brief Get node density. */
-  virtual boost::optional<double>
+  virtual std::optional<double>
   get_node_density(Utils::Vector3i const &node,
                    bool consider_ghosts = false) const = 0;
 
@@ -156,7 +155,7 @@ class LBWalberlaBase : public LatticeModel {
                                  std::vector<double> const &density) = 0;
 
   /** @brief Get node velocity boundary conditions. */
-  virtual boost::optional<Utils::Vector3d>
+  virtual std::optional<Utils::Vector3d>
   get_node_velocity_at_boundary(Utils::Vector3i const &node,
                                 bool consider_ghosts = false) const = 0;
 
@@ -166,24 +165,24 @@ class LBWalberlaBase : public LatticeModel {
                                 Utils::Vector3d const &velocity) = 0;
 
   /** @brief Get slice velocity boundary conditions. */
-  virtual std::vector<boost::optional<Utils::Vector3d>>
+  virtual std::vector<std::optional<Utils::Vector3d>>
   get_slice_velocity_at_boundary(Utils::Vector3i const &lower_corner,
                                  Utils::Vector3i const &upper_corner) const = 0;
 
   /** @brief Set slice velocity boundary conditions. */
   virtual void set_slice_velocity_at_boundary(
       Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner,
-      std::vector<boost::optional<Utils::Vector3d>> const &velocity) = 0;
+      std::vector<std::optional<Utils::Vector3d>> const &velocity) = 0;
 
   /** @brief Get (stored) force applied on node due to boundary condition. */
-  virtual boost::optional<Utils::Vector3d>
+  virtual std::optional<Utils::Vector3d>
   get_node_boundary_force(Utils::Vector3i const &node) const = 0;
 
   /** @brief Remove a node from the boundaries. */
   virtual bool remove_node_from_boundary(Utils::Vector3i const &node) = 0;
 
   /** @brief Check if node has velocity boundary conditions. */
-  virtual boost::optional<bool>
+  virtual std::optional<bool>
   get_node_is_boundary(Utils::Vector3i const &node,
                        bool consider_ghosts = false) const = 0;
 
@@ -214,7 +213,7 @@ class LBWalberlaBase : public LatticeModel {
                           unsigned int shear_plane_normal) const = 0;
 
   /** @brief Get node pressure tensor. */
-  virtual boost::optional<Utils::VectorXd<9>>
+  virtual std::optional<Utils::VectorXd<9>>
   get_node_pressure_tensor(Utils::Vector3i const &node) const = 0;
 
   /** @brief Get slice pressure tensor. */
@@ -247,7 +246,7 @@ class LBWalberlaBase : public LatticeModel {
   virtual double get_kT() const noexcept = 0;
 
   /** @brief Set the RNG counter (if thermalized). */
-  [[nodiscard]] virtual boost::optional<uint64_t> get_rng_state() const = 0;
+  [[nodiscard]] virtual std::optional<uint64_t> get_rng_state() const = 0;
 
   /** @brief Set the rng state of thermalized LBs */
   virtual void set_rng_state(uint64_t counter) = 0;
diff --git a/src/walberla_bridge/include/walberla_bridge/utils/boundary_utils.hpp b/src/walberla_bridge/include/walberla_bridge/utils/boundary_utils.hpp
index 98150003f7..f91a6c833a 100644
--- a/src/walberla_bridge/include/walberla_bridge/utils/boundary_utils.hpp
+++ b/src/walberla_bridge/include/walberla_bridge/utils/boundary_utils.hpp
@@ -23,8 +23,6 @@
 
 #include <utils/Vector.hpp>
 
-#include <boost/multi_array/multi_array_ref.hpp>
-
 #include <cassert>
 #include <cstddef>
 #include <iterator>
@@ -82,32 +80,31 @@ void set_boundary_from_grid(BoundaryModel &boundary,
                             LatticeWalberla const &lattice,
                             std::vector<int> const &raster_flat,
                             std::vector<DataType> const &data_flat) {
-  // reshape grids
+
   auto const grid_size = lattice.get_grid_dimensions();
+  auto const offset = lattice.get_local_grid_range().first;
+  auto const gl = static_cast<int>(lattice.get_ghost_layers());
   assert(raster_flat.size() == Utils::product(grid_size));
-  boost::const_multi_array_ref<DataType, 3> data_grid(data_flat.data(),
-                                                      grid_size);
-  boost::const_multi_array_ref<int, 3> raster(raster_flat.data(), grid_size);
+  auto const n_y = grid_size[1];
+  auto const n_z = grid_size[2];
+  auto const off_i = offset[0];
+  auto const off_j = offset[1];
+  auto const off_k = offset[2];
 
   auto const &blocks = lattice.get_blocks();
   for (auto block = blocks->begin(); block != blocks->end(); ++block) {
     auto const [size_i, size_j, size_k] = boundary.block_dims(*block);
-    auto const offset = lattice.get_local_grid_range().first;
-    auto const off_i = offset[0];
-    auto const off_j = offset[1];
-    auto const off_k = offset[2];
     // Get field data which knows about the indices
     // In the loop, x,y,z are in block-local coordinates
-    auto const n_ghost_layers = lattice.get_ghost_layers();
-    auto const ghosts = static_cast<int>(n_ghost_layers);
-    for (int i = off_i - ghosts; i < size_i + off_i + ghosts; ++i) {
-      for (int j = off_j - ghosts; j < size_j + off_j + ghosts; ++j) {
-        for (int k = off_k - ghosts; k < size_k + off_k + ghosts; ++k) {
+    for (int i = off_i - gl; i < size_i + off_i + gl; ++i) {
+      for (int j = off_j - gl; j < size_j + off_j + gl; ++j) {
+        for (int k = off_k - gl; k < size_k + off_k + gl; ++k) {
           auto const node = Utils::Vector3i{{i, j, k}};
           auto const idx = (node + grid_size) % grid_size;
-          if (raster(idx)) {
+          auto const index = idx[0] * n_y * n_z + idx[1] * n_z + idx[2];
+          if (raster_flat[index]) {
             auto const bc = get_block_and_cell(lattice, node, true);
-            boundary.set_node_value_at_boundary(node, data_grid(idx), *bc);
+            boundary.set_node_value_at_boundary(node, data_flat[index], *bc);
           }
         }
       }
diff --git a/src/walberla_bridge/src/LatticeWalberla.cpp b/src/walberla_bridge/src/LatticeWalberla.cpp
index e48bd17f6e..b3d9059d5d 100644
--- a/src/walberla_bridge/src/LatticeWalberla.cpp
+++ b/src/walberla_bridge/src/LatticeWalberla.cpp
@@ -29,11 +29,10 @@
 
 #include <utils/Vector.hpp>
 
-#include <boost/optional.hpp>
-
 #include <cmath>
 #include <initializer_list>
 #include <limits>
+#include <optional>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -80,11 +79,11 @@ LatticeWalberla::get_local_domain() const {
 [[nodiscard]] bool
 LatticeWalberla::node_in_local_domain(Utils::Vector3i const &node) const {
   // Note: Lattice constant =1, cell centers offset by .5
-  return ::walberla::get_block_and_cell(*this, node, false) != boost::none;
+  return ::walberla::get_block_and_cell(*this, node, false) != std::nullopt;
 }
 [[nodiscard]] bool
 LatticeWalberla::node_in_local_halo(Utils::Vector3i const &node) const {
-  return ::walberla::get_block_and_cell(*this, node, true) != boost::none;
+  return ::walberla::get_block_and_cell(*this, node, true) != std::nullopt;
 }
 [[nodiscard]] bool
 LatticeWalberla::pos_in_local_domain(Utils::Vector3d const &pos) const {
diff --git a/src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp b/src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp
index 5d6c1cb5d7..b60bfa935a 100644
--- a/src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp
+++ b/src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp
@@ -43,6 +43,7 @@
 #include <cstddef>
 #include <iterator>
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <string>
 #include <type_traits>
@@ -120,7 +121,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
   // TODO: kernel for that
   // std::shared_ptr<ResetForce<PdfField, VectorField>> m_reset_force;
 
-  [[nodiscard]] boost::optional<CellInterval>
+  [[nodiscard]] std::optional<CellInterval>
   get_interval(Utils::Vector3i const &lower_corner,
                Utils::Vector3i const &upper_corner) const {
     auto const &lattice = get_lattice();
@@ -129,7 +130,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     auto const lower_bc = get_block_and_cell(lattice, cell_min, true);
     auto const upper_bc = get_block_and_cell(lattice, cell_max, true);
     if (not lower_bc or not upper_bc) {
-      return {};
+      return std::nullopt;
     }
     assert(&(*(lower_bc->block)) == &(*(upper_bc->block)));
     return {CellInterval(lower_bc->cell, upper_bc->cell)};
@@ -378,13 +379,13 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     return true;
   }
 
-  [[nodiscard]] boost::optional<double>
+  [[nodiscard]] std::optional<double>
   get_node_density(Utils::Vector3i const &node,
                    bool consider_ghosts = false) const override {
     auto bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
 
     if (!bc)
-      return {boost::none};
+      return std::nullopt;
 
     auto const density_field =
         bc->block->template getData<DensityField>(m_density_field_id);
@@ -457,12 +458,12 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     return true;
   }
 
-  [[nodiscard]] boost::optional<Utils::Vector3d>
+  [[nodiscard]] std::optional<Utils::Vector3d>
   get_node_flux_at_boundary(Utils::Vector3i const &node,
                             bool consider_ghosts = false) const override {
     auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc or !m_boundary_flux->node_is_boundary(node))
-      return {boost::none};
+      return std::nullopt;
 
     return {m_boundary_flux->get_node_value_at_boundary(node)};
   }
@@ -489,19 +490,19 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     return true;
   }
 
-  [[nodiscard]] boost::optional<double>
+  [[nodiscard]] std::optional<double>
   get_node_density_at_boundary(Utils::Vector3i const &node,
                                bool consider_ghosts = false) const override {
     auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc or !m_boundary_density->node_is_boundary(node))
-      return {boost::none};
+      return std::nullopt;
 
     return {double_c(m_boundary_density->get_node_value_at_boundary(node))};
   }
 
   void set_slice_density_boundary(
       Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner,
-      std::vector<boost::optional<double>> const &density) override {
+      std::vector<std::optional<double>> const &density) override {
     if (auto const ci = get_interval(lower_corner, upper_corner)) {
       auto const &lattice = get_lattice();
       auto const local_offset = std::get<0>(lattice.get_local_grid_range());
@@ -527,11 +528,11 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     }
   }
 
-  [[nodiscard]] std::vector<boost::optional<double>>
+  [[nodiscard]] std::vector<std::optional<double>>
   get_slice_density_at_boundary(
       Utils::Vector3i const &lower_corner,
       Utils::Vector3i const &upper_corner) const override {
-    std::vector<boost::optional<double>> out;
+    std::vector<std::optional<double>> out;
     if (auto const ci = get_interval(lower_corner, upper_corner)) {
       auto const &lattice = get_lattice();
       auto const local_offset = std::get<0>(lattice.get_local_grid_range());
@@ -547,7 +548,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
               out.emplace_back(
                   m_boundary_density->get_node_value_at_boundary(node));
             } else {
-              out.emplace_back(boost::none);
+              out.emplace_back(std::nullopt);
             }
           }
         }
@@ -559,7 +560,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
 
   void set_slice_flux_boundary(
       Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner,
-      std::vector<boost::optional<Utils::Vector3d>> const &flux) override {
+      std::vector<std::optional<Utils::Vector3d>> const &flux) override {
     if (auto const ci = get_interval(lower_corner, upper_corner)) {
       auto const &lattice = get_lattice();
       auto const local_offset = std::get<0>(lattice.get_local_grid_range());
@@ -585,11 +586,11 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     }
   }
 
-  [[nodiscard]] std::vector<boost::optional<Utils::Vector3d>>
+  [[nodiscard]] std::vector<std::optional<Utils::Vector3d>>
   get_slice_flux_at_boundary(
       Utils::Vector3i const &lower_corner,
       Utils::Vector3i const &upper_corner) const override {
-    std::vector<boost::optional<Utils::Vector3d>> out;
+    std::vector<std::optional<Utils::Vector3d>> out;
     if (auto const ci = get_interval(lower_corner, upper_corner)) {
       auto const &lattice = get_lattice();
       auto const local_offset = std::get<0>(lattice.get_local_grid_range());
@@ -605,7 +606,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
               out.emplace_back(
                   m_boundary_flux->get_node_value_at_boundary(node));
             } else {
-              out.emplace_back(boost::none);
+              out.emplace_back(std::nullopt);
             }
           }
         }
@@ -650,32 +651,32 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     return true;
   }
 
-  [[nodiscard]] boost::optional<bool>
+  [[nodiscard]] std::optional<bool>
   get_node_is_flux_boundary(Utils::Vector3i const &node,
                             bool consider_ghosts) const override {
     auto bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc)
-      return {boost::none};
+      return std::nullopt;
 
     return {m_boundary_flux->node_is_boundary(node)};
   }
 
-  [[nodiscard]] boost::optional<bool>
+  [[nodiscard]] std::optional<bool>
   get_node_is_density_boundary(Utils::Vector3i const &node,
                                bool consider_ghosts) const override {
     auto bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc)
-      return {boost::none};
+      return std::nullopt;
 
     return {m_boundary_density->node_is_boundary(node)};
   }
 
-  [[nodiscard]] boost::optional<bool>
+  [[nodiscard]] std::optional<bool>
   get_node_is_boundary(Utils::Vector3i const &node,
                        bool consider_ghosts = false) const override {
     auto bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc)
-      return {boost::none};
+      return std::nullopt;
 
     return {m_boundary_density->node_is_boundary(node) or
             m_boundary_flux->node_is_boundary(node)};
diff --git a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp
index 23b9d21991..90c14b9069 100644
--- a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp
+++ b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp
@@ -30,10 +30,9 @@
 #include <domain_decomposition/IBlock.h>
 #include <field/AddToStorage.h>
 
-#include <boost/optional.hpp>
-
 #include <cstddef>
 #include <memory>
+#include <optional>
 #include <tuple>
 #include <vector>
 
@@ -183,11 +182,11 @@ void EKReactionImplIndexed::set_node_is_boundary(Utils::Vector3i const &node,
   m_pending_changes = true;
 }
 
-boost::optional<bool>
+std::optional<bool>
 EKReactionImplIndexed::get_node_is_boundary(Utils::Vector3i const &node) {
   auto bc = get_block_and_cell(*get_lattice(), node, true);
   if (!bc)
-    return {boost::none};
+    return std::nullopt;
 
   auto [flag_field, boundary_flag] =
       detail::get_flag_field_and_flag<detail::FlagField>(
diff --git a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp
index b62f85d301..65686644f7 100644
--- a/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp
+++ b/src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.hpp
@@ -25,10 +25,9 @@
 
 #include <utils/Vector.hpp>
 
-#include <boost/optional.hpp>
-
 #include <cstddef>
 #include <memory>
+#include <optional>
 #include <vector>
 
 namespace walberla {
@@ -53,7 +52,7 @@ class EKReactionImplIndexed : public EKReactionBase {
   void perform_reaction() override;
 
   void set_node_is_boundary(Utils::Vector3i const &node, bool is_boundary);
-  [[nodiscard]] boost::optional<bool>
+  [[nodiscard]] std::optional<bool>
   get_node_is_boundary(Utils::Vector3i const &node);
 
   [[nodiscard]] auto get_indexvector_id() const noexcept {
diff --git a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp
index 4b6d0c8122..f4a6e2560e 100644
--- a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp
+++ b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp
@@ -61,19 +61,18 @@
 #include <utils/interpolation/bspline_3d.hpp>
 #include <utils/math/make_lin_space.hpp>
 
-#include <boost/optional.hpp>
-#include <boost/variant.hpp>
-
 #include <array>
 #include <cmath>
 #include <cstddef>
 #include <initializer_list>
 #include <limits>
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <string>
 #include <type_traits>
 #include <utility>
+#include <variant>
 #include <vector>
 
 namespace walberla {
@@ -97,7 +96,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
                        typename detail::BoundaryHandlingTrait<
                            FloatType, Architecture>::Dynamic_UBB>;
   using CollisionModel =
-      boost::variant<CollisionModelThermalized, CollisionModelLeesEdwards>;
+      std::variant<CollisionModelThermalized, CollisionModelLeesEdwards>;
 
 public:
   // Type definitions
@@ -131,7 +130,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
   }
 
 private:
-  class : public boost::static_visitor<> {
+  class CollideSweepVisitor {
   public:
     void operator()(CollisionModelThermalized &cm, IBlock *b) { cm(b); }
 
@@ -140,15 +139,16 @@ class LBWalberlaImpl : public LBWalberlaBase {
           m_lees_edwards_callbacks->get_shear_velocity());
       cm(b);
     }
-    void register_lees_edwards_callbacks(
-        std::shared_ptr<LeesEdwardsPack> lees_edwards_callbacks) {
-      m_lees_edwards_callbacks = std::move(lees_edwards_callbacks);
+
+    CollideSweepVisitor() = default;
+    explicit CollideSweepVisitor(std::shared_ptr<LeesEdwardsPack> callbacks) {
+      m_lees_edwards_callbacks = std::move(callbacks);
     }
 
   private:
-    std::shared_ptr<LeesEdwardsPack> m_lees_edwards_callbacks;
-
-  } run_collide_sweep;
+    std::shared_ptr<LeesEdwardsPack> m_lees_edwards_callbacks{};
+  };
+  CollideSweepVisitor m_run_collide_sweep{};
 
   FloatType shear_mode_relaxation_rate() const {
     return FloatType{2} / (FloatType{6} * m_viscosity + FloatType{1});
@@ -262,7 +262,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
   // lattice
   std::shared_ptr<LatticeWalberla> m_lattice;
 
-  [[nodiscard]] boost::optional<CellInterval>
+  [[nodiscard]] std::optional<CellInterval>
   get_interval(Utils::Vector3i const &lower_corner,
                Utils::Vector3i const &upper_corner) const {
     auto const &lattice = get_lattice();
@@ -271,7 +271,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
     auto const lower_bc = get_block_and_cell(lattice, cell_min, true);
     auto const upper_bc = get_block_and_cell(lattice, cell_max, true);
     if (not lower_bc or not upper_bc) {
-      return {};
+      return std::nullopt;
     }
     assert(&(*(lower_bc->block)) == &(*(upper_bc->block)));
     return {CellInterval(lower_bc->cell, upper_bc->cell)};
@@ -392,17 +392,17 @@ class LBWalberlaImpl : public LBWalberlaBase {
   }
 
   void integrate_collide(std::shared_ptr<Lattice_T> const &blocks) {
+    auto &cm_variant = *m_collision_model;
     for (auto b = blocks->begin(); b != blocks->end(); ++b)
-      boost::apply_visitor(run_collide_sweep, *m_collision_model,
-                           boost::variant<IBlock *>(&*b));
-    if (auto *cm = boost::get<CollisionModelThermalized>(&*m_collision_model)) {
+      std::visit(m_run_collide_sweep, cm_variant, std::variant<IBlock *>(&*b));
+    if (auto *cm = std::get_if<CollisionModelThermalized>(&cm_variant)) {
       cm->time_step_++;
     }
   }
 
   auto has_lees_edwards_bc() const {
-    return boost::get<CollisionModelLeesEdwards>(&*m_collision_model) !=
-           nullptr;
+    return std::holds_alternative<CollisionModelLeesEdwards>(
+        *m_collision_model);
   }
 
   void apply_lees_edwards_pdf_interpolation(
@@ -534,7 +534,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
         m_last_applied_force_field_id, m_pdf_field_id, agrid, omega, shear_vel);
     m_collision_model = std::make_shared<CollisionModel>(std::move(obj));
     m_lees_edwards_callbacks = std::move(lees_edwards_pack);
-    run_collide_sweep.register_lees_edwards_callbacks(m_lees_edwards_callbacks);
+    m_run_collide_sweep = CollideSweepVisitor{m_lees_edwards_callbacks};
     m_lees_edwards_pdf_interpol_sweep =
         std::make_shared<InterpolateAndShiftAtBoundary<PdfField, FloatType>>(
             blocks, m_pdf_field_id, m_pdf_tmp_field_id, n_ghost_layers,
@@ -577,7 +577,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
   }
 
   // Velocity
-  boost::optional<Utils::Vector3d>
+  std::optional<Utils::Vector3d>
   get_node_velocity(Utils::Vector3i const &node,
                     bool consider_ghosts = false) const override {
     auto const is_boundary = get_node_is_boundary(node, consider_ghosts);
@@ -586,7 +586,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
         return get_node_velocity_at_boundary(node, consider_ghosts);
     auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc)
-      return {};
+      return std::nullopt;
 
     auto field = bc->block->template uncheckedFastGetData<VectorField>(
         m_velocity_field_id);
@@ -680,13 +680,13 @@ class LBWalberlaImpl : public LBWalberlaBase {
     }
   }
 
-  boost::optional<Utils::Vector3d>
+  std::optional<Utils::Vector3d>
   get_velocity_at_pos(Utils::Vector3d const &pos,
                       bool consider_points_in_halo = false) const override {
     if (!consider_points_in_halo and !m_lattice->pos_in_local_domain(pos))
-      return {};
+      return std::nullopt;
     if (consider_points_in_halo and !m_lattice->pos_in_local_halo(pos))
-      return {};
+      return std::nullopt;
     Utils::Vector3d v{0., 0., 0.};
     interpolate_bspline_at_pos(
         pos, [this, &v, &pos](std::array<int, 3> const node, double weight) {
@@ -703,13 +703,13 @@ class LBWalberlaImpl : public LBWalberlaBase {
     return {std::move(v)};
   }
 
-  boost::optional<double>
+  std::optional<double>
   get_density_at_pos(Utils::Vector3d const &pos,
                      bool consider_points_in_halo = false) const override {
     if (!consider_points_in_halo and !m_lattice->pos_in_local_domain(pos))
-      return {};
+      return std::nullopt;
     if (consider_points_in_halo and !m_lattice->pos_in_local_halo(pos))
-      return {};
+      return std::nullopt;
     double dens = 0.;
     interpolate_bspline_at_pos(
         pos, [this, &dens, &pos](std::array<int, 3> const node, double weight) {
@@ -747,11 +747,11 @@ class LBWalberlaImpl : public LBWalberlaBase {
     return true;
   }
 
-  boost::optional<Utils::Vector3d>
+  std::optional<Utils::Vector3d>
   get_node_force_to_be_applied(Utils::Vector3i const &node) const override {
     auto const bc = get_block_and_cell(get_lattice(), node, true);
     if (!bc)
-      return {};
+      return std::nullopt;
 
     auto field =
         bc->block->template getData<VectorField>(m_force_to_be_applied_id);
@@ -759,12 +759,12 @@ class LBWalberlaImpl : public LBWalberlaBase {
     return to_vector3d(vec);
   }
 
-  boost::optional<Utils::Vector3d>
+  std::optional<Utils::Vector3d>
   get_node_last_applied_force(Utils::Vector3i const &node,
                               bool consider_ghosts = false) const override {
     auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc)
-      return {};
+      return std::nullopt;
 
     auto const field =
         bc->block->template getData<VectorField>(m_last_applied_force_field_id);
@@ -842,12 +842,12 @@ class LBWalberlaImpl : public LBWalberlaBase {
   }
 
   // Population
-  boost::optional<std::vector<double>>
+  std::optional<std::vector<double>>
   get_node_population(Utils::Vector3i const &node,
                       bool consider_ghosts = false) const override {
     auto bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc)
-      return {boost::none};
+      return std::nullopt;
 
     auto pdf_field = bc->block->template getData<PdfField>(m_pdf_field_id);
     auto const pop = lbm::accessor::Population::get(pdf_field, bc->cell);
@@ -909,12 +909,12 @@ class LBWalberlaImpl : public LBWalberlaBase {
   }
 
   // Density
-  boost::optional<double>
+  std::optional<double>
   get_node_density(Utils::Vector3i const &node,
                    bool consider_ghosts = false) const override {
     auto bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc)
-      return {boost::none};
+      return std::nullopt;
 
     auto pdf_field =
         bc->block->template uncheckedFastGetData<PdfField>(m_pdf_field_id);
@@ -966,12 +966,12 @@ class LBWalberlaImpl : public LBWalberlaBase {
     }
   }
 
-  boost::optional<Utils::Vector3d>
+  std::optional<Utils::Vector3d>
   get_node_velocity_at_boundary(Utils::Vector3i const &node,
                                 bool consider_ghosts = false) const override {
     auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc or !m_boundary->node_is_boundary(node))
-      return {boost::none};
+      return std::nullopt;
 
     return {m_boundary->get_node_value_at_boundary(node)};
   }
@@ -987,10 +987,10 @@ class LBWalberlaImpl : public LBWalberlaBase {
     return true;
   }
 
-  std::vector<boost::optional<Utils::Vector3d>> get_slice_velocity_at_boundary(
+  std::vector<std::optional<Utils::Vector3d>> get_slice_velocity_at_boundary(
       Utils::Vector3i const &lower_corner,
       Utils::Vector3i const &upper_corner) const override {
-    std::vector<boost::optional<Utils::Vector3d>> out;
+    std::vector<std::optional<Utils::Vector3d>> out;
     if (auto const ci = get_interval(lower_corner, upper_corner)) {
       auto const &lattice = get_lattice();
       auto const local_offset = std::get<0>(lattice.get_local_grid_range());
@@ -1005,7 +1005,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
             if (m_boundary->node_is_boundary(node)) {
               out.emplace_back(m_boundary->get_node_value_at_boundary(node));
             } else {
-              out.emplace_back(boost::none);
+              out.emplace_back(std::nullopt);
             }
           }
         }
@@ -1017,7 +1017,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
 
   void set_slice_velocity_at_boundary(
       Utils::Vector3i const &lower_corner, Utils::Vector3i const &upper_corner,
-      std::vector<boost::optional<Utils::Vector3d>> const &velocity) override {
+      std::vector<std::optional<Utils::Vector3d>> const &velocity) override {
     if (auto const ci = get_interval(lower_corner, upper_corner)) {
       auto const &lattice = get_lattice();
       auto const local_offset = std::get<0>(lattice.get_local_grid_range());
@@ -1043,11 +1043,11 @@ class LBWalberlaImpl : public LBWalberlaBase {
     }
   }
 
-  boost::optional<Utils::Vector3d>
+  std::optional<Utils::Vector3d>
   get_node_boundary_force(Utils::Vector3i const &node) const override {
     auto const bc = get_block_and_cell(get_lattice(), node, true);
     if (!bc or !m_boundary->node_is_boundary(node))
-      return {boost::none};
+      return std::nullopt;
 
     return get_node_last_applied_force(node, true);
   }
@@ -1062,12 +1062,12 @@ class LBWalberlaImpl : public LBWalberlaBase {
     return true;
   }
 
-  boost::optional<bool>
+  std::optional<bool>
   get_node_is_boundary(Utils::Vector3i const &node,
                        bool consider_ghosts = false) const override {
     auto const bc = get_block_and_cell(get_lattice(), node, consider_ghosts);
     if (!bc)
-      return {boost::none};
+      return std::nullopt;
 
     return {m_boundary->node_is_boundary(node)};
   }
@@ -1109,11 +1109,11 @@ class LBWalberlaImpl : public LBWalberlaBase {
   }
 
   // Pressure tensor
-  boost::optional<Utils::VectorXd<9>>
+  std::optional<Utils::VectorXd<9>>
   get_node_pressure_tensor(Utils::Vector3i const &node) const override {
     auto bc = get_block_and_cell(get_lattice(), node, false);
     if (!bc)
-      return {boost::none};
+      return std::nullopt;
 
     auto pdf_field = bc->block->template getData<PdfField>(m_pdf_field_id);
     auto tensor = lbm::accessor::PressureTensor::get(pdf_field, bc->cell);
@@ -1193,16 +1193,16 @@ class LBWalberlaImpl : public LBWalberlaBase {
     return numeric_cast<double>(m_kT);
   }
 
-  [[nodiscard]] boost::optional<uint64_t> get_rng_state() const override {
-    auto const *cm = boost::get<CollisionModelThermalized>(&*m_collision_model);
+  [[nodiscard]] std::optional<uint64_t> get_rng_state() const override {
+    auto const cm = std::get_if<CollisionModelThermalized>(&*m_collision_model);
     if (!cm or m_kT == 0.) {
-      return {boost::none};
+      return std::nullopt;
     }
     return {static_cast<uint64_t>(cm->time_step_)};
   }
 
   void set_rng_state(uint64_t counter) override {
-    auto *cm = boost::get<CollisionModelThermalized>(&*m_collision_model);
+    auto const cm = std::get_if<CollisionModelThermalized>(&*m_collision_model);
     if (!cm or m_kT == 0.) {
       throw std::runtime_error("This LB instance is unthermalized");
     }
diff --git a/src/walberla_bridge/tests/CMakeLists.txt b/src/walberla_bridge/tests/CMakeLists.txt
index 5abe24f149..581a123f44 100644
--- a/src/walberla_bridge/tests/CMakeLists.txt
+++ b/src/walberla_bridge/tests/CMakeLists.txt
@@ -22,53 +22,43 @@ include(unit_test)
 function(ESPRESSO_WALBERLA_UNIT_TEST)
   cmake_parse_arguments(TEST "" "NAME;NUM_PROC" "SRC;DEPENDS" ${ARGN})
   unit_test(NAME ${TEST_NAME} NUM_PROC ${TEST_NUM_PROC} SRC ${TEST_SRC} DEPENDS
-            ${TEST_DEPENDS})
+            ${TEST_DEPENDS} espresso::walberla espresso::walberla::cpp_flags
+            espresso::utils)
   set_target_properties(${TEST_NAME} PROPERTIES CXX_CLANG_TIDY "")
   target_include_directories(${TEST_NAME} PRIVATE ${WALBERLA_INCLUDE_DIRS}
                                                   ${walberla_BINARY_DIR}/src)
   target_link_libraries(${TEST_NAME} PRIVATE ${WALBERLA_LIBS})
 endfunction()
 
-espresso_walberla_unit_test(
-  NAME ResourceManager_test SRC ResourceManager_test.cpp DEPENDS
-  espresso::walberla espresso::walberla::cpp_flags)
+espresso_walberla_unit_test(NAME ResourceManager_test SRC
+                            ResourceManager_test.cpp)
 
-espresso_walberla_unit_test(
-  NAME kernels_unit_tests SRC kernels_unit_tests.cpp DEPENDS espresso::walberla
-  espresso::walberla::cpp_flags espresso::utils)
+espresso_walberla_unit_test(NAME kernels_unit_tests SRC kernels_unit_tests.cpp)
 
 espresso_walberla_unit_test(
   NAME LatticeWalberla_unit_tests SRC LatticeWalberla_unit_tests.cpp DEPENDS
-  espresso::walberla espresso::walberla::cpp_flags espresso::utils Boost::mpi
-  NUM_PROC 2)
+  Boost::boost Boost::mpi NUM_PROC 2)
 
 espresso_walberla_unit_test(
   NAME LBWalberlaImpl_unit_tests SRC LBWalberlaImpl_unit_tests.cpp DEPENDS
-  espresso::walberla espresso::walberla::cpp_flags espresso::utils Boost::mpi
-  NUM_PROC 2)
+  Boost::boost Boost::mpi NUM_PROC 2)
 
 espresso_walberla_unit_test(
   NAME LBWalberlaImpl_bspline_tests SRC LBWalberlaImpl_bspline_tests.cpp
-  DEPENDS espresso::walberla espresso::walberla::cpp_flags espresso::utils
-  Boost::mpi NUM_PROC 2)
+  DEPENDS Boost::mpi NUM_PROC 2)
 
 if(NOT (ESPRESSO_BUILD_WITH_ASAN OR ESPRESSO_BUILD_WITH_UBSAN))
   espresso_walberla_unit_test(
     NAME LBWalberlaImpl_statistical_tests SRC
-    LBWalberlaImpl_statistical_tests.cpp DEPENDS espresso::walberla
-    espresso::walberla::cpp_flags espresso::utils Boost::mpi)
+    LBWalberlaImpl_statistical_tests.cpp DEPENDS Boost::mpi)
 endif()
 
-espresso_walberla_unit_test(
-  NAME LBWalberlaImpl_flow_tests SRC LBWalberlaImpl_flow_tests.cpp DEPENDS
-  espresso::walberla espresso::walberla::cpp_flags espresso::utils Boost::mpi)
+espresso_walberla_unit_test(NAME LBWalberlaImpl_flow_tests SRC
+                            LBWalberlaImpl_flow_tests.cpp DEPENDS Boost::mpi)
 
-espresso_walberla_unit_test(
-  NAME LBWalberlaImpl_lees_edwards_test SRC LBWalberlaImpl_lees_edwards.cpp
-  DEPENDS espresso::walberla espresso::walberla::cpp_flags espresso::utils
-  Boost::mpi)
+espresso_walberla_unit_test(NAME LBWalberlaImpl_lees_edwards_test SRC
+                            LBWalberlaImpl_lees_edwards.cpp DEPENDS Boost::mpi)
 
 espresso_walberla_unit_test(
   NAME EKinWalberlaImpl_unit_tests SRC EKinWalberlaImpl_unit_tests.cpp DEPENDS
-  espresso::walberla espresso::walberla::cpp_flags espresso::utils Boost::mpi
-  NUM_PROC 2)
+  Boost::boost Boost::mpi NUM_PROC 2)