From b342ab1d2dfadd440fe01efa29962cf284aa6376 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Mon, 16 Oct 2023 08:03:42 +0200 Subject: [PATCH] Rip out concentration change by current. --- arbor/backends/gpu/diffusion.cu | 24 ++------------ arbor/backends/gpu/diffusion_state.hpp | 20 ++---------- arbor/backends/gpu/shared_state.cpp | 3 -- arbor/backends/gpu/shared_state.hpp | 3 +- arbor/backends/multicore/diffusion_solver.hpp | 31 +++---------------- arbor/backends/multicore/shared_state.cpp | 3 -- arbor/backends/multicore/shared_state.hpp | 3 +- arbor/backends/shared_state_base.hpp | 7 +---- arbor/include/arbor/mechanism_abi.h | 3 +- modcc/identifier.hpp | 1 - modcc/module.cpp | 13 -------- modcc/printer/printerutil.cpp | 5 --- 12 files changed, 15 insertions(+), 101 deletions(-) diff --git a/arbor/backends/gpu/diffusion.cu b/arbor/backends/gpu/diffusion.cu index b5cc5c8245..660d770ce8 100644 --- a/arbor/backends/gpu/diffusion.cu +++ b/arbor/backends/gpu/diffusion.cu @@ -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, @@ -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; } } @@ -199,11 +186,6 @@ 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, @@ -211,7 +193,7 @@ ARB_ARBOR_API void assemble_diffusion(arb_value_type* d, launch_1d(n, 128, kernels::assemble_diffusion, - d, rhs, invariant_d, concentration, voltage, current, q, conductivity, area, volume, dt, perm, n); + d, rhs, invariant_d, volume, dt, perm, n); } // Example: diff --git a/arbor/backends/gpu/diffusion_state.hpp b/arbor/backends/gpu/diffusion_state.hpp index dc950d95bf..5171fcc545 100644 --- a/arbor/backends/gpu/diffusion_state.hpp +++ b/arbor/backends/gpu/diffusion_state.hpp @@ -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 @@ -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 invariant_d_tmp(matrix_size, 0); @@ -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); @@ -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(), @@ -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); } diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp index a15a3f9545..ec7f2fdd64 100644 --- a/arbor/backends/gpu/shared_state.cpp +++ b/arbor/backends/gpu/shared_state.cpp @@ -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)), @@ -77,7 +76,6 @@ void ion_state::init_concentration() { } void ion_state::zero_current() { - memory::fill(gX_, 0); memory::fill(iX_, 0); } @@ -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. diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp index 8c9beb3a10..a5ab2d338d 100644 --- a/arbor/backends/gpu/shared_state.hpp +++ b/arbor/backends/gpu/shared_state.hpp @@ -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 @@ -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 diff --git a/arbor/backends/multicore/diffusion_solver.hpp b/arbor/backends/multicore/diffusion_solver.hpp index 75913c5c80..0c702adfb1 100644 --- a/arbor/backends/multicore/diffusion_solver.hpp +++ b/arbor/backends/multicore/diffusion_solver.hpp @@ -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 @@ -36,15 +35,12 @@ struct diffusion_solver { diffusion_solver(const std::vector& p, const std::vector& cell_cv_divs, const std::vector& diff, - const std::vector& area, const std::vector& 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()); @@ -65,21 +61,11 @@ 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 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(); @@ -87,17 +73,10 @@ struct diffusion_solver { // 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); diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp index a32ea19b0f..fcc1041698 100644 --- a/arbor/backends/multicore/shared_state.cpp +++ b/arbor/backends/multicore/shared_state.cpp @@ -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)), @@ -88,7 +87,6 @@ void ion_state::init_concentration() { } void ion_state::zero_current() { - util::zero(gX_); util::zero(iX_); } @@ -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. diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp index 606527823d..b30f3fa5b2 100644 --- a/arbor/backends/multicore/shared_state.hpp +++ b/arbor/backends/multicore/shared_state.hpp @@ -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 @@ -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 diff --git a/arbor/backends/shared_state_base.hpp b/arbor/backends/shared_state_base.hpp index 643cd924f4..87065b6dbd 100644 --- a/arbor/backends/shared_state_base.hpp +++ b/arbor/backends/shared_state_base.hpp @@ -78,7 +78,6 @@ struct shared_state_base { if (data.is_diffusive) solver = std::make_unique(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)); } @@ -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); } } } diff --git a/arbor/include/arbor/mechanism_abi.h b/arbor/include/arbor/mechanism_abi.h index bf0a9b30ba..103595fde9 100644 --- a/arbor/include/arbor/mechanism_abi.h +++ b/arbor/include/arbor/mechanism_abi.h @@ -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) @@ -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; diff --git a/modcc/identifier.hpp b/modcc/identifier.hpp index a2ce9fbbc3..00f78ff89a 100644 --- a/modcc/identifier.hpp +++ b/modcc/identifier.hpp @@ -49,7 +49,6 @@ enum class sourceKind { dt, time, ion_current, - ion_conductivity, ion_current_density, ion_diffusive, ion_revpot, diff --git a/modcc/module.cpp b/modcc/module.cpp index 27b7bae246..d32a1c2f03 100644 --- a/modcc/module.cpp +++ b/modcc/module.cpp @@ -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; @@ -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_) { @@ -79,11 +76,6 @@ class NrnCurrentRewriter: public BlockRewriterBase { conductivity_sum = make_expression( Location{}, std::move(conductivity_sum), cond->clone()); } - if (name != non_specific_current) { - statements_.push_back(make_expression(loc_, - i2ion(name), - cond->clone())); - } } if (current_sum) { statements_.push_back(make_expression(loc_, @@ -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) { diff --git a/modcc/printer/printerutil.cpp b/modcc/printer/printerutil.cpp index abda0e0a94..af89777ffa 100644 --- a/modcc/printer/printerutil.cpp +++ b/modcc/printer/printerutil.cpp @@ -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";