Newton error propagation based convergence check.

src/codegen/codegen_neuron_cpp_visitor.cpp

@@ -1269,7 +1269,7 @@
 
 void CodegenNeuronCppVisitor::print_mechanism_range_var_structure(bool print_initializers) {
     auto const value_initialize = print_initializers ? "{}" : "";
     auto int_type = default_int_data_type();
     printer->add_newline(2);
     printer->add_line("/** all mechanism instance variables and global variables */");
     printer->fmt_push_block("struct {} ", instance_struct());
@@ -1483,9 +1483,21 @@
 
     print_rename_state_vars();
 
+    if (!info.changed_dt.empty()) {
+        printer->fmt_line("double _save_prev_dt = {};",
+                          get_variable_name(naming::NTHREAD_DT_VARIABLE));
+        printer->fmt_line("{} = {};",
+                          get_variable_name(naming::NTHREAD_DT_VARIABLE),
+                          info.changed_dt);
+    }
+
     print_initial_block(info.initial_node);
-    printer->pop_block();
 
+    if (!info.changed_dt.empty()) {
+        printer->fmt_line("{} = _save_prev_dt;", get_variable_name(naming::NTHREAD_DT_VARIABLE));
+    }
+
+    printer->pop_block();
     printer->pop_block();
 }
 
@@ -1576,7 +1588,7 @@
             }
             const auto& var_name = var->get_name();
             auto var_pos = position_of_float_var(var_name);
             double var_value = var->get_value() == nullptr ? 0.0 : *var->get_value();
 
             printer->fmt_line("_lmc.template fpfield<{}>(_iml) = {}; /* {} */",
                               var_pos,

src/solver/newton/newton.hpp

@@ -34,8 +34,22 @@ namespace newton {
  * @{
  */
 
-static constexpr int MAX_ITER = 1e3;
-static constexpr double EPS = 1e-12;
+static constexpr int MAX_ITER = 50;
+static constexpr double EPS = 1e-13;
+
+template <int N>
+bool is_converged(const Eigen::Matrix<double, N, 1>& X,
+                  const Eigen::Matrix<double, N, N>& J,
+                  const Eigen::Matrix<double, N, 1>& F,
+                  double eps) {
+    for (Eigen::Index i = 0; i < N; ++i) {
+        double square_error = J(i, Eigen::all).cwiseAbs2() * (eps * X).cwiseAbs2();
+        if (F(i) * F(i) > square_error) {
+            return false;
+        }
+    }
+    return true;
+}
 
 /**
  * \brief Newton method with user-provided Jacobian
@@ -58,17 +72,14 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
                                     int max_iter = MAX_ITER) {
     // Vector to store result of function F(X):
     Eigen::Matrix<double, N, 1> F;
-    // Matrix to store jacobian of F(X):
+    // Matrix to store Jacobian of F(X):
     Eigen::Matrix<double, N, N> J;
     // Solver iteration count:
     int iter = -1;
     while (++iter < max_iter) {
         // calculate F, J from X using user-supplied functor
         functor(X, F, J);
-        // get error norm: here we use sqrt(|F|^2)
-        double error = F.norm();
-        if (error < eps) {
-            // we have converged: return iteration count
+        if (is_converged(X, J, F, eps)) {
             return iter;
         }
         // In Eigen the default storage order is ColMajor.
@@ -109,8 +120,7 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
     int iter = -1;
     while (++iter < max_iter) {
         functor(X, F, J);
-        double error = F.norm();
-        if (error < eps) {
+        if (is_converged(X, J, F, eps)) {
             return iter;
         }
         // The inverse can be called from within OpenACC regions without any issue, as opposed to

src/visitors/sympy_solver_visitor.cpp

@@ -543,15 +543,15 @@ void SympySolverVisitor::visit_derivative_block(ast::DerivativeBlock& node) {
                     pre_solve_statements.push_back(std::move(expression));
                 }
                 // replace ODE with Euler equation
-                eq = x;
+                eq = "(";
+                eq.append(x);
                 eq.append(x_array_index);
-                eq.append(" = ");
+                eq.append(" - ");
                 eq.append(old_x);
-                eq.append(" + ");
+                eq.append(") / ");
                 eq.append(codegen::naming::NTHREAD_DT_VARIABLE);
-                eq.append(" * (");
+                eq.append(" = ");
                 eq.append(dxdt);
-                eq.append(")");
                 logger->debug("SympySolverVisitor :: -> constructed Euler eq: {}", eq);
             }
         }