Revert "Newton error propagation based convergence check. (#1403)" (#…

…1417) This reverts commit 7984a45.
BlueBrain · Aug 30, 2024 · d244105 · d244105
1 parent 3df436c
commit d244105
Show file tree

Hide file tree

Showing 7 changed files with 130 additions and 282 deletions.
diff --git a/src/codegen/codegen_neuron_cpp_visitor.cpp b/src/codegen/codegen_neuron_cpp_visitor.cpp
@@ -1475,21 +1475,9 @@ void CodegenNeuronCppVisitor::print_nrn_init(bool skip_init_check) {
 
     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);
-
-    if (!info.changed_dt.empty()) {
-        printer->fmt_line("{} = _save_prev_dt;", get_variable_name(naming::NTHREAD_DT_VARIABLE));
-    }
-
     printer->pop_block();
+
     printer->pop_block();
 }
 

diff --git a/src/solver/newton/newton.hpp b/src/solver/newton/newton.hpp
@@ -34,22 +34,8 @@ namespace newton {
  * @{
  */
 
-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;
-}
+static constexpr int MAX_ITER = 1e3;
+static constexpr double EPS = 1e-12;
 
 /**
  * \brief Newton method with user-provided Jacobian
@@ -72,14 +58,17 @@ 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);
-        if (is_converged(X, J, F, eps)) {
+        // get error norm: here we use sqrt(|F|^2)
+        double error = F.norm();
+        if (error < eps) {
+            // we have converged: return iteration count
             return iter;
         }
         // In Eigen the default storage order is ColMajor.
@@ -120,7 +109,8 @@ 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);
-        if (is_converged(X, J, F, eps)) {
+        double error = F.norm();
+        if (error < eps) {
             return iter;
         }
         // The inverse can be called from within OpenACC regions without any issue, as opposed to

diff --git a/src/visitors/sympy_solver_visitor.cpp b/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 = "(";
-                eq.append(x);
+                eq = 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);
             }
         }