Fast membrane current calculation

coreneuron/nrniv/main1.cpp

@@ -37,6 +37,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/engine.h"
 #include "coreneuron/utils/randoms/nrnran123.h"
 #include "coreneuron/nrnconf.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrnmpi/nrnmpi.h"
@@ -126,10 +127,13 @@ char* prepare_args(int& argc, char**& argv, int use_mpi, const char* arg) {
     return first;
 }
 
-int corenrn_embedded_run(int nthread, int have_gaps, int use_mpi, const char* arg) {
+int corenrn_embedded_run(int nthread, int have_gaps, int use_mpi, int use_fast_imem, const char* arg) {
     corenrn_embedded = 1;
     corenrn_embedded_nthread = nthread;
     coreneuron::nrn_have_gaps = have_gaps;
+    if (use_fast_imem) {
+        coreneuron::nrn_use_fast_imem = true;
+    }
 
     set_openmp_threads(nthread);
     int argc = 0;

coreneuron/nrniv/nrn_setup.cpp

@@ -34,6 +34,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/nrnconf.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrniv/nrniv_decl.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrniv/vrecitem.h"
 #include "coreneuron/nrniv/multisend.h"
@@ -50,6 +51,8 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/utils/reports/nrnsection_mapping.h"
 
 // callbacks into nrn/src/nrniv/nrnbbcore_write.cpp
+#include "coreneuron/nrnoc/fast_imem.h"
+#include "coreneuron/nrniv/nrniv_decl.h"
 #include "coreneuron/nrniv/nrn2core_direct.h"
 
 int corenrn_embedded;
@@ -795,6 +798,9 @@ void nrn_setup(const char* filesdat,
     mk_cell_indices();
 #endif
 
+    /// Allocate memory for fast_imem calculation
+    nrn_fast_imem_alloc();
+
     /// Generally, tables depend on a few parameters. And if those parameters change,
     /// then the table needs to be recomputed. This is obviously important in NEURON
     /// since the user can change those parameters at any time. However, there is no
@@ -912,17 +918,28 @@ int nrn_i_layout(int icnt, int cnt, int isz, int sz, int layout) {
     return 0;
 }
 
+// This function is related to nrn_dblpntr2nrncore in Neuron to determine which values should
+// be transferred from CoreNeuron. Types correspond to the value to be transferred based on
+// mech_type enum or non-artificial cell mechanisms.
 // take into account alignment, layout, permutation
-// only voltage or mechanism data index allowed. (mtype 0 means time)
+// only voltage, i_membrane_ or mechanism data index allowed. (mtype 0 means time)
 double* stdindex2ptr(int mtype, int index, NrnThread& nt) {
-    if (mtype == -5) {  // voltage
+    if (mtype == voltage) {  // voltage
         int v0 = nt._actual_v - nt._data;
         int ix = index;  // relative to _actual_v
         nrn_assert((ix >= 0) && (ix < nt.end));
         if (nt._permute) {
             node_permute(&ix, 1, nt._permute);
         }
         return nt._data + (v0 + ix);                // relative to nt._data
+    } else if (mtype == i_membrane_) {              // membrane current from fast_imem calculation
+            int i_mem = nt.nrn_fast_imem->nrn_sav_rhs - nt._data;
+            int ix = index;  // relative to nrn_fast_imem->nrn_sav_rhs
+            nrn_assert((ix >= 0) && (ix < nt.end));
+            if (nt._permute) {
+                node_permute(&ix, 1, nt._permute);
+            }
+            return nt._data + (i_mem + ix);         // relative to nt._data
     } else if (mtype > 0 && mtype < n_memb_func) {  //
         Memb_list* ml = nt._ml_list[mtype];
         nrn_assert(ml);
@@ -1150,6 +1167,10 @@ void nrn_cleanup(bool clean_ion_global_map) {
         }
 
         free_memory(nt->_ml_list);
+
+        if (nt->nrn_fast_imem) {
+            fast_imem_free();
+        }
     }
 
 #if NRN_MULTISEND

coreneuron/nrniv/nrniv_decl.h

@@ -35,6 +35,12 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/utils/endianness.h"
 #include "coreneuron/nrniv/nrnoptarg.h"
 namespace coreneuron {
+
+/// Mechanism type to be used from stdindex2ptr and nrn_dblpntr2nrncore (in Neuron)
+/// Values of the mechanism types should be negative numbers to avoid any conflict with
+/// mechanism types of Memb_list(>0) or time(0) passed from Neuron
+enum mech_type {voltage = -1, i_membrane_ = -2};
+
 extern int cvode_active_;
 /// Vector of maps for negative presyns
 extern std::vector<std::map<int, PreSyn*> > neg_gid2out;

coreneuron/nrnoc/fadvance_core.cpp

@@ -30,6 +30,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include "coreneuron/nrnconf.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrnmpi/nrnmpi.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrniv/nrn_acc_manager.h"
 #include "coreneuron/utils/reports/nrnreport.h"
@@ -190,6 +191,9 @@ void update(NrnThread* _nt) {
         assert(_nt->tml->index == CAP);
         nrn_cur_capacitance(_nt, _nt->tml->ml, _nt->tml->index);
     }
+    if (nrn_use_fast_imem) { 
+        nrn_calc_fast_imem(_nt);
+    }
 }
 
 void nonvint(NrnThread* _nt) {

coreneuron/nrnoc/fast_imem.cpp

@@ -0,0 +1,78 @@
+/*
+Copyright (c) 2019, Blue Brain Project
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+1. Redistributions of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+3. Neither the name of the copyright holder nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#include "coreneuron/nrnconf.h"
+#include "coreneuron/nrnoc/fast_imem.h"
+#include "coreneuron/nrniv/memory.h"
+#include "coreneuron/nrnmpi/nrnmpi.h"
+
+namespace coreneuron {
+
+extern int nrn_nthread;
+extern NrnThread *nrn_threads;
+bool nrn_use_fast_imem;
+
+void fast_imem_free() {
+    for (NrnThread* nt = nrn_threads; nt < nrn_threads + nrn_nthread; ++nt) {
+        if (nt->nrn_fast_imem) {
+            free(nt->nrn_fast_imem->nrn_sav_rhs);
+            free(nt->nrn_fast_imem->nrn_sav_d);
+            free(nt->nrn_fast_imem);
+            nt->nrn_fast_imem = nullptr;
+        }
+    }
+}
+
+void nrn_fast_imem_alloc() {
+    if (nrn_use_fast_imem) {
+        fast_imem_free();
+        for (NrnThread* nt = nrn_threads; nt < nrn_threads + nrn_nthread; ++nt) {
+            int n = nt->end;
+            nt->nrn_fast_imem = (NrnFastImem*)ecalloc(1, sizeof(NrnFastImem));
+            nt->nrn_fast_imem->nrn_sav_rhs = (double*)emalloc_align(n * sizeof(double));
+            nt->nrn_fast_imem->nrn_sav_d = (double*)emalloc_align(n * sizeof(double));
+        }
+    }
+}
+
+void nrn_calc_fast_imem(NrnThread* nt) {
+    int i1 = 0;
+    int i3 = nt->end;
+
+    double* vec_rhs = nt->_actual_rhs;
+    double* vec_area = nt->_actual_area;
+
+    double* fast_imem_d = nt->nrn_fast_imem->nrn_sav_d;
+    double* fast_imem_rhs = nt->nrn_fast_imem->nrn_sav_rhs;
+    for (int i = i1; i < i3 ; ++i) {
+        fast_imem_rhs[i] = (fast_imem_d[i]*vec_rhs[i] + fast_imem_rhs[i])*vec_area[i]*0.01;
+    }
+}
+
+}
+

coreneuron/nrnoc/fast_imem.h

@@ -0,0 +1,63 @@
+/*
+Copyright (c) 2019, Blue Brain Project
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted provided that the following conditions are met:
+1. Redistributions of source code must retain the above copyright notice,
+   this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+3. Neither the name of the copyright holder nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
+THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef fast_imem_h
+#define fast_imem_h
+
+#include "coreneuron/nrnoc/multicore.h"
+
+namespace coreneuron {
+
+/* Bool global variable to define if the fast_imem
+ * calculations should be enabled.
+ */
+extern bool nrn_use_fast_imem;
+
+/* Free memory allocated for the fast current membrane calculation.
+ * Found in src/nrnoc/multicore.c in NEURON.
+ */
+void fast_imem_free();
+
+/* Allocate memory for the rhs and d arrays needed for the fast
+ * current membrane calculation.
+ * Found in src/nrnoc/multicore.c in NEURON.
+ */
+static void fast_imem_alloc();
+
+/* fast_imem_alloc() wrapper.
+ * Found in src/nrnoc/multicore.c in NEURON.
+ */
+void nrn_fast_imem_alloc();
+
+/* Calculate the new values of rhs array at every timestep.
+ * Found in src/nrnoc/fadvance.c in NEURON.
+ */
+void nrn_calc_fast_imem(NrnThread* _nt);
+
+}  // namespace coreneuron
+#endif //fast_imem_h

coreneuron/nrnoc/finitialize.cpp

@@ -27,11 +27,13 @@ THE POSSIBILITY OF SUCH DAMAGE.
 */
 
 #include "coreneuron/nrnconf.h"
+#include "coreneuron/nrnoc/fast_imem.h"
 #include "coreneuron/nrnoc/multicore.h"
 #include "coreneuron/nrnoc/nrnoc_decl.h"
 #include "coreneuron/nrniv/profiler_interface.h"
 
 namespace coreneuron {
+
 void nrn_finitialize(int setv, double v) {
     int i;
     NrnThread* _nt;
@@ -98,6 +100,9 @@ void nrn_finitialize(int setv, double v) {
     }
     for (i = 0; i < nrn_nthread; ++i) {
         setup_tree_matrix_minimal(nrn_threads + i);
+        if (nrn_use_fast_imem) {
+            nrn_calc_fast_imem(nrn_threads + i);
+        }
     }
     for (i = 0; i < nrn_nthread; ++i) {
         nrn_ba(nrn_threads + i, BEFORE_STEP);

coreneuron/nrnoc/multicore.cpp

@@ -137,6 +137,8 @@ void nrn_threads_create(int n) {
                 nt->_watch_types = NULL;
                 nt->mapping = NULL;
                 nt->trajec_requests = NULL;
+
+                nt->nrn_fast_imem = nullptr;
             }
         }
         v_structure_change = 1;

coreneuron/nrnoc/multicore.h

@@ -59,6 +59,11 @@ struct NrnThreadBAList {
     NrnThreadBAList* next;
 };
 
+struct NrnFastImem {
+    double* nrn_sav_rhs;
+    double* nrn_sav_d;
+};
+
 struct TrajectoryRequests {
     void** vpr;       /* PlayRecord Objects known by NEURON */
     double** scatter; /* if bsize == 0, each time step */
@@ -116,6 +121,10 @@ struct NrnThread {
                              compartment */
     double* _shadow_d;    /* Not pointer into _data. Avoid race for multiple POINT_PROCESS in same
                              compartment */
+
+    /* Fast membrane current calculation struct */
+    NrnFastImem* nrn_fast_imem;
+
     int* _v_parent_index;
     int* _permute;
     char* _sp13mat;              /* handle to general sparse matrix */

coreneuron/nrnoc/treeset_core.cpp

@@ -67,6 +67,13 @@ static void nrn_rhs(NrnThread* _nt) {
         vec_d[i] = 0.;
     }
 
+    if (_nt->nrn_fast_imem) {
+        for (i = i1; i < i3; ++i) {
+            _nt->nrn_fast_imem->nrn_sav_rhs[i] = 0.;
+            _nt->nrn_fast_imem->nrn_sav_d[i] = 0.;
+        }
+    }
+
     nrn_ba(_nt, BEFORE_BREAKPOINT);
     /* note that CAP has no current */
     for (tml = _nt->tml; tml; tml = tml->next)
@@ -82,6 +89,16 @@ static void nrn_rhs(NrnThread* _nt) {
             }
 #endif
         }
+
+    if (_nt->nrn_fast_imem) {
+        /* _nrn_save_rhs has only the contribution of electrode current
+           so here we transform so it only has membrane current contribution
+        */
+        double* p = _nt->nrn_fast_imem->nrn_sav_rhs;
+        for (i = i1; i < i3; ++i) {
+            p[i] -= vec_rhs[i];
+        }
+    }
 
 /* now the internal axial currents.
 The extracellular mechanism contribution is already done.
@@ -151,6 +168,16 @@ static void nrn_lhs(NrnThread* _nt) {
     double* vec_b = &(VEC_B(0));
     int* parent_index = _nt->_v_parent_index;
 
+    if (_nt->nrn_fast_imem) {
+        /* _nrn_save_d has only the contribution of electrode current
+           so here we transform so it only has membrane current contribution
+        */
+        double* p = _nt->nrn_fast_imem->nrn_sav_d;
+        for (i = i1; i < i3; ++i) {
+            p[i] += vec_d[i];
+        }
+    }
+
 /* now add the axial currents */
 // clang-format off
     #pragma acc parallel loop present(          \