diff --git a/arbor/backends/gpu/matrix_state_fine.hpp b/arbor/backends/gpu/matrix_state_fine.hpp index 714e1f399..f052db0d9 100644 --- a/arbor/backends/gpu/matrix_state_fine.hpp +++ b/arbor/backends/gpu/matrix_state_fine.hpp @@ -3,7 +3,6 @@ #include #include -#include #include @@ -37,15 +36,11 @@ struct matrix_state_fine { array rhs; // [nA] // Required for matrix assembly - array cv_area; // [μm^2] array cv_capacitance; // [pF] // Invariant part of the matrix diagonal array invariant_d; // [μS] - // Solution in unpacked format - array solution_; - // Maximum number of branches in each level per block unsigned max_branches_per_level; @@ -82,16 +77,13 @@ struct matrix_state_fine { // `solver_format[perm[i]] = external_format[i]` iarray perm; - matrix_state_fine() = default; // constructor for fine-grained matrix. matrix_state_fine(const std::vector& p, - const std::vector& cell_cv_divs, - const std::vector& cap, - const std::vector& face_conductance, - const std::vector& area) - { + const std::vector& cell_cv_divs, + const std::vector& cap, + const std::vector& face_conductance) { using util::make_span; constexpr unsigned npos = unsigned(-1); @@ -360,7 +352,6 @@ struct matrix_state_fine { // cv_capacitance : flat // 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); @@ -386,9 +377,6 @@ struct matrix_state_fine { // transform u_shuffled values into packed u vector. 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); - // the cv_capacitance can be copied directly because it is // to be stored in flat format cv_capacitance = memory::make_const_view(cap); @@ -408,19 +396,18 @@ struct matrix_state_fine { // voltage [mV] // current density [A/m²] // conductivity [kS/m²] - void assemble(const T dt, 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(), - dt, - perm.data(), - size()); + void assemble(const T dt, const_view voltage, const_view current, const_view conductivity, const_view area_um2) { + assemble_matrix_fine(d.data(), + rhs.data(), + invariant_d.data(), + voltage.data(), + current.data(), + conductivity.data(), + cv_capacitance.data(), + area_um2.data(), + dt, + perm.data(), + size()); } void solve(array& to) { @@ -441,8 +428,8 @@ struct matrix_state_fine { void solve(array& voltage, - const T dt, const_view current, const_view conductivity) { - assemble(dt, voltage, current, conductivity); + const T dt, const_view current, const_view conductivity, const_view area_um2) { + assemble(dt, voltage, current, conductivity, area_um2); solve(voltage); } diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp index 928c938ba..c6213b54a 100644 --- a/arbor/backends/gpu/shared_state.cpp +++ b/arbor/backends/gpu/shared_state.cpp @@ -49,23 +49,42 @@ ion_state::ion_state(const fvm_ion_config& ion_data, write_eX_(ion_data.revpot_written), write_Xo_(ion_data.econc_written), write_Xi_(ion_data.iconc_written), + write_Xd_(ion_data.is_diffusive), + read_Xo_(ion_data.econc_written || ion_data.econc_read), // ensure that if we have W access, also R access is flagged + read_Xi_(ion_data.iconc_written || ion_data.iconc_read), node_index_(make_const_view(ion_data.cv)), iX_(ion_data.cv.size(), NAN), eX_(ion_data.init_revpot.begin(), ion_data.init_revpot.end()), - 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)), - reset_Xo_(make_const_view(ion_data.reset_econc)), - init_eX_(make_const_view(ion_data.init_revpot)), charge(1u, static_cast(ion_data.charge)), solver(std::move(ptr)) { - arb_assert(node_index_.size()==init_Xi_.size()); - arb_assert(node_index_.size()==init_Xo_.size()); - arb_assert(node_index_.size()==init_eX_.size()); + // We don't need to allocate these if we never use them... + if (read_Xi_) { + Xi_ = make_const_view(ion_data.init_iconc); + } + if (read_Xo_) { + Xo_ = make_const_view(ion_data.init_econc); + } + if (write_Xi_ || write_Xd_) { + // ... but this is used by Xd and Xi! + reset_Xi_ = make_const_view(ion_data.reset_iconc); + } + if (write_Xi_) { + init_Xi_ = make_const_view(ion_data.init_iconc); + arb_assert(node_index_.size()==init_Xi_.size()); + } + if (write_Xo_) { + init_Xo_ = make_const_view(ion_data.init_econc); + reset_Xo_ = make_const_view(ion_data.reset_econc); + arb_assert(node_index_.size()==init_Xo_.size()); + } + if (write_eX_) { + init_eX_ = make_const_view(ion_data.init_revpot); + arb_assert(node_index_.size()==init_eX_.size()); + } + if (write_Xd_) { + Xd_ = array(ion_data.cv.size(), NAN); + } } void ion_state::init_concentration() { @@ -81,10 +100,13 @@ void ion_state::zero_current() { void ion_state::reset() { zero_current(); - memory::copy(reset_Xi_, Xd_); if (write_Xi_) memory::copy(reset_Xi_, Xi_); if (write_Xo_) memory::copy(reset_Xo_, Xo_); if (write_eX_) memory::copy(init_eX_, eX_); + // This goes _last_ or at least after Xi since we might have removed reset_Xi + // when Xi is constant. Thus conditionally resetting Xi first and then copying + // Xi -> Xd is save in all cases. + if (write_Xd_) memory::copy(reset_Xi_, Xd_); } // istim_state methods: diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp index 030cbcdab..792fc7476 100644 --- a/arbor/backends/gpu/shared_state.hpp +++ b/arbor/backends/gpu/shared_state.hpp @@ -40,9 +40,14 @@ struct ARB_ARBOR_API ion_state { using solver_type = arb::gpu::diffusion_state; using solver_ptr = std::unique_ptr; - bool write_eX_; // is eX written? - bool write_Xo_; // is Xo written? - bool write_Xi_; // is Xi written? + bool write_eX_:1; // is eX written? + bool write_Xo_:1; // is Xo written? + bool write_Xi_:1; // is Xi written? + bool write_Xd_:1; // is Xd written? + bool read_eX_:1; // is eX read? + bool read_Xo_:1; // is Xo read? + bool read_Xi_:1; // is Xi read? + bool read_Xd_:1; // is Xd read? iarray node_index_; // Instance to CV map. array iX_; // (A/m²) current density diff --git a/arbor/backends/multicore/cable_solver.hpp b/arbor/backends/multicore/cable_solver.hpp index 3622ee688..bc748486a 100644 --- a/arbor/backends/multicore/cable_solver.hpp +++ b/arbor/backends/multicore/cable_solver.hpp @@ -23,7 +23,6 @@ struct cable_solver { array d; // [μS] array u; // [μS] array cv_capacitance; // [pF] - array cv_area; // [μm^2] array invariant_d; // [μS] invariant part of matrix diagonal cable_solver() = default; @@ -36,13 +35,11 @@ struct cable_solver { cable_solver(const std::vector& p, const std::vector& cell_cv_divs, const std::vector& cap, - const std::vector& cond, - const std::vector& area): + const std::vector& cond): parent_index(p.begin(), p.end()), cell_cv_divs(cell_cv_divs.begin(), cell_cv_divs.end()), d(size(), 0), u(size(), 0), cv_capacitance(cap.begin(), cap.end()), - cv_area(area.begin(), area.end()), invariant_d(size(), 0) { // Sanity check @@ -67,7 +64,7 @@ struct cable_solver { // * expects the voltage from its first argument // * will likewise overwrite the first argument with the solction template - void solve(T& rhs, const value_type dt, const_view current, const_view conductivity) { + void solve(T& rhs, const value_type dt, const_view current, const_view conductivity, const_view cv_area) { value_type * const ARB_NO_ALIAS d_ = d.data(); value_type * const ARB_NO_ALIAS r_ = rhs.data(); diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp index 17a57789f..e625c5b97 100644 --- a/arbor/backends/multicore/shared_state.cpp +++ b/arbor/backends/multicore/shared_state.cpp @@ -1,6 +1,5 @@ #include #include -#include #include #include #include @@ -58,24 +57,43 @@ ion_state::ion_state(const fvm_ion_config& ion_data, write_eX_(ion_data.revpot_written), write_Xo_(ion_data.econc_written), write_Xi_(ion_data.iconc_written), + write_Xd_(ion_data.is_diffusive), + read_Xo_(ion_data.econc_written || ion_data.econc_read), // ensure that if we have W access, also R access is flagged + read_Xi_(ion_data.iconc_written || ion_data.iconc_read), + read_Xd_(ion_data.is_diffusive), node_index_(ion_data.cv.begin(), ion_data.cv.end(), pad(alignment)), iX_(ion_data.cv.size(), NAN, pad(alignment)), eX_(ion_data.init_revpot.begin(), ion_data.init_revpot.end(), pad(alignment)), - 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)), - reset_Xo_(ion_data.reset_econc.begin(), ion_data.reset_econc.end(), pad(alignment)), - init_eX_(ion_data.init_revpot.begin(), ion_data.init_revpot.end(), pad(alignment)), charge(1u, ion_data.charge, pad(alignment)), solver(std::move(ptr)) { - arb_assert(node_index_.size()==init_Xi_.size()); - arb_assert(node_index_.size()==init_Xo_.size()); - arb_assert(node_index_.size()==eX_.size()); - arb_assert(node_index_.size()==init_eX_.size()); + // We don't need to allocate these if we never use them... + if (read_Xi_) { + Xi_ = {ion_data.init_iconc.begin(), ion_data.init_iconc.end(), pad(alignment)}; + } + if (read_Xo_) { + Xo_ = {ion_data.init_econc.begin(), ion_data.init_econc.end(), pad(alignment)}; + } + if (write_Xi_ || write_Xd_) { + // ... but this is used by Xd and Xi! + reset_Xi_ = {ion_data.reset_iconc.begin(), ion_data.reset_iconc.end(), pad(alignment)}; + } + if (write_Xi_) { + init_Xi_ = {ion_data.init_iconc.begin(), ion_data.init_iconc.end(), pad(alignment)}; + arb_assert(node_index_.size()==init_Xi_.size()); + } + if (write_Xo_) { + init_Xo_ = {ion_data.init_econc.begin(), ion_data.init_econc.end(), pad(alignment)}; + reset_Xo_ = {ion_data.reset_econc.begin(), ion_data.reset_econc.end(), pad(alignment)}; + arb_assert(node_index_.size()==init_Xo_.size()); + } + if (write_eX_) { + init_eX_ = {ion_data.init_revpot.begin(), ion_data.init_revpot.end(), pad(alignment)}; + arb_assert(node_index_.size()==init_eX_.size()); + } + if (read_Xd_) { + Xd_ = {ion_data.reset_iconc.begin(), ion_data.reset_iconc.end(), pad(alignment)}; + } } void ion_state::init_concentration() { @@ -91,10 +109,13 @@ void ion_state::zero_current() { void ion_state::reset() { zero_current(); - std::copy(reset_Xi_.begin(), reset_Xi_.end(), Xd_.begin()); if (write_Xi_) std::copy(reset_Xi_.begin(), reset_Xi_.end(), Xi_.begin()); if (write_Xo_) std::copy(reset_Xo_.begin(), reset_Xo_.end(), Xo_.begin()); if (write_eX_) std::copy(init_eX_.begin(), init_eX_.end(), eX_.begin()); + // This goes _last_ or at least after Xi since we might have removed reset_Xi + // when Xi is constant. Thus conditionally resetting Xi first and then copying + // Xi -> Xd is safe in all cases. + if (write_Xd_) std::copy(Xi_.begin(), Xi_.end(), Xd_.begin()); } // istim_state methods: diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp index 99bbe37b5..81a10bcc8 100644 --- a/arbor/backends/multicore/shared_state.hpp +++ b/arbor/backends/multicore/shared_state.hpp @@ -50,9 +50,14 @@ struct ARB_ARBOR_API ion_state { unsigned alignment = 1; // Alignment and padding multiple. - bool write_eX_; // is eX written? - bool write_Xo_; // is Xo written? - bool write_Xi_; // is Xi written? + bool write_eX_:1; // is eX written? + bool write_Xo_:1; // is Xo written? + bool write_Xi_:1; // is Xi written? + bool write_Xd_:1; // is Xd written? + bool read_eX_:1; // is eX read? + bool read_Xo_:1; // is Xo read? + bool read_Xi_:1; // is Xi read? + bool read_Xd_:1; // is Xd read? iarray node_index_; // Instance to CV map. array iX_; // (A/m²) current density diff --git a/arbor/backends/shared_state_base.hpp b/arbor/backends/shared_state_base.hpp index c20247d36..92b3a5cb5 100644 --- a/arbor/backends/shared_state_base.hpp +++ b/arbor/backends/shared_state_base.hpp @@ -7,8 +7,9 @@ #include "backends/common_types.hpp" #include "fvm_layout.hpp" -#include "event_lane.hpp" +#include "util/rangeutil.hpp" #include "timestep_range.hpp" +#include "event_lane.hpp" namespace arb { @@ -48,7 +49,7 @@ struct shared_state_base { void configure_solver(const fvm_cv_discretization& disc) { auto d = static_cast(this); - d->solver = {disc.geometry.cv_parent, disc.geometry.cell_cv_divs, disc.cv_capacitance, disc.face_conductance, disc.cv_area}; + d->solver = {disc.geometry.cv_parent, disc.geometry.cell_cv_divs, disc.cv_capacitance, disc.face_conductance}; } void add_ion(const std::string& ion_name, @@ -134,7 +135,7 @@ struct shared_state_base { void integrate_cable_state() { auto d = static_cast(this); - d->solver.solve(d->voltage, d->dt, d->current_density, d->conductivity); + d->solver.solve(d->voltage, d->dt, d->current_density, d->conductivity, d->area_um2); for (auto& [ion, data]: d->ion_data) { if (data.solver) { data.solver->solve(data.Xd_, diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp index 04b1ac9ce..c0b049cd5 100644 --- a/arbor/fvm_layout.cpp +++ b/arbor/fvm_layout.cpp @@ -721,9 +721,11 @@ fvm_mechanism_data& append(fvm_mechanism_data& left, const fvm_mechanism_data& r append(L.reset_econc, R.reset_econc); append(L.init_revpot, R.init_revpot); append(L.face_diffusivity, R.face_diffusivity); - L.is_diffusive |= R.is_diffusive; + L.is_diffusive |= R.is_diffusive; L.econc_written |= R.econc_written; L.iconc_written |= R.iconc_written; + L.econc_read |= R.econc_read; + L.iconc_read |= R.iconc_read; L.revpot_written |= R.revpot_written; } @@ -826,14 +828,26 @@ struct fvm_ion_build_data { mcable_map init_econc_mask; bool write_xi = false; bool write_xo = false; + bool read_xi = false; + bool read_xo = false; std::vector support; - void add_to_support(const std::vector& cvs) { + auto& add_to_support(const std::vector& cvs) { arb_assert(util::is_sorted(cvs)); support = unique_union(support, cvs); + return *this; + } + + auto& add_ion_dep(const ion_dependency& dep) { + write_xi |= dep.write_concentration_int; + write_xo |= dep.write_concentration_ext; + read_xi |= dep.read_concentration_int; + read_xo |= dep.read_concentration_ext; + return *this; } }; + using fvm_mechanism_config_map = std::map; using fvm_ion_config_map = std::unordered_map; using fvm_ion_map = std::unordered_map; @@ -913,7 +927,7 @@ make_gj_mechanism_config(const std::unordered_map make_revpot_mechanism_config(const std::unordered_map& method, - const std::unordered_map& ions, + std::unordered_map& ions, const cell_build_data& data, fvm_mechanism_config_map&); @@ -1228,10 +1242,9 @@ make_density_mechanism_config(const region_assignment& assignments, apply_parameters_on_cv(config, data, param_maps, support); for (const auto& [ion, dep]: info.ions) { - auto& build_data = ion_build_data[ion]; - build_data.write_xi |= dep.write_concentration_int; - build_data.write_xo |= dep.write_concentration_ext; - build_data.add_to_support(config.cv); + auto& build_data = ion_build_data[ion] + .add_ion_dep(dep) + .add_to_support(config.cv); auto ok = true; if (dep.write_concentration_int) { @@ -1343,6 +1356,8 @@ make_ion_config(fvm_ion_map build_data, config.econc_written = build_data.write_xo; config.iconc_written = build_data.write_xi; + config.econc_read = build_data.read_xo; + config.iconc_read = build_data.read_xi; if (!config.cv.empty()) result[ion] = std::move(config); } } @@ -1528,10 +1543,9 @@ make_point_mechanism_config(const std::unordered_map make_revpot_mechanism_config(const std::unordered_map& method, - const std::unordered_map& ions, + std::unordered_map& ions, const cell_build_data& data, fvm_mechanism_config_map& result) { std::unordered_map revpot_tbl; @@ -1667,10 +1680,13 @@ make_revpot_mechanism_config(const std::unordered_map(n_cv, val)); } - if (!config.cv.empty()) result[name] = std::move(config); } + ion_conf.econc_read |= dep.read_concentration_ext; + ion_conf.iconc_read |= dep.read_concentration_int; } } diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp index 6c2684ca7..73c5c03b9 100644 --- a/arbor/fvm_layout.hpp +++ b/arbor/fvm_layout.hpp @@ -242,6 +242,8 @@ struct fvm_ion_config { bool revpot_written = false; bool iconc_written = false; bool econc_written = false; + bool iconc_read = false; + bool econc_read = false; // Ordered CV indices where ion must be present. std::vector cv; diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp index 730471d10..27199db20 100644 --- a/arbor/fvm_lowered_cell_impl.hpp +++ b/arbor/fvm_lowered_cell_impl.hpp @@ -288,14 +288,12 @@ void fvm_lowered_cell_impl::update_ion_state() { template void fvm_lowered_cell_impl::assert_voltage_bounded(arb_value_type bound) { - auto v_minmax = state_->voltage_bounds(); - if (v_minmax.first>=-bound && v_minmax.second<=bound) { - return; - } + const auto& [vmin, vmax] = state_->voltage_bounds(); + if (vmin >= -bound && vmax <= bound) return; throw range_check_failure( util::pprintf("voltage solution out of bounds for at t = {}", state_->time), - v_minmax.first<-bound? v_minmax.first: v_minmax.second); + vmin < -bound ? vmin : vmax); } inline @@ -1002,31 +1000,35 @@ void resolve_probe(const cable_probe_ion_current_cell& p, probe_resolution_data< template void resolve_probe(const cable_probe_ion_int_concentration& p, probe_resolution_data& R) { + const auto& ion = p.ion; + const auto& xi = R.state->ion_data.at(ion).Xi_; + if (xi.empty()) return; for (mlocation loc: thingify(p.locations, R.cell.provider())) { - auto opt_i = R.ion_location_index(p.ion, loc); + auto opt_i = R.ion_location_index(ion, loc); if (!opt_i) continue; - - R.result.push_back(fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xi_.data()+*opt_i}, loc}); + R.result.push_back(fvm_probe_scalar{{xi.data() + *opt_i}, loc}); } } template void resolve_probe(const cable_probe_ion_ext_concentration& p, probe_resolution_data& R) { + const auto& ion = p.ion; + const auto& xo = R.state->ion_data.at(ion).Xo_; for (mlocation loc: thingify(p.locations, R.cell.provider())) { - auto opt_i = R.ion_location_index(p.ion, loc); + auto opt_i = R.ion_location_index(ion, loc); if (!opt_i) continue; - - R.result.push_back(fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xo_.data()+*opt_i}, loc}); + R.result.push_back(fvm_probe_scalar{{xo.data() + *opt_i}, loc}); } } template void resolve_probe(const cable_probe_ion_diff_concentration& p, probe_resolution_data& R) { + const auto& ion = p.ion; + const auto& xd = R.state->ion_data.at(ion).Xd_; for (mlocation loc: thingify(p.locations, R.cell.provider())) { - auto opt_i = R.ion_location_index(p.ion, loc); + auto opt_i = R.ion_location_index(ion, loc); if (!opt_i) continue; - - R.result.push_back(fvm_probe_scalar{{R.state->ion_data.at(p.ion).Xd_.data()+*opt_i}, loc}); + R.result.push_back(fvm_probe_scalar{{xd.data() + *opt_i}, loc}); } } @@ -1035,7 +1037,6 @@ template void resolve_ion_conc_common(const std::vector& ion_cvs, const arb_value_type* src, probe_resolution_data& R) { fvm_probe_multi r; mcable_list cables; - for (auto i: util::count_along(ion_cvs)) { for (auto cable: R.D.geometry.cables(ion_cvs[i])) { if (cable.prox_pos!=cable.dist_pos) { @@ -1052,18 +1053,21 @@ void resolve_ion_conc_common(const std::vector& ion_cvs, const a template void resolve_probe(const cable_probe_ion_int_concentration_cell& p, probe_resolution_data& R) { if (!R.state->ion_data.count(p.ion)) return; + if (R.state->ion_data.at(p.ion).Xi_.empty()) return; resolve_ion_conc_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xi_.data(), R); } template void resolve_probe(const cable_probe_ion_ext_concentration_cell& p, probe_resolution_data& R) { if (!R.state->ion_data.count(p.ion)) return; + if (R.state->ion_data.at(p.ion).Xo_.empty()) return; resolve_ion_conc_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xo_.data(), R); } template void resolve_probe(const cable_probe_ion_diff_concentration_cell& p, probe_resolution_data& R) { if (!R.state->ion_data.count(p.ion)) return; + if (R.state->ion_data.at(p.ion).Xd_.empty()) return; resolve_ion_conc_common(R.M.ions.at(p.ion).cv, R.state->ion_data.at(p.ion).Xd_.data(), R); } diff --git a/arbor/include/arbor/mechanism_abi.h b/arbor/include/arbor/mechanism_abi.h index bf0a9b30b..b5bca6f24 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) @@ -193,6 +193,8 @@ typedef struct arb_ion_info { const char* name; bool write_int_concentration; bool write_ext_concentration; + bool read_int_concentration; + bool read_ext_concentration; bool use_diff_concentration; bool write_rev_potential; bool read_rev_potential; diff --git a/arbor/include/arbor/mechinfo.hpp b/arbor/include/arbor/mechinfo.hpp index 6d717a65f..a4527ed99 100644 --- a/arbor/include/arbor/mechinfo.hpp +++ b/arbor/include/arbor/mechinfo.hpp @@ -7,8 +7,6 @@ #include #include #include -#include -#include #include #include @@ -35,6 +33,8 @@ struct mechanism_field_spec { struct ion_dependency { bool write_concentration_int = false; bool write_concentration_ext = false; + bool read_concentration_int = false; + bool read_concentration_ext = false; bool access_concentration_diff = false; diff --git a/arbor/include/arbor/util/expected.hpp b/arbor/include/arbor/util/expected.hpp index 9d860f653..1af88eb42 100644 --- a/arbor/include/arbor/util/expected.hpp +++ b/arbor/include/arbor/util/expected.hpp @@ -486,7 +486,7 @@ struct expected { // Swap ops. void swap(expected& other) { - data_.swap(other.data); + data_.swap(other.data_); } // Accessors. diff --git a/arbor/mechinfo.cpp b/arbor/mechinfo.cpp index 1d5a15205..b96cca7bd 100644 --- a/arbor/mechinfo.cpp +++ b/arbor/mechinfo.cpp @@ -22,14 +22,16 @@ mechanism_info::mechanism_info(const arb_mechanism_type& m) { } for (auto idx: util::make_span(m.n_ions)) { const auto& v = m.ions[idx]; - ions[v.name] = { v.write_int_concentration, - v.write_ext_concentration, - v.use_diff_concentration, - v.read_rev_potential, - v.write_rev_potential, - v.read_valence, - v.verify_valence, - v.expected_valence }; + ions[v.name] = {v.write_int_concentration, + v.write_ext_concentration, + v.read_int_concentration, + v.read_ext_concentration, + v.use_diff_concentration, + v.read_rev_potential, + v.write_rev_potential, + v.read_valence, + v.verify_valence, + v.expected_valence }; } for (auto idx: util::make_span(m.n_random_variables)) { const auto& rv = m.random_variables[idx]; diff --git a/arbor/morph/embed_pwlin.cpp b/arbor/morph/embed_pwlin.cpp index 945dffdf6..02a4da607 100644 --- a/arbor/morph/embed_pwlin.cpp +++ b/arbor/morph/embed_pwlin.cpp @@ -103,7 +103,7 @@ struct embed_pwlin_data { template double interpolate(double pos, const pw_ratpoly& f) { - auto [extent, poly] = f(pos); + const auto& [extent, poly] = f(pos); auto [left, right] = extent; return left==right? poly[0]: poly((pos-left)/(right-left)); diff --git a/modcc/blocks.hpp b/modcc/blocks.hpp index 1e12b5428..448e2f2d4 100644 --- a/modcc/blocks.hpp +++ b/modcc/blocks.hpp @@ -4,10 +4,8 @@ #include #include #include -#include #include "identifier.hpp" -#include "location.hpp" #include "token.hpp" #include @@ -48,6 +46,15 @@ struct IonDep { bool writes_concentration_ext() const { return writes_variable(name + "o"); }; + bool reads_current() const { + return reads_variable("i" + name); + }; + bool reads_concentration_int() const { + return reads_variable(name + "i"); + }; + bool reads_concentration_ext() const { + return reads_variable(name + "o"); + }; bool writes_rev_potential() const { return writes_variable("e" + name); }; diff --git a/modcc/printer/infoprinter.cpp b/modcc/printer/infoprinter.cpp index f59f52f54..f201fd5fb 100644 --- a/modcc/printer/infoprinter.cpp +++ b/modcc/printer/infoprinter.cpp @@ -49,12 +49,13 @@ ARB_LIBMODCC_API std::string build_info_header(const Module& m, const printer_op id.unit_string(), val, lo, hi); }; auto fmt_ion = [&](const auto& ion) { - return fmt::format(FMT_COMPILE("{{ \"{}\", {}, {}, {}, {}, {}, {}, {}, {} }}"), - ion.name, - ion.writes_concentration_int(), ion.writes_concentration_ext(), - ion.uses_concentration_diff(), - ion.writes_rev_potential(), ion.uses_rev_potential(), - ion.uses_valence(), ion.verifies_valence(), ion.expected_valence); + return fmt::format(FMT_COMPILE("{{ \"{}\", {}, {}, {}, {}, {}, {}, {}, {}, {}, {} }}"), + ion.name, + ion.writes_concentration_int(), ion.writes_concentration_ext(), + ion.reads_concentration_int(), ion.reads_concentration_ext(), + ion.uses_concentration_diff(), + ion.writes_rev_potential(), ion.uses_rev_potential(), + ion.uses_valence(), ion.verifies_valence(), ion.expected_valence); }; auto fmt_random_variable = [&](const auto& rv_id) { return fmt::format(FMT_COMPILE("{{ \"{}\", {} }}"), rv_id.first.name(), rv_id.second); diff --git a/python/strprintf.hpp b/python/strprintf.hpp index 638016772..2ccf12c1a 100644 --- a/python/strprintf.hpp +++ b/python/strprintf.hpp @@ -168,7 +168,7 @@ namespace impl { for (auto& x: s.seq_) { if (!first) o << s.sep_; first = false; - o << s.f(x); + o << s.f_(x); } return o; } diff --git a/python/test/unit/test_probes.py b/python/test/unit/test_probes.py index f72bfb4e3..27c69e32d 100644 --- a/python/test/unit/test_probes.py +++ b/python/test/unit/test_probes.py @@ -19,12 +19,15 @@ def __init__(self): st = A.segment_tree() st.append(A.mnpos, (0, 0, 0, 10), (1, 0, 0, 10), 1) - dec = A.decor() - - dec.place("(location 0 0.08)", A.synapse("expsyn"), "syn0") - dec.place("(location 0 0.09)", A.synapse("exp2syn"), "syn1") - dec.place("(location 0 0.1)", A.iclamp(20.0 * U.nA), "iclamp") - dec.paint("(all)", A.density("hh")) + dec = ( + A.decor() + # This ensures we _read_ nai/nao and can probe them, too + .set_ion(ion="na", method="nernst/x=na") + .place("(location 0 0.08)", A.synapse("expsyn"), "syn0") + .place("(location 0 0.09)", A.synapse("exp2syn"), "syn1") + .place("(location 0 0.1)", A.iclamp(20.0 * U.nA), "iclamp") + .paint("(all)", A.density("hh")) + ) self.cell = A.cable_cell(st, dec) @@ -40,7 +43,7 @@ def cell_kind(self, gid): def global_properties(self, kind): return self.props - def probes(self, _): + def probes(self, gid): # Use keyword arguments to check that the wrappers have actually declared keyword arguments correctly. # Place single-location probes at (location 0 0.01*j) where j is the index of the probe address in # the returned list. diff --git a/test/unit/test_math.cpp b/test/unit/test_math.cpp index 755e325c6..22e788b5b 100644 --- a/test/unit/test_math.cpp +++ b/test/unit/test_math.cpp @@ -97,15 +97,13 @@ TEST(math, infinity) { EXPECT_GT(ldinf, 0.0l); // check default value promotes correctly (i.e., acts like INFINITY) - struct { - float f; - double d; - long double ld; - } check = {infinity<>, infinity<>, infinity<>}; - - EXPECT_EQ(std::numeric_limits::infinity(), check.f); - EXPECT_EQ(std::numeric_limits::infinity(), check.d); - EXPECT_EQ(std::numeric_limits::infinity(), check.ld); + float f = infinity<>; + double d = infinity<>; + long double ld = infinity<>; + + EXPECT_EQ(std::numeric_limits::infinity(), f); + EXPECT_EQ(std::numeric_limits::infinity(), d); + EXPECT_EQ(std::numeric_limits::infinity(), ld); } TEST(math, signum) { diff --git a/test/unit/test_matrix.cpp b/test/unit/test_matrix.cpp index 4a7a3b320..c6a42f1c1 100644 --- a/test/unit/test_matrix.cpp +++ b/test/unit/test_matrix.cpp @@ -21,10 +21,9 @@ using value_type = arb_value_type; using vvec = std::vector; -TEST(matrix, construct_from_parent_only) -{ +TEST(matrix, construct_from_parent_only) { std::vector p = {0,0,1}; - solver_type m(p, {0, 3}, vvec(3), vvec(3), vvec(3)); + solver_type m(p, {0, 3}, vvec(3), vvec(3)); EXPECT_EQ(m.num_cells(), 1u); EXPECT_EQ(m.size(), 3u); EXPECT_EQ(p.size(), 3u); @@ -35,14 +34,13 @@ TEST(matrix, construct_from_parent_only) EXPECT_EQ(mp[2], index_type(1)); } -TEST(matrix, solve_host) -{ +TEST(matrix, solve_host) { using util::make_span; using util::fill; // trivial case : 1x1 matrix { - solver_type m({0}, {0,1}, vvec(1), vvec(1), vvec(1)); + solver_type m({0}, {0,1}, vvec(1), vvec(1)); fill(m.d, 2); fill(m.u, -1); array x({1}); @@ -56,7 +54,7 @@ TEST(matrix, solve_host) for(auto n : make_span(2, 1001)) { auto p = std::vector(n); std::iota(p.begin()+1, p.end(), 0); - solver_type m(p, {0, n}, vvec(n), vvec(n), vvec(n)); + solver_type m(p, {0, n}, vvec(n), vvec(n)); EXPECT_EQ(m.size(), (unsigned)n); EXPECT_EQ(m.num_cells(), 1u); @@ -78,8 +76,7 @@ TEST(matrix, solve_host) } } -TEST(matrix, solve_multi_matrix) -{ +TEST(matrix, solve_multi_matrix) { // Use assemble method to construct same zero-diagonal // test case from CV data. @@ -107,7 +104,7 @@ TEST(matrix, solve_multi_matrix) // Initial voltage of zero; currents alone determine rhs. auto v = vvec{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; - vvec area(7, 1.0); + array area(7, 1.0); // (Scaled) membrane conductances contribute to diagonal. array mg = { 1000, 2000, 3000, 4000, 5000, 6000, 7000 }; @@ -121,9 +118,9 @@ TEST(matrix, solve_multi_matrix) // Expected solution: // x = [ 4 5 6 7 8 9 10 ] - solver_type m(p, c, Cm, g, area); + solver_type m(p, c, Cm, g); std::vector expected = {4, 5, 6, 7, 8, 9, 10}; - m.solve(v, dt, i, mg); + m.solve(v, dt, i, mg, area); EXPECT_TRUE(testing::seq_almost_eq(expected, v)); } diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp index 43449d017..52153fac5 100644 --- a/test/unit/test_probe.cpp +++ b/test/unit/test_probe.cpp @@ -475,6 +475,18 @@ void run_expsyn_g_cell_probe_test(context ctx) { } } +template +typename B::array mk_array(size_t n, size_t a) { + return typename B::array(n, 0, util::padded_allocator<>(a)); +} + +#ifdef ARB_GPU_ENABLED +template<> +typename arb::gpu::backend::array mk_array(size_t n, size_t a) { + return arb::gpu::backend::array(n, 0); +} +#endif + template void run_ion_density_probe_test(context ctx) { using fvm_cell = typename backend_access::fvm_cell; @@ -530,15 +542,39 @@ void run_ion_density_probe_test(context ctx) { fvm_cell lcell(*ctx); + using array = typename Backend::array; + auto fvm_info = lcell.initialize({0}, rec); // We skipped FVM layout here, so we need to set these manually auto& state = backend_access::state(lcell); - state.ion_data["ca"].write_Xi_ = true; - state.ion_data["ca"].write_Xo_ = true; - state.ion_data["ca"].init_concentration(); - state.ion_data["na"].write_Xi_ = true; - state.ion_data["na"].write_Xo_ = true; - state.ion_data["na"].init_concentration(); + auto align = state.alignment; + + auto& ca = state.ion_data["ca"]; + auto nca = ca.node_index_.size(); + auto cai = mk_array(nca, align); + ca.write_Xi_ = true; + ca.Xi_ = cai; + ca.init_Xi_ = cai; + ca.reset_Xi_ = cai; + auto cao = mk_array(nca, align); + ca.write_Xo_ = true; + ca.Xo_ = cao; + ca.init_Xo_ = cao; + ca.reset_Xo_ = cao; + + auto& na = state.ion_data["na"]; + auto nna = na.node_index_.size(); + auto nai = mk_array(nna, align); + na.write_Xi_ = true; + na.Xi_ = nai; + na.init_Xi_ = nai; + na.reset_Xi_ = nai; + auto nao = mk_array(nna, align); + na.write_Xo_ = true; + na.Xo_ = nao; + na.init_Xo_ = nao; + na.reset_Xo_ = nao; + // Now, re-init cell lcell.reset(); diff --git a/test/unit/test_sde.cpp b/test/unit/test_sde.cpp index ed9b1a1ef..103046d61 100644 --- a/test/unit/test_sde.cpp +++ b/test/unit/test_sde.cpp @@ -732,8 +732,7 @@ TEST(sde, solver) { // context auto context = make_context({arbenv::default_concurrency(), -1}); - for (unsigned s=0; s