Enable fast_imem on GPU.

Updates mod2c/nmodl submodule commits to include relevant fixes,
BlueBrain/mod2c#64 and
BlueBrain/nmodl#681.

Closes #197.

Author: olupton

coreneuron/gpu/nrn_acc_manager.cpp

@@ -291,6 +291,23 @@ void setup_nrnthreads_on_device(NrnThread* threads, int nthreads) {
             acc_memcpy_to_device(&(d_nt->_shadow_d), &d_shadow_ptr, sizeof(double*));
         }
 
+        /* Fast membrane current calculation struct */
+        if (nt->nrn_fast_imem) {
+            auto* d_fast_imem = reinterpret_cast<NrnFastImem*>(
+                acc_copyin(nt->nrn_fast_imem, sizeof(NrnFastImem)));
+            acc_memcpy_to_device(&(d_nt->nrn_fast_imem), &d_fast_imem, sizeof(NrnFastImem*));
+            {
+                auto* d_ptr = reinterpret_cast<double*>(
+                    acc_copyin(nt->nrn_fast_imem->nrn_sav_rhs, nt->end * sizeof(double)));
+                acc_memcpy_to_device(&(d_fast_imem->nrn_sav_rhs), &d_ptr, sizeof(double*));
+            }
+            {
+                auto* d_ptr = reinterpret_cast<double*>(
+                    acc_copyin(nt->nrn_fast_imem->nrn_sav_d, nt->end * sizeof(double)));
+                acc_memcpy_to_device(&(d_fast_imem->nrn_sav_d), &d_ptr, sizeof(double*));
+            }
+        }
+
         if (nt->n_pntproc) {
             /* copy Point_processes array and fix the pointer to execute net_receive blocks on GPU
              */
@@ -659,6 +676,11 @@ void update_nrnthreads_on_host(NrnThread* threads, int nthreads) {
                 acc_update_self(nt->_shadow_d, pcnt * sizeof(double));
             }
 
+            if (nt->nrn_fast_imem) {
+                acc_update_self(nt->nrn_fast_imem->nrn_sav_rhs, nt->end * sizeof(double));
+                acc_update_self(nt->nrn_fast_imem->nrn_sav_d, nt->end * sizeof(double));
+            }
+
             if (nt->n_pntproc) {
                 acc_update_self(nt->pntprocs, nt->n_pntproc * sizeof(Point_process));
             }
@@ -748,6 +770,11 @@ void update_nrnthreads_on_device(NrnThread* threads, int nthreads) {
                 acc_update_device(nt->_shadow_d, pcnt * sizeof(double));
             }
 
+            if (nt->nrn_fast_imem) {
+                acc_update_device(nt->nrn_fast_imem->nrn_sav_rhs, nt->end * sizeof(double));
+                acc_update_device(nt->nrn_fast_imem->nrn_sav_d, nt->end * sizeof(double));
+            }
+
             if (nt->n_pntproc) {
                 acc_update_device(nt->pntprocs, nt->n_pntproc * sizeof(Point_process));
             }
@@ -787,6 +814,19 @@ void update_voltage_from_gpu(NrnThread* nt) {
     }
 }
 
