fix area units in unit tests for mc models

src/fvm.hpp

@@ -310,9 +310,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
         mechanisms_.push_back(
             helper->new_mechanism(voltage_, current_, node_index)
         );
-        std::cout << "created mech " << mech.first
-                  << " with size " << mechanisms_.back()->size()
-                  << " and indexes " << node_index << "\n";
     }
 
     /////////////////////////////////////////////
@@ -400,9 +397,12 @@ void fvm_cell<T, I>::setup_matrix(T dt)
     //        .     .  .
     //       l[i] . . d[i]
     //
-    d(all) = 1.0;
+    //d(all) = 1.0;
+    d(all) = cv_areas_;
     for(auto i=1u; i<d.size(); ++i) {
-        auto a = dt * face_alpha_[i]; // TODO get this right
+        // TODO get this right
+        // probably requires scaling a by cv_areas_[i] and cv_areas_[p[i]]
+        auto a = 1e7*dt * face_alpha_[i];
 
         d[i] +=  a;
         l[i]  = -a;
@@ -411,12 +411,14 @@ void fvm_cell<T, I>::setup_matrix(T dt)
         // add contribution to the diagonal of parent
         d[p[i]] += a;
     }
+    //std::cout << "d " << d << " l " << l << " u " << u << "\n";
 
     // the RHS of the linear system is
     //      V[i] - dt/cm*(im - ie)
     auto factor = 10.*dt;
     for(auto i=0u; i<d.size(); ++i) {
-        rhs[i] = voltage_[i] - factor/cv_capacitance_[i]*current_[i];
+        //rhs[i] = voltage_[i] - factor/cv_capacitance_[i]*current_[i];
+        rhs[i] = cv_areas_[i]*(voltage_[i] - factor/cv_capacitance_[i]*current_[i]);
     }
 }
 
@@ -442,9 +444,13 @@ void fvm_cell<T, I>::advance(T dt)
         m->nrn_current();
     }
 
-    if(t_>=5. && t_<8.) {
-        current_[0] -= 0.01;
-    }
+    // the factor scales the injected current to 10^2.nA
+    auto ie_factor = 100.;
+    auto ie = 0.1;
+    auto loc = size()-1;
+    //auto loc = 0;
+    if(t_>=5. && t_<8.)
+      current_[loc] -= ie_factor*ie/cv_areas_[loc];
 
     //std::cout << "t " << t_ << " current " << current_;
 

tests/test_run.cpp

@@ -13,8 +13,8 @@ TEST(run, single_compartment)
     // setup global state for the mechanisms
     nest::mc::mechanisms::setup_mechanism_helpers();
 
-    // Soma with diameter 18.8um -> 1.88e-3 cm and HH channel
-    auto soma = cell.add_soma(1.88e-3/2.0);
+    // Soma with diameter 18.8um and HH channel
+    auto soma = cell.add_soma(18.8/2.0);
     soma->mechanism("membrane").set("r_L", 123); // no effect for single compartment cell
     soma->add_mechanism(hh_parameters());
 
@@ -36,8 +36,7 @@ TEST(run, single_compartment)
 
     // run the simulation
     auto dt     =   0.02; // ms
-    //auto tfinal =   0.10; // ms
-    auto tfinal =   10.; // ms
+    auto tfinal =   30.; // ms
     int nt = tfinal/dt;
     std::vector<double> result;
     result.push_back(model.voltage()[0]);
@@ -66,13 +65,13 @@ TEST(run, ball_and_stick)
     // setup global state for the mechanisms
     nest::mc::mechanisms::setup_mechanism_helpers();
 
-    // Soma with diameter 12.6157 and HH channel
+    // Soma with diameter 12.6157 um and HH channel
     auto soma = cell.add_soma(12.6157/2.0);
     soma->add_mechanism(hh_parameters());
 
-    // add dendrite with passive channel and 10 compartments
+    // add dendrite of length 200 um and diameter 1 um with passive channel
     auto dendrite = cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
-    dendrite->set_compartments(5);
+    dendrite->set_compartments(5); // 5 compartments
     dendrite->add_mechanism(pas_parameters());
 
     // make the lowered finite volume cell
@@ -83,12 +82,30 @@ TEST(run, ball_and_stick)
     // set initial conditions
     using memory::all;
     model.voltage()(all) = -65.;
+    model.initialize(); // have to do this _after_ initial conditions are set
 
     // run the simulation
     auto dt = 0.02; // ms
-    model.advance(dt);
+    auto tfinal = 20.; // ms
+    int nt = tfinal/dt;
+    std::vector<double> result;
+    result.push_back(model.voltage()[0]);
+    for(auto i=0; i<nt; ++i) {
+        model.advance(dt);
+        result.push_back(model.voltage()[0]);
+    }
+    std::cout << "took " << nt << " time steps" << std::endl;
 
     std::cout << model.voltage() << "\n";
+
+    {
+        std::ofstream fid("v.dat");
+        auto t = 0.;
+        for(auto v:result) {
+            fid << t << " " << v << "\n";
+            t += dt;
+        }
+    }
 }
 
 TEST(run, cable)

tests/vis.py

@@ -5,8 +5,8 @@ data = np.loadtxt('v.dat')
 t = data[:,0]
 v = data[:,1]
 
-
 pyplot.plot(t, v)
+pyplot.ylim(-80, 50)
 pyplot.xlabel('time (ms)')
 pyplot.ylabel('mV')
 pyplot.show()