CV area exposed to NMODL.

arbor/backends/gpu/shared_state.cpp

@@ -174,6 +174,7 @@ shared_state::shared_state(task_system_handle tp,
                            const std::vector<arb_value_type>& init_membrane_potential,
                            const std::vector<arb_value_type>& temperature_K,
                            const std::vector<arb_value_type>& diam,
+                           const std::vector<arb_value_type>& area,
                            const std::vector<arb_index_type>& src_to_spike_,
                            const fvm_detector_info& detector,
                            unsigned, // align parameter ignored
@@ -188,6 +189,7 @@ shared_state::shared_state(task_system_handle tp,
     init_voltage(make_const_view(init_membrane_potential)),
     temperature_degC(make_const_view(temperature_K)),
     diam_um(make_const_view(diam)),
+    area_um2(make_const_view(area)),
     time_since_spike(n_cell*n_detector),
     src_to_spike(make_const_view(src_to_spike_)),
     cbprng_seed(cbprng_seed_),
@@ -237,6 +239,7 @@ void shared_state::instantiate(mechanism& m,
     m.ppack_.vec_g            = conductivity.data();
     m.ppack_.temperature_degC = temperature_degC.data();
     m.ppack_.diam_um          = diam_um.data();
+    m.ppack_.area_um2         = area_um2.data();
     m.ppack_.time_since_spike = time_since_spike.data();
     m.ppack_.n_detectors      = n_detector;
 
@@ -393,26 +396,25 @@ void shared_state::take_samples() {
 ARB_ARBOR_API std::ostream& operator<<(std::ostream& o, shared_state& s) {
     using io::csv;
 
-    o << " n_cv         " << s.n_cv << "\n";
-    o << " time         " << s.time << "\n";
-    o << " time_to      " << s.time_to << "\n";
-    o << " voltage      " << s.voltage << "\n";
-    o << " init_voltage " << s.init_voltage << "\n";
-    o << " temperature  " << s.temperature_degC << "\n";
-    o << " diameter     " << s.diam_um << "\n";
-    o << " current      " << s.current_density << "\n";
-    o << " conductivity " << s.conductivity << "\n";
-    for (auto& ki: s.ion_data) {
-        auto& kn = ki.first;
-        auto& i = ki.second;
-        o << " " << kn << "/current_density        " << csv(i.iX_) << "\n";
-        o << " " << kn << "/reversal_potential     " << csv(i.eX_) << "\n";
-        o << " " << kn << "/internal_concentration " << csv(i.Xi_) << "\n";
-        o << " " << kn << "/external_concentration " << csv(i.Xo_) << "\n";
-        o << " " << kn << "/intconc_initial        " << csv(i.init_Xi_) << "\n";
-        o << " " << kn << "/extconc_initial        " << csv(i.init_Xo_) << "\n";
-        o << " " << kn << "/revpot_initial         " << csv(i.init_eX_) << "\n";
-        o << " " << kn << "/node_index             " << csv(i.node_index_) << "\n";
+    o << " n_cv         " << s.n_cv << "\n"
+      << " time         " << s.time << "\n"
+      << " time_to      " << s.time_to << "\n"
+      << " voltage      " << s.voltage << "\n"
+      << " init_voltage " << s.init_voltage << "\n"
+      << " temperature  " << s.temperature_degC << "\n"
+      << " diameter     " << s.diam_um << "\n"
+      << " area         " << s.area_um2 << "\n"
+      << " current      " << s.current_density << "\n"
+      << " conductivity " << s.conductivity << "\n";
+    for (auto& [kn, i]: s.ion_data) {
+        o << " " << kn << "/current_density        " << csv(i.iX_) << "\n"
+          << " " << kn << "/reversal_potential     " << csv(i.eX_) << "\n"
+          << " " << kn << "/internal_concentration " << csv(i.Xi_) << "\n"
+          << " " << kn << "/external_concentration " << csv(i.Xo_) << "\n"
+          << " " << kn << "/intconc_initial        " << csv(i.init_Xi_) << "\n"
+          << " " << kn << "/extconc_initial        " << csv(i.init_Xo_) << "\n"
+          << " " << kn << "/revpot_initial         " << csv(i.init_eX_) << "\n"
+          << " " << kn << "/node_index             " << csv(i.node_index_) << "\n";
     }
     return o;
 }

arbor/backends/gpu/shared_state.hpp

@@ -147,22 +147,23 @@ struct ARB_ARBOR_API shared_state: shared_state_base<shared_state, array, ion_st
     arb_size_type n_detector = 0; // Max number of detectors on all cells.
     arb_size_type n_cv = 0;       // Total number of CVs.
 
-    iarray cv_to_cell;       // Maps CV index to cell index.
-    arb_value_type time;     // integration start time [ms].
-    arb_value_type time_to;  // integration end time [ms]
-    arb_value_type dt;       // 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²].
+    iarray cv_to_cell;            // Maps CV index to cell index.
+    arb_value_type time;          // integration start time [ms].
+    arb_value_type time_to;       // integration end time [ms]
+    arb_value_type dt;            // 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 init_voltage;      // Maps CV index to initial membrane voltage [mV].
-    array temperature_degC;  // Maps CV to local temperature (read only) [°C].
-    array diam_um;           // Maps CV to local diameter (read only) [µm].
+    array init_voltage;           // Maps CV index to initial membrane voltage [mV].
+    array temperature_degC;       // Maps CV to local temperature (read only) [°C].
+    array diam_um;                // Maps CV to local diameter (read only) [µm].
+    array area_um2;               // Maps CV to local diameter (read only) [µm²].
 
-    array time_since_spike;   // Stores time since last spike on any detector, organized by cell.
-    iarray src_to_spike;      // Maps spike source index to spike index
+    array time_since_spike;       // Stores time since last spike on any detector, organized by cell.
+    iarray src_to_spike;          // Maps spike source index to spike index
 
-    arb_seed_type cbprng_seed; // random number generator seed
+    arb_seed_type cbprng_seed;    // random number generator seed
 
     sample_event_stream sample_events;
     array sample_time;
@@ -179,13 +180,43 @@ struct ARB_ARBOR_API shared_state: shared_state_base<shared_state, array, ion_st
 
     shared_state() = default;
 
+    shared_state(task_system_handle tp,
+                 arb_size_type n_cell,
+                 const std::vector<arb_index_type>& cv_to_cell_vec,
+                 const fvm_cv_discretization& D,
+                 const std::vector<arb_index_type>& src_to_spike,
+                 const fvm_detector_info& detector,
+                 const std::unordered_map<std::string, fvm_ion_config>& ions,
+                 const fvm_stimulus_config& stims,
+                 unsigned align,
+                 arb_seed_type cbprng_seed_ = 0u)
+        : shared_state{std::move(tp),
+                       n_cell,
+                       (arb_size_type) D.size(),
+                       cv_to_cell_vec,
+                       D.init_membrane_potential,
+                       D.temperature_K,
+                       D.diam_um,
+                       D.cv_area,
+                       src_to_spike,
+                       detector,
+                       align,
+                       cbprng_seed_}
+    {
+        configure_stimulus(stims);
+        configure_solver(D);
+        add_ions(D, ions);
+    }
+
+
     shared_state(task_system_handle tp,
                  arb_size_type n_cell,
                  arb_size_type n_cv,
                  const std::vector<arb_index_type>& cv_to_cell_vec,
                  const std::vector<arb_value_type>& init_membrane_potential,
                  const std::vector<arb_value_type>& temperature_K,
                  const std::vector<arb_value_type>& diam,
+                 const std::vector<arb_value_type>& area,
                  const std::vector<arb_index_type>& src_to_spike,
                  const fvm_detector_info& detector,
                  unsigned, // align parameter ignored

arbor/backends/multicore/shared_state.cpp

@@ -54,11 +54,9 @@ inline unsigned min_alignment(unsigned align) {
 
 using pad = util::padded_allocator<>;
 
-ion_state::ion_state(
-    int charge,
-    const fvm_ion_config& ion_data,
-    unsigned align,
-    solver_ptr ptr):
+ion_state::ion_state(const fvm_ion_config& ion_data,
+                     unsigned align,
+                     solver_ptr ptr):
     alignment(min_alignment(align)),
     write_eX_(ion_data.revpot_written),
     write_Xo_(ion_data.econc_written),
@@ -75,7 +73,7 @@ ion_state::ion_state(
     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, charge, 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());
@@ -200,6 +198,7 @@ shared_state::shared_state(task_system_handle,    // ignored in mc backend
                            const std::vector<arb_value_type>& init_membrane_potential,
                            const std::vector<arb_value_type>& temperature_K,
                            const std::vector<arb_value_type>& diam,
+                           const std::vector<arb_value_type>& area,
                            const std::vector<arb_index_type>& src_to_spike_,
                            const fvm_detector_info& detector,
                            unsigned align,
@@ -215,6 +214,7 @@ shared_state::shared_state(task_system_handle,    // ignored in mc backend
     init_voltage(init_membrane_potential.begin(), init_membrane_potential.end(), pad(alignment)),
     temperature_degC(n_cv_, pad(alignment)),
     diam_um(diam.begin(), diam.end(), pad(alignment)),
+    area_um2(area.begin(), area.end(), pad(alignment)),
     time_since_spike(n_cell*n_detector, pad(alignment)),
     src_to_spike(src_to_spike_.begin(), src_to_spike_.end(), pad(alignment)),
     cbprng_seed(cbprng_seed_),
@@ -278,27 +278,26 @@ void shared_state::take_samples() {
 ARB_ARBOR_API std::ostream& operator<<(std::ostream& out, const shared_state& s) {
     using io::csv;
 
-    out << "n_cv         " << s.n_cv << "\n";
-    out << "time         " << s.time << "\n";
-    out << "time_to      " << s.time_to << "\n";
-    out << "dt           " << s.dt << "\n";
-    out << "voltage      " << csv(s.voltage) << "\n";
-    out << "init_voltage " << csv(s.init_voltage) << "\n";
-    out << "temperature  " << csv(s.temperature_degC) << "\n";
-    out << "diameter     " << csv(s.diam_um) << "\n";
-    out << "current      " << csv(s.current_density) << "\n";
-    out << "conductivity " << csv(s.conductivity) << "\n";
-    for (const auto& ki: s.ion_data) {
-        auto& kn = ki.first;
-        auto& i = ki.second;
-        out << kn << "/current_density        " << csv(i.iX_) << "\n";
-        out << kn << "/reversal_potential     " << csv(i.eX_) << "\n";
-        out << kn << "/internal_concentration " << csv(i.Xi_) << "\n";
-        out << kn << "/external_concentration " << csv(i.Xo_) << "\n";
-        out << kn << "/intconc_initial        " << csv(i.init_Xi_) << "\n";
-        out << kn << "/extconc_initial        " << csv(i.init_Xo_) << "\n";
-        out << kn << "/revpot_initial         " << csv(i.init_eX_) << "\n";
-        out << kn << "/node_index             " << csv(i.node_index_) << "\n";
+    out << "n_cv         " << s.n_cv << "\n"
+        << "time         " << s.time << "\n"
+        << "time_to      " << s.time_to << "\n"
+        << "dt           " << s.dt << "\n"
+        << "voltage      " << csv(s.voltage) << "\n"
+        << "init_voltage " << csv(s.init_voltage) << "\n"
+        << "temperature  " << csv(s.temperature_degC) << "\n"
+        << "diameter     " << csv(s.diam_um) << "\n"
+        << "area         " << csv(s.area_um2) << "\n"
+        << "current      " << csv(s.current_density) << "\n"
+        << "conductivity " << csv(s.conductivity) << "\n";
+    for (const auto& [kn, i]: s.ion_data) {
+        out << kn << "/current_density        " << csv(i.iX_) << "\n"
+            << kn << "/reversal_potential     " << csv(i.eX_) << "\n"
+            << kn << "/internal_concentration " << csv(i.Xi_) << "\n"
+            << kn << "/external_concentration " << csv(i.Xo_) << "\n"
+            << kn << "/intconc_initial        " << csv(i.init_Xi_) << "\n"
+            << kn << "/extconc_initial        " << csv(i.init_Xo_) << "\n"
+            << kn << "/revpot_initial         " << csv(i.init_eX_) << "\n"
+            << kn << "/node_index             " << csv(i.node_index_) << "\n";
     }
 
     return out;
@@ -406,6 +405,7 @@ void shared_state::instantiate(arb::mechanism& m,
     m.ppack_.vec_g            = conductivity.data();
     m.ppack_.temperature_degC = temperature_degC.data();
     m.ppack_.diam_um          = diam_um.data();
+    m.ppack_.area_um2         = area_um2.data();
     m.ppack_.time_since_spike = time_since_spike.data();
     m.ppack_.n_detectors      = n_detector;
     m.ppack_.events           = {};

arbor/backends/multicore/shared_state.hpp

@@ -76,12 +76,7 @@ struct ARB_ARBOR_API ion_state {
 
     ion_state() = default;
 
-    ion_state(
-        int charge,
-        const fvm_ion_config& ion_data,
-        unsigned align,
-        solver_ptr ptr
-    );
+    ion_state(const fvm_ion_config& ion_data, unsigned align, solver_ptr ptr);
 
     // Set ion concentrations to weighted proportion of default concentrations.
     void init_concentration();
@@ -164,6 +159,7 @@ struct ARB_ARBOR_API shared_state: shared_state_base<shared_state, array, ion_st
     array init_voltage;       // Maps CV index to initial membrane voltage [mV].
     array temperature_degC;   // Maps CV to local temperature (read only) [°C].
     array diam_um;            // Maps CV to local diameter (read only) [µm].
+    array area_um2;           // Maps CV to local lateral surface area (read only) [µm²].
 
     array time_since_spike;   // Stores time since last spike on any detector, organized by cell.
     iarray src_to_spike;      // Maps spike source index to spike index
@@ -192,11 +188,40 @@ struct ARB_ARBOR_API shared_state: shared_state_base<shared_state, array, ion_st
                  const std::vector<arb_value_type>& init_membrane_potential,
                  const std::vector<arb_value_type>& temperature_K,
                  const std::vector<arb_value_type>& diam,
+                 const std::vector<arb_value_type>& area,
                  const std::vector<arb_index_type>& src_to_spike,
                  const fvm_detector_info& detector,
                  unsigned align,
                  arb_seed_type cbprng_seed_ = 0u);
 
+    shared_state(task_system_handle tp,
+                 arb_size_type n_cell,
+                 std::vector<arb_index_type> cv_to_cell_vec,
+                 const fvm_cv_discretization& D,
+                 std::vector<arb_index_type> src_to_spike,
+                 const fvm_detector_info& detector,
+                 std::unordered_map<std::string, fvm_ion_config> ions,
+                 const fvm_stimulus_config& stims,
+                 unsigned align,
+                 arb_seed_type cbprng_seed_ = 0u)
+        : shared_state{std::move(tp),
+                       n_cell,
+                       D.size(),
+                       cv_to_cell_vec,
+                       D.init_membrane_potential,
+                       D.temperature_K,
+                       D.diam_um,
+                       D.cv_area,
+                       src_to_spike,
+                       detector,
+                       align,
+                       cbprng_seed_}
+    {
+        configure_stimulus(stims);
+        configure_solver(D);
+        add_ions(D, ions);
+    }
+
     // Setup a mechanism and tie its backing store to this object
     void instantiate(mechanism&,
                      unsigned,

arbor/backends/shared_state_base.hpp

@@ -49,14 +49,31 @@ struct shared_state_base {
         d->watcher.clear_crossings();
     }
 
+    void configure_solver(const fvm_cv_discretization& disc) {
+        auto d = static_cast<D*>(this);
+        d->solver = {disc.geometry.cv_parent, disc.geometry.cell_cv_divs, disc.cv_capacitance, disc.face_conductance, disc.cv_area};
+    }
+
     void add_ion(const std::string& ion_name,
-                 int charge,
                  const fvm_ion_config& ion_info,
                  typename ion_state::solver_ptr ptr=nullptr) {
         auto d = static_cast<D*>(this);
         d->ion_data.emplace(std::piecewise_construct,
                             std::forward_as_tuple(ion_name),
-                            std::forward_as_tuple(charge, ion_info, d->alignment, std::move(ptr)));
+                            std::forward_as_tuple(ion_info, d->alignment, std::move(ptr)));
+    }
+
+    void add_ions(const fvm_cv_discretization& disc,
+                  const std::unordered_map<std::string, fvm_ion_config>& ions) {
+        auto d = static_cast<D*>(this);
+        for (const auto& [ion, data]: ions) {
+            std::unique_ptr<typename ion_state::solver_type> solver = nullptr;
+            if (data.is_diffusive) solver = std::make_unique<typename ion_state::solver_type>(disc.geometry.cv_parent,
+                                                                                              disc.geometry.cell_cv_divs,
+                                                                                              data.face_diffusivity,
+                                                                                              disc.cv_area);
+            d->add_ion(ion, data, std::move(solver));
+        }
     }
 
     arb_value_type* mechanism_state_data(const mechanism& m,
@@ -103,8 +120,10 @@ struct shared_state_base {
     }
 
     void configure_stimulus(const fvm_stimulus_config& stims) {
-        auto d = static_cast<D*>(this);
-        d->stim_data = {stims, d->alignment};
+        if (!stims.cv.empty()) {
+            auto d = static_cast<D*>(this);
+            d->stim_data = {stims, d->alignment};
+        }
     }
 
     void add_stimulus_current() {