diff --git a/dumux/assembly/partialreassembler.hh b/dumux/assembly/partialreassembler.hh
index 52defa5d3ca999867a3b7bc037a026ad8bf2feaa..88c436fa7d00e364938a57c1c87dfccfa0941a71 100644
--- a/dumux/assembly/partialreassembler.hh
+++ b/dumux/assembly/partialreassembler.hh
@@ -30,6 +30,7 @@
 #include <dune/grid/common/gridenums.hh>
 #include <dune/istl/multitypeblockvector.hh>
 
+#include <dumux/io/format.hh>
 #include <dumux/common/typetraits/isvalid.hh>
 #include <dumux/discretization/method.hh>
 #include <dumux/parallel/vectorcommdatahandle.hh>
@@ -486,11 +487,10 @@ public:
         if (comm.size() > 1)
             greenElems_ = comm.sum(greenElems_);
 
-        auto reassembledElems = totalElems_ - greenElems_;
-        auto width = std::to_string(totalElems_).size();
-        outStream << ", reassembled " << std::setw(width)
-            << reassembledElems << " (" << std::setw(3)
-            << 100*reassembledElems/totalElems_ << "%) elements";
+        const auto reassembledElems = totalElems_ - greenElems_;
+        const auto width = Fmt::formatted_size("{}", totalElems_);
+        outStream << Fmt::format(", reassembled {:{}} ({:3}%) elements",
+                                 reassembledElems, width, 100*reassembledElems/totalElems_);
     }
 
     EntityColor elementColor(size_t idx) const
diff --git a/dumux/common/timeloop.hh b/dumux/common/timeloop.hh
index d7c3cdbea94002c5deadc5822ac386b97618cce6..9086c67689d6e91b41012f98a532f0c1433c7c96 100644
--- a/dumux/common/timeloop.hh
+++ b/dumux/common/timeloop.hh
@@ -38,6 +38,7 @@
 #include <dune/common/exceptions.hh>
 
 #include <dumux/common/parameters.hh>
