From 0de836a1b6e622a37300da72b788270c3a347152 Mon Sep 17 00:00:00 2001 From: Ian Lateur Date: Mon, 20 Jan 2025 16:29:36 +0100 Subject: [PATCH] Add adaptive time stepping to elastic sims via enum in Variable --- examples/magnetoelastic.py | 3 -- examples/magnetoelastic_dispersion.py | 4 -- examples/pure_elastic.py | 4 -- mumaxplus/timesolver.py | 66 +++++++++++++++++++++++---- src/bindings/wrap_timesolver.cpp | 7 ++- src/core/rungekutta.cpp | 18 ++++++-- src/core/rungekutta.hpp | 2 + src/core/timesolver.cpp | 14 ++++++ src/core/timesolver.hpp | 15 ++++-- src/core/variable.cpp | 13 ++++-- src/core/variable.hpp | 14 +++++- src/physics/ferromagnet.cpp | 2 +- src/physics/magnet.cpp | 6 ++- src/physics/relaxer.cpp | 12 ++--- 14 files changed, 135 insertions(+), 45 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/magnetoelastic_dispersion.py b/examples/magnetoelastic_dispersion.py index 5d4f6e0e..b3a081b7 100644 --- a/examples/magnetoelastic_dispersion.py +++ b/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)) 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/timesolver.py b/mumaxplus/timesolver.py index d083983a..011fa07e 100644 --- a/mumaxplus/timesolver.py +++ b/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 diff --git a/src/bindings/wrap_timesolver.cpp b/src/bindings/wrap_timesolver.cpp index 61527f4f..fc88a910 100644 --- a/src/bindings/wrap_timesolver.cpp +++ b/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) diff --git a/src/core/rungekutta.cpp b/src/core/rungekutta.cpp index 7ecbb7a1..dcd3d158 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 stepper.solver_->getMaxError(x.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.cpp b/src/core/timesolver.cpp index 603c7bc7..6ab61972 100644 --- a/src/core/timesolver.cpp +++ b/src/core/timesolver.cpp @@ -9,6 +9,7 @@ #include "reduce.hpp" #include "rungekutta.hpp" #include "stepper.hpp" +#include "variable.hpp" std::unique_ptr TimeSolver::Factory::create() { return std::unique_ptr(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 diff --git a/src/core/timesolver.hpp b/src/core/timesolver.hpp index 35dcc308..38cd9ce4 100644 --- a/src/core/timesolver.hpp +++ b/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 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; diff --git a/src/core/variable.cpp b/src/core/variable.cpp index bccfe5df..60641916 100644 --- a/src/core/variable.cpp +++ b/src/core/variable.cpp @@ -9,8 +9,8 @@ #include "system.hpp" Variable::Variable(std::shared_ptr 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 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 diff --git a/src/core/variable.hpp b/src/core/variable.hpp index fd1dc5d9..875bc685 100644 --- a/src/core/variable.hpp +++ b/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 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 system() const; std::string name() const; std::string unit() const; + MaxError maxError() const; Field eval() const; @@ -39,6 +47,7 @@ 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 @@ -46,7 +55,8 @@ class Variable : public FieldQuantity { class NormalizedVariable : public Variable { public: NormalizedVariable(std::shared_ptr 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; diff --git a/src/physics/ferromagnet.cpp b/src/physics/ferromagnet.cpp index 54ad2040..52a585b0 100644 --- a/src/physics/ferromagnet.cpp +++ b/src/physics/ferromagnet.cpp @@ -20,7 +20,7 @@ Ferromagnet::Ferromagnet(std::shared_ptr 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"), diff --git a/src/physics/magnet.cpp b/src/physics/magnet.cpp index 93a9cf59..8b52ee31 100644 --- a/src/physics/magnet.cpp +++ b/src/physics/magnet.cpp @@ -187,10 +187,12 @@ void Magnet::setEnableElastodynamics(bool value) { if (value) { // properly initialize Variables now elasticDisplacement_ = std::make_unique(system(), 3, - name() + ":elastic_displacement", "m"); + name() + ":elastic_displacement", "m", + MaxError::DISPLACEMENT); elasticDisplacement_->set(real3{0,0,0}); elasticVelocity_ = std::make_unique(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 diff --git a/src/physics/relaxer.cpp b/src/physics/relaxer.cpp index 99e6a55c..1df2f9ba 100644 --- a/src/physics/relaxer.cpp +++ b/src/physics/relaxer.cpp @@ -90,7 +90,7 @@ void Relaxer::exec() { real time = timesolver_.time(); real timestep = timesolver_.timestep(); bool adaptive = timesolver_.hasAdaptiveTimeStep(); - real maxerr = timesolver_.maxError(); + real maxerr = timesolver_.magnetizationMaxError(); std::string method = getRungeKuttaNameFromMethod(timesolver_.getRungeKuttaMethod()); auto eqs = timesolver_.equations(); @@ -120,11 +120,11 @@ void Relaxer::exec() { real t0 = 0; real t1 = calcTorque(torque); - real err = timesolver_.maxError(); + real err = timesolver_.magnetizationMaxError(); while (err > tol_) { err /= std::sqrt(2); - timesolver_.setMaxError(err); + timesolver_.setMagnetizationMaxError(err); timesolver_.steps(N); t0 = t1; @@ -144,7 +144,7 @@ void Relaxer::exec() { // If threshold is set by user: relax until torque is smaller than or equal to threshold. else { - real err = timesolver_.maxError(); + real err = timesolver_.magnetizationMaxError(); std::vector torque = getTorque(); while (err > tol_) { @@ -158,7 +158,7 @@ void Relaxer::exec() { if (torqueConverged) { err /= std::sqrt(2); - timesolver_.setMaxError(err); + timesolver_.setMagnetizationMaxError(err); } timesolver_.steps(N); } @@ -166,7 +166,7 @@ void Relaxer::exec() { // Restore solver settings after relaxing timesolver_.setRungeKuttaMethod(method); - timesolver_.setMaxError(maxerr); + timesolver_.setMagnetizationMaxError(maxerr); if (!adaptive) { timesolver_.disableAdaptiveTimeStep(); } timesolver_.setTime(time); timesolver_.setTimeStep(timestep);