From 160e15915edc52d16c5de3dbd45a59dfc6dfec5b Mon Sep 17 00:00:00 2001
From: Ian Lateur <ian.lateur@gmail.com>
Date: Wed, 15 Jan 2025 12:02:22 +0100
Subject: [PATCH] Add proper adaptive timestepping to elastic simulations

---
 examples/magnetoelastic.py        |  3 ---
 examples/pure_elastic.py          |  4 ----
 mumaxplus/ferromagnet.py          | 15 ++++++++++++
 mumaxplus/magnet.py               | 38 +++++++++++++++++++++++++++++++
 mumaxplus/timesolver.py           | 26 ++++-----------------
 src/bindings/wrap_ferromagnet.cpp |  1 +
 src/bindings/wrap_magnet.cpp      |  3 +++
 src/bindings/wrap_timesolver.cpp  |  1 -
 src/core/dynamicequation.cpp      | 15 ++++++++++--
 src/core/dynamicequation.hpp      | 15 +++++++++++-
 src/core/rungekutta.cpp           | 18 +++++++++++----
 src/core/rungekutta.hpp           |  2 ++
 src/core/timesolver.hpp           |  4 +---
 src/physics/ferromagnet.cpp       |  1 +
 src/physics/ferromagnet.hpp       |  1 +
 src/physics/magnet.cpp            |  5 +++-
 src/physics/magnet.hpp            |  3 +++
 src/physics/mumaxworld.cpp        | 12 ++++++----
 src/physics/relaxer.cpp           | 27 ++++++++++------------
 src/physics/relaxer.hpp           |  1 +
 20 files changed, 134 insertions(+), 61 deletions(-)

diff --git a/examples/magnetoelastic.py b/examples/magnetoelastic.py
index f4796d10..7244b29f 100644
--- a/examples/magnetoelastic.py
+++ b/examples/magnetoelastic.py
@@ -72,9 +72,6 @@ def displacement_to_scatter_data(magnet, scale, skip):
 u_scale = 5e4  # amplification of displacement
 u_skip  = 5  # don't show every displacement
 
-world.timesolver.adaptive_timestep = False
-world.timesolver.timestep = 1e-12
-
 steps = 400
 time_max = 0.8e-9
 duration = time_max/steps
diff --git a/examples/pure_elastic.py b/examples/pure_elastic.py
index 10623210..6289db5e 100644
--- a/examples/pure_elastic.py
+++ b/examples/pure_elastic.py
@@ -74,10 +74,6 @@ def displacement_to_scatter_data(magnet, scale, skip):
     u[i, ...] -= u_avg[i]
 magnet.elastic_displacement = u
 
