diff --git a/src/fvm.hpp b/src/fvm.hpp index 7abf6ed4b2c465122f8b30e0b5c97c1d9a135df0..33439dd95ae23431eb8d605c37eb2e441cc3ffe8 100644 --- a/src/fvm.hpp +++ b/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_; diff --git a/tests/test_run.cpp b/tests/test_run.cpp index e4002736c32b1075c3a2e6e01bedd27a4923e1d8..12012f0f8fb885f50152edefdf87e05fd5010ef8 100644 --- a/tests/test_run.cpp +++ b/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) diff --git a/tests/vis.py b/tests/vis.py index 751c1fd26e915cf9c725c2146098311a8420e3c4..9e6290e137c56b552fc25e22b3d28217f5526de5 100644 --- a/tests/vis.py +++ b/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()