Skip to content

Commit

Permalink
Rip out concentration change by current.
Browse files Browse the repository at this point in the history
  • Loading branch information
@thorstenhater
thorstenhater committed Oct 16, 2023
1 parent 124833e commit b342ab1
Show file tree
Hide file tree
Showing 12 changed files with 15 additions and 101 deletions.
24 changes: 3 additions & 21 deletions arbor/backends/gpu/diffusion.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,6 @@ void assemble_diffusion(T* __restrict__ const d,
T* __restrict__ const rhs,
const T* __restrict__ const invariant_d,
const T* __restrict__ const concentration,
const T* __restrict__ const voltage,
const T* __restrict__ const current,
const T q,
const T* __restrict__ const conductivity,
const T* __restrict__ const area,
const T* __restrict__ const volume,
const T dt,
const I* __restrict__ const perm,
Expand All @@ -33,18 +28,10 @@ void assemble_diffusion(T* __restrict__ const d,
if (tid < n) {
auto _1_dt = 1e-3/dt; // 1/µs
auto pid = perm[tid];
auto u = voltage[tid]; // mV
auto g = conductivity[tid]; // µS
auto J = current[tid]; // A/m^2
auto A = 1e-3*area[tid]; // 1e-9·m²
auto V = volume[tid]; // um^3
auto X = concentration[tid]; // mM
// conversion from current density to concentration change
// using Faraday's constant
auto F = A/(q*96.485332);

d[pid] = _1_dt*V + F*g + invariant_d[tid];
rhs[pid] = _1_dt*V*X + F*(u*g - J);
d[pid] = _1_dt*V + invariant_d[tid];
rhs[pid] = _1_dt*V*X;
}
}

Expand Down Expand Up @@ -199,19 +186,14 @@ ARB_ARBOR_API void assemble_diffusion(arb_value_type* d,
arb_value_type* rhs,
const arb_value_type* invariant_d,
const arb_value_type* concentration,
const arb_value_type* voltage,
const arb_value_type* current,
arb_value_type q,
const arb_value_type* conductivity,
const arb_value_type* area,
const arb_value_type* volume,
const arb_value_type dt,
const arb_index_type* perm,
unsigned n) {
launch_1d(n,
128,
kernels::assemble_diffusion<arb_value_type, arb_index_type>,
d, rhs, invariant_d, concentration, voltage, current, q, conductivity, area, volume, dt, perm, n);
d, rhs, invariant_d, volume, dt, perm, n);
}

