From e07f64fa4d2c7d2176b3f90c5eecb72bb3006c8d Mon Sep 17 00:00:00 2001
From: Sam Yates <halfflat@gmail.com>
Date: Tue, 7 May 2019 16:28:14 +0930
Subject: [PATCH] Add conductivity to implicit voltage solve step for
stability. (#735)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
Changes implicit solve step from:
> solve (c/δt + L) v' = c/δt v - J
to
> solve (c/δt + g + L)v' = (c/δt +g) v - J
where _c_ is capacitance, _g_ is membrane conductance, _J_ is membrane current.
* Compute conductivity contribution for mechanisms from symbolic d/dv of current contribution (extracted from linearity test.)
* Add new modcc 'source kind' for conductivity; tie to `vec_g_`.
* Add conductivity field to fvm shared state.
* Include conductivity in matrix assemblies for solution.
Fixes #633.
---
arbor/backends/gpu/matrix_assemble.cu | 28 ++++++++++++-------
arbor/backends/gpu/matrix_fine.cu | 11 +++++---
arbor/backends/gpu/matrix_fine.hpp | 1 +
arbor/backends/gpu/matrix_state_fine.hpp | 10 ++++---
arbor/backends/gpu/matrix_state_flat.hpp | 5 ++--
.../backends/gpu/matrix_state_interleaved.hpp | 8 ++++--
arbor/backends/gpu/mechanism.cpp | 1 +
arbor/backends/gpu/mechanism_ppack_base.hpp | 1 +
arbor/backends/gpu/shared_state.cpp | 4 +++
arbor/backends/gpu/shared_state.hpp | 25 +++++++++--------
arbor/backends/multicore/matrix_state.hpp | 17 ++++++-----
arbor/backends/multicore/mechanism.cpp | 1 +
arbor/backends/multicore/mechanism.hpp | 1 +
arbor/backends/multicore/shared_state.cpp | 4 +++
arbor/backends/multicore/shared_state.hpp | 15 +++++-----
arbor/fvm_lowered_cell_impl.hpp | 2 +-
arbor/matrix.hpp | 4 +--
modcc/identifier.hpp | 1 +
modcc/module.cpp | 23 +++++++++++++--
modcc/printer/printerutil.cpp | 3 ++
test/unit/test_matrix.cpp | 21 ++++++++------
test/unit/test_matrix.cu | 15 ++++++----
test/unit/test_matrix_cpuvsgpu.cpp | 6 ++--
23 files changed, 135 insertions(+), 72 deletions(-)
diff --git a/arbor/backends/gpu/matrix_assemble.cu b/arbor/backends/gpu/matrix_assemble.cu
index bcc589d2..4b998413 100644
--- a/arbor/backends/gpu/matrix_assemble.cu
+++ b/arbor/backends/gpu/matrix_assemble.cu
@@ -21,8 +21,9 @@ void assemble_matrix_flat(
const T* invariant_d,
const T* voltage,
const T* current,
+ const T* conductivity,
const T* cv_capacitance,
- const T* area,
+ const T* cv_area,
const I* cv_to_cell,
const T* dt_intdom,
const I* cell_to_intdom,
@@ -42,11 +43,12 @@ void assemble_matrix_flat(
// The 1e-3 is a constant of proportionality required to ensure that the
// conductance (gi) values have units μS (micro-Siemens).
// See the model documentation in docs/model for more information.
- T factor = 1e-3/dt;
+ T oodt_factor = 1e-3/dt; // [1/μs]
+ T area_factor = 1e-3*cv_area[tid]; // [1e-9·m²]
- auto gi = factor * cv_capacitance[tid];
+ auto gi = oodt_factor * cv_capacitance[tid] + area_factor*conductivity[tid]; // [μS]
d[tid] = gi + invariant_d[tid];
- rhs[tid] = gi*voltage[tid] - T(1e-3)*area[tid]*current[tid];
+ rhs[tid] = gi*voltage[tid] - area_factor*current[tid];
}
else {
d[tid] = 0;
@@ -69,6 +71,7 @@ void assemble_matrix_interleaved(
const T* invariant_d,
const T* voltage,
const T* current,
+ const T* conductivity,
const T* cv_capacitance,
const T* area,
const I* sizes,
@@ -82,6 +85,7 @@ void assemble_matrix_interleaved(
"number of threads must equal number of values to process per block");
__shared__ T buffer_v[Threads];
__shared__ T buffer_i[Threads];
+ __shared__ T buffer_g[Threads];
const unsigned tid = threadIdx.x + blockIdx.x*blockDim.x;
const unsigned lid = threadIdx.x;
@@ -103,7 +107,7 @@ void assemble_matrix_interleaved(
const unsigned max_size = sizes[0];
- T factor = 0;
+ T oodt_factor = 0;
T dt = 0;
const unsigned permuted_cid = blk_id*BlockWidth + blk_lane;
@@ -115,23 +119,25 @@ void assemble_matrix_interleaved(
// conductance (gi) values have units μS (micro-Siemens).
// See the model documentation in docs/model for more information.
- factor = dt>0? 1e-3/dt: 0;
+ oodt_factor = dt>0? T(1e-3)/dt: 0;
}
for (unsigned j=0u; j<max_size; j+=LoadWidth) {
if (do_load && load_pos<end) {
buffer_v[lid] = voltage[load_pos];
buffer_i[lid] = current[load_pos];
+ buffer_g[lid] = conductivity[load_pos];
}
__syncthreads();
if (j+blk_row<padded_size) {
- const auto gi = factor * cv_capacitance[store_pos];
+ T area_factor = T(1e-3)*area[store_pos];
+ const auto gi = oodt_factor*cv_capacitance[store_pos] + area_factor*buffer_g[blk_pos];
if (dt>0) {
d[store_pos] = (gi + invariant_d[store_pos]);
- rhs[store_pos] = (gi*buffer_v[blk_pos] - T(1e-3)*area[store_pos]*buffer_i[blk_pos]);
+ rhs[store_pos] = (gi*buffer_v[blk_pos] - area_factor*buffer_i[blk_pos]);
}
else {
d[store_pos] = 0;
@@ -154,6 +160,7 @@ void assemble_matrix_flat(
const fvm_value_type* invariant_d,
const fvm_value_type* voltage,
const fvm_value_type* current,
+ const fvm_value_type* conductivity,
const fvm_value_type* cv_capacitance,
const fvm_value_type* area,
const fvm_index_type* cv_to_cell,
@@ -167,7 +174,7 @@ void assemble_matrix_flat(
kernels::assemble_matrix_flat
<fvm_value_type, fvm_index_type>
<<<grid_dim, block_dim>>>
- (d, rhs, invariant_d, voltage, current, cv_capacitance,
+ (d, rhs, invariant_d, voltage, current, conductivity, cv_capacitance,
area, cv_to_cell, dt_intdom, cell_to_intdom, n);
}
@@ -178,6 +185,7 @@ void assemble_matrix_interleaved(
const fvm_value_type* invariant_d,
const fvm_value_type* voltage,
const fvm_value_type* current,
+ const fvm_value_type* conductivity,
const fvm_value_type* cv_capacitance,
const fvm_value_type* area,
const fvm_index_type* sizes,
@@ -197,7 +205,7 @@ void assemble_matrix_interleaved(
kernels::assemble_matrix_interleaved
<fvm_value_type, fvm_index_type, bd, lw, block_dim>
<<<grid_dim, block_dim>>>
- (d, rhs, invariant_d, voltage, current, cv_capacitance, area,
+ (d, rhs, invariant_d, voltage, current, conductivity, cv_capacitance, area,
sizes, starts, matrix_to_cell,
dt_intdom, cell_to_intdom, padded_size, num_mtx);
}
diff --git a/arbor/backends/gpu/matrix_fine.cu b/arbor/backends/gpu/matrix_fine.cu
index f3ad5372..6fe430df 100644
--- a/arbor/backends/gpu/matrix_fine.cu
+++ b/arbor/backends/gpu/matrix_fine.cu
@@ -50,6 +50,7 @@ void assemble_matrix_fine(
const T* invariant_d,
const T* voltage,
const T* current,
+ const T* conductivity,
const T* cv_capacitance,
const T* area,
const I* cv_to_cell,
@@ -68,12 +69,13 @@ void assemble_matrix_fine(
// The 1e-3 is a constant of proportionality required to ensure that the
// conductance (gi) values have units μS (micro-Siemens).
// See the model documentation in docs/model for more information.
- T factor = 1e-3/dt;
+ T oodt_factor = T(1e-3)/dt;
+ T area_factor = T(1e-3)*area[tid];
- const auto gi = factor * cv_capacitance[tid];
+ const auto gi = oodt_factor*cv_capacitance[tid] + area_factor*conductivity[tid];
const auto pid = perm[tid];
d[pid] = gi + invariant_d[tid];
- rhs[pid] = gi*voltage[tid] - T(1e-3)*area[tid]*current[tid];
+ rhs[pid] = gi*voltage[tid] - area_factor*current[tid];
}
else {
const auto pid = perm[tid];
@@ -262,6 +264,7 @@ void assemble_matrix_fine(
const fvm_value_type* invariant_d,
const fvm_value_type* voltage,
const fvm_value_type* current,
+ const fvm_value_type* conductivity,
const fvm_value_type* cv_capacitance,
const fvm_value_type* area,
const fvm_index_type* cv_to_cell,
@@ -274,7 +277,7 @@ void assemble_matrix_fine(
const unsigned num_blocks = impl::block_count(n, block_dim);
kernels::assemble_matrix_fine<<<num_blocks, block_dim>>>(
- d, rhs, invariant_d, voltage, current, cv_capacitance, area,
+ d, rhs, invariant_d, voltage, current, conductivity, cv_capacitance, area,
cv_to_cell, dt_intdom, cell_to_intdom,
perm, n);
}
diff --git a/arbor/backends/gpu/matrix_fine.hpp b/arbor/backends/gpu/matrix_fine.hpp
index 896bb28b..1a15bd7f 100644
--- a/arbor/backends/gpu/matrix_fine.hpp
+++ b/arbor/backends/gpu/matrix_fine.hpp
@@ -55,6 +55,7 @@ void assemble_matrix_fine(
const fvm_value_type* invariant_d,
const fvm_value_type* voltage,
const fvm_value_type* current,
+ const fvm_value_type* conductivity,
const fvm_value_type* cv_capacitance,
const fvm_value_type* area,
const fvm_index_type* cv_to_cell,
diff --git a/arbor/backends/gpu/matrix_state_fine.hpp b/arbor/backends/gpu/matrix_state_fine.hpp
index f5c44cbe..87ebf2ae 100644
--- a/arbor/backends/gpu/matrix_state_fine.hpp
+++ b/arbor/backends/gpu/matrix_state_fine.hpp
@@ -437,17 +437,19 @@ public:
}
// Assemble the matrix
- // Afterwards the diagonal and RHS will have been set given dt, voltage and current
- // dt_intdom [ms] (per cell)
+ // Afterwards the diagonal and RHS will have been set given dt, voltage, current, and conductivity.
+ // dt_intdom [ms] (per integration domain)
// voltage [mV]
- // current [nA]
- void assemble(const_view dt_intdom, const_view voltage, const_view current) {
+ // current density [A/m²]
+ // conductivity [kS/m²]
+ void assemble(const_view dt_intdom, const_view voltage, const_view current, const_view conductivity) {
assemble_matrix_fine(
d.data(),
rhs.data(),
invariant_d.data(),
voltage.data(),
current.data(),
+ conductivity.data(),
cv_capacitance.data(),
cv_area.data(),
cv_to_cell.data(),
diff --git a/arbor/backends/gpu/matrix_state_flat.hpp b/arbor/backends/gpu/matrix_state_flat.hpp
index 6f4cbd01..df1c8c84 100644
--- a/arbor/backends/gpu/matrix_state_flat.hpp
+++ b/arbor/backends/gpu/matrix_state_flat.hpp
@@ -27,6 +27,7 @@ void assemble_matrix_flat(
const fvm_value_type* invariant_d,
const fvm_value_type* voltage,
const fvm_value_type* current,
+ const fvm_value_type* conductivity,
const fvm_value_type* cv_capacitance,
const fvm_value_type* cv_area,
const fvm_index_type* cv_to_cell,
@@ -123,11 +124,11 @@ struct matrix_state_flat {
// dt_intdom [ms] (per integration domain)
// voltage [mV]
// current [nA]
- void assemble(const_view dt_intdom, const_view voltage, const_view current) {
+ void assemble(const_view dt_intdom, const_view voltage, const_view current, const_view conductance) {
// perform assembly on the gpu
assemble_matrix_flat(
d.data(), rhs.data(), invariant_d.data(), voltage.data(),
- current.data(), cv_capacitance.data(), cv_area.data(),
+ current.data(), conductance.data(), cv_capacitance.data(), cv_area.data(),
cv_to_cell.data(), dt_intdom.data(), cell_to_intdom.data(), size());
}
diff --git a/arbor/backends/gpu/matrix_state_interleaved.hpp b/arbor/backends/gpu/matrix_state_interleaved.hpp
index 864b77eb..89e63cb7 100644
--- a/arbor/backends/gpu/matrix_state_interleaved.hpp
+++ b/arbor/backends/gpu/matrix_state_interleaved.hpp
@@ -23,6 +23,7 @@ void assemble_matrix_interleaved(
const fvm_value_type* invariant_d,
const fvm_value_type* voltage,
const fvm_value_type* current,
+ const fvm_value_type* conductivity,
const fvm_value_type* cv_capacitance,
const fvm_value_type* area,
const fvm_index_type* sizes,
@@ -266,11 +267,12 @@ struct matrix_state_interleaved {
// Afterwards the diagonal and RHS will have been set given dt, voltage and current.
// dt_intdom [ms] (per integration domain)
// voltage [mV] (per compartment)
- // current density [A.m^-2] (per compartment)
- void assemble(const_view dt_intdom, const_view voltage, const_view current) {
+ // current density [A/m²] (per compartment)
+ // conductivity [kS/m²] (per compartment)
+ void assemble(const_view dt_intdom, const_view voltage, const_view current, const_view conductivity) {
assemble_matrix_interleaved
(d.data(), rhs.data(), invariant_d.data(),
- voltage.data(), current.data(), cv_capacitance.data(), cv_area.data(),
+ voltage.data(), current.data(), conductivity.data(), cv_capacitance.data(), cv_area.data(),
matrix_sizes.data(), matrix_index.data(),
matrix_to_cell_index.data(),
dt_intdom.data(), cell_to_intdom.data(), padded_matrix_size(), num_matrices());
diff --git a/arbor/backends/gpu/mechanism.cpp b/arbor/backends/gpu/mechanism.cpp
index d73e46e9..f1c74225 100644
--- a/arbor/backends/gpu/mechanism.cpp
+++ b/arbor/backends/gpu/mechanism.cpp
@@ -76,6 +76,7 @@ void mechanism::instantiate(unsigned id,
pp->vec_v_ = shared.voltage.data();
pp->vec_i_ = shared.current_density.data();
+ pp->vec_g_ = shared.conductivity.data();
pp->temperature_degC_ = shared.temperature_degC.data();
diff --git a/arbor/backends/gpu/mechanism_ppack_base.hpp b/arbor/backends/gpu/mechanism_ppack_base.hpp
index 16518c09..457f3eeb 100644
--- a/arbor/backends/gpu/mechanism_ppack_base.hpp
+++ b/arbor/backends/gpu/mechanism_ppack_base.hpp
@@ -35,6 +35,7 @@ struct mechanism_ppack_base {
const value_type* vec_dt_;
const value_type* vec_v_;
value_type* vec_i_;
+ value_type* vec_g_;
const value_type* temperature_degC_;
const index_type* node_index_;
diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp
index fe499023..1f8c7331 100644
--- a/arbor/backends/gpu/shared_state.cpp
+++ b/arbor/backends/gpu/shared_state.cpp
@@ -126,6 +126,7 @@ shared_state::shared_state(
dt_cv(n_cv),
voltage(n_cv),
current_density(n_cv),
+ conductivity(n_cv),
temperature_degC(1),
deliverable_events(n_intdom)
{}
@@ -144,6 +145,7 @@ void shared_state::add_ion(
void shared_state::reset(fvm_value_type initial_voltage, fvm_value_type temperature_K) {
memory::fill(voltage, initial_voltage);
memory::fill(current_density, 0);
+ memory::fill(conductivity, 0);
memory::fill(time, 0);
memory::fill(time_to, 0);
memory::fill(temperature_degC, temperature_K - 273.15);
@@ -155,6 +157,7 @@ void shared_state::reset(fvm_value_type initial_voltage, fvm_value_type temperat
void shared_state::zero_currents() {
memory::fill(current_density, 0);
+ memory::fill(conductivity, 0);
for (auto& i: ion_data) {
i.second.zero_current();
}
@@ -205,6 +208,7 @@ std::ostream& operator<<(std::ostream& o, shared_state& s) {
o << " dt_cv " << s.dt_cv << "\n";
o << " voltage " << s.voltage << "\n";
o << " current " << s.current_density << "\n";
+ o << " conductivity " << s.conductivity << "\n";
for (auto& ki: s.ion_data) {
auto kn = to_string(ki.first);
auto& i = const_cast<ion_state&>(ki.second);
diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp
index 778d82b2..fe662720 100644
--- a/arbor/backends/gpu/shared_state.hpp
+++ b/arbor/backends/gpu/shared_state.hpp
@@ -66,18 +66,19 @@ struct ion_state {
struct shared_state {
fvm_size_type n_intdom = 0; // Number of distinct integration domains.
- fvm_size_type n_cv = 0; // Total number of CVs.
- fvm_size_type n_gj = 0; // Total number of GJs.
-
- iarray cv_to_intdom; // Maps CV index to intdom index.
- gjarray gap_junctions; // Stores gap_junction info.
- array time; // Maps intdom index to integration start time [ms].
- array time_to; // Maps intdom index to integration stop time [ms].
- array dt_intdom; // Maps intdom index to (stop time) - (start time) [ms].
- array dt_cv; // Maps CV index to dt [ms].
- array voltage; // Maps CV index to membrane voltage [mV].
- array current_density; // Maps CV index to current density [A/m²].
- array temperature_degC; // Global temperature [°C] (length 1 array).
+ fvm_size_type n_cv = 0; // Total number of CVs.
+ fvm_size_type n_gj = 0; // Total number of GJs.
+
+ iarray cv_to_intdom; // Maps CV index to intdom index.
+ gjarray gap_junctions; // Stores gap_junction info.
+ array time; // Maps intdom index to integration start time [ms].
+ array time_to; // Maps intdom index to integration stop time [ms].
+ array dt_intdom; // Maps intdom index to (stop time) - (start time) [ms].
+ array dt_cv; // Maps CV index to dt [ms].
+ array voltage; // Maps CV index to membrane voltage [mV].
+ array current_density; // Maps CV index to current density [A/m²].
+ array conductivity; // Maps CV index to membrane conductivity [kS/m²].
+ array temperature_degC; // Global temperature [°C] (length 1 array).
std::unordered_map<ionKind, ion_state> ion_data;
diff --git a/arbor/backends/multicore/matrix_state.hpp b/arbor/backends/multicore/matrix_state.hpp
index 9e54b8eb..7cba8e96 100644
--- a/arbor/backends/multicore/matrix_state.hpp
+++ b/arbor/backends/multicore/matrix_state.hpp
@@ -74,10 +74,11 @@ public:
// Assemble the matrix
// Afterwards the diagonal and RHS will have been set given dt, voltage and current.
- // dt_intdom [ms] (per integration domain)
- // voltage [mV] (per compartment)
- // current density [A.m^-2] (per compartment)
- void assemble(const_view dt_intdom, const_view voltage, const_view current) {
+ // dt_intdom [ms] (per integration domain)
+ // voltage [mV] (per control volume)
+ // current density [A.m^-2] (per control volume)
+ // conductivity [kS.m^-2] (per control volume)
+ void assemble(const_view dt_intdom, const_view voltage, const_view current, const_view conductivity) {
auto cell_cv_part = util::partition_view(cell_cv_divs);
const index_type ncells = cell_cv_part.size();
@@ -86,13 +87,15 @@ public:
auto dt = dt_intdom[cell_to_intdom[m]];
if (dt>0) {
- value_type factor = 1e-3/dt;
+ value_type oodt_factor = 1e-3/dt; // [1/µs]
for (auto i: util::make_span(cell_cv_part[m])) {
- auto gi = factor*cv_capacitance[i];
+ auto area_factor = 1e-3*cv_area[i]; // [1e-9·m²]
+
+ auto gi = oodt_factor*cv_capacitance[i] + area_factor*conductivity[i]; // [μS]
d[i] = gi + invariant_d[i];
// convert current to units nA
- rhs[i] = gi*voltage[i] - 1e-3*cv_area[i]*current[i];
+ rhs[i] = gi*voltage[i] - area_factor*current[i];
}
}
else {
diff --git a/arbor/backends/multicore/mechanism.cpp b/arbor/backends/multicore/mechanism.cpp
index 97aa457e..fde15d3e 100644
--- a/arbor/backends/multicore/mechanism.cpp
+++ b/arbor/backends/multicore/mechanism.cpp
@@ -78,6 +78,7 @@ void mechanism::instantiate(unsigned id, backend::shared_state& shared, const la
vec_v_ = shared.voltage.data();
vec_i_ = shared.current_density.data();
+ vec_g_ = shared.conductivity.data();
temperature_degC_ = &shared.temperature_degC;
diff --git a/arbor/backends/multicore/mechanism.hpp b/arbor/backends/multicore/mechanism.hpp
index 21493e94..f3819ab1 100644
--- a/arbor/backends/multicore/mechanism.hpp
+++ b/arbor/backends/multicore/mechanism.hpp
@@ -79,6 +79,7 @@ protected:
const value_type* vec_dt_; // CV to integration time step.
const value_type* vec_v_; // CV to cell membrane voltage.
value_type* vec_i_; // CV to cell membrane current density.
+ value_type* vec_g_; // CV to cell membrane conductivity.
const value_type* temperature_degC_; // Pointer to global temperature scalar.
deliverable_event_stream* event_stream_ptr_;
diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp
index 6313d429..6af8064d 100644
--- a/arbor/backends/multicore/shared_state.cpp
+++ b/arbor/backends/multicore/shared_state.cpp
@@ -133,6 +133,7 @@ shared_state::shared_state(
dt_cv(n_cv, pad(alignment)),
voltage(n_cv, pad(alignment)),
current_density(n_cv, pad(alignment)),
+ conductivity(n_cv, pad(alignment)),
temperature_degC(NAN),
deliverable_events(n_intdom)
{
@@ -161,6 +162,7 @@ void shared_state::add_ion(
void shared_state::reset(fvm_value_type initial_voltage, fvm_value_type temperature_K) {
util::fill(voltage, initial_voltage);
util::fill(current_density, 0);
+ util::fill(conductivity, 0);
util::fill(time, 0);
util::fill(time_to, 0);
temperature_degC = temperature_K - 273.15;
@@ -172,6 +174,7 @@ void shared_state::reset(fvm_value_type initial_voltage, fvm_value_type temperat
void shared_state::zero_currents() {
util::fill(current_density, 0);
+ util::fill(conductivity, 0);
for (auto& i: ion_data) {
i.second.zero_current();
}
@@ -262,6 +265,7 @@ std::ostream& operator<<(std::ostream& out, const shared_state& s) {
out << "dt_cv " << csv(s.dt_cv) << "\n";
out << "voltage " << csv(s.voltage) << "\n";
out << "current " << csv(s.current_density) << "\n";
+ out << "conductivity " << csv(s.conductivity) << "\n";
for (auto& ki: s.ion_data) {
auto kn = to_string(ki.first);
auto& i = const_cast<ion_state&>(ki.second);
diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp
index fc2c0e7f..be3b5725 100644
--- a/arbor/backends/multicore/shared_state.hpp
+++ b/arbor/backends/multicore/shared_state.hpp
@@ -88,14 +88,15 @@ struct shared_state {
fvm_size_type n_cv = 0; // Total number of CVs.
fvm_size_type n_gj = 0; // Total number of GJs.
- iarray cv_to_intdom; // Maps CV index to integration domain index.
+ iarray cv_to_intdom; // Maps CV index to integration domain index.
gjarray gap_junctions; // Stores gap_junction info.
- array time; // Maps intdom index to integration start time [ms].
- array time_to; // Maps intdom index to integration stop time [ms].
- array dt_intdom; // Maps index to (stop time) - (start time) [ms].
- array dt_cv; // Maps CV index to dt [ms].
- array voltage; // Maps CV index to membrane voltage [mV].
- array current_density; // Maps CV index to current density [A/m²].
+ array time; // Maps intdom index to integration start time [ms].
+ array time_to; // Maps intdom index to integration stop time [ms].
+ array dt_intdom; // Maps index to (stop time) - (start time) [ms].
+ array dt_cv; // Maps CV index to dt [ms].
+ array voltage; // Maps CV index to membrane voltage [mV].
+ array current_density; // Maps CV index to membrane current density contributions [A/m²].
+ array conductivity; // Maps CV index to membrane conductivity [kS/m²].
fvm_value_type temperature_degC; // Global temperature [°C].
std::unordered_map<ionKind, ion_state> ion_data;
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 1bd3889d..c3e1d8da 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -237,7 +237,7 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
// Integrate voltage by matrix solve.
PE(advance_integrate_matrix_build);
- matrix_.assemble(state_->dt_intdom, state_->voltage, state_->current_density);
+ matrix_.assemble(state_->dt_intdom, state_->voltage, state_->current_density, state_->conductivity);
PL();
PE(advance_integrate_matrix_solve);
matrix_.solve();
diff --git a/arbor/matrix.hpp b/arbor/matrix.hpp
index 94958c00..8609cb7c 100644
--- a/arbor/matrix.hpp
+++ b/arbor/matrix.hpp
@@ -67,8 +67,8 @@ public:
}
/// Assemble the matrix for given dt
- void assemble(const array& dt_cell, const array& voltage, const array& current) {
- state_.assemble(dt_cell, voltage, current);
+ void assemble(const array& dt_cell, const array& voltage, const array& current, const array& conductivity) {
+ state_.assemble(dt_cell, voltage, current, conductivity);
}
/// Get a view of the solution
diff --git a/modcc/identifier.hpp b/modcc/identifier.hpp
index 3ddd3c56..891e6f49 100644
--- a/modcc/identifier.hpp
+++ b/modcc/identifier.hpp
@@ -46,6 +46,7 @@ enum class ionKind {
enum class sourceKind {
voltage,
current,
+ conductivity,
dt,
ion_current,
ion_revpot,
diff --git a/modcc/module.cpp b/modcc/module.cpp
index 93e5ac7c..c78d9c72 100644
--- a/modcc/module.cpp
+++ b/modcc/module.cpp
@@ -45,16 +45,25 @@ public:
virtual void finalize() override {
if (has_current_update_) {
- // Initialize current_ as first statement.
+ // Initialize current_ and conductivity_ as first statements.
+ statements_.push_front(make_expression<AssignmentExpression>(loc_,
+ id("conductivity_"),
+ make_expression<NumberExpression>(loc_, 0.0)));
statements_.push_front(make_expression<AssignmentExpression>(loc_,
id("current_"),
make_expression<NumberExpression>(loc_, 0.0)));
+ // Scale current and conductivity contributions by weight.
statements_.push_back(make_expression<AssignmentExpression>(loc_,
id("current_"),
make_expression<MulBinaryExpression>(loc_,
id("weight_"),
id("current_"))));
+ statements_.push_back(make_expression<AssignmentExpression>(loc_,
+ id("conductivity_"),
+ make_expression<MulBinaryExpression>(loc_,
+ id("weight_"),
+ id("conductivity_"))));
for (auto& v: ion_current_vars_) {
statements_.push_back(make_expression<AssignmentExpression>(loc_,
@@ -79,7 +88,8 @@ public:
}
has_current_update_ = true;
- if (!linear_test(e->rhs(), {"v"}).is_linear) {
+ linear_test_result L = linear_test(e->rhs(), {"v"});
+ if (!L.is_linear) {
error({"current update expressions must be linear in v: "+e->rhs()->to_string(),
e->location()});
return;
@@ -90,6 +100,13 @@ public:
make_expression<AddBinaryExpression>(loc,
id("current_", loc),
e->lhs()->clone())));
+ if (L.coef.count("v")) {
+ statements_.push_back(make_expression<AssignmentExpression>(loc,
+ id("conductivity_", loc),
+ make_expression<AddBinaryExpression>(loc,
+ id("conductivity_", loc),
+ L.coef.at("v")->clone())));
+ }
}
}
}
@@ -492,6 +509,8 @@ void Module::add_variables_to_symbols() {
create_indexed_variable("current_", "vec_i", sourceKind::current, tok::plus,
accessKind::write, ionKind::none, Location());
+ create_indexed_variable("conductivity_", "vec_g", sourceKind::conductivity, tok::plus,
+ accessKind::write, ionKind::none, Location());
create_indexed_variable("v", "vec_v", sourceKind::voltage, tok::eq,
accessKind::read, ionKind::none, Location());
create_indexed_variable("dt", "vec_dt", sourceKind::dt, tok::eq,
diff --git a/modcc/printer/printerutil.cpp b/modcc/printer/printerutil.cpp
index 09bc4c73..f5f95cf5 100644
--- a/modcc/printer/printerutil.cpp
+++ b/modcc/printer/printerutil.cpp
@@ -129,6 +129,9 @@ indexed_variable_info decode_indexed_variable(IndexedVariable* sym) {
case sourceKind::current:
data_var="vec_i_";
break;
+ case sourceKind::conductivity:
+ data_var="vec_g_";
+ break;
case sourceKind::dt:
data_var="vec_dt_";
break;
diff --git a/test/unit/test_matrix.cpp b/test/unit/test_matrix.cpp
index eb72dacf..6ead914f 100644
--- a/test/unit/test_matrix.cpp
+++ b/test/unit/test_matrix.cpp
@@ -144,18 +144,21 @@ TEST(matrix, zero_diagonal_assembled)
// Intial voltage of zero; currents alone determine rhs.
array v(7, 0.0);
vvec area(7, 1.0);
- array i = {-3000, -5000, -7000, -6000, -9000, -16000, -32000};
+
+ // (Scaled) membrane conductances contribute to diagonal.
+ array mg = { 1000, 2000, 3000, 4000, 5000, 6000, 7000 };
+ array i = {-7000, -15000, -25000, -34000, -49000, -70000, -102000};
// Expected matrix and rhs:
// u = [ 0 -1 -1 0 -1 0 -2]
- // d = [ 2 3 2 2 2 4 5]
- // rhs = [ 3 5 7 2 4 16 32]
+ // d = [ 3 5 5 6 7 10 12]
+ // rhs = [ 7 15 25 34 49 70 102 ]
//
// Expected solution:
- // x = [ 4 5 6 7 8 9 10]
+ // x = [ 4 5 6 7 8 9 10 ]
matrix_type m(p, c, Cm, g, area, s);
- m.assemble(dt, v, i);
+ m.assemble(dt, v, i, mg);
m.solve();
vvec x;
@@ -168,13 +171,13 @@ TEST(matrix, zero_diagonal_assembled)
// should then return voltage values for that submatrix.
dt[1] = 0;
- v[3] = 20;
- v[4] = 30;
- m.assemble(dt, v, i);
+ v[3] = -20;
+ v[4] = -30;
+ m.assemble(dt, v, i, mg);
m.solve();
assign(x, m.solution());
- expected = {4, 5, 6, 20, 30, 9, 10};
+ expected = {4, 5, 6, -20, -30, 9, 10};
EXPECT_TRUE(testing::seq_almost_eq<double>(expected, x));
}
diff --git a/test/unit/test_matrix.cu b/test/unit/test_matrix.cu
index 8d01632d..706f3a34 100644
--- a/test/unit/test_matrix.cu
+++ b/test/unit/test_matrix.cu
@@ -270,8 +270,8 @@ TEST(matrix, backends)
auto group_size = cell_cv_divs.back();
- // Build the capacitance and conductance vectors and
- // populate with nonzero random values
+ // Build the capacitance, (axial) conductance, voltage, current density,
+ // and membrane conductance vectors. Populate them with nonzero random values.
auto gen = std::mt19937();
gen.seed(100);
auto dist = std::uniform_real_distribution<T>(1, 200);
@@ -280,12 +280,14 @@ TEST(matrix, backends)
std::vector<T> g(group_size);
std::vector<T> v(group_size);
std::vector<T> i(group_size);
+ std::vector<T> mg(group_size);
std::vector<T> area(group_size, 1e3);
std::generate(Cm.begin(), Cm.end(), [&](){return dist(gen);});
std::generate(g.begin(), g.end(), [&](){return dist(gen);});
std::generate(v.begin(), v.end(), [&](){return dist(gen);});
std::generate(i.begin(), i.end(), [&](){return dist(gen);});
+ std::generate(mg.begin(), mg.end(), [&](){return dist(gen);});
// Make the reference matrix and the gpu matrix
auto flat = state_flat(p, cell_cv_divs, Cm, g, area, cell_to_intdom); // flat
@@ -298,14 +300,15 @@ TEST(matrix, backends)
auto dt_dist = std::uniform_real_distribution<T>(0.01, 0.02);
std::generate(dt.begin(), dt.end(), [&](){return dt_dist(gen);});
- // Voltage and current values.
+ // Voltage, current, and membrane conductance values.
auto gpu_dt = on_gpu(dt);
auto gpu_v = on_gpu(v);
auto gpu_i = on_gpu(i);
+ auto gpu_mg = on_gpu(mg);
- flat.assemble(gpu_dt, gpu_v, gpu_i);
- intl.assemble(gpu_dt, gpu_v, gpu_i);
- fine.assemble(gpu_dt, gpu_v, gpu_i);
+ flat.assemble(gpu_dt, gpu_v, gpu_i, gpu_mg);
+ intl.assemble(gpu_dt, gpu_v, gpu_i, gpu_mg);
+ fine.assemble(gpu_dt, gpu_v, gpu_i, gpu_mg);
flat.solve();
intl.solve();
diff --git a/test/unit/test_matrix_cpuvsgpu.cpp b/test/unit/test_matrix_cpuvsgpu.cpp
index 9d90e062..f89bb900 100644
--- a/test/unit/test_matrix_cpuvsgpu.cpp
+++ b/test/unit/test_matrix_cpuvsgpu.cpp
@@ -123,10 +123,10 @@ TEST(matrix, assemble)
auto dt_dist = std::uniform_real_distribution<T>(0.1, 0.2);
std::generate(dt.begin(), dt.end(), [&](){return dt_dist(gen);});
- // Voltage and current values
- m_mc.assemble(host_array(dt.begin(), dt.end()), host_array(group_size, -64), host_array(group_size, 10));
+ // Voltage, current, and conductance values
+ m_mc.assemble(host_array(dt.begin(), dt.end()), host_array(group_size, -64), host_array(group_size, 10), host_array(group_size, 3));
m_mc.solve();
- m_gpu.assemble(on_gpu(dt), gpu_array(group_size, -64), gpu_array(group_size, 10));
+ m_gpu.assemble(on_gpu(dt), gpu_array(group_size, -64), gpu_array(group_size, 10), gpu_array(group_size, 3));
m_gpu.solve();
// Compare the GPU and CPU results.
--
GitLab