Implement diam and area.

src/codegen/codegen_neuron_cpp_visitor.cpp

@@ -1113,11 +1113,6 @@
         printer->add_newline();
     }
 
-    if (info.diam_used) {
-        printer->add_line("_morphology_sym = hoc_lookup(\"morphology\");");
-        printer->add_newline();
-    }
-
     for (const auto& ion: info.ions) {
         printer->fmt_line("_{0}_sym = hoc_lookup(\"{0}_ion\");", ion.name);
     }
@@ -1616,15 +1611,24 @@
         printer->add_line("_nrn_mechanism_access_dparam(_prop) = _ppvar;");
     }
 
-    if (info.diam_used) {
-        throw std::runtime_error("Diam allocation not implemented.");
-    }
+    const auto codegen_int_variables_size = codegen_int_variables.size();
 
-    if (info.area_used) {
-        throw std::runtime_error("Area allocation not implemented.");
-    }
+    if (info.diam_used || info.area_used) {
+        for (size_t i = 0; i < codegen_int_variables.size(); ++i) {
+            auto var_info = codegen_int_variables[i];
+            if (var_info.symbol->get_name() == "diam") {
+                printer->add_line("Symbol * morphology_sym = hoc_lookup(\"morphology\");");
+                printer->fmt_line("Prop * morphology_prop = need_memb(morphology_sym);");
 
-    const auto codegen_int_variables_size = codegen_int_variables.size();
+                printer->fmt_line(
+                    "_ppvar[{}] = _nrn_mechanism_get_param_handle(morphology_prop, 0);", i);
+            }
+            if (var_info.symbol->get_name() == "area") {
+                printer->fmt_line("_ppvar[{}] = _nrn_mechanism_get_area_handle(nrn_alloc_node_);",
+                                  i);
+            }
+        }
+    }
 
     for (const auto& ion: info.ions) {
         printer->fmt_line("Symbol * {}_sym = hoc_lookup(\"{}_ion\");", ion.name, ion.name);
@@ -1937,6 +1941,8 @@
 void CodegenNeuronCppVisitor::print_macro_definitions() {
     print_global_macros();
     print_mechanism_variables_macros();
+
+    printer->add_line("extern Node* nrn_alloc_node_;");
 }
 
 

test/usecases/CMakeLists.txt

@@ -5,6 +5,7 @@ set(NMODL_USECASE_DIRS
     global
     hodgkin_huxley
     kinetic
+    morphology
     nonspecific_current
     neuron_variables
     net_event

test/usecases/morphology/test_morphology.py

@@ -0,0 +1,60 @@
+import numpy as np
+from neuron import gui
+from neuron import h
+
+# This checks the use of `diam` and `area` in MOD files.
+
+
+def check_morphology(sections):
+    d2_expected = np.array([s(0.5).diam ** 2 for s in sections])
+    a2_expected = np.array([s(0.5).area() ** 2 for s in sections])
+
+    d2_actual = np.array([s(0.5).two_radii.square_diam() for s in sections])
+    a2_actual = np.array([s(0.5).two_radii.square_area() for s in sections])
+
+    assert np.all(d2_actual == d2_expected), f"delta: {d2_actual - d2_expected}"
+    assert np.all(a2_actual == a2_expected), f"delta: {a2_actual - a2_expected}"
+
+
+def check_voltage(t, v, sections):
+    v0 = -65.0
+    erev = 20.0
+    rates = [s(0.5).two_radii.square_diam() + s(0.5).area() for s in sections]
+
+    v_exact = np.array([erev + (v0 - erev) * np.exp(-c * t) for c in rates])
+    assert np.all(
+        np.abs(v - v_exact) < 0.01 * np.abs(v0)
+    ), f"{np.max(np.abs(v - v_exact))=}"
+
+
+def simulate():
+    nsec = 10
+
+    # Independent sections are needed to create independent ODEs.
+    diams = 1.0 + 0.1 * np.arange(1, nsec + 1)
+    sections = [h.Section() for _ in range(nsec)]
+    for d, s in zip(diams, sections):
+        s.nseg = 1
+        s.insert("two_radii")
+        s.pt3dadd(0.0, 0.0, 0.0, d)
+        s.pt3dadd(1.0, 2.0, 0.0, d)
+
+    t_hoc = h.Vector().record(h._ref_t)
+    v_hoc = [h.Vector().record(s(0.5)._ref_v) for s in sections]
+
+    check_morphology(sections)
+
+    h.stdinit()
+    h.dt = 0.001
+    h.run()
+
+    check_morphology(sections)
+
+    t = np.array(t_hoc.as_numpy())
+    v = np.array([vv.as_numpy() for vv in v_hoc])
+
+    check_voltage(t, v, sections)
+
+
+if __name__ == "__main__":
+    simulate()

test/usecases/morphology/two_radii.mod

@@ -0,0 +1,28 @@
+NEURON {
+  SUFFIX two_radii
+  NONSPECIFIC_CURRENT il
+  RANGE inv
+}
+
+ASSIGNED {
+  il
+  inv
+  diam
+  area
+}
+
+INITIAL {
+  inv = 1.0 / (square_diam() + area)
+}
+
+BREAKPOINT {
+  il = (square_diam() + area) * 0.001 * (v - 20.0)
+}
+
+FUNCTION square_diam() {
+  square_diam = diam * diam
+}
+
+FUNCTION square_area() {
+  square_area = area * area
+}