// Example:
Expand Down
20 changes: 3 additions & 17 deletions arbor/backends/gpu/diffusion_state.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,6 @@ struct diffusion_state {
array rhs; // [nA]

// Required for matrix assembly
array cv_area; // [μm^2]
array cv_volume; // [μm^3]

// Invariant part of the matrix diagonal
Expand Down Expand Up @@ -356,7 +355,6 @@ struct diffusion_state {
// d, u, rhs : packed
// invariant_d : flat
// cv_to_cell : flat
// area : flat

// the invariant part of d is stored in in flat form
std::vector<value_type> invariant_d_tmp(matrix_size, 0);
Expand All @@ -383,7 +381,6 @@ struct diffusion_state {
flat_to_packed(u_shuffled, u);

// the invariant part of d and cv_area are in flat form
cv_area = memory::make_const_view(area);
cv_volume = memory::make_const_view(volume);
invariant_d = memory::make_const_view(invariant_d_tmp);

Expand All @@ -400,18 +397,11 @@ struct diffusion_state {
// Afterwards the diagonal and RHS will have been set given dt, voltage, current, and conductivity.
// dt [ms] (scalar)
// voltage [mV]
// current density [A/m²]
// conductivity [kS/m²]
void assemble(const value_type dt, const_view concentration, const_view voltage, const_view current, const_view conductivity, arb_value_type q) {
void assemble(const value_type dt, const_view concentration) {
assemble_diffusion(d.data(),
rhs.data(),
invariant_d.data(),
concentration.data(),
voltage.data(),
current.data(),
q,
conductivity.data(),
cv_area.data(),
cv_volume.data(),
dt,
perm.data(),
Expand All @@ -435,12 +425,8 @@ struct diffusion_state {
}

void solve(array& concentration,
const value_type dt,
const_view voltage,
const_view current,
const_view conductivity,
arb_value_type q) {
assemble(dt, concentration, voltage, current, conductivity, q);
const value_type dt) {
assemble(dt, concentration);
solve(concentration);
}

Expand Down
3 changes: 0 additions & 3 deletions arbor/backends/gpu/shared_state.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@ ion_state::ion_state(const fvm_ion_config& ion_data,
Xi_(ion_data.init_iconc.begin(), ion_data.init_iconc.end()),
Xd_(ion_data.cv.size(), NAN),
Xo_(ion_data.init_econc.begin(), ion_data.init_econc.end()),
gX_(ion_data.cv.size(), NAN),
init_Xi_(make_const_view(ion_data.init_iconc)),
init_Xo_(make_const_view(ion_data.init_econc)),
reset_Xi_(make_const_view(ion_data.reset_iconc)),
Expand All @@ -77,7 +76,6 @@ void ion_state::init_concentration() {
}

void ion_state::zero_current() {
memory::fill(gX_, 0);
memory::fill(iX_, 0);
}

Expand Down Expand Up @@ -263,7 +261,6 @@ void shared_state::instantiate(mechanism& m,
ion_state.external_concentration = oion->Xo_.data();
ion_state.diffusive_concentration = oion->Xd_.data();
ion_state.ionic_charge = oion->charge.data();
ion_state.conductivity = oion->gX_.data();
}

// If there are no sites (is this ever meaningful?) there is nothing more to do.
Expand Down
3 changes: 1 addition & 2 deletions arbor/backends/gpu/shared_state.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ struct ARB_ARBOR_API ion_state {
array Xi_; // (mM) internal concentration
array Xd_; // (mM) diffusive concentration
array Xo_; // (mM) external concentration
array gX_; // (kS/m²) per-species conductivity

array init_Xi_; // (mM) area-weighted initial internal concentration
array init_Xo_; // (mM) area-weighted initial external concentration
Expand Down Expand Up @@ -246,7 +245,7 @@ ARB_ARBOR_API std::ostream& operator<<(std::ostream& o, shared_state& s);

} // namespace gpu

ARB_SERDES_ENABLE_EXT(gpu::ion_state, Xd_, gX_);
ARB_SERDES_ENABLE_EXT(gpu::ion_state, Xd_);
ARB_SERDES_ENABLE_EXT(gpu::mech_storage,
data_,
// NOTE(serdes) ion_states_, this is just a bunch of pointers
Expand Down
31 changes: 5 additions & 26 deletions arbor/backends/multicore/diffusion_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@ struct diffusion_solver {

array d; // [μS]
array u; // [μS]
array cv_area; // [μm^2]
array cv_volume; // [μm^3]
array invariant_d; // [μS] invariant part of matrix diagonal

Expand All @@ -36,15 +35,12 @@ struct diffusion_solver {
diffusion_solver(const std::vector<index_type>& p,
const std::vector<index_type>& cell_cv_divs,
const std::vector<value_type>& diff,
const std::vector<value_type>& area,
const std::vector<value_type>& volume):
parent_index(p.begin(), p.end()),
cell_cv_divs(cell_cv_divs.begin(), cell_cv_divs.end()),
d(size(), 0), u(size(), 0),
cv_area(area.begin(), area.end()),
cv_volume(volume.begin(), volume.end()),
invariant_d(size(), 0)
{
invariant_d(size(), 0) {
// Sanity check
arb_assert(diff.size() == size());
arb_assert(cell_cv_divs.back() == (index_type)size());
Expand All @@ -65,39 +61,22 @@ struct diffusion_solver {

// Assemble and solve the matrix
// Assemble the matrix
// concentration [mM] (per control volume)
// dt [ms] (scalar)
// internal conc [mM] (per control volume)
// voltage [mV] (per control volume)
// current density [A.m^-2] (per control volume and species)
// diffusivity [m^2/s] (per control volume)
// charge [e]
// diff. concentration
// * will be overwritten by the solution
template<typename T>
void solve(T& concentration,
value_type dt,
const_view voltage,
const_view current,
const_view conductivity,
arb_value_type q) {
value_type dt) {
auto cell_cv_part = util::partition_view(cell_cv_divs);
index_type ncells = cell_cv_part.size();

value_type _1_dt = 1e-3/dt; // 1/µs
// loop over submatrices
for (auto m: util::make_span(0, ncells)) {
for (auto i: util::make_span(cell_cv_part[m])) {
auto u = voltage[i]; // mV
auto g = conductivity[i]; // µS
auto J = current[i]; // A/m^2
auto A = 1e-3*cv_area[i]; // 1e-9·m²
auto X = concentration[i]; // mM
auto V = cv_volume[i]; // m^3
// conversion from current density to concentration change
// using Faraday's constant
auto F = A/(q*96.485332);
d[i] = _1_dt*V + F*g + invariant_d[i];
concentration[i] = _1_dt*V*X + F*(u*g - J);
d[i] = _1_dt*V + invariant_d[i];
concentration[i] = _1_dt*V*X;
}
}
solve(concentration);
Expand Down
3 changes: 0 additions & 3 deletions arbor/backends/multicore/shared_state.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,6 @@ ion_state::ion_state(const fvm_ion_config& ion_data,
Xi_(ion_data.init_iconc.begin(), ion_data.init_iconc.end(), pad(alignment)),
Xd_(ion_data.reset_iconc.begin(), ion_data.reset_iconc.end(), pad(alignment)),
Xo_(ion_data.init_econc.begin(), ion_data.init_econc.end(), pad(alignment)),
gX_(ion_data.cv.size(), NAN, pad(alignment)),
init_Xi_(ion_data.init_iconc.begin(), ion_data.init_iconc.end(), pad(alignment)),
init_Xo_(ion_data.init_econc.begin(), ion_data.init_econc.end(), pad(alignment)),
reset_Xi_(ion_data.reset_iconc.begin(), ion_data.reset_iconc.end(), pad(alignment)),
Expand All @@ -88,7 +87,6 @@ void ion_state::init_concentration() {
}

void ion_state::zero_current() {
util::zero(gX_);
util::zero(iX_);
}

Expand Down Expand Up @@ -432,7 +430,6 @@ void shared_state::instantiate(arb::mechanism& m,
ion_state.external_concentration = oion->Xo_.data();
ion_state.diffusive_concentration = oion->Xd_.data();
ion_state.ionic_charge = oion->charge.data();
ion_state.conductivity = oion->gX_.data();
}

// Initialize state and parameter vectors with default values.
Expand Down
3 changes: 1 addition & 2 deletions arbor/backends/multicore/shared_state.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,6 @@ struct ARB_ARBOR_API ion_state {
array Xi_; // (mM) internal concentration
array Xd_; // (mM) diffusive internal concentration
array Xo_; // (mM) external concentration
array gX_; // (kS/m²) per-species conductivity

array init_Xi_; // (mM) area-weighted initial internal concentration
array init_Xo_; // (mM) area-weighted initial external concentration
Expand Down Expand Up @@ -249,7 +248,7 @@ ARB_ARBOR_API std::ostream& operator<<(std::ostream& o, const shared_state& s);
} // namespace multicore

// Xd and gX are the only things that persist
ARB_SERDES_ENABLE_EXT(multicore::ion_state, Xd_, gX_);
ARB_SERDES_ENABLE_EXT(multicore::ion_state, Xd_);
ARB_SERDES_ENABLE_EXT(multicore::mech_storage,
data_,
// NOTE(serdes) ion_states_, this is just a bunch of pointers
Expand Down
7 changes: 1 addition & 6 deletions arbor/backends/shared_state_base.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,6 @@ struct shared_state_base {
if (data.is_diffusive) solver = std::make_unique<diff_solver>(disc.geometry.cv_parent,
disc.geometry.cell_cv_divs,
data.face_diffusivity,
disc.cv_area,
disc.cv_volume);
d->add_ion(ion, data, std::move(solver));
}
Expand Down Expand Up @@ -152,11 +151,7 @@ struct shared_state_base {
for (auto& [ion, data]: d->ion_data) {
if (data.solver) {
data.solver->solve(data.Xd_,
d->dt,
d->voltage,
data.iX_,
data.gX_,
data.charge[0]);
d->dt);
}
}
}
Expand Down
3 changes: 1 addition & 2 deletions arbor/include/arbor/mechanism_abi.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ extern "C" {

// Version
#define ARB_MECH_ABI_VERSION_MAJOR 0
#define ARB_MECH_ABI_VERSION_MINOR 6
#define ARB_MECH_ABI_VERSION_MINOR 7
#define ARB_MECH_ABI_VERSION_PATCH 0
#define ARB_MECH_ABI_VERSION ((ARB_MECH_ABI_VERSION_MAJOR * 10000L * 10000L) + (ARB_MECH_ABI_VERSION_MAJOR * 10000L) + ARB_MECH_ABI_VERSION_PATCH)

Expand Down Expand Up @@ -54,7 +54,6 @@ inline const char* arb_mechanism_kind_str(const arb_mechanism_kind& mech) {
// Ion state variables; view into shared_state
typedef struct arb_ion_state {
arb_value_type* current_density;
arb_value_type* conductivity;
arb_value_type* reversal_potential;
arb_value_type* internal_concentration;
arb_value_type* external_concentration;
Expand Down
1 change: 0 additions & 1 deletion modcc/identifier.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,6 @@ enum class sourceKind {
dt,
time,
ion_current,
ion_conductivity,
ion_current_density,
ion_diffusive,
ion_revpot,
Expand Down
13 changes: 0 additions & 13 deletions modcc/module.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,6 @@ class NrnCurrentRewriter: public BlockRewriterBase {
if (src==sourceKind::current_density ||
src==sourceKind::current ||
src==sourceKind::ion_current_density ||
src==sourceKind::ion_conductivity ||
src==sourceKind::ion_current)
{
return src;
Expand All @@ -59,8 +58,6 @@ class NrnCurrentRewriter: public BlockRewriterBase {
using BlockRewriterBase::visit;

virtual void finalize() override {
// Take current name 'iX' and strip off leading 'i' to get ion name.
auto i2ion = [this](const auto& name) { return id("conductivity_" + name.substr(1) + "_"); };
if (has_current_update_) {
expression_ptr current_sum, conductivity_sum;
for (auto& curr: current_vars_) {
Expand All @@ -79,11 +76,6 @@ class NrnCurrentRewriter: public BlockRewriterBase {
conductivity_sum = make_expression<AddBinaryExpression>(
Location{}, std::move(conductivity_sum), cond->clone());
}
if (name != non_specific_current) {
statements_.push_back(make_expression<AssignmentExpression>(loc_,
i2ion(name),
cond->clone()));
}
}
if (current_sum) {
statements_.push_back(make_expression<AssignmentExpression>(loc_,
Expand Down Expand Up @@ -784,11 +776,6 @@ void Module::add_variables_to_symbols() {
for(auto const& ion: neuron_block_.ions) {
for(auto const& var: ion.write) {
update_ion_symbols(var, accessKind::write, ion.name);
auto name = "conductivity_" + ion.name + "_";
if (cond.find(name) == cond.end()) {
create_indexed_variable(name, sourceKind::ion_conductivity, accessKind::write, ion.name, var.location);
cond.insert(name);
}
}

for(auto const& var: ion.read) {
Expand Down
5 changes: 0 additions & 5 deletions modcc/printer/printerutil.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,11 +195,6 @@ ARB_LIBMODCC_API indexed_variable_info decode_indexed_variable(IndexedVariable*
v.scale = 0.1;
v.readonly = false;
break;
case sourceKind::ion_conductivity:
v.data_var = ion_pfx+".conductivity";
v.scale = 0.1;
v.readonly = false;
break;
case sourceKind::ion_current:
// unit scale; sourceKind for point processes updating an ionic current variable.
v.data_var = ion_pfx+".current_density";
Expand Down

0 comments on commit b342ab1

Please sign in to comment.