From d7778e2bb64729e3d6779685f0c133b7450e5bec Mon Sep 17 00:00:00 2001 From: Michael Hines Date: Thu, 2 Jan 2025 18:54:31 -0500 Subject: [PATCH 1/7] Lvardt minimize size of CvMembList.ml Memb_list vector. --- src/nrncvode/netcvode.cpp | 112 +++++++++++++++++++++++++++----------- 1 file changed, 81 insertions(+), 31 deletions(-) diff --git a/src/nrncvode/netcvode.cpp b/src/nrncvode/netcvode.cpp index 94df2cdf8b..2470c67d08 100644 --- a/src/nrncvode/netcvode.cpp +++ b/src/nrncvode/netcvode.cpp @@ -1615,6 +1615,9 @@ bool NetCvode::init_global() { } } + // Modified to also count the nodes and set the offsets for + // each CvMembList.ml[contig_region] + // See the long comment after this loop. for (NrnThreadMembList* tml = _nt->tml; tml; tml = tml->next) { i = tml->index; const Memb_func& mf = memb_func[i]; @@ -1625,42 +1628,77 @@ bool NetCvode::init_global() { // singly linked list built below int j; for (j = 0; j < ml->nodecount; ++j) { + auto offset = ml->get_storage_offset() + j; + // for each Memb_list instance constructed, keep + // track of its initial storage offset (i.e. offset) + int inode = ml->nodelist[j]->v_node_index; + int icell = cellnum[inode]; Cvode& cv = d.lcv_[cellnum[inode]]; CvodeThreadData& z = cv.ctd_[0]; - if (!z.cv_memb_list_) { + + // Circumstances for creating a new CvMembList + // or (due to non-contiguity of a cell), + // appending a Memb_list instance to cml->ml + if (!z.cv_memb_list_) { // initialize the first cml = new CvMembList{i}; - cml->next = nullptr; - assert(cml->ml.size() == 1); - cml->ml[0].nodecount = 0; z.cv_memb_list_ = cml; + cml->next = nullptr; last[cellnum[inode]] = cml; - } - if (last[cellnum[inode]]->index == i) { - assert(last[cellnum[inode]]->ml.size() == 1); - ++last[cellnum[inode]]->ml[0].nodecount; - } else { + assert(cml->ml.size() == 1); + assert(cml->ml[0].nodecount == 0); + } else if (last[cellnum[inode]]->index != i) { // initialize next cml = new CvMembList{i}; last[cellnum[inode]]->next = cml; cml->next = nullptr; last[cellnum[inode]] = cml; assert(cml->ml.size() == 1); - cml->ml[0].nodecount = 1; + assert(cml->ml[0].nodecount == 0); + } else { // if non-contiguous, append Memb_list + cml = last[cellnum[inode]]; + auto& cvml = cml->ml.back(); + auto cvml_offset = cvml.get_storage_offset() + cvml.nodecount; + if (cvml_offset != offset) { + // not contiguous, add another Memb_list + // instance to cml->ml + cml->ml.emplace_back(cml->index); + assert(cml->ml.back().nodecount == 0); + } } + + auto& cvml = cml->ml.back(); + if (cvml.nodecount == 0) { //first time for this Memb_List + cvml.set_storage_offset(offset); + } + // Increment count of last Memb_list in cml->ml. + ++cvml.nodecount; } } } - // allocate and re-initialize count + // each NetCvode.p[id].lcv_[nt->ncell].ctd_[0].cv_memb_list_->ml[0].nodecount + // is setup and ml.size() == 1. + + // if ml[0] data is not contiguous, need to count number + // of contiguous data regions and reserve that number of Memb_list + // in the ml vector. The sum of the ml[i].nodecount must equal + // the present ml[0].nodecount and each ml[i] data must be + // contiguous. + // Ideally the node permutation would be such that each cell + // is contiguous. That is almost the case with the default + // permutation. But the cell root nodes are all at the + // beginning, and thereafter the rest of the cell is contiguous. + // For all density mechanims (except extracellular) this is fine + // since those mechanisms do not exist at the root node. And + // most of the time, there are no POINT_PROCESSes in the + // root nodes. So, we expect that most of the time, we have + // the default node permutation and that will result in + // almost all CvMembList.ml.size() == 1 with a very occasional + // size() == 2 + + std::vector cvml(d.nlcv_); for (i = 0; i < d.nlcv_; ++i) { - cvml[i] = d.lcv_[i].ctd_[0].cv_memb_list_; - for (cml = cvml[i]; cml; cml = cml->next) { - // non-contiguous mode, so we're going to create a lot of 1-element Memb_list - // inside cml->ml - cml->ml.reserve(cml->ml[0].nodecount); - // remove the single entry from contiguous mode - cml->ml.clear(); - } + cvml[i] = d.lcv_[i].ctd_[0].cv_memb_list_; // whole cell in thread } // fill pointers (and nodecount) // now list order is from 0 to n_memb_func @@ -1674,20 +1712,32 @@ bool NetCvode::init_global() { int icell = cellnum[ml->nodelist[j]->v_node_index]; if (cvml[icell]->index != i) { cvml[icell] = cvml[icell]->next; - assert(cvml[icell] && cvml[icell]->index); + assert(cvml[icell] && cvml[icell]->index == i); } - cml = cvml[icell]; - auto& newml = cml->ml.emplace_back(cml->index /* mechanism type */); - newml.nodecount = 1; - newml.nodelist = new Node*[1]; - newml.nodelist[0] = ml->nodelist[j]; - newml.nodeindices = new int[1]{ml->nodeindices[j]}; - newml.prop = new Prop* [1] { ml->prop[j] }; - if (!mf.hoc_mech) { - newml.set_storage_offset(ml->get_storage_offset() + j); - newml.pdata = new Datum* [1] { ml->pdata[j] }; + auto& cml = cvml[icell]; + int kk = j; + for (auto& newml: cml->ml) { + auto nodecount = newml.nodecount; + if (! newml.nodelist) { + newml.nodelist = new Node*[nodecount]; + newml.nodeindices = new int[nodecount]; + newml.prop = new Prop* [nodecount]; + if (!mf.hoc_mech) { + newml.pdata = new Datum* [nodecount]; + } + for (int k = 0; k < nodecount; ++k) { + newml.nodelist[k] = ml->nodelist[kk+k]; + newml.nodeindices[k] = ml->nodeindices[kk+k]; + assert (cellnum[newml.nodeindices[k]] == cellnum[ml->nodeindices[j]]); + newml.prop[k] = ml->prop[kk+k]; + if (!mf.hoc_mech) { + newml.pdata[k] = ml->pdata[kk+k]; + } + } + kk += nodecount; + newml._thread = ml->_thread; + } } - newml._thread = ml->_thread; } } } From 5c1537b4a479865782671755a7c906f2207ba05b Mon Sep 17 00:00:00 2001 From: Michael Hines Date: Fri, 3 Jan 2025 12:58:41 -0500 Subject: [PATCH 2/7] For lvardt, allow CvMembList.ml.size() > 1 and nodecount > 1 --- src/nrncvode/cvodeobj.h | 3 +- src/nrncvode/netcvode.cpp | 70 +++++++++++++++++++-------------------- src/nrncvode/occvode.cpp | 18 ++++++---- 3 files changed, 47 insertions(+), 44 deletions(-) diff --git a/src/nrncvode/cvodeobj.h b/src/nrncvode/cvodeobj.h index eb8b52a665..3b687bece6 100644 --- a/src/nrncvode/cvodeobj.h +++ b/src/nrncvode/cvodeobj.h @@ -30,7 +30,8 @@ struct model_sorted_token; * contiguous * - with ml.size() >= 1 and ml[i].nodecount == 1 when non-contiguous instances need to be processed * - * generic configurations with ml.size() and ml[i].nodecount both larger than one are not supported. + * generic configurations with ml.size() and ml[i].nodecount both larger than one are only + * supported for the local variable time step method. */ struct CvMembList { CvMembList(int type) diff --git a/src/nrncvode/netcvode.cpp b/src/nrncvode/netcvode.cpp index 2470c67d08..59139f91b1 100644 --- a/src/nrncvode/netcvode.cpp +++ b/src/nrncvode/netcvode.cpp @@ -1616,8 +1616,19 @@ bool NetCvode::init_global() { } // Modified to also count the nodes and set the offsets for - // each CvMembList.ml[contig_region] - // See the long comment after this loop. + // each CvMembList.ml[contig_region]. + // The sum of the ml[i].nodecount must equal the mechanism + // nodecount for the cell and each ml[i] data must be contiguous. + // Ideally the node permutation would be such that each cell + // is contiguous. So only needing a ml[0]. That is sadly not + // the case with the default permutation. The cell root nodes are + // all at the beginning, and thereafter only Section nodes are + // contiguous. It would be easy to permute nodes so that each cell + // is contiguous (except root node). This would result in a + // CvMembList.ml.size() == 1 almost always with an exception of + // size() == 2 only for extracellular and for POINT_PROCESSes + // located both in the root node and other cell nodes. + for (NrnThreadMembList* tml = _nt->tml; tml; tml = tml->next) { i = tml->index; const Memb_func& mf = memb_func[i]; @@ -1638,23 +1649,23 @@ bool NetCvode::init_global() { CvodeThreadData& z = cv.ctd_[0]; // Circumstances for creating a new CvMembList - // or (due to non-contiguity of a cell), + // or (due to non-contiguity of a cell), // appending a Memb_list instance to cml->ml - if (!z.cv_memb_list_) { // initialize the first + if (!z.cv_memb_list_) { // initialize the first cml = new CvMembList{i}; z.cv_memb_list_ = cml; cml->next = nullptr; last[cellnum[inode]] = cml; assert(cml->ml.size() == 1); assert(cml->ml[0].nodecount == 0); - } else if (last[cellnum[inode]]->index != i) { // initialize next + } else if (last[cellnum[inode]]->index != i) { // initialize next cml = new CvMembList{i}; last[cellnum[inode]]->next = cml; cml->next = nullptr; last[cellnum[inode]] = cml; assert(cml->ml.size() == 1); assert(cml->ml[0].nodecount == 0); - } else { // if non-contiguous, append Memb_list + } else { // if non-contiguous, append Memb_list cml = last[cellnum[inode]]; auto& cvml = cml->ml.back(); auto cvml_offset = cvml.get_storage_offset() + cvml.nodecount; @@ -1667,7 +1678,7 @@ bool NetCvode::init_global() { } auto& cvml = cml->ml.back(); - if (cvml.nodecount == 0) { //first time for this Memb_List + if (cvml.nodecount == 0) { // first time for this Memb_List cvml.set_storage_offset(offset); } // Increment count of last Memb_list in cml->ml. @@ -1675,30 +1686,10 @@ bool NetCvode::init_global() { } } } - // each NetCvode.p[id].lcv_[nt->ncell].ctd_[0].cv_memb_list_->ml[0].nodecount - // is setup and ml.size() == 1. - - // if ml[0] data is not contiguous, need to count number - // of contiguous data regions and reserve that number of Memb_list - // in the ml vector. The sum of the ml[i].nodecount must equal - // the present ml[0].nodecount and each ml[i] data must be - // contiguous. - // Ideally the node permutation would be such that each cell - // is contiguous. That is almost the case with the default - // permutation. But the cell root nodes are all at the - // beginning, and thereafter the rest of the cell is contiguous. - // For all density mechanims (except extracellular) this is fine - // since those mechanisms do not exist at the root node. And - // most of the time, there are no POINT_PROCESSes in the - // root nodes. So, we expect that most of the time, we have - // the default node permutation and that will result in - // almost all CvMembList.ml.size() == 1 with a very occasional - // size() == 2 - std::vector cvml(d.nlcv_); for (i = 0; i < d.nlcv_; ++i) { - cvml[i] = d.lcv_[i].ctd_[0].cv_memb_list_; // whole cell in thread + cvml[i] = d.lcv_[i].ctd_[0].cv_memb_list_; // whole cell in thread } // fill pointers (and nodecount) // now list order is from 0 to n_memb_func @@ -1718,24 +1709,31 @@ bool NetCvode::init_global() { int kk = j; for (auto& newml: cml->ml) { auto nodecount = newml.nodecount; - if (! newml.nodelist) { + if (!newml.nodelist) { + // do nodecount of these for ml and then + // skip forward by nodecount in the outer + // ml->nodecount j loop (i.e. a contiguity + // region) newml.nodelist = new Node*[nodecount]; newml.nodeindices = new int[nodecount]; - newml.prop = new Prop* [nodecount]; + newml.prop = new Prop*[nodecount]; if (!mf.hoc_mech) { - newml.pdata = new Datum* [nodecount]; + newml.pdata = new Datum*[nodecount]; } for (int k = 0; k < nodecount; ++k) { - newml.nodelist[k] = ml->nodelist[kk+k]; - newml.nodeindices[k] = ml->nodeindices[kk+k]; - assert (cellnum[newml.nodeindices[k]] == cellnum[ml->nodeindices[j]]); - newml.prop[k] = ml->prop[kk+k]; + newml.nodelist[k] = ml->nodelist[kk + k]; + newml.nodeindices[k] = ml->nodeindices[kk + k]; + assert(cellnum[newml.nodeindices[k]] == + cellnum[ml->nodeindices[j]]); + newml.prop[k] = ml->prop[kk + k]; if (!mf.hoc_mech) { - newml.pdata[k] = ml->pdata[kk+k]; + newml.pdata[k] = ml->pdata[kk + k]; } } kk += nodecount; newml._thread = ml->_thread; + j += nodecount - 1; + break; } } } diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 630e1d733e..6589c0ec6c 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -176,7 +176,8 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize= if (z.cmlcap_) { for (auto& ml: z.cmlcap_->ml) { // support `1 x n` and `n x 1` but not `n x m` - assert(z.cmlcap_->ml.size() == 1 || ml.nodecount == 1); + // why not? (I'm trying to improve lvardt performance) + // assert(z.cmlcap_->ml.size() == 1 || ml.nodecount == 1); zneq_cap_v += ml.nodecount; } } @@ -207,12 +208,15 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize= // sentinal values for determining no_cap NODERHS(z.v_node_[i]) = 1.; } - for (i = 0; i < zneq_cap_v; ++i) { - auto* const node = z.cmlcap_->ml.size() == 1 ? z.cmlcap_->ml[0].nodelist[i] - : z.cmlcap_->ml[i].nodelist[0]; - z.pv_[i] = node->v_handle(); - z.pvdot_[i] = node->rhs_handle(); - *z.pvdot_[i] = 0.; // only ones = 1 are no_cap + i = 0; + for (auto& ml: z.cmlcap_->ml) { + for (int j = 0; j < ml.nodecount; ++j) { + auto* const node = ml.nodelist[j]; + z.pv_[i] = node->v_handle(); + z.pvdot_[i] = node->rhs_handle(); + *z.pvdot_[i] = 0.; // only ones = 1 are no_cap + ++i; + } } // the remainder are no_cap nodes From 01b41496a28cef28d842c2b3786d096a9e79b9ae Mon Sep 17 00:00:00 2001 From: Michael Hines Date: Sat, 4 Jan 2025 17:15:34 -0500 Subject: [PATCH 3/7] fix lvardt segfault when a cvode instance is empty. --- src/nrncvode/occvode.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 6589c0ec6c..1ab5c08ce0 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -209,7 +209,7 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize= NODERHS(z.v_node_[i]) = 1.; } i = 0; - for (auto& ml: z.cmlcap_->ml) { + if (zneq_cap_v) for (auto& ml: z.cmlcap_->ml) { for (int j = 0; j < ml.nodecount; ++j) { auto* const node = ml.nodelist[j]; z.pv_[i] = node->v_handle(); From d73896e341ad3ad580efa303b0b6827e6f83b3f3 Mon Sep 17 00:00:00 2001 From: Michael Hines Date: Sun, 5 Jan 2025 07:42:03 -0500 Subject: [PATCH 4/7] clang format --- src/nrncvode/occvode.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index 1ab5c08ce0..b88c42ca64 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -209,15 +209,16 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize= NODERHS(z.v_node_[i]) = 1.; } i = 0; - if (zneq_cap_v) for (auto& ml: z.cmlcap_->ml) { - for (int j = 0; j < ml.nodecount; ++j) { - auto* const node = ml.nodelist[j]; - z.pv_[i] = node->v_handle(); - z.pvdot_[i] = node->rhs_handle(); - *z.pvdot_[i] = 0.; // only ones = 1 are no_cap - ++i; + if (zneq_cap_v) + for (auto& ml: z.cmlcap_->ml) { + for (int j = 0; j < ml.nodecount; ++j) { + auto* const node = ml.nodelist[j]; + z.pv_[i] = node->v_handle(); + z.pvdot_[i] = node->rhs_handle(); + *z.pvdot_[i] = 0.; // only ones = 1 are no_cap + ++i; + } } - } // the remainder are no_cap nodes if (z.no_cap_node_) { From 11ba674a7408edd054f9ca04aa883689a0c3f2f8 Mon Sep 17 00:00:00 2001 From: Michael Hines Date: Tue, 7 Jan 2025 13:15:29 -0500 Subject: [PATCH 5/7] incorporate reviewer suggestions --- src/nrncvode/netcvode.cpp | 18 +++++++++--------- src/nrncvode/occvode.cpp | 3 --- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/nrncvode/netcvode.cpp b/src/nrncvode/netcvode.cpp index 59139f91b1..b3a88385b6 100644 --- a/src/nrncvode/netcvode.cpp +++ b/src/nrncvode/netcvode.cpp @@ -1699,21 +1699,23 @@ bool NetCvode::init_global() { Memb_list* ml = tml->ml; if (ml->nodecount && (mf.current || mf.ode_count || mf.ode_matsol || mf.ode_spec || mf.state || i == CAP || ba_candidate.count(i) == 1)) { - for (int j = 0; j < ml->nodecount; ++j) { + int increment = 1; // newml.nodecount is handled in the newml loop below + for (int j = 0; j < ml->nodecount; j += increment) { int icell = cellnum[ml->nodelist[j]->v_node_index]; if (cvml[icell]->index != i) { cvml[icell] = cvml[icell]->next; assert(cvml[icell] && cvml[icell]->index == i); } auto& cml = cvml[icell]; - int kk = j; + increment = 1; for (auto& newml: cml->ml) { - auto nodecount = newml.nodecount; if (!newml.nodelist) { + auto nodecount = newml.nodecount; // do nodecount of these for ml and then // skip forward by nodecount in the outer // ml->nodecount j loop (i.e. a contiguity // region) + increment = nodecount; newml.nodelist = new Node*[nodecount]; newml.nodeindices = new int[nodecount]; newml.prop = new Prop*[nodecount]; @@ -1721,18 +1723,16 @@ bool NetCvode::init_global() { newml.pdata = new Datum*[nodecount]; } for (int k = 0; k < nodecount; ++k) { - newml.nodelist[k] = ml->nodelist[kk + k]; - newml.nodeindices[k] = ml->nodeindices[kk + k]; + newml.nodelist[k] = ml->nodelist[j + k]; + newml.nodeindices[k] = ml->nodeindices[j + k]; assert(cellnum[newml.nodeindices[k]] == cellnum[ml->nodeindices[j]]); - newml.prop[k] = ml->prop[kk + k]; + newml.prop[k] = ml->prop[j + k]; if (!mf.hoc_mech) { - newml.pdata[k] = ml->pdata[kk + k]; + newml.pdata[k] = ml->pdata[j + k]; } } - kk += nodecount; newml._thread = ml->_thread; - j += nodecount - 1; break; } } diff --git a/src/nrncvode/occvode.cpp b/src/nrncvode/occvode.cpp index b88c42ca64..8cb2befbd6 100644 --- a/src/nrncvode/occvode.cpp +++ b/src/nrncvode/occvode.cpp @@ -175,9 +175,6 @@ printf("%d Cvode::init_eqn id=%d neq_v_=%d #nonvint=%d #nonvint_extra=%d nvsize= zneq_cap_v = 0; if (z.cmlcap_) { for (auto& ml: z.cmlcap_->ml) { - // support `1 x n` and `n x 1` but not `n x m` - // why not? (I'm trying to improve lvardt performance) - // assert(z.cmlcap_->ml.size() == 1 || ml.nodecount == 1); zneq_cap_v += ml.nodecount; } } From ac1c6a9f2bc7808565e496a4c60d950438e9b5d0 Mon Sep 17 00:00:00 2001 From: Michael Hines Date: Tue, 7 Jan 2025 14:51:23 -0500 Subject: [PATCH 6/7] Full coverage for this PR. --- test/cover/test_netcvode.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/test/cover/test_netcvode.py b/test/cover/test_netcvode.py index 19ac3f2612..2ca2138b1b 100644 --- a/test/cover/test_netcvode.py +++ b/test/cover/test_netcvode.py @@ -12,6 +12,7 @@ cv = h.CVode() pc = h.ParallelContext() + # remove address info from cv.debug_event output def debug_event_filter(s): s = re.sub(r"cvode_0x[0-9abcdef]* ", "cvode_0x... ", s) @@ -762,6 +763,21 @@ def nc_event_before_init(): expect_err("nc.event(0)") # nrn_assert triggered if outside of finitialize +def contiguous(): + # cover a couple of lines in the lvardt part of init_global(). + # Need same POINT_PROCESS type in rootnode and somewhere else on cell + net = Net(8) + syns = [h.ExpSyn(seg) for cell in net.cells for seg in cell.soma.allseg()] + cv.active(1) + cv.use_local_dt(1) + h.finitialize(-65) + pc.nthread(2) + h.finitialize(-65) + pc.nthread(1) + cv.use_local_dt(0) + cv.active(0) + + def test_netcvode_cover(): nrn_use_daspk() node() @@ -774,6 +790,7 @@ def test_netcvode_cover(): scatter_gather() playrecord() interthread() + contiguous() nc_event_before_init() From 32921b2836bd7a50734a3bc5dabab44f615cdee9 Mon Sep 17 00:00:00 2001 From: Michael Hines Date: Wed, 8 Jan 2025 09:50:07 -0500 Subject: [PATCH 7/7] accept review suggestion about icell --- src/nrncvode/netcvode.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/nrncvode/netcvode.cpp b/src/nrncvode/netcvode.cpp index b3a88385b6..6e4522c852 100644 --- a/src/nrncvode/netcvode.cpp +++ b/src/nrncvode/netcvode.cpp @@ -1645,7 +1645,7 @@ bool NetCvode::init_global() { int inode = ml->nodelist[j]->v_node_index; int icell = cellnum[inode]; - Cvode& cv = d.lcv_[cellnum[inode]]; + Cvode& cv = d.lcv_[icell]; CvodeThreadData& z = cv.ctd_[0]; // Circumstances for creating a new CvMembList @@ -1655,18 +1655,18 @@ bool NetCvode::init_global() { cml = new CvMembList{i}; z.cv_memb_list_ = cml; cml->next = nullptr; - last[cellnum[inode]] = cml; + last[icell] = cml; assert(cml->ml.size() == 1); assert(cml->ml[0].nodecount == 0); - } else if (last[cellnum[inode]]->index != i) { // initialize next + } else if (last[icell]->index != i) { // initialize next cml = new CvMembList{i}; - last[cellnum[inode]]->next = cml; + last[icell]->next = cml; cml->next = nullptr; - last[cellnum[inode]] = cml; + last[icell] = cml; assert(cml->ml.size() == 1); assert(cml->ml[0].nodecount == 0); } else { // if non-contiguous, append Memb_list - cml = last[cellnum[inode]]; + cml = last[icell]; auto& cvml = cml->ml.back(); auto cvml_offset = cvml.get_storage_offset() + cvml.nodecount; if (cvml_offset != offset) {