From 0c81a3f6d20e64311f013c362f49fa1ae8e9aa01 Mon Sep 17 00:00:00 2001
From: bcumming <bcumming@cscs.ch>
Date: Fri, 15 Apr 2016 16:32:38 +0200
Subject: [PATCH] generate full lin. system and RHS
---
.ycm_extra_conf.py | 4 +-
docs/formulation.tex | 37 ++++++++-----
src/fvm.hpp | 115 +++++++++++++++++++++++++++++++----------
src/matrix.hpp | 2 +-
src/parameter_list.cpp | 6 +--
src/parameter_list.hpp | 36 +++++++------
tests/test_matrix.cpp | 11 ++--
tests/test_run.cpp | 57 ++++++++++++++++++--
8 files changed, 199 insertions(+), 69 deletions(-)
diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 09a57a74..b8378fb8 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -40,7 +40,9 @@ flags = [
'-x',
'c++',
'-I',
- './src',
+ 'src',
+ '-I',
+ '.',
# '-isystem',
# '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1',
# '-I',
diff --git a/docs/formulation.tex b/docs/formulation.tex
index 88e8e0a5..0676ff4f 100644
--- a/docs/formulation.tex
+++ b/docs/formulation.tex
@@ -193,7 +193,7 @@ where $a_{i,\ell}$ and $a_{i,r}$ are the radii of at the left and right end of t
%-------------------------------------------------------------------------------
\subsubsection{Putting it all together}
%-------------------------------------------------------------------------------
-By substituting the volume averaging of the temporal derivative in~\eq{eq:dvdt} approximations for the flux over the surfaces in~\eq{eq:q_ij} and~\eq{eq:cv_volume} respectively into the conservation equation~\eq{eq:cable_balance} we get the following ODE defined for each node in the cell
+By substituting the volume averaging of the temporal derivative in~\eq{eq:dvdt} approximations for the flux over the surfaces in~\eq{eq:J_ij} and~\eq{eq:cv_volume} respectively into the conservation equation~\eq{eq:cable_balance} we get the following ODE defined for each node in the cell
\begin{equation}
\sigma_i c_m \dder{V_i}{t}
= -\sum_{j\in\mathcal{N}_i} {\frac{\sigma_{i,j}}{r_L \Delta x_{i,j}} (V_i-V_j)} - \sigma_i\cdot(i_m(V_i) - i_e(x_i)),
@@ -228,18 +228,25 @@ The current $i_m$ is often a nonlinear function of voltage, so if it was formula
The equations can be rearranged to have all unknown voltage values on the lhs, and values that can be calculated directly on the rhs:
\begin{align}
- & V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\frac{\alpha_{ij}}{\sigma_i} (V_i^{k+1}-V_j^{k+1})}
+ & \sigma_i V_i^{k+1} + \sum_{j\in\mathcal{N}_i} {\Delta t \alpha_{ij} (V_i^{k+1}-V_j^{k+1})}
\nonumber \\
- = & V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e),
+ = & \sigma_i \left( V_i^k - \frac{\Delta t}{c_m}(i_m^{k} - i_e) \right),
\label{eq:ode_linsys}
\end{align}
where the value
\begin{equation}
- \alpha_{ij} = \alpha_{ji} = \frac{\Delta t \sigma_{ij}}{ c_m r_L \Delta x_{ij}}
+ \alpha_{ij} = \alpha_{ji} = \frac{\sigma_{ij}}{ c_m r_L \Delta x_{ij}}
\label{eq:alpha_linsys}
\end{equation}
is a constant that can be computed for each interface between adjacent compartments during set up.
+The left hand side of \eq{eq:ode_linsys} can be rearranged
+\begin{equation}
+ \left[ \sigma_i + \sum_{j\in\mathcal{N}_i} {\Delta t \alpha_{ij}} \right] V_i^{k+1}
+ - \sum_{j\in\mathcal{N}_i} { \Delta t \alpha_{ij} V_j^{k+1}},
+\end{equation}
+which gives the coefficients for the linear system.
+
%-------------------------------------------------------------------------------
\subsubsection{Example: unbranched uniform cable}
%-------------------------------------------------------------------------------
@@ -247,18 +254,22 @@ For an unrbanched uniform cable of constant radius $a$, with length $L$ and $n$
\begin{align}
\Delta x_{ij} &= \Delta x = \frac{L}{n-1}, \nonumber \\
\sigma_{ij} &= \pi a^2, \nonumber \\
- \sigma_{i} &= 2 \pi a \Delta x, \nonumber \\
- \alpha_{ij} &= \frac{\pi a^2\Delta t}{c_m r_L\Delta x}, \nonumber \\
- \frac{\alpha_{ij}}{\sigma_i}
- &= \frac{a\Delta t}{2c_m r_L\Delta x^2}. \nonumber
+ \sigma_{i} &= \sigma_s = 2 \pi a \Delta x, \nonumber \\
+ \alpha_{ij} &= \alpha = \frac{\pi a^2}{c_m r_L\Delta x}, \nonumber
\end{align}
-With these simplifications, the lhs of the linear system is
+With these simplifications, the LHS of the linear system is
\begin{align}
- & V_i^{k+1} + \beta (V_i^{k+1}-V_{i+1}^{k+1}) + \beta (V_i^{k+1}-V_{i-1}^{k+1})
- \nonumber \\
- = & (1+2\beta)V_i^{k+1} - \beta V_{i+1}^{k+1} - \beta V_{i-1}^{k+1}.
+ \left[\sigma_s + 2\Delta t\alpha\right] & V_{i}^{k+1}
+ - \Delta t \alpha V_{i+1}^{k+1}
+ - \Delta t \alpha V_{i-1}^{k+1}
+ \nonumber \\
+ \left[\sigma_s/2 + \Delta t\alpha\right] & V_{1}^{k+1}
+ - \Delta t \alpha V_{2}^{k+1}
+ \nonumber \\
+ \left[\sigma_s/2 + \Delta t\alpha\right] & V_{n}^{k+1}
+ - \Delta t \alpha V_{n-1}^{k+1}
+ \nonumber
\end{align}
-where $\beta=\frac{a\Delta t}{2c_m r_L\Delta x^2}$.
The end points of the cable, i.e. the compartments for $x_1$ and $x_n$, have to be handled differently.
If we assume that a no-flux boundary condition, i.e. $\vv{J}\cdot\vv{n}=0$, is imposed at the end of the cable, the lhs of the linear system are
diff --git a/src/fvm.hpp b/src/fvm.hpp
index db19aaae..9e84119f 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -40,11 +40,28 @@ class fvm_cell {
fvm_cell(nest::mc::cell const& cell);
/// build the matrix for a given time step
- void setup_matrx(value_type dt);
-
- matrix_type& matrix()
+ void setup_matrix(value_type dt);
+
+ /// TODO this should be const
+ /// which requires const_view in the vector library
+ matrix_type& jacobian();
+
+ /// TODO this should be const
+ /// return list of CV areas in :
+ /// um^2
+ /// 1e-6.mm^2
+ /// 1e-8.cm^2
+ vector_view cv_areas();
+
+ /// TODO this should be const
+ /// return the capacitance of each CV surface
+ /// note that this is the total capacitance, not per unit area
+ /// which is equivalent to sigma_i * c_m
+ vector_view cv_capacitance();
+
+ std::size_t size() const
{
- return matrix_;
+ return matrix_.size();
}
private:
@@ -57,8 +74,21 @@ class fvm_cell {
/// alpha_[i] is the following value at the CV face between
/// CV i and its parent, required when constructing linear system
- /// alpha_[i] = area_face / (c_m * r_L * delta_x);
- vector_type alpha_;
+ /// face_alpha_[i] = area_face / (c_m * r_L * delta_x);
+ vector_type face_alpha_;
+
+ /// cv_capacitance_[i] is the capacitance of CV i per unit area (i.e. c_m)
+ vector_type cv_capacitance_;
+
+ /// the average current over the surface of each CV
+ /// current_ = i_m - i_e
+ /// so the total current over the surface of CV i is
+ /// current_[i] * cv_areas_
+ vector_type current_;
+
+ /// the potential in mV in each CV
+ vector_type voltage_;
+
};
////////////////////////////////////////////////////////////
@@ -67,9 +97,12 @@ class fvm_cell {
template <typename T, typename I>
fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
-: matrix_(cell.parent_index())
-, cv_areas_(matrix_.size())
-, alpha_(matrix_.size())
+: matrix_ {cell.parent_index()}
+, cv_areas_ {size(), T(0)}
+, face_alpha_ {size(), T(0)}
+, cv_capacitance_{size(), T(0)}
+, current_ {size(), T(0)}
+, voltage_ {size(), T(0)}
{
using util::left;
using util::right;
@@ -77,14 +110,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
auto parent_index = matrix_.p();
auto const& segment_index = cell.segment_index();
- // Use the membrane parameters for the first segment everywhere
- // in the cell to start with.
- // This has to be extended to use compartment/segment specific
- // membrane properties.
- auto membrane_params = cell.segments()[0]->mechanism("membrane");
- auto c_m = membrane_params.get("c_m").value;
- auto r_L = membrane_params.get("r_L").value;
-
auto seg_idx = 0;
for(auto const& s : cell.segments()) {
if(auto soma = s->as_soma()) {
@@ -94,7 +119,9 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
"FVM lowering encountered soma with non-zero index"
);
}
- cv_areas_[0] += math::area_sphere(soma->radius());
+ auto area = math::area_sphere(soma->radius());
+ cv_areas_[0] += area;
+ cv_capacitance_[0] += area * soma->mechanism("membrane").get("c_m").value;
}
else if(auto cable = s->as_cable()) {
// loop over each compartment in the cable
@@ -111,13 +138,18 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// (i.e. it follows the minimal degree ordering)
// The face is at the center, marked C.
// The full control volume to the left (marked with .)
+ auto c_m = cable->mechanism("membrane").get("c_m").value;
+ auto r_L = cable->mechanism("membrane").get("r_L").value;
for(auto c : cable->compartments()) {
auto i = segment_index[seg_idx] + c.index;
auto j = parent_index[i];
auto radius_center = math::mean(c.radius);
auto area_face = math::area_circle( radius_center );
- alpha_[i] = area_face / (c_m * r_L * c.length);
+ face_alpha_[i] = area_face / (c_m * r_L * c.length);
+ cv_capacitance_[i] = c_m;
+
+ std::cout << "radius " << radius_center << ", c_m " << c_m << ", r_L " << r_L << ", dx " << c.length << "\n";
auto halflen = c.length/2;
@@ -125,6 +157,8 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
cv_areas_[j] += al;
cv_areas_[i] += ar;
+ cv_capacitance_[j] += al * c_m;
+ cv_capacitance_[i] += ar * c_m;
}
}
else {
@@ -132,10 +166,36 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
}
++seg_idx;
}
+
+ // normalize the capacitance by cv_area
+ for(auto i=0u; i<size(); ++i) {
+ cv_capacitance_[i] /= cv_areas_[i];
+ }
+}
+
+template <typename T, typename I>
+typename fvm_cell<T,I>::matrix_type&
+fvm_cell<T,I>::jacobian()
+{
+ return matrix_;
}
template <typename T, typename I>
-void fvm_cell<T, I>::setup_matrx(T dt)
+typename fvm_cell<T,I>::vector_view
+fvm_cell<T,I>::cv_areas()
+{
+ return cv_areas_;
+}
+
+template <typename T, typename I>
+typename fvm_cell<T,I>::vector_view
+fvm_cell<T,I>::cv_capacitance()
+{
+ return cv_capacitance_;
+}
+
+template <typename T, typename I>
+void fvm_cell<T, I>::setup_matrix(T dt)
{
using memory::all;
@@ -144,6 +204,7 @@ void fvm_cell<T, I>::setup_matrx(T dt)
auto d = matrix_.d();
auto u = matrix_.u();
auto p = matrix_.p();
+ auto rhs = matrix_.rhs();
// The matrix has the following layout in memory
// where j is the parent index of i, i.e. i<j
@@ -157,13 +218,9 @@ void fvm_cell<T, I>::setup_matrx(T dt)
// . . .
// l[i] . . d[i]
//
-
- //d(all) = cv_areas_ + dt*(alpha_ + alpha_(p));
- //d[0] = cv_areas_[0];
-
d(all) = cv_areas_;
- for(auto i=1; i<d.size(); ++i) {
- auto a = dt * alpha_[i];
+ for(auto i=1u; i<d.size(); ++i) {
+ auto a = dt * face_alpha_[i];
d[i] += a;
l[i] = -a;
@@ -172,6 +229,12 @@ void fvm_cell<T, I>::setup_matrx(T dt)
// add contribution to the diagonal of parent
d[p[i]] += a;
}
+
+ // the RHS of the linear system is
+ // sigma_i * (V[i] - dt/cm*(im - ie))
+ for(auto i=0u; i<d.size(); ++i) {
+ rhs[i] = cv_areas_[i] * (voltage_[i] - dt/cv_capacitance_[i]*current_[i]);
+ }
}
} // namespace fvm
diff --git a/src/matrix.hpp b/src/matrix.hpp
index 4fbd3346..3e5f4442 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -57,7 +57,7 @@ class matrix {
}
/// the dimension of the matrix (i.e. the number of rows or colums)
- size_type size() const
+ std::size_t size() const
{
return parent_index_.size();
}
diff --git a/src/parameter_list.cpp b/src/parameter_list.cpp
index 65155100..86c641aa 100644
--- a/src/parameter_list.cpp
+++ b/src/parameter_list.cpp
@@ -85,8 +85,7 @@ namespace mc {
} // namespace mc
} // namespace nest
-/*
-static std::ostream&
+std::ostream&
operator<<(std::ostream& o, nest::mc::parameter const& p)
{
return o
@@ -97,7 +96,7 @@ operator<<(std::ostream& o, nest::mc::parameter const& p)
<< ")";
}
-static std::ostream&
+std::ostream&
operator<<(std::ostream& o, nest::mc::parameter_list const& l)
{
o << "parameters \"" << l.name() << "\" :\n";
@@ -106,4 +105,3 @@ operator<<(std::ostream& o, nest::mc::parameter_list const& l)
}
return o;
}
-*/
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index 74b9b1b3..e7dd1bd1 100644
--- a/src/parameter_list.hpp
+++ b/src/parameter_list.hpp
@@ -54,15 +54,6 @@ namespace mc {
value_type max;
};
- template <typename T>
- std::ostream& operator<<(std::ostream& o, value_range<T> const& r)
- {
- return
- o << "[ "
- << (r.has_lower_bound() ? std::to_string(r.min) : "-inf") << ", "
- << (r.has_upper_bound() ? std::to_string(r.max) : "inf") << "]";
- }
-
struct parameter {
using value_type = double;
using range_type = value_range<value_type>;
@@ -91,8 +82,6 @@ namespace mc {
range_type range;
};
- std::ostream& operator<<(std::ostream& o, parameter const& p);
-
// Use a dumb container class for now
// might have to use a more sophisticated interface in the future if need be
class parameter_list {
@@ -134,8 +123,6 @@ namespace mc {
};
- std::ostream& operator<<(std::ostream& o, parameter_list const& l);
-
///////////////////////////////////////////////////////////////////////////
// predefined parameter sets
///////////////////////////////////////////////////////////////////////////
@@ -159,8 +146,8 @@ namespace mc {
membrane_parameters()
: base("membrane")
{
- base::add_parameter({"r_L", 0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m2
- base::add_parameter({"c_m", 180.00, {0., 1e9}}); // Ohm.cm
+ base::add_parameter({"c_m", 0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m2
+ base::add_parameter({"r_L", 180.00, {0., 1e9}}); // Ohm.cm
}
};
@@ -192,3 +179,22 @@ namespace mc {
} // namespace mc
} // namespace nest
+template <typename T>
+std::ostream& operator<<(std::ostream& o, nest::mc::value_range<T> const& r)
+{
+ o << "[";
+ if(r.has_lower_bound())
+ o << r.min;
+ else
+ o<< "-inf";
+ o << ", ";
+ if(r.has_upper_bound())
+ o << r.max;
+ else
+ o<< "inf";
+ return o << "]";
+}
+
+std::ostream& operator<<(std::ostream& o, nest::mc::parameter const& p);
+std::ostream& operator<<(std::ostream& o, nest::mc::parameter_list const& l);
+
diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp
index 19e6702a..f0271bd4 100644
--- a/tests/test_matrix.cpp
+++ b/tests/test_matrix.cpp
@@ -16,7 +16,7 @@ TEST(matrix, construct_from_parent_only)
std::vector<int> p = {0,0,1};
matrix_type m{p};
EXPECT_EQ(m.num_cells(), 1);
- EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 3u);
auto mp = m.p();
@@ -31,9 +31,8 @@ TEST(matrix, construct_from_parent_only)
std::vector<int> p = {0,0,1};
matrix_type m{std::move(p)};
EXPECT_EQ(m.num_cells(), 1);
- EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 3u);
- EXPECT_EQ(m.size(), 3);
auto mp = m.p();
EXPECT_EQ(mp[0], 0);
@@ -48,7 +47,7 @@ TEST(matrix, construct_from_parent_only)
std::iota(p.begin()+1, p.end(), 0);
matrix_type m{p};
EXPECT_EQ(m.num_cells(), 1);
- EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 3u);
auto mp = m.p();
@@ -63,7 +62,7 @@ TEST(matrix, construct_from_parent_only)
std::iota(p.begin()+1, p.end(), 0);
matrix_type m{std::move(p)};
EXPECT_EQ(m.num_cells(), 1);
- EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(m.size(), 3u);
EXPECT_EQ(p.size(), 0u); // 0 implies moved from
auto mp = m.p();
@@ -98,7 +97,7 @@ TEST(matrix, solve)
std::iota(p.begin()+1, p.end(), 0);
matrix_type m{p};
- EXPECT_EQ(m.size(), (int)n);
+ EXPECT_EQ(m.size(), n);
EXPECT_EQ(m.num_cells(), 1);
m.d()(memory::all) = 2;
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 73e79292..712fca43 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -3,6 +3,56 @@
#include "../src/cell.hpp"
#include "../src/fvm.hpp"
+TEST(run, cable)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ cell.add_soma(6e-4); // 6um in cm
+
+ // 1um radius and 4mm long, all in cm
+ cell.add_cable(0, segmentKind::dendrite, 1e-4, 1e-4, 4e-1);
+
+ std::cout << cell.segment(1)->area() << " is the area\n";
+ EXPECT_EQ(cell.tree().num_segments(), 2u);
+
+ cell.soma()->add_mechanism(hh_parameters());
+
+ auto& soma_hh = cell.soma()->mechanism("hh");
+
+ soma_hh.set("gnabar", 0.12);
+ soma_hh.set("gkbar", 0.036);
+ soma_hh.set("gl", 0.0003);
+ soma_hh.set("el", -54.387);
+
+ cell.segment(1)->set_compartments(4);
+
+ using fvm_cell = fvm::fvm_cell<double, int>;
+ fvm_cell fvcell(cell);
+ auto& J = fvcell.jacobian();
+ EXPECT_EQ(J.size(), 5u);
+
+ fvcell.setup_matrix(0.02);
+ EXPECT_EQ(fvcell.cv_areas().size(), J.size());
+
+ auto& cable_parms = cell.segment(1)->mechanism("membrane");
+ std::cout << soma_hh << std::endl;
+ std::cout << cable_parms << std::endl;
+
+ std::cout << "l " << J.l() << "\n";
+ std::cout << "d " << J.d() << "\n";
+ std::cout << "u " << J.u() << "\n";
+ std::cout << "p " << J.p() << "\n";
+
+ J.rhs()(memory::all) = 1.;
+ J.rhs()[0] = 10.;
+
+ J.solve();
+
+ //std::cout << "x" << J.rhs() << "\n";
+}
+
TEST(run, init)
{
using namespace nest::mc;
@@ -38,11 +88,12 @@ TEST(run, init)
using fvm_cell = fvm::fvm_cell<double, int>;
fvm_cell fvcell(cell);
- EXPECT_EQ(fvcell.matrix().size(), 11);
+ auto& J = fvcell.jacobian();
+ EXPECT_EQ(J.size(), 11u);
- fvcell.setup_matrx(0.01);
+ fvcell.setup_matrix(0.01);
+ std::cout << "areas " << fvcell.cv_areas() << "\n";
- auto& J = fvcell.matrix();
//std::cout << "l" << J.l() << "\n";
//std::cout << "d" << J.d() << "\n";
//std::cout << "u" << J.u() << "\n";
--
GitLab