Implement ELECTRODE_CURRENT. (#1425)

src/codegen/codegen_neuron_cpp_visitor.cpp

@@ -1549,7 +1549,8 @@ void CodegenNeuronCppVisitor::print_nrn_jacob() {
 
     if (breakpoint_exist()) {
         printer->add_line("int node_id = node_data.nodeindices[id];");
-        printer->fmt_line("node_data.node_diagonal[node_id] += inst.{}[id];",
+        printer->fmt_line("node_data.node_diagonal[node_id] {} inst.{}[id];",
+                          operator_for_d(),
                           info.vectorize ? naming::CONDUCTANCE_UNUSED_VARIABLE
                                          : naming::CONDUCTANCE_VARIABLE);
     }
@@ -1933,7 +1934,7 @@ void CodegenNeuronCppVisitor::print_nrn_cur() {
     // }
 
 
-    printer->add_line("node_data.node_rhs[node_id] -= rhs;");
+    printer->fmt_line("node_data.node_rhs[node_id] {} rhs;", operator_for_rhs());
 
     if (breakpoint_exist()) {
         printer->fmt_line("inst.{}[id] = g;",

test/usecases/CMakeLists.txt

@@ -1,6 +1,7 @@
 set(NMODL_USECASE_DIRS
     solve
     constant
+    electrode_current
     function
     procedure
     global

test/usecases/electrode_current/leonhard.mod

@@ -0,0 +1,21 @@
+UNITS {
+    (mA) = (milliamp)
+}
+
+NEURON {
+    SUFFIX leonhard
+    ELECTRODE_CURRENT il
+    RANGE c
+}
+
+ASSIGNED {
+    il (mA/cm2)
+}
+
+PARAMETER {
+    c = 0.005
+}
+
+BREAKPOINT {
+    il = - c * (v - 1.5)
+}

test/usecases/electrode_current/simulate.py

@@ -0,0 +1,60 @@
+import numpy as np
+
+from neuron import gui, h
+from neuron.units import ms
+
+
+# Test the accuracy of the simulation output to the analytical result
+nseg = 1
+s = h.Section()
+s.insert("leonhard")
+s.nseg = nseg
+s.c_leonhard = 0.005
+
+v_hoc = h.Vector().record(s(0.5)._ref_v)
+t_hoc = h.Vector().record(h._ref_t)
+
+h.stdinit()
+h.tstop = 5.0 * ms
+h.run()
+
+v = np.array(v_hoc.as_numpy())
+t = np.array(t_hoc.as_numpy())
+
+erev = 1.5
+rate = s.c_leonhard / 1e-3
+v0 = -65.0
+v_exact = erev + (v0 - erev) * np.exp(-rate * t)
+rel_err = np.abs(v - v_exact) / np.max(np.abs(v_exact))
+
+assert np.all(rel_err < 1e-1), f"rel_err = {rel_err}"
+
+
+# Test the stability of the simulation at final time using a single timestep
+nseg = 1
+s = h.Section()
+s.insert("leonhard")
+s.nseg = nseg
+s.c_leonhard = 0.5
+
+v_hoc = h.Vector().record(s(0.5)._ref_v)
+t_hoc = h.Vector().record(h._ref_t)
+
+h.stdinit()
+h.tstop = 1000.0 * ms
+h.dt = 1000.0 * ms
+h.run()
+
+v = np.array(v_hoc.as_numpy())
+t = np.array(t_hoc.as_numpy())
+
+erev = 1.5
+rate = s.c_leonhard / 1e-3
+v0 = -65.0
+v_exact = erev + (v0 - erev) * np.exp(-rate * t)
+rel_err = np.abs(v - v_exact) / np.max(np.abs(v_exact))
+
+assert np.allclose(v[-1], v_exact[-1], atol=0.0), f"rel_err = {rel_err}"
+
+
+print("leonhard: success")