diff --git a/include/bout/bout_enum_class.hxx b/include/bout/bout_enum_class.hxx
index ce088a4628..9d4249213d 100644
--- a/include/bout/bout_enum_class.hxx
+++ b/include/bout/bout_enum_class.hxx
@@ -64,12 +64,12 @@
 
 #define BOUT_MAKE_FROMSTRING_NAME(enumname) enumname##FromString
 
-/// Create an enum class with toString and <enum name>FromString functions, and an
+/// Create an enum class with format_as and <enum name>FromString functions, and an
 /// Options::as<enum> overload to read the enum
 #define BOUT_ENUM_CLASS(enumname, ...)                                         \
   enum class enumname : std::uint8_t { __VA_ARGS__ };                          \
                                                                                \
-  inline std::string toString(enumname e) {                                    \
+  inline std::string format_as(enumname e) {                                   \
                                                                                \
     const static std::map<enumname, std::string> toString_map = {              \
         BOUT_ENUM_CLASS_MAP_ARGS(BOUT_ENUM_CLASS_STR, enumname, __VA_ARGS__)}; \
@@ -102,51 +102,49 @@
   }                                                                            \
                                                                                \
   inline std::ostream& operator<<(std::ostream& out, const enumname& e) {      \
-    return out << toString(e);                                                 \
+    return out << fmt::to_string(e);                                           \
   }
 
-/// Create an enum class with toString and <enum name>FromString functions, and an
+/// Create an enum class with format_as and <enum name>FromString functions, and an
 /// Options::as<enum> overload to read the enum
-#define BOUT_ENUM_CLASS_NS(ns, enumname, ...)                                      \
-  namespace ns {                                                                   \
-  enum class enumname : std::uint8_t { __VA_ARGS__ };                              \
-  }                                                                                \
-                                                                                   \
-  inline std::string toString(ns::enumname e) {                                    \
-                                                                                   \
-    const static std::map<ns::enumname, std::string> toString_map = {              \
-        BOUT_ENUM_CLASS_MAP_ARGS(BOUT_ENUM_CLASS_STR, ns::enumname, __VA_ARGS__)}; \
-    auto found = toString_map.find(e);                                             \
-    if (found == toString_map.end()) {                                             \
-      throw BoutException("Did not find enum {:d}", static_cast<int>(e));          \
-    }                                                                              \
-    return found->second;                                                          \
-  }                                                                                \
-  namespace ns {                                                                   \
-  inline enumname BOUT_MAKE_FROMSTRING_NAME(enumname)(const std::string& s) {      \
-    const static std::map<std::string, enumname> fromString_map = {                \
-        BOUT_ENUM_CLASS_MAP_ARGS(BOUT_STR_ENUM_CLASS, enumname, __VA_ARGS__)};     \
-    auto found = fromString_map.find(s);                                           \
-    if (found == fromString_map.end()) {                                           \
-      std::string valid_values{};                                                  \
-      for (auto const& entry : fromString_map) {                                   \
-        valid_values += std::string(" ") + entry.first;                            \
-      }                                                                            \
-      throw BoutException("Did not find enum {:s}. Valid values: {:s}", s,         \
-                          valid_values);                                           \
-    }                                                                              \
-    return found->second;                                                          \
-  }                                                                                \
-  }                                                                                \
-                                                                                   \
-  template <>                                                                      \
-  inline ns::enumname Options::as<ns::enumname>(const ns::enumname&) const {       \
-    return ns::BOUT_MAKE_FROMSTRING_NAME(enumname)(this->as<std::string>());       \
-  }                                                                                \
-  namespace ns {                                                                   \
-  inline std::ostream& operator<<(std::ostream& out, const enumname& e) {          \
-    return out << toString(e);                                                     \
-  }                                                                                \
+#define BOUT_ENUM_CLASS_NS(ns, enumname, ...)                                  \
+  namespace ns {                                                               \
+  enum class enumname : std::uint8_t { __VA_ARGS__ };                          \
+                                                                               \
+  inline std::string format_as(const enumname& e) {                            \
+                                                                               \
+    const static std::map<enumname, std::string> toString_map = {              \
+        BOUT_ENUM_CLASS_MAP_ARGS(BOUT_ENUM_CLASS_STR, enumname, __VA_ARGS__)}; \
+    auto found = toString_map.find(e);                                         \
+    if (found == toString_map.end()) {                                         \
+      throw BoutException("Did not find enum {:d}", static_cast<int>(e));      \
+    }                                                                          \
+    return found->second;                                                      \
+  }                                                                            \
+                                                                               \
+  inline enumname BOUT_MAKE_FROMSTRING_NAME(enumname)(const std::string& s) {  \
+    const static std::map<std::string, enumname> fromString_map = {            \
+        BOUT_ENUM_CLASS_MAP_ARGS(BOUT_STR_ENUM_CLASS, enumname, __VA_ARGS__)}; \
+    auto found = fromString_map.find(s);                                       \
+    if (found == fromString_map.end()) {                                       \
+      std::string valid_values{};                                              \
+      for (auto const& entry : fromString_map) {                               \
+        valid_values += std::string(" ") + entry.first;                        \
+      }                                                                        \
+      throw BoutException("Did not find enum {:s}. Valid values: {:s}", s,     \
+                          valid_values);                                       \
+    }                                                                          \
+    return found->second;                                                      \
+  }                                                                            \
+                                                                               \
+  inline std::ostream& operator<<(std::ostream& out, const enumname& e) {      \
+    return out << fmt::to_string(e);                                           \
+  }                                                                            \
+  } /* end namespace ns */                                                     \
+                                                                               \
+  template <>                                                                  \
+  inline ns::enumname Options::as<ns::enumname>(const ns::enumname&) const {   \
+    return ns::BOUT_MAKE_FROMSTRING_NAME(enumname)(this->as<std::string>());   \
   }
 
 #endif // BOUT_ENUM_CLASS_H
diff --git a/include/bout/bout_types.hxx b/include/bout/bout_types.hxx
index 41ca37c99d..4cb51cbe37 100644
--- a/include/bout/bout_types.hxx
+++ b/include/bout/bout_types.hxx
@@ -49,9 +49,8 @@ constexpr CELL_LOC CELL_YLOW = CELL_LOC::ylow;
 constexpr CELL_LOC CELL_ZLOW = CELL_LOC::zlow;
 constexpr CELL_LOC CELL_VSHIFT = CELL_LOC::vshift;
 
-std::string toString(CELL_LOC location);
+std::string format_as(const CELL_LOC& location);
 CELL_LOC CELL_LOCFromString(const std::string& location_string);
-inline std::string format_as(const CELL_LOC& loc) { return toString(loc); }
 
 /// Differential methods. Both central and upwind
 enum class DIFF_METHOD { deflt, u1, u2, c2, w2, w3, c4, u3, fft, split, s2 };
@@ -67,8 +66,7 @@ constexpr DIFF_METHOD DIFF_FFT = DIFF_METHOD::fft;
 constexpr DIFF_METHOD DIFF_SPLIT = DIFF_METHOD::split;
 constexpr DIFF_METHOD DIFF_S2 = DIFF_METHOD::s2;
 
-std::string toString(DIFF_METHOD location);
-inline std::string format_as(const DIFF_METHOD& method) { return toString(method); }
+std::string format_as(const DIFF_METHOD& method);
 
 /// Specify grid region for looping
 enum class REGION { all, nobndry, nox, noy, noz };
@@ -78,8 +76,7 @@ constexpr REGION RGN_NOX = REGION::nox;
 constexpr REGION RGN_NOY = REGION::noy;
 constexpr REGION RGN_NOZ = REGION::noz;
 
-std::string toString(REGION region);
-inline std::string format_as(const REGION& region) { return toString(region); }
+std::string format_as(const REGION& region);
 
 /// To identify particular directions (in index space):
 ///   - X, Y, Z are the coordinate directions
@@ -90,8 +87,7 @@ inline std::string format_as(const REGION& region) { return toString(region); }
 ///     field-aligned
 enum class DIRECTION { X, Y, Z, YAligned, YOrthogonal };
 
-std::string toString(DIRECTION direction);
-inline std::string format_as(const DIRECTION& dir) { return toString(dir); }
+std::string format_as(const DIRECTION& dir);
 
 /// Identify kind of a field's y-direction
 /// - Standard is the default for the Mesh/Coordinates/ParallelTransform
@@ -99,9 +95,8 @@ inline std::string format_as(const DIRECTION& dir) { return toString(dir); }
 ///   coordinates
 enum class YDirectionType { Standard, Aligned };
 
-std::string toString(YDirectionType d);
+std::string format_as(const YDirectionType& d);
 YDirectionType YDirectionTypeFromString(const std::string& y_direction_string);
-inline std::string format_as(const YDirectionType& d) { return toString(d); }
 
 /// Identify kind of a field's z-direction
 /// - Standard is the default
@@ -110,9 +105,8 @@ inline std::string format_as(const YDirectionType& d) { return toString(d); }
 ///   Field2D)
 enum class ZDirectionType { Standard, Average };
 
-std::string toString(ZDirectionType d);
+std::string format_as(const ZDirectionType& d);
 ZDirectionType ZDirectionTypeFromString(const std::string& z_direction_string);
-inline std::string format_as(const ZDirectionType& d) { return toString(d); }
 
 /// Container for direction types
 struct DirectionTypes {
@@ -120,8 +114,7 @@ struct DirectionTypes {
   ZDirectionType z;
 };
 
-std::string toString(const DirectionTypes& dir);
-inline std::string format_as(const DirectionTypes& dir) { return toString(dir); }
+std::string format_as(const DirectionTypes& dir);
 
 /// Check whether direction types are compatible, so two fields with attributes
 /// d1 and d2 respectively can be added, subtracted, etc.
@@ -132,14 +125,12 @@ void swap(const DirectionTypes& first, const DirectionTypes& second);
 /// To identify valid staggering combinations
 enum class STAGGER { None, C2L, L2C };
 
-std::string toString(STAGGER stagger);
-inline std::string format_as(const STAGGER& stagger) { return toString(stagger); }
+std::string format_as(const STAGGER& stagger);
 
 /// To identify types of derivative method combinations
 enum class DERIV { Standard, StandardSecond, StandardFourth, Upwind, Flux };
 
-std::string toString(DERIV deriv);
-inline std::string format_as(const DERIV& deriv) { return toString(deriv); }
+std::string format_as(const DERIV& deriv);
 
 // A small struct that can be used to wrap a specific enum value, giving
 // it a unique type that can be passed as a valid type to templates and
diff --git a/include/bout/deriv_store.hxx b/include/bout/deriv_store.hxx
index 36412ec531..51d3de65f4 100644
--- a/include/bout/deriv_store.hxx
+++ b/include/bout/deriv_store.hxx
@@ -29,6 +29,7 @@
 #ifndef __DERIV_STORE_HXX__
 #define __DERIV_STORE_HXX__
 
+#include "fmt/format.h"
 #include <cstddef>
 #include <functional>
 #include <map>
@@ -87,7 +88,7 @@ struct DerivativeStore {
                STAGGER stagger = STAGGER::None) const {
 
     // Get the key
-    auto key = getKey(direction, stagger, toString(derivType));
+    auto key = getKey(direction, stagger, fmt::to_string(derivType));
     return isEmpty(key);
   }
 
@@ -97,7 +98,7 @@ struct DerivativeStore {
                                             STAGGER stagger = STAGGER::None) const {
 
     // Get the key
-    auto key = getKey(direction, stagger, toString(derivType));
+    auto key = getKey(direction, stagger, fmt::to_string(derivType));
     if (isEmpty(key)) {
       return std::set<std::string>{};
     } else {
@@ -112,11 +113,11 @@ struct DerivativeStore {
 
     // Introductory information
     output_info << "Available methods for derivative type '";
-    output_info << toString(derivType);
+    output_info << fmt::to_string(derivType);
     output_info << "' in direction ";
-    output_info << toString(direction);
+    output_info << fmt::to_string(direction);
     output_info << " ( Staggering : ";
-    output_info << toString(stagger) << " ) : \n";
+    output_info << fmt::to_string(stagger) << " ) : \n";
 
     for (const auto& i : getAvailableMethods(derivType, direction, stagger)) {
       output_info << "\t" << i << "\n";
@@ -162,7 +163,8 @@ struct DerivativeStore {
     };
 
     // Register this method name in lookup of known methods
-    registeredMethods[getKey(direction, stagger, toString(derivType))].insert(methodName);
+    registeredMethods[getKey(direction, stagger, fmt::to_string(derivType))].insert(
+        methodName);
   };
 
   /// Register a function with upwindFunc/fluxFunc interface. Which map is used
@@ -196,7 +198,8 @@ struct DerivativeStore {
     };
 
     // Register this method name in lookup of known methods
-    registeredMethods[getKey(direction, stagger, toString(derivType))].insert(methodName);
+    registeredMethods[getKey(direction, stagger, fmt::to_string(derivType))].insert(
+        methodName);
   };
 
   /// Templated versions of the above registration routines.
@@ -225,7 +228,7 @@ struct DerivativeStore {
                                      DERIV derivType = DERIV::Standard) const {
 
     const auto realName = nameLookup(
-        name, defaultMethods.at(getKey(direction, stagger, toString(derivType))));
+        name, defaultMethods.at(getKey(direction, stagger, fmt::to_string(derivType))));
     const auto key = getKey(direction, stagger, realName);
 
     const storageType<std::size_t, standardFunc>* theMap = nullptr;
@@ -269,7 +272,7 @@ struct DerivativeStore {
                              DERIV derivType = DERIV::Upwind) const {
 
     const auto realName = nameLookup(
-        name, defaultMethods.at(getKey(direction, stagger, toString(derivType))));
+        name, defaultMethods.at(getKey(direction, stagger, fmt::to_string(derivType))));
     const auto key = getKey(direction, stagger, realName);
 
     const storageType<std::size_t, flowFunc>* theMap = nullptr;
@@ -331,7 +334,7 @@ struct DerivativeStore {
 
         const auto theDirection = direction.first;
         const auto theDerivType = deriv.first;
-        const auto theDerivTypeString = toString(theDerivType);
+        const auto theDerivTypeString = fmt::to_string(theDerivType);
 
         // Note both staggered and unstaggered have the same fallback default currently
         std::string theDefault{
@@ -360,7 +363,7 @@ struct DerivativeStore {
         defaultMethods[getKey(theDirection, STAGGER::None, theDerivTypeString)] =
             theDefault;
         output_verbose << "The default method for derivative type " << theDerivTypeString
-                       << " in direction " << toString(theDirection) << " is "
+                       << " in direction " << fmt::to_string(theDirection) << " is "
                        << theDefault << "\n";
 
         //-------------------------------------------------------------
@@ -384,8 +387,8 @@ struct DerivativeStore {
         defaultMethods[getKey(theDirection, STAGGER::C2L, theDerivTypeString)] =
             theDefault;
         output_verbose << "The default method for staggered derivative type "
-                       << theDerivTypeString << " in direction " << toString(theDirection)
-                       << " is " << theDefault << "\n";
+                       << theDerivTypeString << " in direction "
+                       << fmt::to_string(theDirection) << " is " << theDefault << "\n";
       }
     }
   }
@@ -393,7 +396,7 @@ struct DerivativeStore {
   /// Provide a method to override/force a specific default method
   void forceDefaultMethod(std::string methodName, DERIV deriv, DIRECTION direction,
                           STAGGER stagger = STAGGER::None) {
-    const auto key = getKey(direction, stagger, toString(deriv));
+    const auto key = getKey(direction, stagger, fmt::to_string(deriv));
     defaultMethods[key] = uppercase(methodName);
   }
 
@@ -461,7 +464,7 @@ private:
     for (const auto& direction : directions) {
       for (const auto& deriv : derivTypes) {
         const auto theDirection = direction.first;
-        const auto theDerivTypeString = toString(deriv.first);
+        const auto theDerivTypeString = fmt::to_string(deriv.first);
         const std::string theDefault{uppercase(initialDefaultMethods[deriv.first])};
 
         //-------------------------------------------------------------
@@ -482,11 +485,11 @@ private:
   std::string getMethodName(std::string name, DIRECTION direction,
                             STAGGER stagger = STAGGER::None) const {
 
-    return name + " (" + toString(direction) + ", " + toString(stagger) + ")";
+    return name + " (" + fmt::to_string(direction) + ", " + fmt::to_string(stagger) + ")";
   };
 
   std::string nameLookup(const std::string name, const std::string defaultName) const {
-    return name != toString(DIFF_DEFAULT) ? name : defaultName;
+    return name != fmt::to_string(DIFF_DEFAULT) ? name : defaultName;
   }
 
   /// Provides a routine to produce a unique key given information
@@ -503,8 +506,8 @@ private:
     // maps to store the different field types as the signature is
     // different.
     std::size_t result;
-    result = std::hash<std::string>{}(toString(direction));
-    result = result ^ std::hash<std::string>{}(toString(stagger));
+    result = std::hash<std::string>{}(fmt::to_string(direction));
+    result = result ^ std::hash<std::string>{}(fmt::to_string(stagger));
     result = result ^ std::hash<std::string>{}(key);
     return result;
   }
diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx
index e20aef2459..dfd1587694 100644
--- a/include/bout/difops.hxx
+++ b/include/bout/difops.hxx
@@ -61,7 +61,7 @@ bout::FieldMetric Grad_par(const Field2D& var, CELL_LOC outloc = CELL_DEFAULT,
                            const std::string& method = "DEFAULT");
 inline bout::FieldMetric Grad_par(const Field2D& var, CELL_LOC outloc,
                                   DIFF_METHOD method) {
-  return Grad_par(var, outloc, toString(method));
+  return Grad_par(var, outloc, fmt::to_string(method));
 }
 
 Field3D Grad_par(const Field3DParallel& var, CELL_LOC outloc = CELL_DEFAULT,
@@ -99,7 +99,7 @@ bout::FieldMetric Vpar_Grad_par(const Field2D& v, const Field2D& f,
 
 inline bout::FieldMetric Vpar_Grad_par(const Field2D& v, const Field2D& f,
                                        CELL_LOC outloc, DIFF_METHOD method) {
-  return Vpar_Grad_par(v, f, outloc, toString(method));
+  return Vpar_Grad_par(v, f, outloc, fmt::to_string(method));
 }
 
 Field3D Vpar_Grad_par(const Field3DParallel& v, const Field3DParallel& f,
@@ -108,7 +108,7 @@ Field3D Vpar_Grad_par(const Field3DParallel& v, const Field3DParallel& f,
 
 inline Field3D Vpar_Grad_par(const Field3DParallel& v, const Field3DParallel& f,
                              CELL_LOC outloc, DIFF_METHOD method) {
-  return Vpar_Grad_par(v, f, outloc, toString(method));
+  return Vpar_Grad_par(v, f, outloc, fmt::to_string(method));
 }
 
 /*!
@@ -127,14 +127,14 @@ bout::FieldMetric Div_par(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
                           const std::string& method = "DEFAULT");
 
 inline bout::FieldMetric Div_par(const Field2D& f, CELL_LOC outloc, DIFF_METHOD method) {
-  return Div_par(f, outloc, toString(method));
+  return Div_par(f, outloc, fmt::to_string(method));
 }
 
 Field3D Div_par(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
                 const std::string& method = "DEFAULT");
 
 inline Field3D Div_par(const Field3DParallel& f, CELL_LOC outloc, DIFF_METHOD method) {
-  return Div_par(f, outloc, toString(method));
+  return Div_par(f, outloc, fmt::to_string(method));
 }
 
 // Divergence of a parallel flow: Div(f*v)
@@ -151,7 +151,7 @@ Field3D Div_par_flux(const Field3DParallel& v, const Field3DParallel& f,
 
 inline Field3D Div_par_flux(const Field3DParallel& v, const Field3DParallel& f,
                             CELL_LOC outloc, DIFF_METHOD method) {
-  return Div_par_flux(v, f, outloc, toString(method));
+  return Div_par_flux(v, f, outloc, fmt::to_string(method));
 }
 
 /*!
@@ -169,13 +169,13 @@ bout::FieldMetric Grad2_par2(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
                              const std::string& method = "DEFAULT");
 inline bout::FieldMetric Grad2_par2(const Field2D& f, CELL_LOC outloc,
                                     DIFF_METHOD method) {
-  return Grad2_par2(f, outloc, toString(method));
+  return Grad2_par2(f, outloc, fmt::to_string(method));
 }
 
 Field3D Grad2_par2(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT,
                    const std::string& method = "DEFAULT");
 inline Field3D Grad2_par2(const Field3DParallel& f, CELL_LOC outloc, DIFF_METHOD method) {
-  return Grad2_par2(f, outloc, toString(method));
+  return Grad2_par2(f, outloc, fmt::to_string(method));
 }
 
 /*!
@@ -219,10 +219,6 @@ Field3D Div_par_K_Grad_par(const Field3DParallel& kY, const Field3DParallel& f,
 /// different results when coefficients or cell sizes vary strongly.
 BOUT_ENUM_CLASS_NS(bout, ConductionMethod, Original, ProductJK, Harmonic);
 
-namespace bout {
-inline std::string format_as(const ConductionMethod& method) { return toString(method); }
-} // namespace bout
-
 /// Version with energy flow diagnostic
 /// For FCI fields, `flow_ylow` is currently returned as zero.
 Field3D
diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx
index 6e89275ce8..02a2ebe620 100644
--- a/include/bout/field2d.hxx
+++ b/include/bout/field2d.hxx
@@ -24,11 +24,14 @@
  * along with BOUT++.  If not, see <http://www.gnu.org/licenses/>.
  *
  */
+
+#pragma once
+
 #include "bout/utils.hxx"
 #include <optional>
+
 class Field2D;
 
-#pragma once
 #ifndef BOUT_FIELD2D_H
 #define BOUT_FIELD2D_H
 
@@ -43,6 +46,11 @@ class Field2D;
 #include "bout/fieldperp.hxx"
 #include "bout/region.hxx"
 
+#include "fmt/format.h"
+#include <fmt/base.h>
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
 #include <cstddef>
 #include <iostream>
 #include <ostream>
@@ -617,4 +625,8 @@ bool operator==(const Field2D& a, const Field2D& b);
 
 std::ostream& operator<<(std::ostream& out, const Field2D& value);
 
+// Tell {fmt} not to treat these as ranges so it uses format_as() instead
+template <>
+struct fmt::is_range<Field2D, char> : std::false_type {};
+
 #endif /* BOUT_FIELD2D_H */
diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx
index 04a566678e..c5af15ebdb 100644
--- a/include/bout/field3d.hxx
+++ b/include/bout/field3d.hxx
@@ -19,12 +19,14 @@
  * along with BOUT++.  If not, see <http://www.gnu.org/licenses/>.
  *
  **************************************************************************/
+#pragma once
 
+#include "fmt/base.h"
+#include "fmt/format.h"
 #include "bout/utils.hxx"
 #include <cstddef>
 class Field3D;
 
-#pragma once
 #ifndef BOUT_FIELD3D_H
 #define BOUT_FIELD3D_H
 
@@ -950,7 +952,7 @@ Field3D lowPass(const Field3D& var, int zmax, bool keep_zonal,
 [[deprecated("Please use a bool for `keep_zonal`")]] inline Field3D
 lowPass(const Field3D& var, int zmax, int keep_zonal, REGION rgn = RGN_ALL) {
   ASSERT0(static_cast<bool>(keep_zonal) == keep_zonal);
-  return lowPass(var, zmax, static_cast<bool>(keep_zonal), toString(rgn));
+  return lowPass(var, zmax, static_cast<bool>(keep_zonal), fmt::to_string(rgn));
 }
 
 /// Fourier low pass filtering. Removes modes higher than \p zmax
@@ -1237,4 +1239,8 @@ filledFrom(const Field3DParallel& f,
   return result;
 }
 
+// Tell {fmt} not to treat these as ranges so it uses format_as() instead
+template <>
+struct fmt::is_range<Field3D, char> : std::false_type {};
+
 #endif /* BOUT_FIELD3D_H */
diff --git a/include/bout/fieldperp.hxx b/include/bout/fieldperp.hxx
index 2f8408d82b..d4364d36d9 100644
--- a/include/bout/fieldperp.hxx
+++ b/include/bout/fieldperp.hxx
@@ -38,6 +38,9 @@ class FieldPerp;
 #include "bout/unused.hxx"
 #include "bout/utils.hxx"
 
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
 #include <cstddef>
 #include <optional>
 #include <ostream>
@@ -482,4 +485,8 @@ struct is_expr_fieldperp<BinaryExpr<ResT, L, R, Fun>>
                                        || (is_expr_fieldperp_v<std::decay_t<L>>
                                            && is_expr_constant_v<std::decay_t<R>>)> {};
 
+// Tell {fmt} not to treat these as ranges so it uses format_as() instead
+template <>
+struct fmt::is_range<FieldPerp, char> : std::false_type {};
+
 #endif
diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx
index d235746250..797f6f4ac3 100644
--- a/include/bout/hypre_interface.hxx
+++ b/include/bout/hypre_interface.hxx
@@ -1143,7 +1143,7 @@ public:
       break;
     }
     default: {
-      throw BoutException("Unsupported hypre_solver_type {}", toString(solver_type));
+      throw BoutException("Unsupported hypre_solver_type {}", solver_type);
     }
     }
 
@@ -1178,7 +1178,7 @@ public:
       break;
     }
     default: {
-      throw BoutException("Unsupported hypre_solver_type {}", toString(solver_type));
+      throw BoutException("Unsupported hypre_solver_type {}", solver_type);
     }
     }
 
diff --git a/include/bout/index_derivs.hxx b/include/bout/index_derivs.hxx
index a8d1bbb1d3..2127452f4e 100644
--- a/include/bout/index_derivs.hxx
+++ b/include/bout/index_derivs.hxx
@@ -68,7 +68,7 @@ inline std::ostream& operator<<(std::ostream& out, const metaData& meta) {
   out << ", ";
   out << "nGuards : " << meta.nGuards;
   out << ", ";
-  out << "type : " << toString(meta.derivType);
+  out << "type : " << fmt::to_string(meta.derivType);
   return out;
 }
 
diff --git a/include/bout/interpolation.hxx b/include/bout/interpolation.hxx
index bf9b1de466..3ec77994a2 100644
--- a/include/bout/interpolation.hxx
+++ b/include/bout/interpolation.hxx
@@ -88,7 +88,7 @@ interp_to(const T& var, CELL_LOC loc, const std::string& region = "RGN_ALL") {
   T result{emptyFrom(var).setLocation(loc)};
 
   // Staggered grids enabled, and need to perform interpolation
-  TRACE("Interpolating {} -> {}", toString(var.getLocation()), toString(loc));
+  TRACE("Interpolating {} -> {}", fmt::to_string(var.getLocation()), fmt::to_string(loc));
 
   if (region != "RGN_NOBNDRY") {
     // result is requested in some boundary region(s)
diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index 3d630b5ce3..dd7ccbfc0b 100644
--- a/include/bout/mesh.hxx
+++ b/include/bout/mesh.hxx
@@ -88,8 +88,6 @@ public:
 BOUT_ENUM_CLASS(BoundaryParType, all, xin, xout, fwd, bwd, xin_fwd, xout_fwd, xin_bwd,
                 xout_bwd, SIZE);
 
-inline auto format_as(const BoundaryParType& type) { return toString(type); }
-
 template <class DerivedType>
 using RegisterMesh = MeshFactory::RegisterInFactory<DerivedType>;
 
diff --git a/include/bout/options.hxx b/include/bout/options.hxx
index b7eb9eb29c..bac56ff113 100644
--- a/include/bout/options.hxx
+++ b/include/bout/options.hxx
@@ -608,9 +608,8 @@ public:
     if (bout::utils::variantEqualTo(attributes.at("source"), DEFAULT_SOURCE)) {
       // Check that the default values are the same
       if (!similar(val, def)) {
-        throw BoutException("Inconsistent default values for '{:s}': '{:s}' then '{:s}'",
-                            full_name, bout::utils::variantToString(value),
-                            toString(def));
+        throw BoutException("Inconsistent default values for '{}': '{}' then '{}'",
+                            full_name, bout::utils::variantToString(value), def);
       }
     }
     return val;
@@ -665,9 +664,8 @@ public:
     if (bout::utils::variantEqualTo(attributes.at("source"), DEFAULT_SOURCE)) {
       // Check that the default values are the same
       if (!similar(val, def)) {
-        throw BoutException("Inconsistent default values for '{:s}': '{:s}' then '{:s}'",
-                            full_name, bout::utils::variantToString(value),
-                            toString(def));
+        throw BoutException("Inconsistent default values for '{}': '{}' then '{}'",
+                            full_name, bout::utils::variantToString(value), def);
       }
     }
     return val;
@@ -913,13 +911,13 @@ private:
                       << bout::utils::variantToString(value) << " ("
                       << bout::utils::variantToString(attributes["source"])
                       << _(") overwritten with:") << "\n"
-                      << "\t\t" << full_name << " = " << toString(val) << " (" << source
-                      << ")\n";
+                      << "\t\t" << full_name << " = " << fmt::to_string(val) << " ("
+                      << source << ")\n";
         } else {
           throw BoutException(
-              _f("Options: Setting a value from same source ({:s}) to new value "
-                 "'{:s}' - old value was '{:s}'."),
-              source, toString(val), bout::utils::variantToString(value));
+              _f("Options: Setting a value from same source ({}) to new value "
+                 "'{}' - old value was '{:s}'."),
+              source, val, bout::utils::variantToString(value));
         }
       }
     }
@@ -1018,19 +1016,16 @@ auto Options::as(const Tensor<BoutReal>& similar_to) const -> Tensor<BoutReal>;
 template <>
 auto Options::as(const Tensor<int>& similar_to) const -> Tensor<int>;
 
-/// Convert \p value to string
-std::string toString(const Options& value);
-
 /// Save the parallel fields
 void saveParallel(Options& opt, const std::string& name, const Field3D& tosave);
 
 /// Output a stringified \p value to a stream
 ///
-/// This is templated to avoid implict casting: anything is
+/// This is templated to avoid implicit casting: anything is
 /// convertible to an `Options`, and we want _exactly_ `Options`
 template <class T, typename = bout::utils::EnableIfOptions<T>>
 inline std::ostream& operator<<(std::ostream& out, const T& value) {
-  return out << toString(value);
+  return out << fmt::to_string(value);
 }
 
 namespace bout {
diff --git a/include/bout/output_bout_types.hxx b/include/bout/output_bout_types.hxx
index 77ea8ebf2d..2bdef607bb 100644
--- a/include/bout/output_bout_types.hxx
+++ b/include/bout/output_bout_types.hxx
@@ -116,9 +116,10 @@ constexpr auto parse_nonnegative_int(const char*& begin, const char* end,
 /// - ``e<N>``: Number of elements at each edge to show
 /// - ``f``: Show full field
 template <class T>
-struct fmt::formatter<T, std::enable_if_t<bout::utils::is_Field_v<T>, char>> {
+struct fmt::formatter<T, char,
+                      std::enable_if_t<bout::utils::is_Field_v<std::remove_cv_t<T>>>> {
 private:
-  fmt::formatter<BoutReal> underlying;
+  formatter<BoutReal> underlying;
 
   static constexpr auto default_region = "RGN_ALL";
   std::string_view region = default_region;
diff --git a/include/bout/region.hxx b/include/bout/region.hxx
index 802b3758ba..88abec0b3f 100644
--- a/include/bout/region.hxx
+++ b/include/bout/region.hxx
@@ -364,16 +364,16 @@ using IndPerp = SpecificInd<IND_TYPE::IND_PERP>;
 using IndG3D = SpecificInd<IND_TYPE::IND_GLOBAL_3D>;
 
 /// Get string representation of Ind3D
-inline std::string toString(const Ind3D& i) {
+inline std::string format_as(const Ind3D& i) {
   return "(" + std::to_string(i.x()) + ", " + std::to_string(i.y()) + ", "
          + std::to_string(i.z()) + ")";
 }
 /// Get string representation of Ind2D
-inline std::string toString(const Ind2D& i) {
+inline std::string format_as(const Ind2D& i) {
   return "(" + std::to_string(i.x()) + ", " + std::to_string(i.y()) + ")";
 }
 /// Get string representation of IndPerp
-inline std::string toString(const IndPerp& i) {
+inline std::string format_as(const IndPerp& i) {
   return "(" + std::to_string(i.x()) + ", " + std::to_string(i.z()) + ")";
 }
 
diff --git a/include/bout/sys/variant.hxx b/include/bout/sys/variant.hxx
index 1b541f096d..f43c9b5926 100644
--- a/include/bout/sys/variant.hxx
+++ b/include/bout/sys/variant.hxx
@@ -15,6 +15,7 @@
 
 #pragma once
 
+#include <string>
 #ifndef __VARIANT_HXX__
 #define __VARIANT_HXX__
 
diff --git a/include/bout/utils.hxx b/include/bout/utils.hxx
index 420ded10d6..b80399561a 100644
--- a/include/bout/utils.hxx
+++ b/include/bout/utils.hxx
@@ -37,6 +37,8 @@
 #include "bout/region.hxx"
 #include "bout/unused.hxx"
 
+#include <fmt/ranges.h>
+
 #include <algorithm>
 #include <cmath>
 #include <ctime>
@@ -530,11 +532,6 @@ std::string toString(const T& val) {
   return ss.str();
 }
 
-/// Simple case where input is already a string
-/// This is so that toString can be used in templates
-/// where the type may be std::string.
-inline std::string toString(const std::string& val) { return val; }
-
 template <typename T>
 inline std::string toString(const Array<T>& UNUSED(val)) {
   return "<Array>";
@@ -550,6 +547,22 @@ inline std::string toString(const Tensor<T>& UNUSED(val)) {
   return "<Tensor>";
 }
 
+// Tell {fmt} not to treat these as ranges so it uses format_as() instead
+template <typename T, typename Char>
+struct fmt::range_format_kind<Array<T>, Char> {
+  static constexpr auto value = fmt::range_format::disabled;
+};
+
+template <typename T, typename Char>
+struct fmt::range_format_kind<Matrix<T>, Char> {
+  static constexpr auto value = fmt::range_format::disabled;
+};
+
+template <typename T, typename Char>
+struct fmt::range_format_kind<Tensor<T>, Char> {
+  static constexpr auto value = fmt::range_format::disabled;
+};
+
 /// Convert a bool to "true" or "false"
 inline std::string toString(const bool& val) {
   if (val) {
@@ -558,10 +571,6 @@ inline std::string toString(const bool& val) {
   return "false";
 }
 
-/// Convert a time stamp to a string
-/// This uses std::localtime and std::put_time
-std::string toString(const time_t& time);
-
 /*!
  * Convert a string to lower case
  */
diff --git a/include/bout/vecops.hxx b/include/bout/vecops.hxx
index 9166503855..60985ec80e 100644
--- a/include/bout/vecops.hxx
+++ b/include/bout/vecops.hxx
@@ -64,7 +64,7 @@ Vector3D Grad(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
 Vector3D Grad_perp(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
                    const std::string& method = "DEFAULT");
 inline Vector3D Grad_perp(const Field3D& f, CELL_LOC outloc, DIFF_METHOD method) {
-  return Grad_perp(f, outloc, toString(method));
+  return Grad_perp(f, outloc, fmt::to_string(method));
 }
 Vector2D Grad_perp(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
                    const std::string& method = "DEFAULT");
@@ -92,7 +92,7 @@ Field3D Div(const Vector3D& v, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
 
 inline Field3D Div(const Vector3D& v, const Field3D& f, CELL_LOC outloc,
                    DIFF_METHOD method = DIFF_DEFAULT) {
-  return Div(v, f, outloc, toString(method));
+  return Div(v, f, outloc, fmt::to_string(method));
 }
 
 /// Curl of a vector
diff --git a/src/field/field2d.cxx b/src/field/field2d.cxx
index 6025909550..dca47b8dc6 100644
--- a/src/field/field2d.cxx
+++ b/src/field/field2d.cxx
@@ -25,6 +25,7 @@
  *
  **************************************************************************/
 
+#include "fmt/format.h"
 #include "bout/bout_types.hxx"
 #include "bout/build_defines.hxx"
 
@@ -43,6 +44,7 @@
 #include <cmath>
 #include <cstddef>
 #include <optional>
+#include <ostream>
 
 Field2D::Field2D(Mesh* localmesh, CELL_LOC location_in, DirectionTypes directions_in,
                  [[maybe_unused]] std::optional<size_t> regionID)
@@ -119,7 +121,7 @@ Field2D* Field2D::timeDeriv() {
 ////////////// Indexing ///////////////////
 
 const Region<Ind2D>& Field2D::getRegion(REGION region) const {
-  return fieldmesh->getRegion2D(toString(region));
+  return fieldmesh->getRegion2D(fmt::to_string(region));
 }
 const Region<Ind2D>& Field2D::getRegion(const std::string& region_name) const {
   return fieldmesh->getRegion2D(region_name);
diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx
index 16b9e7a50b..34a8594bc5 100644
--- a/src/field/field3d.cxx
+++ b/src/field/field3d.cxx
@@ -40,6 +40,7 @@
 #include <cstddef>
 #include <memory>
 #include <optional>
+#include <ostream>
 #include <string>
 #include <utility>
 
@@ -264,14 +265,14 @@ const BoutReal& Field3D::operator()(const Ind2D& d, int jz) const {
 }
 
 const Region<Ind3D>& Field3D::getRegion(REGION region) const {
-  return fieldmesh->getRegion3D(toString(region));
+  return fieldmesh->getRegion3D(fmt::to_string(region));
 }
 const Region<Ind3D>& Field3D::getRegion(const std::string& region_name) const {
   return fieldmesh->getRegion3D(region_name);
 }
 
 const Region<Ind2D>& Field3D::getRegion2D(REGION region) const {
-  return fieldmesh->getRegion2D(toString(region));
+  return fieldmesh->getRegion2D(fmt::to_string(region));
 }
 const Region<Ind2D>& Field3D::getRegion2D(const std::string& region_name) const {
   return fieldmesh->getRegion2D(region_name);
@@ -715,7 +716,7 @@ Field3D filter(const Field3D& var, int N0, const std::string& rgn) {
 
   Field3D result{emptyFrom(var)};
 
-  const auto region_str = toString(rgn);
+  const auto region_str = fmt::to_string(rgn);
 
   // Only allow a whitelist of regions for now
   ASSERT2(region_str == "RGN_ALL" || region_str == "RGN_NOBNDRY"
@@ -765,7 +766,7 @@ Field3D lowPass(const Field3D& var, int zmax, bool keep_zonal, const std::string
 
   Field3D result{emptyFrom(var)};
 
-  const auto region_str = toString(rgn);
+  const auto region_str = fmt::to_string(rgn);
 
   // Only allow a whitelist of regions for now
   ASSERT2(region_str == "RGN_ALL" || region_str == "RGN_NOBNDRY"
@@ -829,7 +830,7 @@ void shiftZ(Field3D& var, int jx, int jy, double zangle) {
 }
 
 void shiftZ(Field3D& var, double zangle, const std::string& rgn) {
-  const auto region_str = toString(rgn);
+  const auto region_str = fmt::to_string(rgn);
 
   // Only allow a whitelist of regions for now
   ASSERT2(region_str == "RGN_ALL" || region_str == "RGN_NOBNDRY"
diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx
index bca551c29d..cb0bfe6883 100644
--- a/src/field/fieldperp.cxx
+++ b/src/field/fieldperp.cxx
@@ -107,7 +107,7 @@ FieldPerp& FieldPerp::operator=(const BoutReal rhs) {
 }
 
 const Region<IndPerp>& FieldPerp::getRegion(REGION region) const {
-  return fieldmesh->getRegionPerp(toString(region));
+  return fieldmesh->getRegionPerp(fmt::to_string(region));
 }
 const Region<IndPerp>& FieldPerp::getRegion(const std::string& region_name) const {
   return fieldmesh->getRegionPerp(region_name);
diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx
index 0b558a15ec..713ffd7549 100644
--- a/src/mesh/coordinates.cxx
+++ b/src/mesh/coordinates.cxx
@@ -624,7 +624,7 @@ void Coordinates::readParallelMetricComponents() {
 void Coordinates::outputVars(Options& output_options) {
   const Timer time("io");
   const std::string loc_string =
-      (location == CELL_CENTRE) ? "" : "_" + toString(location);
+      (location == CELL_CENTRE) ? "" : "_" + fmt::to_string(location);
 
   output_options["dx" + loc_string].force(dx(), "Coordinates");
   output_options["dy" + loc_string].force(dy(), "Coordinates");
diff --git a/src/mesh/interpolation_xz.cxx b/src/mesh/interpolation_xz.cxx
index bbc859361e..b7bf2862c4 100644
--- a/src/mesh/interpolation_xz.cxx
+++ b/src/mesh/interpolation_xz.cxx
@@ -35,8 +35,8 @@
 #include <bout/unused.hxx>
 #include <bout/utils.hxx>
 
-void printLocation(const Field3D& var) { output << toString(var.getLocation()); }
-void printLocation(const Field2D& var) { output << toString(var.getLocation()); }
+void printLocation(const Field3D& var) { output << fmt::to_string(var.getLocation()); }
+void printLocation(const Field2D& var) { output << fmt::to_string(var.getLocation()); }
 
 const Field3D interpolate(const Field3D& f, const Field3D& delta_x,
                           const Field3D& delta_z) {
diff --git a/src/mesh/parallel/shiftedmetric.cxx b/src/mesh/parallel/shiftedmetric.cxx
index c98333edc3..18e993e285 100644
--- a/src/mesh/parallel/shiftedmetric.cxx
+++ b/src/mesh/parallel/shiftedmetric.cxx
@@ -48,7 +48,7 @@ void ShiftedMetric::checkInputGrid() {
 void ShiftedMetric::outputVars(Options& output_options) {
   Timer time("io");
   const std::string loc_string =
-      (location == CELL_CENTRE) ? "" : "_" + toString(location);
+      (location == CELL_CENTRE) ? "" : "_" + fmt::to_string(location);
 
   output_options["zShift" + loc_string].force(zShift, "ShiftedMetric");
 }
diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx
index c5dab45dd9..c2fd745988 100644
--- a/src/solver/impls/arkode/arkode.cxx
+++ b/src/solver/impls/arkode/arkode.cxx
@@ -232,7 +232,7 @@ int ArkodeSolver::init() {
         callWithSUNContext(ARKStepCreate, suncontext, nullptr, arkode_rhs, simtime, uvec);
     break;
   default:
-    throw BoutException("Invalid treatment: {}\n", toString(treatment));
+    throw BoutException("Invalid treatment: {}\n", treatment);
   }
   if (arkode_mem == nullptr) {
     throw BoutException("ARKStepCreate failed\n");
diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx
index df0013dd3b..cceb5b2794 100644
--- a/src/solver/solver.cxx
+++ b/src/solver/solver.cxx
@@ -281,7 +281,7 @@ void Solver::add(Field3D& v, const std::string& name, const std::string& descrip
 
   if (mesh->StaggerGrids && (v.getLocation() != CELL_CENTRE)) {
     output_info.write("\tVariable {:s} shifted to {:s}\n", name,
-                      toString(v.getLocation()));
+                      fmt::to_string(v.getLocation()));
     ddt(v).setLocation(v.getLocation()); // Make sure both at the same location
   }
 
@@ -605,7 +605,7 @@ int Solver::solve(int nout, BoutReal timestep) {
   }
 
   time_t start_time = time(nullptr);
-  output_progress.write(_f("\nRun started at  : {:s}\n"), toString(start_time));
+  output_progress.write(_f("\nRun started at  : {:s}\n"), fmt::to_string(start_time));
 
   Timer timer("run"); // Start timer
 
@@ -651,7 +651,7 @@ int Solver::solve(int nout, BoutReal timestep) {
     status = run();
 
     time_t end_time = time(nullptr);
-    output_progress.write(_f("\nRun finished at  : {:s}\n"), toString(end_time));
+    output_progress.write(_f("\nRun finished at  : {:s}\n"), fmt::to_string(end_time));
     output_progress.write(_("Run time : "));
 
     int dt = end_time - start_time;
@@ -1274,9 +1274,9 @@ std::vector<Solver::JacobianVariableMetadata> Solver::getJacobianMetadata2D() co
   metadata.reserve(f2d.size());
 
   for (int i = 0; i < static_cast<int>(f2d.size()); ++i) {
-    metadata.push_back(JacobianVariableMetadata{i, f2d[i].name, toString(f2d[i].location),
-                                                f2d[i].evolve_bndry, f2d[i].constraint,
-                                                f2d[i].description});
+    metadata.push_back(JacobianVariableMetadata{
+        i, f2d[i].name, fmt::to_string(f2d[i].location), f2d[i].evolve_bndry,
+        f2d[i].constraint, f2d[i].description});
   }
 
   return metadata;
@@ -1287,9 +1287,9 @@ std::vector<Solver::JacobianVariableMetadata> Solver::getJacobianMetadata3D() co
   metadata.reserve(f3d.size());
 
   for (int i = 0; i < static_cast<int>(f3d.size()); ++i) {
-    metadata.push_back(JacobianVariableMetadata{i, f3d[i].name, toString(f3d[i].location),
-                                                f3d[i].evolve_bndry, f3d[i].constraint,
-                                                f3d[i].description});
+    metadata.push_back(JacobianVariableMetadata{
+        i, f3d[i].name, fmt::to_string(f3d[i].location), f3d[i].evolve_bndry,
+        f3d[i].constraint, f3d[i].description});
   }
 
   return metadata;
@@ -1389,8 +1389,9 @@ std::string Solver::getJacobianMatrixFilename(const std::string& jacobian_export
                                               bout::PetscMatrixExportFormat format) {
   // The directory the output data is stored in
   const std::string datadir = Options::root()["datadir"];
-  const std::string stem = fmt::format("{}/{}_{}_{:06d}", datadir, jacobian_export_prefix,
-                                       toString(kind), this->jacobian_export_counter++);
+  const std::string stem =
+      fmt::format("{}/{}_{}_{:06d}", datadir, jacobian_export_prefix,
+                  fmt::to_string(kind), this->jacobian_export_counter++);
   return stem + (format == bout::PetscMatrixExportFormat::binary ? ".dat" : ".txt");
 }
 
diff --git a/src/sys/bout_types.cxx b/src/sys/bout_types.cxx
index cbcf19e996..817068de76 100644
--- a/src/sys/bout_types.cxx
+++ b/src/sys/bout_types.cxx
@@ -1,3 +1,4 @@
+#include "fmt/format.h"
 #include <bout/bout_types.hxx>
 #include <bout/boutexception.hxx>
 
@@ -25,7 +26,7 @@ const T& safeAt(const std::map<std::string, T>& mymap, const std::string& s) {
 }
 } // namespace
 
-std::string toString(CELL_LOC location) {
+std::string format_as(const CELL_LOC& location) {
 
   const static std::map<CELL_LOC, std::string> CELL_LOCtoString = {
       ENUMSTR(CELL_DEFAULT), ENUMSTR(CELL_CENTRE), ENUMSTR(CELL_XLOW),
@@ -43,17 +44,17 @@ CELL_LOC CELL_LOCFromString(const std::string& location_string) {
   return safeAt(stringtoCELL_LOC, location_string);
 }
 
-std::string toString(DIFF_METHOD location) {
+std::string format_as(const DIFF_METHOD& method) {
 
   const static std::map<DIFF_METHOD, std::string> DIFF_METHODtoString = {
       {DIFF_DEFAULT, "DEFAULT"}, {DIFF_U1, "U1"},   {DIFF_U2, "U2"},      {DIFF_U3, "U3"},
       {DIFF_C2, "C2"},           {DIFF_C4, "C4"},   {DIFF_S2, "S2"},      {DIFF_W2, "W2"},
       {DIFF_W3, "W3"},           {DIFF_FFT, "FFT"}, {DIFF_SPLIT, "SPLIT"}};
 
-  return safeAt(DIFF_METHODtoString, location);
+  return safeAt(DIFF_METHODtoString, method);
 }
 
-std::string toString(REGION region) {
+std::string format_as(const REGION& region) {
 
   const static std::map<REGION, std::string> REGIONtoString = {
       ENUMSTR(RGN_ALL), ENUMSTR(RGN_NOBNDRY), ENUMSTR(RGN_NOX), ENUMSTR(RGN_NOY),
@@ -61,7 +62,7 @@ std::string toString(REGION region) {
   return safeAt(REGIONtoString, region);
 }
 
-std::string toString(DIRECTION direction) {
+std::string format_as(const DIRECTION& direction) {
 
   const static std::map<DIRECTION, std::string> DIRECTIONtoString = {
       {DIRECTION::X, "X"},
@@ -73,8 +74,8 @@ std::string toString(DIRECTION direction) {
   return safeAt(DIRECTIONtoString, direction);
 }
 
-std::string toString(const DirectionTypes& dir) {
-  return "{" + toString(dir.y) + ", " + toString(dir.z) + "}";
+std::string format_as(const DirectionTypes& dir) {
+  return "{" + fmt::to_string(dir.y) + ", " + fmt::to_string(dir.z) + "}";
 }
 
 void swap(DirectionTypes& first, DirectionTypes& second) {
@@ -119,7 +120,7 @@ bool areDirectionsCompatible(const DirectionTypes& d1, const DirectionTypes& d2)
   return false;
 }
 
-std::string toString(STAGGER stagger) {
+std::string format_as(const STAGGER& stagger) {
 
   const static std::map<STAGGER, std::string> STAGGERtoString = {
       {STAGGER::None, "No staggering"},
@@ -129,7 +130,7 @@ std::string toString(STAGGER stagger) {
   return safeAt(STAGGERtoString, stagger);
 }
 
-std::string toString(DERIV deriv) {
+std::string format_as(const DERIV& deriv) {
 
   const static std::map<DERIV, std::string> DERIVtoString = {
       {DERIV::Standard, "Standard"},
@@ -141,7 +142,7 @@ std::string toString(DERIV deriv) {
   return safeAt(DERIVtoString, deriv);
 }
 
-std::string toString(YDirectionType d) {
+std::string format_as(const YDirectionType& d) {
 
   const static std::map<YDirectionType, std::string> YDirectionTypeToString = {
       {YDirectionType::Standard, "Standard"}, {YDirectionType::Aligned, "Aligned"}};
@@ -157,7 +158,7 @@ YDirectionType YDirectionTypeFromString(const std::string& y_direction_string) {
   return safeAt(stringToYDirectionType, y_direction_string);
 }
 
-std::string toString(ZDirectionType d) {
+std::string format_as(const ZDirectionType& d) {
 
   const static std::map<ZDirectionType, std::string> ZDirectionTypeToString = {
       {ZDirectionType::Standard, "Standard"}, {ZDirectionType::Average, "Average"}};
diff --git a/src/sys/options.cxx b/src/sys/options.cxx
index 09e5ba1ba4..1503d6200d 100644
--- a/src/sys/options.cxx
+++ b/src/sys/options.cxx
@@ -312,25 +312,25 @@ bool Options::isSection(const std::string& name) const {
 
 template <>
 Options& Options::assign<>(Field2D val, std::string source) {
-  attributes["cell_location"] = toString(val.getLocation());
-  attributes["direction_y"] = toString(val.getDirectionY());
-  attributes["direction_z"] = toString(val.getDirectionZ());
+  attributes["cell_location"] = fmt::to_string(val.getLocation());
+  attributes["direction_y"] = fmt::to_string(val.getDirectionY());
+  attributes["direction_z"] = fmt::to_string(val.getDirectionZ());
   _set_no_check(std::move(val), std::move(source));
   return *this;
 }
 template <>
 Options& Options::assign<>(Field3D val, std::string source) {
-  attributes["cell_location"] = toString(val.getLocation());
-  attributes["direction_y"] = toString(val.getDirectionY());
-  attributes["direction_z"] = toString(val.getDirectionZ());
+  attributes["cell_location"] = fmt::to_string(val.getLocation());
+  attributes["direction_y"] = fmt::to_string(val.getDirectionY());
+  attributes["direction_z"] = fmt::to_string(val.getDirectionZ());
   _set_no_check(std::move(val), std::move(source));
   return *this;
 }
 template <>
 Options& Options::assign<>(FieldPerp val, std::string source) {
-  attributes["cell_location"] = toString(val.getLocation());
-  attributes["direction_y"] = toString(val.getDirectionY());
-  attributes["direction_z"] = toString(val.getDirectionZ());
+  attributes["cell_location"] = fmt::to_string(val.getLocation());
+  attributes["direction_y"] = fmt::to_string(val.getDirectionY());
+  attributes["direction_z"] = fmt::to_string(val.getDirectionZ());
   attributes["yindex_global"] = val.getGlobalIndex();
   _set_no_check(std::move(val), std::move(source));
   return *this;
@@ -542,7 +542,7 @@ bool Options::as<bool>(const bool& UNUSED(similar_to)) const {
 
   value_used = true;
 
-  printNameValueSourceLine(*this, toString(result));
+  printNameValueSourceLine(*this, fmt::to_string(result));
 
   return result;
 }
@@ -1197,8 +1197,6 @@ bout::details::OptionsFormatterBase::format(const Options& options,
   return ctx.out();
 }
 
-std::string toString(const Options& value) { return fmt::format("{}", value); }
-
 namespace bout {
 void checkForUnusedOptions() {
   auto& options = Options::root();
diff --git a/src/sys/utils.cxx b/src/sys/utils.cxx
index e15f9248bb..04ab16a2b5 100644
--- a/src/sys/utils.cxx
+++ b/src/sys/utils.cxx
@@ -35,6 +35,7 @@
 #include <vector>
 
 #include "fmt/chrono.h"
+#include "fmt/format.h"
 
 /**************************************************************************
  * String routines
diff --git a/tests/integrated/test-fci-boundary/get_par_bndry.cxx b/tests/integrated/test-fci-boundary/get_par_bndry.cxx
index fd3bc2889f..1954496e5b 100644
--- a/tests/integrated/test-fci-boundary/get_par_bndry.cxx
+++ b/tests/integrated/test-fci-boundary/get_par_bndry.cxx
@@ -22,7 +22,7 @@ int main(int argc, char** argv) {
   for (int i = 0; i < fields.size(); i++) {
     fields[i].allocate();
     const auto boundary = static_cast<BoundaryParType>(i);
-    const auto boundary_name = toString(boundary);
+    const auto boundary_name = fmt::to_string(boundary);
     mesh->communicate(fields[i]);
     for (auto& bndry_par : mesh->getBoundariesPar(boundary)) {
       output.write("{:s} region\n", boundary_name);
diff --git a/tests/integrated/test-invpar/test_invpar.cxx b/tests/integrated/test-invpar/test_invpar.cxx
index 99b4e99c5c..ce50110c0e 100644
--- a/tests/integrated/test-invpar/test_invpar.cxx
+++ b/tests/integrated/test-invpar/test_invpar.cxx
@@ -16,7 +16,7 @@ int test(const std::string& acoef, const std::string& bcoef, const std::string&
          BoutReal tol, FieldFactory& f, CELL_LOC location) {
   output_error.write(
       "acoef={:s} bcoef={:s} ccoef={:s} dcoef={:s} ecoef{:s} with {:s} at {:s} ...",
-      acoef, bcoef, ccoef, dcoef, ecoef, func, toString(location));
+      acoef, bcoef, ccoef, dcoef, ecoef, func, fmt::to_string(location));
   auto inv = InvertPar::create(nullptr, location, mesh);
 
   Field2D A = f.create2D(acoef, nullptr, nullptr, location);
diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx
index cbffccabc5..58370b5cff 100644
--- a/tests/unit/field/test_field_factory.cxx
+++ b/tests/unit/field/test_field_factory.cxx
@@ -32,9 +32,10 @@ using bout::generator::Context;
 // Reuse the "standard" fixture for FakeMesh
 template <typename T>
 class FieldFactoryCreationTest : public FakeMeshFixture {
+  Options options;
+
 public:
   FieldFactoryCreationTest() : FakeMeshFixture() {
-    Options options;
     options["input"]["transform_from_field_aligned"] = false;
     factory = FieldFactory{mesh, &options};
   }
diff --git a/tests/unit/field/test_fieldperp.cxx b/tests/unit/field/test_fieldperp.cxx
index 754316e077..24a506d8c5 100644
--- a/tests/unit/field/test_fieldperp.cxx
+++ b/tests/unit/field/test_fieldperp.cxx
@@ -3,6 +3,7 @@
 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
 #include "gtest/gtest.h"
 
+#include "fmt/format.h"
 #include "test_extras.hxx"
 #include "bout/array.hxx"
 #include "bout/boutexception.hxx"
diff --git a/tests/unit/include/bout/test_bout_enum_class.cxx b/tests/unit/include/bout/test_bout_enum_class.cxx
index 83508a522e..7b22388c67 100644
--- a/tests/unit/include/bout/test_bout_enum_class.cxx
+++ b/tests/unit/include/bout/test_bout_enum_class.cxx
@@ -1,3 +1,4 @@
+#include "fmt/format.h"
 #include <gtest/gtest.h>
 
 #include "bout/bout_enum_class.hxx"
@@ -10,13 +11,13 @@ BOUT_ENUM_CLASS_NS(bout, TestEnum, foo, bar);
 BOUT_ENUM_CLASS_NS(bout::test, TestEnum, foo, bar);
 
 TEST(BoutEnumClassNS, toString) {
-  EXPECT_EQ(toString(bout::TestEnum::foo), "foo");
-  EXPECT_EQ(toString(bout::TestEnum::bar), "bar");
+  EXPECT_EQ(fmt::to_string(bout::TestEnum::foo), "foo");
+  EXPECT_EQ(fmt::to_string(bout::TestEnum::bar), "bar");
 }
 
 TEST(BoutEnumClass, toString) {
-  EXPECT_EQ(toString(TestEnum::foo), "foo");
-  EXPECT_EQ(toString(TestEnum::bar), "bar");
+  EXPECT_EQ(fmt::to_string(TestEnum::foo), "foo");
+  EXPECT_EQ(fmt::to_string(TestEnum::bar), "bar");
 }
 
 TEST(BoutEnumClassNS, fromString) {
diff --git a/tests/unit/sys/test_options_adios2.cxx b/tests/unit/sys/test_options_adios2.cxx
index 1e299e21ab..5b336766d5 100644
--- a/tests/unit/sys/test_options_adios2.cxx
+++ b/tests/unit/sys/test_options_adios2.cxx
@@ -236,7 +236,7 @@ TEST_F(OptionsAdios2Test, Field2DWriteCellCentre) {
   Options data = OptionsIO::create(file_options)->read();
 
   EXPECT_EQ(data["f2d"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_CENTRE));
+            fmt::to_string(CELL_CENTRE));
 }
 
 TEST_F(OptionsAdios2Test, Field2DWriteCellYLow) {
@@ -252,7 +252,7 @@ TEST_F(OptionsAdios2Test, Field2DWriteCellYLow) {
   Options data = OptionsIO::create(file_options)->read();
 
   EXPECT_EQ(data["f2d"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_YLOW));
+            fmt::to_string(CELL_YLOW));
 }
 
 TEST_F(OptionsAdios2Test, Field3DWriteCellCentre) {
@@ -268,7 +268,7 @@ TEST_F(OptionsAdios2Test, Field3DWriteCellCentre) {
   Options data = OptionsIO::create(file_options)->read();
 
   EXPECT_EQ(data["f3d"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_CENTRE));
+            fmt::to_string(CELL_CENTRE));
 }
 
 TEST_F(OptionsAdios2Test, Field3DWriteCellYLow) {
@@ -284,7 +284,7 @@ TEST_F(OptionsAdios2Test, Field3DWriteCellYLow) {
   Options data = OptionsIO::create(file_options)->read();
 
   EXPECT_EQ(data["f3d"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_YLOW));
+            fmt::to_string(CELL_YLOW));
 }
 
 TEST_F(OptionsAdios2Test, FieldPerpWriteCellCentre) {
@@ -306,7 +306,7 @@ TEST_F(OptionsAdios2Test, FieldPerpWriteCellCentre) {
   Options data = OptionsIO::create(file_options)->read();
 
   EXPECT_EQ(data["fperp"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_CENTRE));
+            fmt::to_string(CELL_CENTRE));
   EXPECT_EQ(data["fperp"].attributes["yindex_global"].as<int>(), 2);
 }
 
diff --git a/tests/unit/sys/test_options_netcdf.cxx b/tests/unit/sys/test_options_netcdf.cxx
index b16d838b54..e30f286e10 100644
--- a/tests/unit/sys/test_options_netcdf.cxx
+++ b/tests/unit/sys/test_options_netcdf.cxx
@@ -189,7 +189,7 @@ TEST_F(OptionsNetCDFTest, Field2DWriteCellCentre) {
   Options data = OptionsIO::create(filename)->read();
 
   EXPECT_EQ(data["f2d"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_CENTRE));
+            fmt::to_string(CELL_CENTRE));
 }
 
 TEST_F(OptionsNetCDFTest, Field2DWriteCellYLow) {
@@ -205,7 +205,7 @@ TEST_F(OptionsNetCDFTest, Field2DWriteCellYLow) {
   Options data = OptionsIO::create(filename)->read();
 
   EXPECT_EQ(data["f2d"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_YLOW));
+            fmt::to_string(CELL_YLOW));
 }
 
 TEST_F(OptionsNetCDFTest, Field3DWriteCellCentre) {
@@ -221,7 +221,7 @@ TEST_F(OptionsNetCDFTest, Field3DWriteCellCentre) {
   Options data = OptionsIO::create(filename)->read();
 
   EXPECT_EQ(data["f3d"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_CENTRE));
+            fmt::to_string(CELL_CENTRE));
 }
 
 TEST_F(OptionsNetCDFTest, Field3DWriteCellYLow) {
@@ -237,7 +237,7 @@ TEST_F(OptionsNetCDFTest, Field3DWriteCellYLow) {
   Options data = OptionsIO::create(filename)->read();
 
   EXPECT_EQ(data["f3d"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_YLOW));
+            fmt::to_string(CELL_YLOW));
 }
 
 TEST_F(OptionsNetCDFTest, FieldPerpWriteCellCentre) {
@@ -259,7 +259,7 @@ TEST_F(OptionsNetCDFTest, FieldPerpWriteCellCentre) {
   Options data = OptionsIO::create(filename)->read();
 
   EXPECT_EQ(data["fperp"].attributes["cell_location"].as<std::string>(),
-            toString(CELL_CENTRE));
+            fmt::to_string(CELL_CENTRE));
   EXPECT_EQ(data["fperp"].attributes["yindex_global"].as<int>(), 2);
 }
 
diff --git a/tests/unit/sys/test_utils.cxx b/tests/unit/sys/test_utils.cxx
index 6d84813c48..7bcb48e663 100644
--- a/tests/unit/sys/test_utils.cxx
+++ b/tests/unit/sys/test_utils.cxx
@@ -667,7 +667,7 @@ TEST(FunctionTraitsTest, SecondArg) {
 #ifndef __cpp_lib_erase_if
 TEST(StdLibBackPorts, EraseIfMultiset) {
   std::multiset<int> data{3, 3, 4, 5, 5, 6, 6, 7, 2, 1, 0};
-  auto divisible_by_3 = [](auto const& x) { return (x % 3) == 0; };
+  auto divisible_by_3 = [](const auto& x) { return (x % 3) == 0; };
 
   bout::utils::erase_if(data, divisible_by_3);
 
@@ -678,7 +678,7 @@ TEST(StdLibBackPorts, EraseIfMultiset) {
 
 TEST(StdLibBackPorts, EraseIfSet) {
   std::set<int> data{3, 4, 5, 6, 7, 2, 1, 0};
-  auto divisible_by_3 = [](auto const& x) { return (x % 3) == 0; };
+  auto divisible_by_3 = [](const auto& x) { return (x % 3) == 0; };
 
   bout::utils::erase_if(data, divisible_by_3);
 
@@ -696,7 +696,7 @@ TEST(StdLibBackPorts, EraseVector) {
 
 TEST(StdLibBackPorts, EraseIfVector) {
   std::vector<int> data{3, 3, 4, 5, 5, 6, 6, 7, 2, 1, 0};
-  auto divisible_by_3 = [](auto const& x) { return (x % 3) == 0; };
+  auto divisible_by_3 = [](const auto& x) { return (x % 3) == 0; };
 
   bout::utils::erase_if(data, divisible_by_3);