Added support for euler (similar to recent thread safe support in NEURON) (#68)

* Added support for euler (similar to thread safe support in NEURON)
    - updated kinderiv.py to generate
    - added euler_thread from neuron (updated for aos/soa layout)
    - mod2c module updated
* Renaming and minor changes for euler support:
 - _derivimplic changed to _derivimplicit
 - _cb changed to newton
 - _NRN_EULER_CASES switch case added
* Added some notes and todos

Author: pramodk

coreneuron/kinderiv.py

@@ -1,6 +1,9 @@
 #!/usr/bin/python
 
 # read the translated mod files and construct _kinderiv.h
+# this is implemented to avoid function callbacks using
+# function pointers while using methods like euler, newton
+# or kinetic schemes. This is required for GPU implementation.
 
 # if _kinderiv.h already exists and is the same as what this script generates
 # then do not update it (to avoid re-compiling all the mech c files if
@@ -12,17 +15,21 @@
 import os
 
 fnames = [f.replace('.mod', '.c') for f in os.listdir('.') if f.endswith('.mod')]
+euler = []
 deriv = []
 kin = []
+
 for fname in fnames:
   f = open(fname, "r")
   for line in f:
     word = line.split()
     if len(word) > 3:
-      if word[0] == '/*' and word[1] == '_derivimplic_':
+      if word[0] == '/*' and word[1] == '_derivimplicit_':
         deriv.append([word[2], word[3], fname, word[1]])
       if word[0] == '/*' and word[1] == '_kinetic_':
         kin.append([word[2], word[3], fname, word[1]])
+      if word[0] == '/*' and word[1] == '_euler_':
+        euler.append([word[2], word[3], fname, word[1]])
   f.close()
 
 fout = open(kftmp, "w")
@@ -32,7 +39,7 @@
 ''')
 
 fout.write("\n/* data used to construct this file */\n")
-for l in [deriv, kin]:
+for l in [deriv, kin,  euler]:
   for item in l:
     fout.write('/*')
     for word in item:
@@ -44,7 +51,7 @@
   fout.write('#pragma acc routine seq\n')
   fout.write('extern int %s%s(_threadargsproto_);\n' % (item[0], item[1]))
   fout.write('#pragma acc routine seq\n')
-  fout.write('extern int _cb_%s%s(_threadargsproto_);\n' % (item[0], item[1]))
+  fout.write('extern int _newton_%s%s(_threadargsproto_);\n' % (item[0], item[1]))
 
 for item in kin:
   fout.write('#pragma acc routine seq\n')
@@ -53,27 +60,35 @@
 fout.write("\n/* callback indices */\n")
 derivoffset = 1
 kinoffset = 1
+
 for i, item in enumerate(deriv):
-  fout.write('#define _derivimplic_%s%s %d\n' % (item[0], item[1], i + derivoffset))
+  fout.write('#define _derivimplicit_%s%s %d\n' % (item[0], item[1], i + derivoffset))
 for i, item in enumerate(kin):
   fout.write('#define _kinetic_%s%s %d\n' % (item[0], item[1], i + kinoffset))
+for i, item in enumerate(euler):
+  fout.write('#define _euler_%s%s %d\n' % (item[0], item[1], i + derivoffset))
 
 fout.write("\n/* switch cases */\n")
-fout.write("\n#define _NRN_DERIVIMPLIC_CASES \\\n")
-for item in deriv:
-  fout.write("  case _derivimplic_%s%s: %s%s(_threadargs_); break; \\\n" % (item[0], item[1], item[0], item[1]))
+fout.write("\n#define _NRN_DERIVIMPLICIT_CASES \\\n")
+for item in (deriv):
+  fout.write("  case _derivimplicit_%s%s: %s%s(_threadargs_); break; \\\n" % (item[0], item[1], item[0], item[1]))
 fout.write("\n")
 
-fout.write("\n#define _NRN_DERIVIMPLIC_CB_CASES \\\n")
+fout.write("\n#define _NRN_DERIVIMPLICIT_NEWTON_CASES \\\n")
 for item in deriv:
-  fout.write("  case _derivimplic_%s%s: _cb_%s%s(_threadargs_); break; \\\n" % (item[0], item[1], item[0], item[1]))
+  fout.write("  case _derivimplicit_%s%s: _newton_%s%s(_threadargs_); break; \\\n" % (item[0], item[1], item[0], item[1]))
 fout.write("\n")
 
 fout.write("\n#define _NRN_KINETIC_CASES \\\n")
 for item in kin:
   fout.write("  case _kinetic_%s%s: %s%s(so, rhs, _threadargs_); break; \\\n" % (item[0], item[1], item[0], item[1]))
 fout.write("\n")
 
+fout.write("\n#define _NRN_EULER_CASES \\\n")
+for item in euler:
+  fout.write("  case _euler_%s%s: %s%s(_threadargs_); break; \\\n" % (item[0], item[1], item[0], item[1]))
+fout.write("\n")
+
 fout.write('\n#endif\n')
 fout.close()
 

coreneuron/mech/mod2c_core_thread.h

@@ -37,15 +37,25 @@ typedef int (*SPFUN)(struct SparseObj*, double*, _threadargsproto_);
 
 #else
 
+/**
+ * \todo: typedefs like DIFUN can be removed
+ * \todo: macros for difun, newtfun, eulerfun are not necessary
+ *        and need to be refactored.
+ */
+
 typedef int DIFUN;
 typedef int NEWTFUN;
 typedef int SPFUN;
+typedef int EULFUN;
 #pragma acc routine seq
-extern int nrn_derivimplic_steer(int, _threadargsproto_);
-#define difun(arg) nrn_derivimplic_steer(arg, _threadargs_);
+extern int nrn_derivimplicit_steer(int, _threadargsproto_);
+#define difun(arg) nrn_derivimplicit_steer(arg, _threadargs_);
 #pragma acc routine seq
 extern int nrn_newton_steer(int, _threadargsproto_);
 #define newtfun(arg) nrn_newton_steer(arg, _threadargs_);
+#pragma acc routine seq
+extern int nrn_euler_steer(int, _threadargsproto_);
+#define eulerfun(arg) nrn_euler_steer(arg, _threadargs_);
 
 #endif
 
@@ -97,6 +107,8 @@ typedef struct SparseObj {  /* all the state information */
 extern int nrn_kinetic_steer(int, SparseObj*, double*, _threadargsproto_);
 #define spfun(arg1, arg2, arg3) nrn_kinetic_steer(arg1, arg2, arg3, _threadargs_);
 
+#pragma acc routine seq
+extern int euler_thread(int, int*, int*, EULFUN, _threadargsproto_);
 #pragma acc routine seq
 extern int derivimplicit_thread(int, int*, int*, DIFUN, _threadargsproto_);
 #pragma acc routine seq

coreneuron/scopmath_core/dimplic.c

@@ -1,6 +1,15 @@
-/*
-  hopefully a temporary expedient to work around the inability to
-  pass function pointers as arguments
+/** When we use solve methods like euler, newton or kinetic schemes,
+ *  the state/current updates function need to call solver methods
+ *  defined in coreneuron. This is typically done via function pointers.
+ *  But for GPU implementation using OpenACC, we can't pass function
+ *  pointers.The temporary "workaround" for this was to generate
+ *  switch case to select the proper callback function. This is implemented
+ *  using python script that look into translated file and generate
+ *  _kinderiv.h which has cases for steer functions defined below.
+ *  This allows OpenACC to select gpu implementations at compile
+ *  time.
+ *  \todo: eulerfun/difun are legacy macros and can be replaced with
+ *         actual steer function for euler/derivimplicit methods.
 */
 
 #include "coreneuron/mech/cfile/scoplib.h"
@@ -12,17 +21,39 @@ int derivimplicit_thread(int n, int* slist, int* dlist, DIFUN fun, _threadargspr
     return 0;
 }
 
-int nrn_derivimplic_steer(int fun, _threadargsproto_) {
-    switch (fun) { _NRN_DERIVIMPLIC_CASES }
+int nrn_derivimplicit_steer(int fun, _threadargsproto_) {
+    switch (fun) { _NRN_DERIVIMPLICIT_CASES }
+    return 0;
+}
+
+int nrn_euler_steer(int fun, _threadargsproto_) {
+    switch (fun) { _NRN_EULER_CASES }
     return 0;
 }
 
 int nrn_newton_steer(int fun, _threadargsproto_) {
-    switch (fun) { _NRN_DERIVIMPLIC_CB_CASES }
+    switch (fun) { _NRN_DERIVIMPLICIT_NEWTON_CASES }
     return 0;
 }
 
 int nrn_kinetic_steer(int fun, SparseObj* so, double* rhs, _threadargsproto_) {
     switch (fun) { _NRN_KINETIC_CASES }
     return 0;
 }
+
+// derived from nrn/src/scopmath/euler.c
+// updated for aos/soa layout index
+int euler_thread(int neqn, int* var, int* der, DIFUN fun, _threadargsproto_) {
+
+    double dt = _nt->_dt;
+    int i;
+
+    /* calculate the derivatives */
+     eulerfun(fun);
+
+    /* update dependent variables */
+    for (i = 0; i < neqn; i++)
+        _p[var[i]*_STRIDE] += dt * (_p[der[i]*_STRIDE]);
+
+    return 0;
+}