Add adaptive time stepping to elastic sims via enum in Variable

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/magnetoelastic_dispersion.py

@@ -81,10 +81,6 @@ def simulation(theta):
     magnet.alpha = 0
     magnet.eta = 0
 
-    # time stepping
-    world.timesolver.adaptive_timestep = False
-    world.timesolver.timestep = 1e-13
-
     # parameters to save
     m = np.zeros(shape=(nt, 3, nz, ny, nx))
     u = np.zeros(shape=(nt, 3, nz, ny, nx))

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/timesolver.py

@@ -159,23 +159,65 @@ def time(self, time):
         self._impl.time = time
 
     @property
-    def max_error(self):
-        """Return the maximum error per step the solver can tollerate.
+    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.
 
         See Also
         --------
         headroom, lower_bound, sensible_factor, upper_bound
+        displacement_max_error, velocity_max_error
         """
 
-        return self._impl.max_error
+        return self._impl.torque_max_error
 
-    @max_error.setter
-    def max_error(self, error):
+    @magnetization_max_error.setter
+    def magnetization_max_error(self, error):
         assert error > 0, "The maximum error should be bigger than 0."
-        self._impl.max_error = error
+        self._impl.magnetization_max_error = error
 
+    @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
+        --------
+        headroom, lower_bound, sensible_factor, upper_bound
+        magnetization_max_error, 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
+        --------
+        headroom, lower_bound, sensible_factor, upper_bound
+        displacement_max_error, magnetization_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
+
     @property
     def headroom(self):
         """Return the solver headroom.
@@ -184,7 +226,8 @@ def headroom(self):
 
         See Also
         --------
-        lower_bound, max_error, sensible_factor, upper_bound
+        lower_bound, sensible_factor, upper_bound
+        displacement_max_error, magnetization_max_error, velocity_max_error
         """
         return self._impl.headroom
 
@@ -202,7 +245,8 @@ def lower_bound(self):
 
         See Also
         --------
-        headroom, max_error, sensible_factor, upper_bound
+        headroom, sensible_factor, upper_bound
+        displacement_max_error, magnetization_max_error, velocity_max_error
         """
         return self._impl.lower_bound
 
@@ -220,7 +264,8 @@ def upper_bound(self):
 
         See Also
         --------
-        headroom, lower_bound, max_error, sensible_factor
+        headroom, lower_bound, sensible_factor
+        displacement_max_error, magnetization_max_error, velocity_max_error
         """
         return self._impl.upper_bound
 
@@ -238,7 +283,8 @@ def sensible_factor(self):
 
         See Also
         --------
