Reading global ion variables from NEURON and setting them in nrn_ion_init

Author: iomaganaris

coreneuron/apps/main1.cpp

@@ -17,6 +17,7 @@
 #include <dlfcn.h>
 #include <memory>
 #include <vector>
+#include "coreneuron/permute/data_layout.hpp"
 
 #include "coreneuron/config/config.h"
 #include "coreneuron/utils/randoms/nrnran123.h"
@@ -218,9 +219,7 @@ void nrn_init_and_load_data(int argc,
     // full path of files.dat file
     std::string filesdat(corenrn_param.datpath + "/" + corenrn_param.filesdat);
 
-    // read the global variable names and set their values from globals.dat
-    set_globals(corenrn_param.datpath.c_str(), (corenrn_param.seed >= 0), corenrn_param.seed);
-
+    
     // set global variables for start time, timestep and temperature
     if (!corenrn_embedded) {
         t = checkPoints.restore_time();
@@ -514,8 +513,10 @@ extern "C" void mk_mech_init(int argc, char** argv) {
         out.close();
     }
 
-    // reads mechanism information from bbcore_mech.dat
+    // read the global variable names and set their values from globals.dat
+    set_globals(corenrn_param.datpath.c_str(), (corenrn_param.seed >= 0), corenrn_param.seed);
 
+    // reads mechanism information from bbcore_mech.dat
     mk_mech((corenrn_param.datpath).c_str());
 }
 
@@ -578,6 +579,21 @@ extern "C" int run_solve_core(int argc, char** argv) {
             abort();
         }
 
+        // for (int i = 0; i < nrn_nthread; ++i) {
+        //     NrnThread* nt = nrn_threads + i;
+        //     for (auto tml = nt->tml; tml; tml = tml->next) {
+        //         auto ml = tml->ml;
+        //         auto type = tml->index;
+        //         int n = ml->nodecount;
+        //         int szdp = corenrn.get_prop_dparam_size()[type];
+        //         std::cout << "ml type " << type << std::endl;
+        //         int pcnt = nrn_soa_padded_size(n, SOA_LAYOUT) * szdp;
+        //         for(int i = 0 ; i < pcnt; i++) {
+        //             printf("ml->pdata[%d] = %d\n", i, ml->pdata[i]);
+        //         }
+        //     }
+        // }
+
         // TODO : if some ranks are empty then restore will go in deadlock
         // phase (as some ranks won't have restored anything and hence return
         // false in checkpoint_initialize

coreneuron/io/global_vars.cpp

@@ -11,6 +11,7 @@
 #include <map>
 #include <string>
 #include <algorithm>
+#include <unordered_map>
 
 #include "coreneuron/utils/randoms/nrnran123.h"
 #include "coreneuron/nrnconf.h"
@@ -28,6 +29,8 @@ using N2V = std::map<std::string, PSD>;
 
 static N2V* n2v;
 
+extern std::unordered_map<std::string, double> nrn_ion_init;
+
 void hoc_register_var(DoubScal* ds, DoubVec* dv, VoidFunc*) {
     if (!n2v) {
         n2v = new N2V();
@@ -109,10 +112,14 @@ void set_globals(const char* path, bool cli_global_seed, int cli_global_seed_val
                 if (strcmp(name, "0") == 0) {
                     break;
                 }
+                printf("name %s value %lf\n", name, val);
                 it = n2v->find(name);
                 if (it != n2v->end()) {
                     nrn_assert(it->second.first == 0);
                     *(it->second.second) = val;
+                } else if (static_cast<std::string>(name).find("ion") != std::string::npos) {
+                    printf("ion %s value %lf\n", name, val);
+                    nrn_ion_init[name] = val;
                 }
             } else if (sscanf(line, "%[^[][%d]\n", name, &n) == 2) {
                 if (strcmp(name, "0") == 0) {

coreneuron/mechanism/eion.cpp

@@ -60,6 +60,15 @@ double nrn_ion_charge(int type) {
     return global_charge(type);
 }
 
+std::unordered_map<std::string, double> nrn_ion_init = {
+    {"nai0_na_ion", DEF_nai},
+    {"nao0_na_ion", DEF_nao},
+    {"ki0_k_ion", DEF_ki},
+    {"ko0_k_ion", DEF_ko},
+    {"cai0_ca_ion", DEF_cai},
+    {"cao0_ca_ion", DEF_cao},
+};
+
 void ion_reg(const char* name, double valence) {
     char buf[7][50];
 #define VAL_SENTINAL -10000.
@@ -109,22 +118,30 @@ void ion_reg(const char* name, double valence) {
         sprintf(buf[1], "%so0_%s", name, buf[0]);
         if (strcmp("na", name) == 0) {
             na_ion = mechtype;
-            global_conci(mechtype) = DEF_nai;
-            global_conco(mechtype) = DEF_nao;
+            global_conci(mechtype) = nrn_ion_init["nai0_na_ion"];
+            global_conco(mechtype) = nrn_ion_init["nao0_na_ion"];
             global_charge(mechtype) = 1.;
         } else if (strcmp("k", name) == 0) {
             k_ion = mechtype;
-            global_conci(mechtype) = DEF_ki;
-            global_conco(mechtype) = DEF_ko;
+            global_conci(mechtype) = nrn_ion_init["ki0_k_ion"];
+            global_conco(mechtype) = nrn_ion_init["ko0_k_ion"];
             global_charge(mechtype) = 1.;
         } else if (strcmp("ca", name) == 0) {
             ca_ion = mechtype;
-            global_conci(mechtype) = DEF_cai;
-            global_conco(mechtype) = DEF_cao;
+            global_conci(mechtype) = nrn_ion_init["cai0_ca_ion"];
+            global_conco(mechtype) = nrn_ion_init["cao0_ca_ion"];
             global_charge(mechtype) = 2.;
         } else {
-            global_conci(mechtype) = DEF_ioni;
-            global_conco(mechtype) = DEF_iono;
+            if (nrn_ion_init[static_cast<std::string>(buf[0])]) {
+                global_conci(mechtype) = nrn_ion_init[static_cast<std::string>(buf[0])];
+            } else {
+                global_conci(mechtype) = DEF_ioni;
+            }
+            if (nrn_ion_init[static_cast<std::string>(buf[1])]) {
+                global_conco(mechtype) = nrn_ion_init[static_cast<std::string>(buf[0])];
+            } else {
+                global_conco(mechtype) = DEF_iono;
+            }
             global_charge(mechtype) = VAL_SENTINAL;
         }
     }
@@ -275,6 +292,7 @@ void nrn_cur_ion(NrnThread* nt, Memb_list* ml, int type) {
     for (int _iml = 0; _iml < _cntml_actual; ++_iml) {
         dcurdv = 0.;
         cur = 0.;
+        // printf("nrn_cur_ion %d iontype[%d] %d\n", type, _iml, iontype);
         if (iontype & 0100) {
             erev = nrn_nernst(conci, conco, charge, celsius);
         }
@@ -295,10 +313,10 @@ void nrn_init_ion(NrnThread* nt, Memb_list* ml, int type) {
         return;
     }
 
-    /*printf("ion_init %s\n", memb_func[type].sym->name);*/
     int _cntml_padded = ml->_nodecount_padded;
     pd = ml->data;
     ppd = ml->pdata;
+    printf("nrn_init_ion %d %s conci %lf conco %lf charge %lf celsius %lf\n", type, nrn_get_mechname(type), conci0, conco0, charge, celsius);
     // There was no async(...) clause in the initial OpenACC implementation, so
     // no `nowait` clause has been added to the OpenMP implementation. TODO:
     // verify if this can be made asynchronous or if there is a strong reason it