-# adaptive time stepping does not work for magnetoelastics
-world.timesolver.adaptive_timestep = False
-world.timesolver.timestep = 1e-13
-
 # simulation
 time = np.linspace(0, 1e-10, 500)
 energies = {"E_kin": lambda : magnet.kinetic_energy.eval()*J_to_eV,
diff --git a/mumaxplus/ferromagnet.py b/mumaxplus/ferromagnet.py
index 6e1eb8cd..00edb84a 100644
--- a/mumaxplus/ferromagnet.py
+++ b/mumaxplus/ferromagnet.py
@@ -208,6 +208,21 @@ def RelaxTorqueThreshold(self, value):
         assert value != 0, "The relax threshold should not be zero."
         self._impl.RelaxTorqueThreshold = value
 
+    @property
+    def magnetization_max_error(self):
+        """Return the maximum error per step the solver can tollerate for the
+        magnetization-torque equations of motion (rad).
+        
+        The default value is 1e-5.
+        """
+
+        return self._impl.torque_max_error
+
+    @magnetization_max_error.setter
+    def magnetization_max_error(self, error):
+        assert error > 0, "The maximum error should be bigger than 0."
+        self._impl.magnetization_max_error = error
+
     # ----- MATERIAL PARAMETERS -----------
 
     @property
diff --git a/mumaxplus/magnet.py b/mumaxplus/magnet.py
index e7e6645d..2f74bfba 100644
--- a/mumaxplus/magnet.py
+++ b/mumaxplus/magnet.py
@@ -303,6 +303,44 @@ def rho(self):
     def rho(self, value):
         self.rho.set(value)
 
+    @property
+    def displacement_max_error(self):
+        """Return the maximum error per step the solver can tollerate for the
+        displacement-velocity equations of motion (m).
+        
+        The default value is 1e-18.
+
+        See Also
+        --------
+        velocity_max_error
+        """
+
+        return self._impl.displacement_max_error
+
+    @displacement_max_error.setter
+    def displacement_max_error(self, error):
+        assert error > 0, "The maximum error should be bigger than 0."
+        self._impl.displacement_max_error = error
+
+    @property
+    def velocity_max_error(self):
+        """Return the maximum error per step the solver can tollerate for the
+        velocity-acceleration equations of motion (m/s).
+        
+        The default value is 1e-7.
+
+        See Also
+        --------
+        displacement_max_error
+        """
+
+        return self._impl.velocity_max_error
+
+    @velocity_max_error.setter
+    def velocity_max_error(self, error):
+        assert error > 0, "The maximum error should be bigger than 0."
+        self._impl.velocity_max_error = error
+
     # ----- ELASTIC QUANTITIES -------
 
     @property
diff --git a/mumaxplus/timesolver.py b/mumaxplus/timesolver.py
index d083983a..1724c9b4 100644
--- a/mumaxplus/timesolver.py
+++ b/mumaxplus/timesolver.py
@@ -158,24 +158,6 @@ def time(self):
     def time(self, time):
         self._impl.time = time
 
-    @property
-    def max_error(self):
-        """Return the maximum error per step the solver can tollerate.
-        
-        The default value is 1e-5.
-
-        See Also
-        --------
-        headroom, lower_bound, sensible_factor, upper_bound
-        """
-
-        return self._impl.max_error
-
-    @max_error.setter
-    def max_error(self, error):
-        assert error > 0, "The maximum error should be bigger than 0."
-        self._impl.max_error = error
-    
     @property
     def headroom(self):
         """Return the solver headroom.
@@ -184,7 +166,7 @@ def headroom(self):
 
         See Also
         --------
-        lower_bound, max_error, sensible_factor, upper_bound
+        lower_bound, sensible_factor, upper_bound
         """
         return self._impl.headroom
 
@@ -202,7 +184,7 @@ def lower_bound(self):
 
         See Also
         --------
-        headroom, max_error, sensible_factor, upper_bound
+        headroom, sensible_factor, upper_bound
         """
         return self._impl.lower_bound
 
@@ -220,7 +202,7 @@ def upper_bound(self):
 
         See Also
         --------
-        headroom, lower_bound, max_error, sensible_factor
+        headroom, lower_bound, sensible_factor
         """
         return self._impl.upper_bound
 
@@ -238,7 +220,7 @@ def sensible_factor(self):
 
         See Also
         --------
-        headroom, lower_bound, max_error, upper_bound
+        headroom, lower_bound, upper_bound
         """
         return self._impl.sensible_factor
 
diff --git a/src/bindings/wrap_ferromagnet.cpp b/src/bindings/wrap_ferromagnet.cpp
index da332410..dcccba85 100644
--- a/src/bindings/wrap_ferromagnet.cpp
+++ b/src/bindings/wrap_ferromagnet.cpp
@@ -63,6 +63,7 @@ void wrap_ferromagnet(py::module& m) {
       .def_readonly("conductivity", &Ferromagnet::conductivity)
       .def_readonly("amr_ratio", &Ferromagnet::amrRatio)
       .def_readwrite("RelaxTorqueThreshold", &Ferromagnet::RelaxTorqueThreshold)
+      .def_readwrite("magnetization_max_error", &Ferromagnet::magnetizationMaxError)
       .def_readonly("poisson_system", &Ferromagnet::poissonSystem)
       .def_readonly("B1", &Ferromagnet::B1)
       .def_readonly("B2", &Ferromagnet::B2)
diff --git a/src/bindings/wrap_magnet.cpp b/src/bindings/wrap_magnet.cpp
index 6d21483b..117c2363 100644
--- a/src/bindings/wrap_magnet.cpp
+++ b/src/bindings/wrap_magnet.cpp
@@ -40,6 +40,9 @@ void wrap_magnet(py::module& m) {
       .def_readonly("eta", &Magnet::eta)
       .def_readonly("rho", &Magnet::rho)
 
+      .def_readwrite("displacement_max_error", &Magnet::displacementMaxError)
+      .def_readwrite("velocity_max_error", &Magnet::velocityMaxError)
+
       .def("stray_field_from_magnet",
           [](const Magnet* m, Magnet* magnet) {
             const StrayField* strayField = m->getStrayField(magnet);
diff --git a/src/bindings/wrap_timesolver.cpp b/src/bindings/wrap_timesolver.cpp
index 61527f4f..6f5e625b 100644
--- a/src/bindings/wrap_timesolver.cpp
+++ b/src/bindings/wrap_timesolver.cpp
@@ -23,7 +23,6 @@ void wrap_timesolver(py::module& m) {
            })
       .def_property("headroom", &TimeSolver::headroom, &TimeSolver::setHeadroom)
       .def_property("lower_bound", &TimeSolver::lowerBound, &TimeSolver::setLowerBound)
-      .def_property("max_error", &TimeSolver::maxError, &TimeSolver::setMaxError)
       .def_property("sensible_factor", &TimeSolver::sensibleFactor, &TimeSolver::setSensibleFactor)
       .def_property("upper_bound", &TimeSolver::upperBound, &TimeSolver::setUpperBound)
       .def("step", &TimeSolver::step)
diff --git a/src/core/dynamicequation.cpp b/src/core/dynamicequation.cpp
index 992cbbb3..409b7ec5 100644
--- a/src/core/dynamicequation.cpp
+++ b/src/core/dynamicequation.cpp
@@ -9,8 +9,9 @@
 
 DynamicEquation::DynamicEquation(const Variable* x,
                                  std::shared_ptr<FieldQuantity> rhs,
-                                 std::shared_ptr<FieldQuantity> noiseTerm)
-    : x(x), rhs(rhs), noiseTerm(noiseTerm) {
+                                 std::shared_ptr<FieldQuantity> noiseTerm,
+                                 const real* maxError)
+    : x(x), rhs(rhs), noiseTerm(noiseTerm), maxError_(maxError) {
   if (x->system() != rhs->system()) {
     throw std::runtime_error(
         "The variable and the r.h.s. of a dynamic equation should have the "
@@ -36,6 +37,11 @@ DynamicEquation::DynamicEquation(const Variable* x,
   }
 }
 
+DynamicEquation::DynamicEquation(const Variable* x,
+                                 std::shared_ptr<FieldQuantity> rhs,
+                                 const real* maxError)
+    : DynamicEquation(x, rhs, nullptr, maxError) {}
+
 int DynamicEquation::ncomp() const {
   return x->ncomp();
 }
@@ -47,3 +53,8 @@ Grid DynamicEquation::grid() const {
 std::shared_ptr<const System> DynamicEquation::system() const {
   return x->system();
 }
+
+real DynamicEquation::maxError() const {
+  if (maxError_) { return *maxError_; }
+  return 1e-5;  // default value, usually good for magnetic systems
+}
diff --git a/src/core/dynamicequation.hpp b/src/core/dynamicequation.hpp
index 0b08d268..69ed88ce 100644
--- a/src/core/dynamicequation.hpp
+++ b/src/core/dynamicequation.hpp
@@ -2,16 +2,28 @@
 
 #include <memory>
 
+#include "datatypes.hpp"
+
 class Variable;
 class FieldQuantity;
 class Grid;
 class System;
 
 class DynamicEquation {
+ private:
+  /** Max error for adaptive time stepping, pointing to set value of a magnet.
+   *  If not set, the default is 1e-5.
+   */
+  const real* maxError_;
+
  public:
   DynamicEquation(const Variable* x,
                   std::shared_ptr<FieldQuantity> rhs,
-                  std::shared_ptr<FieldQuantity> noiseTerm = nullptr);
+                  std::shared_ptr<FieldQuantity> noiseTerm = nullptr,
+                  const real* maxError = nullptr);
+  DynamicEquation(const Variable* x,
+                  std::shared_ptr<FieldQuantity> rhs,
+                  const real* maxError);  // option without thermal noise
 
   const Variable* x;
   std::shared_ptr<FieldQuantity> rhs;
@@ -20,4 +32,5 @@ class DynamicEquation {
   int ncomp() const;
   Grid grid() const;
   std::shared_ptr<const System> system() const;
+  real maxError() const;
 };
diff --git a/src/core/rungekutta.cpp b/src/core/rungekutta.cpp
index 7ecbb7a1..d7f850f3 100644
--- a/src/core/rungekutta.cpp
+++ b/src/core/rungekutta.cpp
@@ -57,20 +57,20 @@ void RungeKuttaStepper::step() {
     if (!solver_->hasAdaptiveTimeStep())
       break;
 
-    // loop over equations and get the largest error
+    // loop over equations and get the largest scaled error
     real error = 0.0;
     for (auto& eq : equations)
-      if (real e = eq.getError(); e > error)
+      if (real e = eq.getScaledError(); e > error)
         error = e;
 
-    success = error < solver_->maxError();
+    success = error < 1.;
 
     // update the timestep
     real corrFactor;
     if (success) {
-      corrFactor = std::pow(solver_->maxError() / error, 1. / (butcher_.order2 + 1));
+      corrFactor = std::pow(1. / error, 1. / (butcher_.order2 + 1));
     } else {
-      corrFactor = std::pow(solver_->maxError() / error, 1. / (butcher_.order1 + 1));
+      corrFactor = std::pow(1. / error, 1. / (butcher_.order1 + 1));
     }
     solver_->adaptTimeStep(corrFactor);
 
@@ -149,3 +149,11 @@ real RungeKuttaStepper::RungeKuttaStageExecutor::getError() const {
 
   return maxVecNorm(err);
 }
+
+real RungeKuttaStepper::RungeKuttaStageExecutor::maxError() const {
+  return eq_.maxError();
+}
+
+real RungeKuttaStepper::RungeKuttaStageExecutor::getScaledError() const {
+  return getError() / maxError();
+}
diff --git a/src/core/rungekutta.hpp b/src/core/rungekutta.hpp
index 36006cfa..8b59ad05 100644
--- a/src/core/rungekutta.hpp
+++ b/src/core/rungekutta.hpp
@@ -31,6 +31,8 @@ class RungeKuttaStepper::RungeKuttaStageExecutor {
   void setFinalX();
   void resetX();
   real getError() const;
+  real maxError() const;
+  real getScaledError() const;
 
  private:
   Field x0;
diff --git a/src/core/timesolver.hpp b/src/core/timesolver.hpp
index 35dcc308..3074058c 100644
--- a/src/core/timesolver.hpp
+++ b/src/core/timesolver.hpp
@@ -32,7 +32,6 @@ class TimeSolver {
   RKmethod getRungeKuttaMethod();
   real headroom() const { return headroom_; }
   real lowerBound() const { return lowerBound_; }
-  real maxError() const { return maxError_; }
   real sensibleFactor() const { return sensibleFactor_; }
   real time() const { return time_; }
   real timestep() const { return timestep_; }
@@ -46,7 +45,6 @@ class TimeSolver {
   void setEquations(std::vector<DynamicEquation> eq);
   void setHeadroom(real headroom) { headroom_ = headroom; }
   void setLowerBound(real lowerBound) { lowerBound_ = lowerBound; }
-  void setMaxError(real maxError) { maxError_ = maxError; }
   void setSensibleFactor(real factor) { sensibleFactor_ = factor; }
   void setTime(real time) { time_ = time; }
   void setTimeStep(real dt) { timestep_ = dt; }
@@ -63,6 +61,7 @@ class TimeSolver {
 
   //------------- HELPER FUNCTIONS FOR ADAPTIVE TIMESTEPPING -------------------
 
+  // TODO: take a look at sensible timestep
   real sensibleTimeStep() const; /** Computes a sensible timestep */
   void adaptTimeStep(real corr);
 
@@ -71,7 +70,6 @@ class TimeSolver {
 
   real headroom_ = 0.8;
   real lowerBound_ = 0.5;
-  real maxError_ = 1e-5;
   real sensibleFactor_ = 0.01;
   real time_ = 0.0;
   real timestep_ = 0.0;
diff --git a/src/physics/ferromagnet.cpp b/src/physics/ferromagnet.cpp
index 54ad2040..cd1d13c6 100644
--- a/src/physics/ferromagnet.cpp
+++ b/src/physics/ferromagnet.cpp
@@ -54,6 +54,7 @@ Ferromagnet::Ferromagnet(std::shared_ptr<System> system_ptr,
       conductivity(system(), 0.0, name + ":conductivity", "S/m"),
       amrRatio(system(), 0.0, name + ":amr_ratio", ""),
       RelaxTorqueThreshold(-1.0),
+      magnetizationMaxError(1e-5),
       poissonSystem(this), 
       // magnetoelasticity
       B1(system(), 0.0, name + ":B1", "J/m3"),
diff --git a/src/physics/ferromagnet.hpp b/src/physics/ferromagnet.hpp
index b7c85ee4..363cda33 100644
--- a/src/physics/ferromagnet.hpp
+++ b/src/physics/ferromagnet.hpp
@@ -87,6 +87,7 @@ class Ferromagnet : public Magnet {
   Parameter conductivity;
   Parameter amrRatio;
   real RelaxTorqueThreshold;
+  real magnetizationMaxError;
   
   curandGenerator_t randomGenerator;
 
diff --git a/src/physics/magnet.cpp b/src/physics/magnet.cpp
index 93a9cf59..d2991286 100644
--- a/src/physics/magnet.cpp
+++ b/src/physics/magnet.cpp
@@ -29,7 +29,10 @@ Magnet::Magnet(std::shared_ptr<System> system_ptr,
       C12(system(), 0.0, name + ":C12", "N/m2"),
       C44(system(), 0.0, name + ":C44", "N/m2"),
       eta(system(), 0.0, name + ":eta", "kg/m3s"),
-      rho(system(), 1.0, name + ":rho", "kg/m3") {
+      rho(system(), 1.0, name + ":rho", "kg/m3"),
+      // TODO: may need to change default values
+      displacementMaxError(1e-18),
+      velocityMaxError(1e-7) {
   // Check that the system has at least size 1
   int3 size = system_->grid().size();
   if (size.x < 1 || size.y < 1 || size.z < 1)
diff --git a/src/physics/magnet.hpp b/src/physics/magnet.hpp
index ba08cbb2..d31febcd 100644
--- a/src/physics/magnet.hpp
+++ b/src/physics/magnet.hpp
@@ -82,6 +82,9 @@ class Magnet {
   Parameter eta;  // Phenomenological elastic damping constant
   Parameter rho;  // Mass density
 
+  real displacementMaxError;
+  real velocityMaxError;
+
 
   // Delete copy constructor and copy assignment operator to prevent shallow copies
   Magnet(const Magnet&) = delete;
diff --git a/src/physics/mumaxworld.cpp b/src/physics/mumaxworld.cpp
index c8226500..5cf3ea50 100644
--- a/src/physics/mumaxworld.cpp
+++ b/src/physics/mumaxworld.cpp
@@ -162,7 +162,8 @@ void MumaxWorld::resetTimeSolverEquations(FM_Field torque) const {
     DynamicEquation eq(
         magnet->magnetization(),
         std::shared_ptr<FieldQuantity>(torque(magnet).clone()),
-        std::shared_ptr<FieldQuantity>(thermalNoiseQuantity(magnet).clone()));
+        std::shared_ptr<FieldQuantity>(thermalNoiseQuantity(magnet).clone()),
+        &magnet->magnetizationMaxError);
     equations.push_back(eq);
   }
 
@@ -172,7 +173,8 @@ void MumaxWorld::resetTimeSolverEquations(FM_Field torque) const {
       DynamicEquation eq(
         sub->magnetization(),
         std::shared_ptr<FieldQuantity>(torque(sub).clone()),
-        std::shared_ptr<FieldQuantity>(thermalNoiseQuantity(sub).clone()));
+        std::shared_ptr<FieldQuantity>(thermalNoiseQuantity(sub).clone()),
+        &sub->magnetizationMaxError);
       equations.push_back(eq);
     }
   }
@@ -187,15 +189,17 @@ void MumaxWorld::resetTimeSolverEquations(FM_Field torque) const {
       // change in displacement = velocity
       DynamicEquation dvEq(
           magnet->elasticDisplacement(),
-          std::shared_ptr<FieldQuantity>(elasticVelocityQuantity(magnet).clone()));
+          std::shared_ptr<FieldQuantity>(elasticVelocityQuantity(magnet).clone()),
           // No thermal noise
+          &magnet->displacementMaxError);
       equations.push_back(dvEq);
 
       // change in velocity = acceleration
       DynamicEquation vaEq(
           magnet->elasticVelocity(),
-          std::shared_ptr<FieldQuantity>(elasticAccelerationQuantity(magnet).clone()));
+          std::shared_ptr<FieldQuantity>(elasticAccelerationQuantity(magnet).clone()),
           // No thermal noise
+          &magnet->velocityMaxError);
       equations.push_back(vaEq);
     }
   }
diff --git a/src/physics/relaxer.cpp b/src/physics/relaxer.cpp
index 99e6a55c..b3be3ca8 100644
--- a/src/physics/relaxer.cpp
+++ b/src/physics/relaxer.cpp
@@ -18,12 +18,14 @@ Relaxer::Relaxer(const Magnet* magnet, std::vector<real> RelaxTorqueThreshold, r
       timesolver_(magnet->world()->timesolver()),
       world_(magnet->mumaxWorld()),
       tol_(tol),
-      threshold_(RelaxTorqueThreshold) {}
+      threshold_(RelaxTorqueThreshold),
+      maxerr_(1e-5) {}  // TODO: make user assignable?
 
 Relaxer::Relaxer(const MumaxWorld* world, real RelaxTorqueThreshold, real tol)
     : timesolver_(world->timesolver()),
       world_(world),
-      tol_(tol) {
+      tol_(tol),
+      maxerr_(1e-5) {
         for (const auto& pair : world->magnets()) {
           magnets_.push_back(pair.second);
           threshold_.push_back(RelaxTorqueThreshold);
@@ -39,7 +41,8 @@ std::vector<DynamicEquation> Relaxer::getEquation(const Magnet* magnet) {
       DynamicEquation eq(
           mag->magnetization(),
           std::shared_ptr<FieldQuantity>(relaxTorqueQuantity(mag).clone()),
-          std::shared_ptr<FieldQuantity>(thermalNoiseQuantity(mag).clone()));
+          std::shared_ptr<FieldQuantity>(thermalNoiseQuantity(mag).clone()),
+          &maxerr_);  // adaptive equations should listen to maxerr_ of Relaxer
       eqs.push_back(eq);
     }
     else if (const Antiferromagnet* mag = magnet->asAFM()) {
@@ -47,7 +50,8 @@ std::vector<DynamicEquation> Relaxer::getEquation(const Magnet* magnet) {
         DynamicEquation eq(
           sub->magnetization(),
           std::shared_ptr<FieldQuantity>(relaxTorqueQuantity(sub).clone()),
-          std::shared_ptr<FieldQuantity>(thermalNoiseQuantity(sub).clone()));
+          std::shared_ptr<FieldQuantity>(thermalNoiseQuantity(sub).clone()),
+          &maxerr_);
         eqs.push_back(eq);
       }
     }
@@ -90,7 +94,6 @@ void Relaxer::exec() {
   real time = timesolver_.time();
   real timestep = timesolver_.timestep();
   bool adaptive = timesolver_.hasAdaptiveTimeStep();
-  real maxerr = timesolver_.maxError();
   std::string method = getRungeKuttaNameFromMethod(timesolver_.getRungeKuttaMethod());
   auto eqs = timesolver_.equations();
 
@@ -119,12 +122,9 @@ void Relaxer::exec() {
     std::vector<FM_FieldQuantity> torque = getTorque();
     real t0 = 0;
     real t1 = calcTorque(torque);
-  
-    real err = timesolver_.maxError();
 
-    while (err > tol_) {
-      err /= std::sqrt(2);
-      timesolver_.setMaxError(err);
+    while (maxerr_ > tol_) {
+      maxerr_ /= std::sqrt(2);  // reduce maxError to which the new DynamicEquations listen
 
       timesolver_.steps(N);
       t0 = t1;
@@ -144,10 +144,9 @@ void Relaxer::exec() {
   // If threshold is set by user: relax until torque is smaller than or equal to threshold.
   else {
 
-    real err = timesolver_.maxError();
     std::vector<FM_FieldQuantity> torque = getTorque();
 
-    while (err > tol_) {
+    while (maxerr_ > tol_) {
       bool torqueConverged = true;
       for (size_t i = 0; i < torque.size(); i++) {
         if (maxVecNorm(torque[i].eval()) > threshold_[i]) {
@@ -157,8 +156,7 @@ void Relaxer::exec() {
       }
 
       if (torqueConverged) {    
-        err /= std::sqrt(2);
-        timesolver_.setMaxError(err);
+        maxerr_ /= std::sqrt(2);
       }
       timesolver_.steps(N);
     }
@@ -166,7 +164,6 @@ void Relaxer::exec() {
 
   // Restore solver settings after relaxing
   timesolver_.setRungeKuttaMethod(method);
-  timesolver_.setMaxError(maxerr);
   if (!adaptive) { timesolver_.disableAdaptiveTimeStep(); }
   timesolver_.setTime(time);
   timesolver_.setTimeStep(timestep); 
diff --git a/src/physics/relaxer.hpp b/src/physics/relaxer.hpp
index ad6dcf47..875c7753 100644
--- a/src/physics/relaxer.hpp
+++ b/src/physics/relaxer.hpp
@@ -31,4 +31,5 @@ class Relaxer {
   TimeSolver &timesolver_;
   const MumaxWorld* world_;
   real tol_;
+  real maxerr_;
 };
\ No newline at end of file