-        headroom, lower_bound, max_error, upper_bound
+        headroom, lower_bound, upper_bound
+        displacement_max_error, magnetization_max_error, velocity_max_error
         """
         return self._impl.sensible_factor
 

src/bindings/wrap_timesolver.cpp

@@ -23,7 +23,12 @@ 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("magnetization_max_error", &TimeSolver::magnetizationMaxError,
+                                               &TimeSolver::setMagnetizationMaxError)
+      .def_property("displacement_max_error", &TimeSolver::displacementMaxError,
+                                              &TimeSolver::setDisplacementMaxError)
+      .def_property("velocity_max_error", &TimeSolver::velocityMaxError,
+                                          &TimeSolver::setVelocityMaxError)
       .def_property("sensible_factor", &TimeSolver::sensibleFactor, &TimeSolver::setSensibleFactor)
       .def_property("upper_bound", &TimeSolver::upperBound, &TimeSolver::setUpperBound)
       .def("step", &TimeSolver::step)

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 stepper.solver_->getMaxError(x.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.cpp

@@ -9,6 +9,7 @@
 #include "reduce.hpp"
 #include "rungekutta.hpp"
 #include "stepper.hpp"
+#include "variable.hpp"
 
 std::unique_ptr<TimeSolver> TimeSolver::Factory::create() {
   return std::unique_ptr<TimeSolver>(new TimeSolver());
@@ -34,6 +35,19 @@ RKmethod TimeSolver::getRungeKuttaMethod() {
   return method_;
 }
 
+real TimeSolver::getMaxError(MaxError maxError) const {
+  switch (maxError) {
+    case MaxError::MAGNETIZATION:
+      return this->magnetizationMaxError();
+    case MaxError::DISPLACEMENT:
+      return this->displacementMaxError();
+    case MaxError::VELOCITY:
+      return this->velocityMaxError();
+    default:
+      return this->magnetizationMaxError();
+  }
+}
+
 real TimeSolver::sensibleTimeStep() const {
   if (eqs_.empty())
     return 0.0;  // Timestep is irrelevant if there are no equations to solve

src/core/timesolver.hpp

@@ -10,6 +10,7 @@
 
 class Stepper;
 enum class RKmethod;
+enum class MaxError;
 
 class TimeSolver {
   //------------- CONSTRUCTORS -------------------------------------------------
@@ -32,7 +33,10 @@ class TimeSolver {
   RKmethod getRungeKuttaMethod();
   real headroom() const { return headroom_; }
   real lowerBound() const { return lowerBound_; }
-  real maxError() const { return maxError_; }
+  real magnetizationMaxError() const { return magnetizationMaxError_; }
+  real displacementMaxError() const { return displacementMaxError_; }
+  real velocityMaxError() const { return velocityMaxError_; }
+  real getMaxError(MaxError) const;  // get approprate max error, given enum of Variable
   real sensibleFactor() const { return sensibleFactor_; }
   real time() const { return time_; }
   real timestep() const { return timestep_; }
@@ -46,7 +50,9 @@ 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 setMagnetizationMaxError(real maxError) { magnetizationMaxError_ = maxError; }
+  void setDisplacementMaxError(real maxError) { displacementMaxError_ = maxError; }
+  void setVelocityMaxError(real maxError) { velocityMaxError_ = maxError; }
   void setSensibleFactor(real factor) { sensibleFactor_ = factor; }
   void setTime(real time) { time_ = time; }
   void setTimeStep(real dt) { timestep_ = dt; }
@@ -63,6 +69,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 +78,9 @@ class TimeSolver {
 
   real headroom_ = 0.8;
   real lowerBound_ = 0.5;
-  real maxError_ = 1e-5;
+  real magnetizationMaxError_ = 1e-5;
+  real displacementMaxError_ = 1e-18;  // TODO: may need to change default value
+  real velocityMaxError_ = 1e-7;  // TODO: may need to change default value
   real sensibleFactor_ = 0.01;
   real time_ = 0.0;
   real timestep_ = 0.0;

src/core/variable.cpp

@@ -9,8 +9,8 @@
 #include "system.hpp"
 
 Variable::Variable(std::shared_ptr<const System> system, int ncomp,
-                   std::string name, std::string unit)
-    : name_(name), unit_(unit) {
+                   std::string name, std::string unit, MaxError maxError)
+    : name_(name), unit_(unit), maxError_(maxError) {
   field_ = new Field(system, ncomp);
 }
 
@@ -34,6 +34,10 @@ std::string Variable::unit() const {
   return unit_;
 }
 
+MaxError Variable::maxError() const {
+  return maxError_;
+}
+
 Field Variable::eval() const {
   return field_->eval();
 }
@@ -70,8 +74,9 @@ void Variable::set(real3 value) const {
 NormalizedVariable::NormalizedVariable(std::shared_ptr<const System> system,
                                        int ncomp,
                                        std::string name,
-                                       std::string unit)
-    : Variable(system, ncomp, name, unit) {}
+                                       std::string unit,
+                                       MaxError maxError)
+    : Variable(system, ncomp, name, unit, maxError) {}
 
 void NormalizedVariable::set(const Field& src) const {
   // TODO: check if this is possible without the extra copy

src/core/variable.hpp

@@ -9,16 +9,24 @@
 class CuField;
 class System;
 
+enum class MaxError {
+  MAGNETIZATION,
+  DISPLACEMENT,
+  VELOCITY
+};
+
 class Variable : public FieldQuantity {
  public:
   Variable(std::shared_ptr<const System> system, int ncomp,
-           std::string name = "", std::string unit = "");
+           std::string name = "", std::string unit = "",
+           MaxError maxError = MaxError::MAGNETIZATION);
   ~Variable();
 
   int ncomp() const;
   std::shared_ptr<const System> system() const;
   std::string name() const;
   std::string unit() const;
+  MaxError maxError() const;
 
   Field eval() const;
 
@@ -39,14 +47,16 @@ class Variable : public FieldQuantity {
  private:
   std::string name_;
   std::string unit_;
+  MaxError maxError_;  // enum to select appropriate max error for adaptive time stepping
 };
 
 // Exactly the same as variable, but when values are set, the values are
 // normalized
 class NormalizedVariable : public Variable {
  public:
   NormalizedVariable(std::shared_ptr<const System> system, int ncomp,
-                     std::string name = "", std::string unit = "");
+                     std::string name = "", std::string unit = "",
+                     MaxError maxError = MaxError::MAGNETIZATION);
   void set(const Field&) const;
   void set(real) const;
   void set(real3) const;

src/physics/ferromagnet.cpp

@@ -20,7 +20,7 @@ Ferromagnet::Ferromagnet(std::shared_ptr<System> system_ptr,
                          Antiferromagnet* hostMagnet)
     : Magnet(system_ptr, name),
       hostMagnet_(hostMagnet),
-      magnetization_(system(), 3, name + ":magnetization", ""),
+      magnetization_(system(), 3, name + ":magnetization", "", MaxError::MAGNETIZATION),
       msat(system(), 1.0, name + ":msat", "A/m"),
       aex(system(), 0.0, name + ":aex", "J/m"),
       interExch(system(), 0.0, name + ":inter_exchange", "J/m"),

src/physics/magnet.cpp

@@ -187,10 +187,12 @@ void Magnet::setEnableElastodynamics(bool value) {
     if (value) {
       // properly initialize Variables now
       elasticDisplacement_ = std::make_unique<Variable>(system(), 3,
-                                        name() + ":elastic_displacement", "m");
+                                        name() + ":elastic_displacement", "m",
+                                        MaxError::DISPLACEMENT);
       elasticDisplacement_->set(real3{0,0,0});
       elasticVelocity_ = std::make_unique<Variable>(system(), 3,
-                                        name() + ":elastic_velocity", "m/s");
+                                        name() + ":elastic_velocity", "m/s",
+                                        MaxError::VELOCITY);
       elasticVelocity_->set(real3{0,0,0});
     } else {
       // free memory of unnecessary Variables