Add proper adaptive timestepping to elastic simulations

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

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,

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

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

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
 

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)

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);

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)

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
+}

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;
 };

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();
+}

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;

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;

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"),

src/physics/ferromagnet.hpp

@@ -87,6 +87,7 @@ class Ferromagnet : public Magnet {
   Parameter conductivity;
   Parameter amrRatio;
   real RelaxTorqueThreshold;
+  real magnetizationMaxError;
 
   curandGenerator_t randomGenerator;
 

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)

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;