+#include <dumux/io/format.hh>
 
 namespace Dumux {
 
@@ -233,7 +234,7 @@ public:
     {
         endTime_ = t;
         if (verbose_)
-            std::cout << "Set new end time to t = " << t << " seconds." << std::endl;
+            std::cout << Fmt::format("Set new end time to t = {:.5g} seconds.\n", t);
     }
 
     /*!
@@ -257,11 +258,8 @@ public:
         using std::min;
         timeStepSize_ = min(dt, maxTimeStepSize());
         if (!finished() && Dune::FloatCmp::le(timeStepSize_, 0.0, 1e-14*endTime_))
-        {
-            std::cerr << "You have set a very small timestep size (dt = "
-                      << timeStepSize_ << "). This might lead to numerical problems!"
-                      << std::endl;
-        }
+            std::cerr << Fmt::format("You have set a very small timestep size (dt = {:.5g}).", timeStepSize_)
+                      << " This might lead to numerical problems!\n";
     }
 
     /*!
@@ -351,14 +349,9 @@ public:
         {
             const auto cpuTime = wallClockTime();
             const auto percent = std::round( time_ / endTime_ * 100 );
-            std::cout << "[" << std::fixed << std::setw( 3 ) << std::setfill( ' ' )
-                      << std::setprecision( 0 )  << percent << "%] "
-                      << "Time step " << timeStepIdx_ << " done in "
-                      << std::setprecision( 6 ) << timeStepWallClockTime_ << " seconds. "
-                      << "Wall clock time: " << std::setprecision( 3 ) << cpuTime
-                      << ", time: " << std::setprecision( 5 ) << time_
-                      << ", time step size: " << std::setprecision( 8 ) << previousTimeStepSize_
-                      << std::endl;
+            std::cout << Fmt::format("[{:3.0f}%] ", percent)
+                      << Fmt::format("Time step {} done in {:.2g} seconds. ", timeStepIdx_, timeStepWallClockTime_)
+                      << Fmt::format("Wall clock time: {:.2g}, time: {:.5g}, time step size: {:.5g}\n", cpuTime, time_, previousTimeStepSize_);
         }
     }
 
@@ -371,18 +364,13 @@ public:
         auto cpuTime = timer_.stop();
 
         if (verbose_)
-        {
-            std::cout << "Simulation took " << cpuTime << " seconds on "
-                      << comm.size() << " processes.\n";
-        }
+            std::cout << Fmt::format("Simulation took {:.2g} seconds on {} processes.\n", cpuTime, comm.size());
 
         if (comm.size() > 1)
             cpuTime = comm.sum(cpuTime);
 
         if (verbose_)
-        {
-            std::cout << "The cumulative CPU time was " << cpuTime << " seconds.\n";
-        }
+            std::cout << Fmt::format("The cumulative CPU time was {:.2g} seconds.\n", cpuTime);
     }
 
     //! If the time loop has verbose output
@@ -524,8 +512,8 @@ public:
             lastPeriodicCheckPoint_ += interval;
 
         if (this->verbose())
-            std::cout << "Enabled periodic check points every " << interval
-                      << " seconds with the next check point at " << lastPeriodicCheckPoint_ + interval << " seconds." << std::endl;
+            std::cout << Fmt::format("Enabled periodic check points every {:.5g} seconds ", interval)
+                      << Fmt::format("with the next check point at {:.5g} seconds.\n", lastPeriodicCheckPoint_ + interval);
 
         // check if the current time point is a check point
         if (Dune::FloatCmp::eq(this->time()-lastPeriodicCheckPoint_, 0.0, 1e-7))
@@ -603,8 +591,8 @@ private:
         if (Dune::FloatCmp::le(t - this->time(), 0.0, this->timeStepSize()*1e-7))
         {
             if (this->verbose())
-                std::cerr << "Couldn't insert checkpoint at t = " << t
-                          << " because that's in the past! (current simulation time is " << this->time() << ")" << std::endl;
+                std::cerr << Fmt::format("Couldn't insert checkpoint at t = {:.5g} ", t)
+                          << Fmt::format("because that's in the past! (current simulation time is {:.5g})\n", this->time());
             return;
         }
 
@@ -621,7 +609,7 @@ private:
 
         checkPoints_.push(t);
         if (this->verbose())
-            std::cout << "Set check point at t = " << t << " seconds." << std::endl;
+            std::cout << Fmt::format("Set check point at t = {:.5g} seconds.\n", t);
     }
 
     /*!
diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh
index 84410c4c883ca7f6a0f481f9a82079ec46409dba..81d10683442c74c862c64d314ef5f3883da8dcb2 100644
--- a/dumux/io/vtkoutputmodule.hh
+++ b/dumux/io/vtkoutputmodule.hh
@@ -43,6 +43,7 @@
 
 #include <dumux/common/parameters.hh>
 #include <dumux/common/typetraits/typetraits.hh>
+#include <dumux/io/format.hh>
 #include <dumux/discretization/method.hh>
 
 #include "vtkfunction.hh"
@@ -86,7 +87,7 @@ public:
         writer_ = std::make_shared<Dune::VTKWriter<GridView>>(gridGeometry.gridView(), dm, coordPrecision);
         sequenceWriter_ = std::make_unique<Dune::VTKSequenceWriter<GridView>>(writer_, name);
     }
-    
+
     virtual ~VtkOutputModuleBase() = default;
 
     //! the parameter group for getting parameter from the parameter tree
@@ -182,9 +183,7 @@ public:
         //! output
         timer.stop();
         if (verbose_)
-        {
-            std::cout << "Writing output for problem \"" << name_ << "\". Took " << timer.elapsed() << " seconds." << std::endl;
-        }
+            std::cout << Fmt::format("Writing output for problem \"{}\". Took {:.2g} seconds.\n", name_, timer.elapsed());
     }
 
 protected:
diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index 0d2d242d6bae97d80b832f9f541c0c18db107bae..b33f113a46d5ad5f16f7c117120f02b620dfeda5 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -48,6 +48,7 @@
 #include <dumux/common/typetraits/isvalid.hh>
 #include <dumux/common/timeloop.hh>
 #include <dumux/common/pdesolver.hh>
+#include <dumux/io/format.hh>
 #include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
 #include <dumux/linear/matrixconverter.hh>
 #include <dumux/assembly/partialreassembler.hh>
@@ -311,9 +312,11 @@ public:
                 this->assembler().resetTimeStep(uCurrentIter);
 
                 if (verbosity_ >= 1)
-                    std::cout << "Newton solver did not converge with dt = "
-                              << timeLoop.timeStepSize() << " seconds. Retrying with time step of "
-                              << timeLoop.timeStepSize() * retryTimeStepReductionFactor_ << " seconds\n";
+                {
+                    const auto dt = timeLoop.timeStepSize();
+                    std::cout << Fmt::format("Newton solver did not converge with dt = {} seconds. ", dt)
+                              << Fmt::format("Retrying with time step of dt = {} seconds.\n", dt*retryTimeStepReductionFactor_);
+                }
 
                 // try again with dt = dt * retryTimeStepReductionFactor_
                 timeLoop.setTimeStepSize(timeLoop.timeStepSize() * retryTimeStepReductionFactor_);
@@ -321,9 +324,9 @@ public:
 
             else
             {
-                DUNE_THROW(NumericalProblem, "Newton solver didn't converge after "
-                                             << maxTimeStepDivisions_ << " time-step divisions. dt="
-                                             << timeLoop.timeStepSize() << '\n');
+                DUNE_THROW(NumericalProblem,
+                    Fmt::format("Newton solver didn't converge after {} time-step divisions; dt = {}.\n",
+                                maxTimeStepDivisions_, timeLoop.timeStepSize()));
             }
         }
     }
@@ -336,8 +339,8 @@ public:
     {
         const bool converged = solve_(uCurrentIter);
         if (!converged)
-            DUNE_THROW(NumericalProblem, "Newton solver didn't converge after "
-                                         << numSteps_ << " iterations.\n");
+            DUNE_THROW(NumericalProblem,
+                Fmt::format("Newton solver didn't converge after {} iterations.\n", numSteps_));
     }
 
     /*!
@@ -463,8 +466,7 @@ public:
                 ++numLinearSolverBreakdowns_;
             }
             else if (!convergedRemote) {
-                DUNE_THROW(NumericalProblem,
-                           "Linear solver did not converge on a remote process");
+                DUNE_THROW(NumericalProblem, "Linear solver did not converge on a remote process");
                 ++numLinearSolverBreakdowns_; // we keep correct count for process 0
             }
         }
@@ -574,22 +576,15 @@ public:
             if (enableDynamicOutput_)
                 std::cout << '\r'; // move cursor to beginning of line
 
-            auto width = std::to_string(maxSteps_).size();
-            std::cout << "Newton iteration " << std::setw(width) << numSteps_ << " done";
-
-            auto formatFlags = std::cout.flags();
-            auto prec = std::cout.precision();
-            std::cout << std::scientific << std::setprecision(3);
+            const auto width = Fmt::formatted_size("{}", maxSteps_);
+            std::cout << Fmt::format("Newton iteration {:{}} done", numSteps_, width);
 
             if (enableShiftCriterion_)
-                std::cout << ", maximum relative shift = " << shift_;
+                std::cout << Fmt::format(", maximum relative shift = {:.4e}", shift_);
             if (enableResidualCriterion_ && enableAbsoluteResidualCriterion_)
-                std::cout << ", residual = " << residualNorm_;
+                std::cout << Fmt::format(", residual = {:.4e}", residualNorm_);
             else if (enableResidualCriterion_)
-                std::cout << ", residual reduction = " << reduction_;
-
-            std::cout.flags(formatFlags);
-            std::cout.precision(prec);
+                std::cout << Fmt::format(", residual reduction = {:.4e}", reduction_);
 
             std::cout << endIterMsgStream_.str() << "\n";
         }
@@ -996,11 +991,10 @@ private:
 
             if (verbosity_ >= 1) {
                 const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
-                std::cout << "Assemble/solve/update time: "
-                          <<  assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
-                          <<  solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
-                          <<  updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
-                          << "\n";
+                std::cout << Fmt::format("Assemble/solve/update time: {:.2g}({:.2f}%)/{:.2g}({:.2f}%)/{:.2g}({:.2f}%)\n",
+                                         assembleTimer.elapsed(), 100*assembleTimer.elapsed()/elapsedTot,
+                                         solveTimer.elapsed(), 100*solveTimer.elapsed()/elapsedTot,
+                                         updateTimer.elapsed(), 100*updateTimer.elapsed()/elapsedTot);
             }
             return true;
 
@@ -1066,7 +1060,7 @@ private:
             computeResidualReduction_(uCurrentIter);
 
             if (reduction_ < lastReduction_ || lambda <= 0.125) {
-                endIterMsgStream_ << ", residual reduction " << lastReduction_ << "->"  << reduction_ << "@lambda=" << lambda;
+                endIterMsgStream_ << Fmt::format(", residual reduction {:.4e}->{:.4e}@lambda={:.4f}", lastReduction_, reduction_, lambda);
                 return;
             }