+/**
+ * @brief Copy fast_imem vectors from GPU to CPU.
+ *
+ */
+void update_fast_imem_from_gpu(NrnThread* nt) {
+    if (nt->compute_gpu && nt->end > 0 && nt->nrn_fast_imem) {
+        int num_fast_imem = nt->end;
+        double* fast_imem_d = nt->nrn_fast_imem->nrn_sav_d;
+        double* fast_imem_rhs = nt->nrn_fast_imem->nrn_sav_rhs;
+#pragma acc update host(fast_imem_d [0:num_fast_imem], fast_imem_rhs [0:num_fast_imem])
+    }
+}
+
 /**
  * Copy weights from GPU to CPU
  *
@@ -940,6 +980,11 @@ void delete_nrnthreads_on_device(NrnThread* threads, int nthreads) {
             acc_delete(nt->pntprocs, nt->n_pntproc * sizeof(Point_process));
         }
 
+        if (nt->nrn_fast_imem) {
+            acc_delete(nt->nrn_fast_imem->nrn_sav_d, nt->end * sizeof(double));
+            acc_delete(nt->nrn_fast_imem->nrn_sav_rhs, nt->end * sizeof(double));
+        }
+
         if (nt->shadow_rhs_cnt) {
             int pcnt = nrn_soa_padded_size(nt->shadow_rhs_cnt, 0);
             acc_delete(nt->_shadow_d, pcnt * sizeof(double));

coreneuron/gpu/nrn_acc_manager.hpp

@@ -29,6 +29,7 @@ void update_net_receive_buffer(NrnThread* _nt);
 void realloc_net_receive_buffer(NrnThread* nt, Memb_list* ml);
 void update_net_send_buffer_on_host(NrnThread* nt, NetSendBuffer_t* nsb);
 void update_voltage_from_gpu(NrnThread* nt);
+void update_fast_imem_from_gpu(NrnThread* nt);
 void update_weights_from_gpu(NrnThread* threads, int nthreads);
 void init_gpu();
 

coreneuron/sim/fadvance_core.cpp

@@ -298,6 +298,9 @@ void nrncore2nrn_send_values(NrnThread* nth) {
         //       Currently we are updating voltages if there is any trajectory
         //       requested by NEURON.
         update_voltage_from_gpu(nth);
+        // \todo Check if this information has been requested by the user for
+        //       this NrnThread object.
+        update_fast_imem_from_gpu(nth);
 
         if (tr->varrays) {  // full trajectories into Vector data
             double** va = tr->varrays;

coreneuron/sim/fast_imem.cpp

@@ -50,6 +50,10 @@ void nrn_calc_fast_imem(NrnThread* nt) {
 
     double* fast_imem_d = nt->nrn_fast_imem->nrn_sav_d;
     double* fast_imem_rhs = nt->nrn_fast_imem->nrn_sav_rhs;
+#pragma acc parallel loop present(vec_rhs,     \
+                                  vec_area,    \
+                                  fast_imem_d, \
+                                  fast_imem_rhs) if (nt->compute_gpu) async(nt->stream_id)
     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/sim/treeset_core.cpp

@@ -44,9 +44,14 @@ static void nrn_rhs(NrnThread* _nt) {
     }
 
     if (_nt->nrn_fast_imem) {
+        double* fast_imem_d = _nt->nrn_fast_imem->nrn_sav_d;
+        double* fast_imem_rhs = _nt->nrn_fast_imem->nrn_sav_rhs;
+#pragma acc parallel loop present(fast_imem_d [i1:i3],                         \
+                                  fast_imem_rhs [i1:i3]) if (_nt->compute_gpu) \
+    async(_nt->stream_id)
         for (int i = i1; i < i3; ++i) {
-            _nt->nrn_fast_imem->nrn_sav_rhs[i] = 0.;
-            _nt->nrn_fast_imem->nrn_sav_d[i] = 0.;
+            fast_imem_d[i] = 0.;
+            fast_imem_rhs[i] = 0.;
         }
     }
 
@@ -71,6 +76,7 @@ static void nrn_rhs(NrnThread* _nt) {
            so here we transform so it only has membrane current contribution
         */
         double* p = _nt->nrn_fast_imem->nrn_sav_rhs;
+#pragma acc parallel loop present(p, vec_rhs) if (_nt->compute_gpu) async(_nt->stream_id)
         for (int i = i1; i < i3; ++i) {
             p[i] -= vec_rhs[i];
         }
@@ -144,6 +150,7 @@ static void nrn_lhs(NrnThread* _nt) {
            so here we transform so it only has membrane current contribution
         */
         double* p = _nt->nrn_fast_imem->nrn_sav_d;
+#pragma acc parallel loop present(p, vec_d) if (_nt->compute_gpu) async(_nt->stream_id)
         for (int i = i1; i < i3; ++i) {
             p[i] += vec_d[i];
         }