MULTI_ION_IF.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
// TODO remove the "get_sv_address()" and "get_sv_size()" calls by moving
// the functions that use them to the IonIf class
 
#include "MULTI_ION_IF.h"
#include "limpet_types.h"
#include "mpi_utils.h"
#include <stdio.h>
#ifdef HAS_CUDA_MODEL
#include <cuda_runtime.h>
#endif
#if defined HAS_ROCM_MODEL && defined __HIP__
#include <hip/hip_runtime.h>
#endif
 
#include "SF_init.h" // for SF::init_xxx()
#include "petsc_utils.h" // TODO for EXIT
 
namespace limpet {
 
using ::opencarp::base_timer;
using ::opencarp::dupstr;
using ::opencarp::f_open;
using ::opencarp::FILE_SPEC;
using ::opencarp::get_global;
using ::opencarp::get_rank;
using ::opencarp::get_size;
using ::opencarp::get_time;
using ::opencarp::log_msg;
using ::opencarp::printable_int;
using ::opencarp::mesh_t;
using ::opencarp::num_msh;
using ::opencarp::Salt_list;
using ::opencarp::sf_mesh;
using ::opencarp::sf_vec;
using ::opencarp::timing;
using ::opencarp::unset_msh;
#ifdef USE_FMEM_WRAPPER
using ::opencarp::fmemopen_;
#endif
 
/*
   PROTOTYPES
 */
void    initialize_params_MIIF(MULTI_IF *pMIIF);
void    initialize_ionic_IF(MULTI_IF *pMIIF);
void    CreateIIFLocalNodeLsts(MULTI_IF *pMIIF);
void    CreateIIFNodeLsts_(int, IIF_Mask_t *, node_count_t **, node_index_t ***, node_count_t);
void    alloc_MIIF(MULTI_IF *pMIIF);
void    initializeIMPData(MULTI_IF *pMIIF);
void    freeIMPData(MULTI_IF *pMIIF);
void    allocate_shared_data(MULTI_IF *);
node_index_t getGlobalNodalIndex(IonIfBase& pIF, node_index_t relIdx);
void    CreateIIFGlobalNodeLsts(MULTI_IF *pMIIF);
bool    isIMPdata(const char *);
 
static MULTI_IF *gMIIF_Error_Recovery;
static int current_IIF_index;
static float current_time = 0;
static float start_time = 0;
 
static void prepare_error_recovery(MULTI_IF *MIIF, int IIF_index, float time) {
  gMIIF_Error_Recovery = MIIF;
  current_IIF_index    = IIF_index;
  current_time         = start_time+time;
}
 
node_index_t current_global_node(node_index_t local_node) {
  return gMIIF_Error_Recovery->NodeLists[current_IIF_index][local_node];
}
 
float current_global_time() {
  return current_time;
}
 
static int g_print_bounds_exceeded = 1;
 
int should_print_bounds_exceeded_messages() {
  return g_print_bounds_exceeded;
}
 
int set_print_bounds_exceeded_messages(int newval) {
  int oldval = g_print_bounds_exceeded;
 
  g_print_bounds_exceeded = newval;
  return oldval;
}
 
void CreateIIFLocalNodeLsts(MULTI_IF *pMIIF) {
  CreateIIFNodeLsts_(pMIIF->N_IIF, pMIIF->IIFmask, &pMIIF->N_Nodes,
                     &pMIIF->NodeLists, pMIIF->numNode);
}
 
void CreateIIFNodeLsts_(int N_IIF, IIF_Mask_t *IIF_Mask, node_count_t **N_Nodes,
                        node_index_t ***NodeLists, node_count_t numNode)
{
  /*  create node lists for each IIF
      NodeNum: store number of nodes for each IIF
      NodeLst: store node indices of each IIF
   */
  node_count_t *NodeNum  = static_cast<node_count_t *>(calloc(N_IIF, sizeof(node_count_t)));
  node_index_t **NodeLst = static_cast<node_index_t **>(calloc(N_IIF, sizeof(node_index_t *)));
 
  // determine number of nodes of each IIF
  for (node_index_t i = 0; i < numNode; i++)
    NodeNum[static_cast<int>(IIF_Mask[i])]++;
 
  /* now we know the number of nodes per IIF and store the indices of
     nodes belonging to a particular IIF into a list.
   */
  for (int i = 0; i < N_IIF; i++) {
    if (NodeNum[i] > 0)
      NodeLst[i] = static_cast<node_index_t *>(malloc(sizeof(node_index_t) * NodeNum[i]));
    else
      NodeLst[i] = NULL;
  }
 
  // just store indices relative to low
  for (int j = 0; j < N_IIF; j++) {
    node_count_t hcount = 0;
    for (node_index_t i = 0; i < numNode; i++)
      if (IIF_Mask[i] == j)
        NodeLst[j][hcount++] = i;
  }
 
  for(int j=0; j<N_IIF; j++)
    std::sort(NodeLst[j], NodeLst[j]+NodeNum[j]);
 
  *N_Nodes   = NodeNum;
  *NodeLists = NodeLst;
 
}  // CreateIIFNodeLsts_
 
 
node_index_t getGlobalNodalIndex(IonIfBase& pIF, node_index_t rIdx) {
  // get pointer to parent MIIF
  MULTI_IF *pMIIF = (MULTI_IF *) pIF.parent(); // TODO check?!?!?!?!?!
 
  if (pIF.get_type().is_plugin())
    pMIIF = (MULTI_IF *) pIF.parent()->parent();
 
  return pMIIF->NodeLists[pIF.miifIdx][rIdx];
}
 
void alloc_MIIF(MULTI_IF *pMIIF) {
  int N_IIF = pMIIF->N_IIF;
 
  pMIIF->contiguous = static_cast<bool *>(calloc(N_IIF, sizeof(bool)));
#if 0
  pMIIF->ldata = (GlobalData_t***) build_matrix_ns<GlobalData_t>(N_IIF, NUM_IMP_DATA_TYPES, sizeof(GlobalData_t **), Target::CPU);
#else
  pMIIF->ldata = allocate_on_target<GlobalData_t**>(Target::CPU, N_IIF);
  for (std::size_t i = 0; i < pMIIF->N_IIF; i++) {
    pMIIF->ldata[i] = allocate_on_target<GlobalData_t*>(pMIIF->iontypes[i].get().select_target(pMIIF->targets[i]), NUM_IMP_DATA_TYPES);
  }
#endif
 
  for (int i = 0; i < pMIIF->N_IIF; i++) {
    // Create an IonIf object on the chosen target
    pMIIF->IIF.push_back(pMIIF->iontypes[i].get().make_ion_if(pMIIF->targets[i], pMIIF->N_Nodes[i], pMIIF->plugtypes[i]));
    // this is questionable. A IIF does not have a parent, parent should be NULL
    // It might make sense in some situtation to be able to refer to the MIIF
    // structure, but pIF->parent is of type IonIf and not MULTI_IF!!!!
    pMIIF->IIF[i]->set_parent((IonIfBase *)pMIIF);
 
    // set unique miifIdx for imp and plugins to enable refering back
    // to global entities such as global node index from within imp
    pMIIF->IIF[i]->for_each([&](IonIfBase& IF) { IF.miifIdx = i; });
  }
}  // alloc_MIIF
 
void initialize_params_MIIF(MULTI_IF *pMIIF) {
  for (auto& imp : pMIIF->IIF) {
    imp->initialize_params();
  }
}
 
void initialize_ionic_IF(MULTI_IF *pMIIF) {
  pMIIF->getRealData();         // needed for initializing Vm, etc.
  for (int i = 0; i < pMIIF->IIF.size(); i++)
    pMIIF->IIF[i]->initialize(pMIIF->dt, pMIIF->ldata[i]);
 
  pMIIF->releaseRealDataDuringInit();
}
 
void MULTI_IF::releaseRealDataDuringInit() {
  // we need to copy all global data which has been set, even if it is not
  // listed as moddat. This is the case with Vm which may not be modified by
  // LIMPET but is initialized by LIMPET
  unsigned int *moddat = static_cast<unsigned int *>(calloc(this->N_IIF, sizeof(unsigned int)));
 
  for (int i = 0; i < this->N_IIF; i++) {
    moddat[i]        = this->IIF[i]->get_moddat();
    this->IIF[i]->set_moddat(moddat[i] | this->IIF[i]->get_reqdat());  // add in the reqdat
  }
  this->releaseRealData();
  for (int i = 0; i < this->N_IIF; i++) {
    this->IIF[i]->set_moddat(moddat[i]);  // remove the reqdat
  }
  free(moddat);
}
 
void MULTI_IF::initialize_MIIF()
{
  CreateIIFLocalNodeLsts(this);
  alloc_MIIF(this);
  initialize_params_MIIF(this);
 
  memset(&svd, 0, sizeof(SV_DUMP) );
  svd.active = 0;
  svd.intv   = 1.0;         // 1. ms by default
}
 
#define FILENAME_BUF 1024
 
int node_index_cmp(const void *a, const void *b) {
  const node_index_t lhs = *static_cast<const node_index_t *>(a);
  const node_index_t rhs = *static_cast<const node_index_t *>(b);
  return (lhs > rhs) - (lhs < rhs);
}
 
void open_trace(MULTI_IF *MIIF, int n_traceNodes, int *traceNodes, int *label, opencarp::sf_mesh* imesh)
{
  if (!n_traceNodes) return;
 
  if (n_traceNodes > 1000)
    log_msg(0, 4, 0, "%s warning: %d trace nodes may impact performance", __func__, n_traceNodes);
 
  MIIF->trace_info                      = (Trace_Info *)calloc(n_traceNodes+1, sizeof(Trace_Info));
  MIIF->trace_info[n_traceNodes].region = -1;  // end of list marker
 
  // bijective index mapping between set A (local petsc indexing) and set B (local nodal indexing)
  SF::index_mapping<mesh_int_t> petsc2nod;
  if(imesh)
    SF::local_petsc_to_nodal_mapping(*imesh, petsc2nod);
 
  for (int iTrace = 0; iTrace < n_traceNodes; iTrace++) {
    mesh_int_t lnode = imesh ? imesh->pl.localize(traceNodes[iTrace]) : traceNodes[iTrace];
 
    // here we check if lnode is in set B (i.e. local nodal indexing)
    if(imesh) {
      // here we check if lnode is in set B (i.e. local nodal indexing)
      if(petsc2nod.in_b(lnode) == false) continue;
 
      // lnode is local nodal, we map it to local petsc
      lnode = petsc2nod.backward_map(lnode);
    }
 
    for (int iRegion = 0; iRegion < MIIF->N_IIF; iRegion++)
    {
      node_index_t local_node = static_cast<node_index_t>(lnode);
      auto target = static_cast<node_index_t *>(bsearch(&local_node, MIIF->NodeLists[iRegion],
                                   MIIF->N_Nodes[iRegion], sizeof(node_index_t), node_index_cmp));
 
      if (target != NULL ) {
        node_index_t idx = target - MIIF->NodeLists[iRegion];
        Trace_Info *trace_info = MIIF->trace_info+iTrace;
        trace_info->found    = true;
        trace_info->region   = iRegion;
        trace_info->node_idx = idx;
      }
    }
  }
 
  for (int iTrace = 0; iTrace < n_traceNodes; iTrace++) {
    if(get_global(int(MIIF->trace_info[iTrace].found), MPI_SUM) == 0) {
      MIIF->trace_info[iTrace].ignored = true;
      log_msg(0,4,0, "trace node %d not found", traceNodes[iTrace]);
      continue;
    }
 
    char traceName[FILENAME_BUF];
    snprintf(traceName, sizeof traceName, "Trace_%d.dat", label ? label[iTrace] : traceNodes[iTrace]);
    MIIF->trace_info[iTrace].file = f_open(traceName, "w");
  }
}  // open_trace
 
void dump_trace(MULTI_IF *MIIF, limpet::Real time) {
  if (MIIF->trace_info == NULL)
    return;
 
//  struct exception_type e;
 
#define MAX_TRACE_LINE_LEN 8196
 
  std::vector<IonIfBase*>& IIF = MIIF->IIF;
  FILE   *fs;
 
  for (int iTrace = 0; MIIF->trace_info[iTrace].region >= 0; iTrace++) {
    Trace_Info *ctrace = MIIF->trace_info+iTrace;
 
    if (ctrace->ignored)
      continue;
 
    if (ctrace->found) {
      // The buffer is printed with %s, so it has to be zeroed for every line: the
      // BSD fmemopen_ fallback does not append a terminating null byte on close.
      char trace_buf[MAX_TRACE_LINE_LEN] = {0};
 
      fs = fmemopen_(trace_buf, MAX_TRACE_LINE_LEN, "w");
      fprintf(fs, "%4.10f\t", time);
      if (IIF[ctrace->region]->get_type().has_trace()) {
        IIF[ctrace->region]->get_type().trace(*IIF[ctrace->region], ctrace->node_idx,
                                        fs, MIIF->ldata[ctrace->region]);
      }
      for (auto& plugin : IIF[ctrace->region]->plugins()) {
        if (plugin->get_type().has_trace()) {
          fprintf(fs, "\t");
          plugin->get_type().trace(
            *plugin, ctrace->node_idx,
            fs, MIIF->ldata[ctrace->region]);
        }
      }
      fprintf(fs, "\n");
      fclose(fs);
 
      fprintf(ctrace->file->fd, "%s", trace_buf);
    }
    fflush(ctrace->file->fd);
  }
}  // dump_trace
 
void close_trace(MULTI_IF *MIIF) {
  if (!MIIF->trace_info) return;
 
  for (int iTrace = 0; MIIF->trace_info[iTrace].region >= 0; iTrace++)
    f_close(MIIF->trace_info[iTrace].file);
 
  free(MIIF->trace_info);
  MIIF->trace_info = NULL;
}
 
void MULTI_IF::initialize_currents(double idt, int subDt) {
  numSubDt = subDt;
  dt       = idt/this->numSubDt;
 
  allocate_shared_data(this);
  initializeIMPData(this);
  initialize_ionic_IF(this);
}
 
void MULTI_IF::dump_luts_MIIF(bool zipped) {
  // if one of the IFs in MIIF does not live in partion 0
  // we won't get a LUT dumped.
  if (get_rank()) return;
 
  std::vector<IonIfBase*>& pIF = this->IIF;
  for (auto& IF : pIF) {
    int ndmps = IF->dump_luts(zipped);
    if (ndmps < IF->tables().size()) {
      log_msg(logger, 4, 0, "LUT dump error %s: only %d out of %d LUTs dumped.\n",
              IF->get_type().get_name().c_str(), ndmps, IF->tables().size());
    }
    for (auto& plugin : IF->plugins()) {
      ndmps = plugin->dump_luts(zipped);
      if (ndmps < plugin->tables().size()) {
        log_msg(logger, 4, 0, "LUT dump error %s: only %d out of %d LUTs dumped.\n",
                plugin->get_type().get_name().c_str(), ndmps, IF->tables().size());
      }
    }
  }
}
 
void MULTI_IF::close_svs_dumps() {
  if (get_rank()) return;
 
  // close files
  for (int i = 0; i < svd.n; i++)
#ifndef USE_HDF5
    if (svd.hdls[i] != NULL)
#endif  // ifndef USE_HDF5
    f_close(svd.hdls[i]);
}
 
char *get_sv(void *tab, int offset, node_count_t n, int svSize, int size, int dlo_vector_size) {
  char *buf = static_cast<char *>(malloc(static_cast<size_t>(n)*size));
  char *bp  = buf;
  char *p   = static_cast<char *>(tab) + offset;
 
  for (node_index_t i = 0; i < n; i += dlo_vector_size) {
    node_count_t dlo_array_size = min(static_cast<node_count_t>(dlo_vector_size), n - i);
    memcpy(bp, p, size * dlo_array_size);
    bp += size * dlo_array_size;
    p  += svSize;
  }
  return buf;
}
 
size_t MULTI_IF::dump_svs(base_timer *iot) {
  size_t nwr = 0;
 
  if (iot->triggered) {
    for (int i = 0; i < svd.n; i++) {
      nwr = 0;
 
      FILE* fd = svd.hdls[i] ? svd.hdls[i]->fd : NULL;
      char *buf = get_sv(svd.svtab[i], svd.offset[i], svd.num[i], svd.svsize[i], svd.size[i], svd.dlo_vs[i]);
      nwr += SF::root_write<char>(fd, (char*) buf,
                                  static_cast<size_t>(svd.size[i]) * static_cast<size_t>(svd.num[i]),
                                  PETSC_COMM_WORLD);
      free(buf);
    }
    svd.nwr += static_cast<long>(nwr);
    svd.n_dumps++;
  }
 
  return nwr;
}
 
#ifdef HAS_GPU_MODEL
__global__
void update_vm(node_index_t start, node_index_t end, double *vm, double *ion, double dt)
{
  node_index_t i = blockIdx.x*blockDim.x + threadIdx.x;
  if (i < end)
     vm[i] = vm[i] + (ion[i] * (-dt));
}
#endif
 
// * compute what has to be computed, current or otherwise
void MULTI_IF::compute_ionic_current(bool flag_send, bool flag_receive)
{
  gdata[Iion]->set(0.0);
  if (flag_send == 1) {
    this->getRealData();
  }
 
  for (int j = 0; j < this->numSubDt; j++)
  {
    for (int i = 0; i < this->N_IIF; i++)
    {
      IonIfBase* pIF = this->IIF[i];
      if (!this->N_Nodes[i]) continue;
 
      prepare_error_recovery(this, i, pIF->get_tstp().cnt * this->dt);
      update_ts(&pIF->get_tstp());
 
      node_index_t current = 0;
      do {
        try {
          pIF->compute(current, this->N_Nodes[i], this->ldata[i]);
          current = this->N_Nodes[i];
        }
        catch(int e) { //FIXME
          if(e==-1) {
            // CHECK this->NodeLists[i][e.node], e.node);
            fprintf(stderr, "LIMPET compute fail in %s at node %jd (local %jd)! Aborting!\n",
                    this->name.c_str(), printable_int(this->NodeLists[i][j]), printable_int(j));
            exit(1);
          }
          current++;
        }
      } while (current < this->N_Nodes[i]);
 
      for (auto& plugin : pIF->plugins()) {
        current = 0;
        update_ts(&plugin->get_tstp());
 
        do {
          try {
            plugin->compute(current, this->N_Nodes[i], this->ldata[i]);
            current = this->N_Nodes[i];
 
            if (plugin == nullptr)
              throw -1;
          }
          catch(int e) { //FIXME
            if(e==-1) {
              // CHECK this->NodeLists[i][e.node], e.node);
              fprintf(stderr, "LIMPET plugin fail in %s at node %jd (local %jd)! Aborting!\n",
                      this->name.c_str(), printable_int(this->NodeLists[i][j]), printable_int(j));
              exit(1);
            }
            current++;
          }
        } while (current < this->N_Nodes[i]);
      }
    }
 
    #ifdef HAS_GPU_MODEL
    //TODO: extUpdateVm can be true outside of the bench executable ! But this should
    // only run in bench! This should be done for each model according to its target
    if(!extUpdateVm && is_gpu(this->targets[0])) {
 
      #if defined __CUDA__ || defined __HIP__
      update_vm<<<(this->N_Nodes[0]/64)+1,64>>>(0, this->N_Nodes[0], this->ldata[0][Vm], this->ldata[0][Iion], this->dt);
        #ifdef __CUDA__
      cudaDeviceSynchronize();
        #elif defined __HIP__
      hipDeviceSynchronize();
#endif
      #else
                fprintf(stderr, "GPU/CUDA not found");
      #endif
  }
    #endif
 
    if (flag_receive == 1) {
      this->releaseRealData();
    }
 
    // TODO: Target should be checked for each region
    if(!extUpdateVm && !is_gpu(this->targets[0]))
      gdata[Vm]->add_scaled(*gdata[Iion], SF_real(-dt));
  }
}
 
void MULTI_IF::free_MIIF() {
  assert(!this->doppel);
 
  // Free data first since this function needs to access memory in the IIFs
  freeIMPData(this);
  for (auto& pIF : this->IIF) {
    pIF->get_type().destroy_ion_if(pIF);
  }
}
 
void initializeIMPData(MULTI_IF *pMIIF) {
  for (int i = 0; i < pMIIF->N_IIF; i++) {
    if (!pMIIF->N_Nodes[i])
      continue;
 
    // if memory contiguous, we don't need to copy data, merely pass pointers
    pMIIF->contiguous[i] = pMIIF->N_Nodes[i]-1 ==
      pMIIF->NodeLists[i][pMIIF->N_Nodes[i]-1]-pMIIF->NodeLists[i][0];
 
    for (int j = 0; j < NUM_IMP_DATA_TYPES; j++) {
      // check if it is used by the IMPs
      if (!USED_DAT(pMIIF->IIF[i], imp_data_flag[j]) )
        continue;
 
      // allocate mem buffer only if not contiguous or if data is on GPU
      if (!pMIIF->contiguous[i] || is_gpu(pMIIF->IIF[i]->get_target()))
        pMIIF->ldata[i][j] =
          allocate_on_target<GlobalData_t>(pMIIF->IIF[i]->get_target(),
                                           pMIIF->N_Nodes[i]);
      // check if data supplied
      if (pMIIF->gdata[j] == NULL) {
        log_msg(pMIIF->logger, 5, LOCAL, "IMP data type %s not supplied for region %d", imp_data_names[j], i);
        exit(1);
      }
    }
  }
}  // initializeIMPData
 
void MULTI_IF::getRealData() {
  // set rdata in the parent vector to point to the
  // local data to be computed on the local processor
  for (int j = 0; j < NUM_IMP_DATA_TYPES; j++) {
    SF_real* rdata;
    if (this->gdata[j] != NULL) {
      rdata             = this->gdata[j]->ptr();
      this->procdata[j] = rdata;
    }
    else continue;
 
    // now map local data so that each ionic model can use it
    // if noncontiguous in the parent vector, copy to contiguous local vector
    for (int i = 0; i < this->N_IIF; i++) {
      if (!this->N_Nodes[i] || !USED_DAT(this->IIF[i], imp_data_flag[j]) )
        continue;
 
      if (this->contiguous[i] && !is_gpu(this->IIF[i]->get_target()))
        this->ldata[i][j] = static_cast<GlobalData_t *>(rdata) + this->NodeLists[i][0];
       else {
        const node_index_t* ip = this->NodeLists[i];
        for (node_index_t k = 0; k < this->N_Nodes[i]; k++)
          this->ldata[i][j][k] = rdata[ip[k]];
       }
    }
  }
}  // getRealData
 
void MULTI_IF::releaseRealData() {
 
  for (int i = 0; i < this->N_IIF; i++) {
    if (!this->N_Nodes[i]) continue;
 
    // Every external variable used by the model must be copied on GPU
    if (is_gpu(this->IIF[i]->get_target())) {
    for (int j = 0; j < NUM_IMP_DATA_TYPES; j++) {
        if (!USED_DAT(this->IIF[i], imp_data_flag[j]) )
          continue;
        if (this->contiguous[i])
          memcpy(&this->procdata[j][this->NodeLists[i][0]], this->ldata[i][j], sizeof(this->ldata[i][j][0])*this->N_Nodes[i]);
        else {
          for (node_index_t k = 0; k < this->N_Nodes[i]; k++)
            this->procdata[j][this->NodeLists[i][k]] = this->ldata[i][j][k];
        }
      }
    }
    else {
      for (int j = 0; j < NUM_IMP_DATA_TYPES; j++) {
        if (this->IIF[i]->get_moddat() & imp_data_flag[j]) {
          if (!this->contiguous[i]) {
            for (node_index_t k = 0; k < this->N_Nodes[i]; k++)
              this->procdata[j][this->NodeLists[i][k]] = this->ldata[i][j][k];
          }
        }
      }
    }
  }
 
  for (int j = 0; j < NUM_IMP_DATA_TYPES; j++)
    if (this->gdata[j] != NULL)
      this->gdata[j]->release_ptr(this->procdata[j]);
}
 
// * free the local storage vectors
void freeIMPData(MULTI_IF *pMIIF) {
  for (int j = 0; j < NUM_IMP_DATA_TYPES; j++) {
    if (pMIIF->gdata[j] != NULL) {
      for (int i = 0; i < pMIIF->N_IIF; i++)
        if ((!pMIIF->contiguous[i] || is_gpu(pMIIF->IIF[i]->get_target())) && USED_DAT(pMIIF->IIF[i], imp_data_flag[j]) && pMIIF->N_Nodes[i] > 0) {
          deallocate_on_target<GlobalData_t>(pMIIF->IIF[i]->get_target(),
                                             pMIIF->ldata[i][j]);
        }
      delete pMIIF->gdata[j];
    }
  }
  free(pMIIF->contiguous);
  deallocate_on_target(Target::CPU, pMIIF->ldata);
}
 
int get_plug_flag(char *plgstr, int *out_num_plugins, IonTypeList& out_plugins) {
  if (plgstr == NULL) {
    *out_num_plugins = 0;
    out_plugins.clear();
    return 1;
  }
 
  // maximum number of plugins is 200 apparently
  char *plugspec = dupstr(plgstr);
  char *saveptr;
  char *token = tokstr_r(plugspec, ":", &saveptr);
 
  *out_num_plugins = 0;
  while (token) {
    IonType* type = get_ion_type(std::string(token));
    if (type == NULL) {
      free(plugspec);
      return 0;
    } else {
      out_plugins.push_back(*type);
    }
    token = tokstr_r(NULL, ":", &saveptr);
  }
 
  *out_num_plugins = out_plugins.size();
 
  free(plugspec);
  return 1;
}  // get_plug_flag
 
int MULTI_IF::adjust_MIIF_variables(const char* variable,
                                    const SF::vector<SF_int> & indices,
                                    const SF::vector<SF_real> & values)
{
  int num_changed = 0;
 
  // determine if we're dealing with an external variable or a state variable.
  if (index(variable, '.') == NULL) {
    // external variable, the easy case.
    // which external variable?
    int data_id = -1;
    for (int ii = 0; ii < NUM_IMP_DATA_TYPES; ii++) {
      if (strcmp(variable, imp_data_names[ii]) == 0) {
        data_id = ii;
        break;
      }
    }
    if (data_id == -1) {
      log_msg(logger, 5, FLUSH, "Error!  No external variable named %s in this build of openCARP.", variable);
      exit(EXIT_FAILURE);
    }
 
    // make sure that external variable is being used.
    if (this->gdata[data_id] == NULL) {
      log_msg(logger, 4, 0,
              "External variable %s is not being used this processor.\n"
              "Is this really what you meant to do?\n"
              "Perhaps you don't have the correct Ionic models selected.",
              imp_data_names[data_id]);
    }
    else {
      // fill in the external variable with the data from the file.
      SF_real *raw_data = this->gdata[data_id]->ptr();
 
      for (size_t i = 0; i < indices.size(); i++) {
        raw_data[indices[i]] = values[i];
        num_changed++;
      }
 
      this->gdata[data_id]->release_ptr(raw_data);
    }
  }
  else {
    // state variable.
    // extract the IMP name and state variable name from the string.
    char *saveptr;
    char *my_variable = dupstr(variable);
    char *IIF_name    = tokstr_r(my_variable, ".", &saveptr);
    char *sv_name     = tokstr_r(NULL, ".", &saveptr);
 
    IonType* type = get_ion_type(std::string(IIF_name));
    if (type == NULL) {
      log_msg(logger, 5, 0, "%s error: %s is not a valid IMP name.", __func__, IIF_name);
      exit(EXIT_FAILURE);
    }
 
    int sv_offset;
    int sv_size;
    SVgetfcn sv_get = type->get_sv_offset(sv_name,  &sv_offset, &sv_size);
    if (sv_get == NULL) {
      log_msg(logger, 5, 0, "%s error: %s is not a valid state variable for the %s model.",
              __func__, sv_name, IIF_name);
      exit(EXIT_FAILURE);
    }
    SVputfcn sv_put = getPutSV(sv_get);
 
    // Ok, go through the ionic models
    for (int i_iif = 0; i_iif < this->N_IIF; i_iif++) {
      IonIfBase *lIIF = NULL;
      if (this->IIF[i_iif]->get_type() == *type) {
        lIIF = this->IIF[i_iif];
      }
      else {
        // Go through the plugins
        for (auto& plugin : this->IIF[i_iif]->plugins()) {
          if (plugin->get_type() == *type) {
            lIIF = plugin;
            break;
          }
        }
      }
 
      if (lIIF == NULL) {
        continue;
      }
 
      for (size_t ii = 0; ii < indices.size(); ii++) {
        GlobalData_t file_value = values[ii];
 
        node_index_t local_node = static_cast<node_index_t>(indices[ii]);
        node_index_t *target = static_cast<node_index_t*>(bsearch(&local_node, this->NodeLists[i_iif],
               this->N_Nodes[i_iif], sizeof(node_index_t), node_index_cmp));
        if (target) {
          // We found a point to adjust! Do the adjustment.
          num_changed++;
          sv_put(*lIIF, target - this->NodeLists[i_iif], sv_offset, file_value);
        }
      }
 
    }
    free(my_variable);
  }
 
  MPI_Allreduce(MPI_IN_PLACE, &num_changed, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
 
  return num_changed;
}  // adjust_MIIF_variables
 
int determine_write_ranges(int N, size_t *offset, size_t bufsize, int **ranges) {
  int nitems;
 
  if (!get_rank() ) {
    Salt_list r; STRUCT_ZERO(r); r.chunk = 10;
    int temp    = 0;
    SLIST_APPEND(&r, temp);
    long loff = offset[0];
    for (int i = 0; i < N; i++)
      if (offset[i]-loff > bufsize) {
        SLIST_APPEND(&r, i);
        loff = offset[i];
      }
 
    nitems = r.nitems;
    SLIST_APPEND(&r, N);
    *ranges = (int *)r.data;
  }
 
  MPI_Bcast(&nitems, 1,        MPI_INT, 0, PETSC_COMM_WORLD);
  if (get_rank() )
    *ranges = static_cast<int *>(malloc( (nitems+1)*sizeof(int) ));
  MPI_Bcast(*ranges, nitems+1, MPI_INT, 0, PETSC_COMM_WORLD);
  return nitems;
}  // determine_write_ranges
 
void MULTI_IF::dump_state(char *fname, float simtime, mesh_t gid,
                          bool append, unsigned int revision)
{
  float t0, t1;
  t0 = get_time();
 
  FILE_SPEC out = NULL;
  int rank = get_rank();
  global_node_index_t miif_node_gsize = get_global(static_cast<global_node_index_t>(this->numNode), MPI_SUM);
  int error = 0;
 
  // EMI ionic state lives on mesh elements (membrane faces), electrics on nodes.
  // Pick the matching canonical permutation and numbering from the data layout.
  const bool elemwise   = this->gdata[Vm] && this->gdata[Vm]->layout == sf_vec::elemwise;
  const int  canon_perm = elemwise ? ELEM_PETSC_TO_CANONICAL : PETSC_TO_CANONICAL;
 
  log_msg(logger, 0, 0, "Saving state at time %f in file: %s", simtime, fname);
 
  if (rank == 0) {
    out = f_open(fname, append ? "a" : "w");
 
    if(out) {
      fseek(out->fd, 0, SEEK_END);  // NOP call to sync disk data for switch to write mode
      write_bin_string(out, Magic_MIIF_ID);
 
      fwrite(&MIIF_Format, sizeof(unsigned int), 1, out->fd);
      fwrite(&revision, sizeof(unsigned int), 1, out->fd);
      time_t tm = time(NULL);
      fwrite(&tm, sizeof(time_t), 1, out->fd);
      fwrite(&simtime, sizeof(float), 1, out->fd);
      fwrite(&miif_node_gsize, sizeof(global_node_index_t), 1, out->fd);
 
      int num_gdata = 0;
      for (int i = 0; i < NUM_IMP_DATA_TYPES; i++)
        if (this->gdata[i]) num_gdata++;
 
      fwrite(&num_gdata, sizeof(int), 1, out->fd);
    }
    else
      error++;
  }
 
// #define CHATTY
 
  if(get_global(error, MPI_SUM))
    EXIT(EXIT_FAILURE);
 
  FILE* fd = rank == 0 ? out->fd : NULL;
  sf_vec* outVec;
 
  for (int i = 0; i < NUM_IMP_DATA_TYPES; i++) {
    if (this->gdata[i]) {
      if(rank == 0) {
        write_bin_string(out, imp_data_names[i]);
#ifdef CHATTY
        log_msg(NULL, 0, 0, "\tDumping %s at %d", imp_data_names[i], ftell(fd) );
#endif
      }
 
      // If we have a proper set up parallel context, we pass a gid != unset_msh. Then
      // we can query petsc-to-canonical scattering. If we are in a simpler context
      // like bench, we assume that no scattering is needed and just shallow-copy the vector.
      if(gid != unset_msh) {
        SF::init_vector(&outVec, this->gdata[Vm]);
        SF::scattering* sc = get_permutation(gid, canon_perm, 1);
        assert(sc != NULL);
 
        sc->forward(*this->gdata[i], *outVec);
        outVec->write_binary<SF_real>(fd);
        delete outVec;
      }
      else {
        this->gdata[i]->write_binary<SF_real>(fd);
      }
    }
  }
 
  // output IMP region info
  // #regions and the IM and plugins for each region
  // determine memory requirements for each region
  int*    imp_mem = new int[this->N_IIF];
  size_t *offset = NULL;
  long filepos   = 0.;
 
  if (rank == 0) {
#ifdef CHATTY
    log_msg(NULL, 0, 0, "\tDumping IMP sizes at %d", ftell(fd) );
#endif
 
    int max = 1;
    fwrite(&this->N_IIF, sizeof(int), 1, fd);  // #IIF's
 
    for (int i = 0; i < this->N_IIF; i++) {
      IonIfBase *imp = this->IIF[i];
      write_bin_string(out, imp->get_type().get_name().c_str());
 
      // record the de-interleaved per-node size (format 3), not the whole SoA block
      int sv_size = (int) imp->get_sv_per_node_size();
      fwrite(&sv_size, sizeof(int), 1, fd);
      // fingerprint of the SV layout, so a restart against a changed model revision
      // is detected instead of silently misread field-for-field
      uint64_t sv_fp = imp->sv_fingerprint();
      fwrite(&sv_fp, sizeof(uint64_t), 1, fd);
      unsigned long n_plugins = (unsigned long) imp->plugins().size();
      fwrite(&n_plugins, sizeof(int), 1, fd);
      imp_mem[i] = sv_size;
 
      for (auto& plug : imp->plugins()) {
        write_bin_string(out, plug->get_type().get_name().c_str());
 
        int plug_size = (int) plug->get_sv_per_node_size();
        fwrite(&plug_size, sizeof(int), 1, fd);
        uint64_t plug_fp = plug->sv_fingerprint();
        fwrite(&plug_fp, sizeof(uint64_t), 1, fd);
        imp_mem[i] += plug_size;
      }
    }
    filepos = ftell(fd);
  }
 
  MPI_Bcast(imp_mem, this->N_IIF, MPI_INT, 0, PETSC_COMM_WORLD);
 
  // If we have a proper set up parallel context, we pass a gid != unset_msh. Then
  // we can query a mesh and its parallel numbering. If we are in a simpler context
  // like bench, we build the numbering ourselves. Then, we assume that no
  // renumbering has taken place.
  SF::vector<global_node_index_t> loc2canon;
  if(gid != unset_msh) {
    const sf_mesh & mesh = get_mesh(gid);
    if(elemwise) {
      const SF::vector<mesh_int_t> & canon_nbr = mesh.get_numbering(SF::NBR_ELEM_SUBMESH);
      loc2canon.resize(canon_nbr.size());
 
      for(size_t i=0; i<canon_nbr.size(); i++) loc2canon[i] = canon_nbr[i];
    }
    else {
      const SF::vector<mesh_int_t> & canon_nbr = mesh.get_numbering(SF::NBR_SUBMESH);
      const SF::vector<mesh_int_t> & alg_idx   = mesh.pl.algebraic_nodes();
      loc2canon.resize(alg_idx.size());
 
      for(size_t i=0; i<alg_idx.size(); i++) loc2canon[i] = canon_nbr[alg_idx[i]];
    }
  }
  else {
    int rank = get_rank();
    SF::vector<global_node_index_t> layout;
    SF::layout_from_count<global_node_index_t>(static_cast<global_node_index_t>(this->numNode),
                                               layout, PETSC_COMM_WORLD);
    loc2canon.resize(this->numNode);
 
    for (node_index_t i = 0; i < this->numNode; i++) loc2canon[i] = layout[rank] + i;
  }
 
  // loc2canon must hold one canonical index per local node; a shorter array
  // would read out of bounds below and write a corrupt checkpoint
  assert(loc2canon.size() == static_cast<size_t>(this->numNode));
 
  // write out IIF_Mask canonically ordered
  global_node_index_t *canord = new global_node_index_t[this->numNode];
  for (node_index_t i = 0; i < this->numNode; i++)
    canord[i] = loc2canon[i];
 
#ifdef CHATTY
  log_msg(NULL, 0, 0, "\tDumping IMP masks at %d", filepos);
#endif  // ifdef CHATTY
 
  SF::root_write_ordered<global_node_index_t, IIF_Mask_t>(fd, canord, IIFmask, numNode, PETSC_COMM_WORLD);
 
  // for every local node, we compute the size of its ionic model
  // and plugins, and the associated displacments between the nodes
  SF::vector<global_node_index_t> impsize(numNode), impdsp;
  for(node_index_t i=0; i<numNode; i++)
    impsize[i] = imp_mem[(int)IIFmask[i]];
  SF::dsp_from_cnt(impsize, impdsp);
 
  // this also gives us the size of the total buffer
  size_t num_entr = static_cast<size_t>(SF::sum(impsize));
  SF::vector<char> impdata(num_entr);
 
  // De-interleave each node's state into a contiguous per-node record (format 3).
  // The state lives as data-layout-optimized SoA blocks: node imp_nod_idx is lane
  // (imp_nod_idx % vec_size) of block (imp_nod_idx / vec_size). We gather that lane
  // of every state-variable field, keyed by the node's own canonical id (via loc).
  // This makes the dump independent of the partitioning. Without DLO vec_size is 1
  // and this reduces to copying each node's contiguous state.
  for (int imp_idx = 0; imp_idx < this->N_IIF; imp_idx++)
  {
    IonIfBase* iif = this->IIF[imp_idx];
 
    // field layouts are model properties: gather them once, not per node
    std::vector<std::pair<int, int>> main_fields;
    iif->get_sv_layout(main_fields);
 
    std::vector<std::vector<std::pair<int, int>>> plug_fields;
    for (auto& plug : iif->plugins()) {
      plug_fields.emplace_back();
      plug->get_sv_layout(plug_fields.back());
    }
 
    for (node_index_t imp_nod_idx = 0; imp_nod_idx < this->N_Nodes[imp_idx]; imp_nod_idx++) {
      node_index_t loc = this->NodeLists[imp_idx][imp_nod_idx];
 
      char* write = impdata.data() + impdsp[loc];
      auto gather = [&](char* src, std::size_t sz) { memcpy(write, src, sz); write += sz; };
 
      iif->for_each_sv_field(imp_nod_idx, main_fields, gather);
      for (size_t p = 0; p < plug_fields.size(); p++)
        iif->plugins()[p]->for_each_sv_field(imp_nod_idx, plug_fields[p], gather);
    }
  }
 
  // finally we can write out the IonIf data to disk
  SF::root_write_ordered<global_node_index_t, char>(fd, canord, impsize.data(), impdata.data(),
                                    numNode, num_entr, PETSC_COMM_WORLD);
  delete [] canord;
  delete [] imp_mem;
 
  if(fd) fclose(fd);
 
  double dump_time = timing(t1, t0);
  log_msg(logger, 0, 0, " in %.3f seconds.\n", dump_time);
}
 
float MULTI_IF::restore_state(const char *fname, mesh_t gid, bool close)
{
  static FILE_SPEC in = NULL;
  static char *last_fn = NULL;
  double t0, t1;
 
  t0 = get_time();
  int error    = 0;
  int my_rank  = get_rank();
  int mpi_size = get_size();
 
  // EMI ionic state lives on mesh elements (membrane faces), electrics on nodes.
  // Pick the matching canonical permutation and numbering from the data layout.
  const bool elemwise   = this->gdata[Vm] && this->gdata[Vm]->layout == sf_vec::elemwise;
  const int  canon_perm = elemwise ? ELEM_PETSC_TO_CANONICAL : PETSC_TO_CANONICAL;
 
  if (my_rank == 0) {
    if (fname) in = f_open(fname, "r");
 
    if (!in) {
      if (fname) log_msg(logger, 5, 0, "Error: cannot open file: %s\n", fname);
      else       log_msg(logger, 5, 0, "Error: file stream not open yet");
      error++;
    }
    else {
      if (strcmp(read_bin_string(in), Magic_MIIF_ID) ) {
        log_msg(logger, 5, 0, "%s is not a recognized MIIF dump file", fname);
        error++;
      }
    }
  }
 
  if (get_global(error, MPI_SUM))
    EXIT(1);
 
  if (fname) {
    free(last_fn);
    last_fn = strdup(fname);
  }
  else {
    fname = last_fn;
  }
 
  unsigned int format, version;
  float time = 0.0f;
  time_t save_date;
  f_read_par(&format,    sizeof(unsigned int), 1, in);
  f_read_par(&version,   sizeof(unsigned int), 1, in);
  f_read_par(&save_date, sizeof(time_t),       1, in);
  f_read_par(&time,      sizeof(float),        1, in);
  log_msg(logger, 0, 0, "Restoring time %f from %s (format v%d) generated\n\tby calling "
                          "program r%d on %s", time, fname, format, version, ctime(&save_date) );
 
  global_node_index_t savedNum = 0;
  global_node_index_t glob_numNode = get_global(static_cast<global_node_index_t>(this->numNode),
                                                MPI_SUM, PETSC_COMM_WORLD);
 
  if (format >= 2) {
    f_read_par(&savedNum, sizeof(global_node_index_t), 1, in);
  } else {
    int savedNum32 = 0;
    f_read_par(&savedNum32, sizeof(int), 1, in);
    savedNum = savedNum32;
  }
  if (savedNum != glob_numNode) {
    log_msg(logger, 5, 0, "expecting %jd nodes but read %jd nodes",
            printable_int(glob_numNode), printable_int(savedNum));
    EXIT(1);
  }
 
  int num_gdata;
  f_read_par(&num_gdata, sizeof(int), 1, in);
  SF::scattering & petsc2canon = *get_permutation(gid, canon_perm, 1);
 
  for (int g = 0; g < num_gdata; g++) {
    char *datatype = read_bin_string_par(in);
    int   i;
    for (i = 0; i < NUM_IMP_DATA_TYPES; i++) {
      if (!strcmp(datatype, imp_data_names[i]) ) {
        if (this->gdata[i]) {
          log_msg(logger, 0, 0, "\tRestoring global data %s", imp_data_names[i], datatype);
 
          FILE* fd = my_rank == 0 ? in->fd : NULL;
          this->gdata[i]->read_binary<SF_real>(fd);
          bool fwd = false;
          petsc2canon(*this->gdata[i], fwd);
        } else {
          log_msg(logger, 0, 0, "\tGlobal data %s not used", datatype);
        }
        break;
      }
    }
    if (i == NUM_IMP_DATA_TYPES)
      log_msg(logger, 3, 0, "\tSaved global data %s not recognized", datatype);
 
    // skip the saved block for a global-data type we cannot restore. The skip
    // length must come from the saved size (one SF_real per global node), not
    // from this->gdata[i], which is exactly what is NULL / out of range here.
    if ( (my_rank == 0) && ((i == NUM_IMP_DATA_TYPES) || !this->gdata[i]) )
      fseek(in->fd, static_cast<long>(glob_numNode) * sizeof(SF_real), SEEK_CUR);
 
    free(datatype);
  }
 
  int N_IIF;
  f_read_par(&N_IIF, sizeof(int), 1, in);
 
  // read in the IMPs for each region
  IMPinfo *imps = static_cast<IMPinfo *>(calloc(N_IIF, sizeof(IMPinfo)));
  int *IMPsz    = static_cast<int *>(malloc(N_IIF * sizeof(int)));
 
  for (int i = 0; i < N_IIF; i++) {
    imps[i].name = read_bin_string_par(in);
    f_read_par(&imps[i].sz,    sizeof(int), 1, in);
    if (format >= 3) f_read_par(&imps[i].fingerprint, sizeof(uint64_t), 1, in);
    f_read_par(&imps[i].nplug, sizeof(int), 1, in);
 
    imps[i].offset = 0;
    IMPsz[i]       = imps[i].sz;
    imps[i].plug   = static_cast<IMPinfo *>(calloc(sizeof(IMPinfo), imps[i].nplug));
 
    for (int j = 0; j < imps[i].nplug; j++) {
      imps[i].plug[j].offset = IMPsz[i];
      imps[i].plug[j].name   = read_bin_string_par(in);
      f_read_par(&imps[i].plug[j].sz,  sizeof(int), 1, in);
      if (format >= 3) f_read_par(&imps[i].plug[j].fingerprint, sizeof(uint64_t), 1, in);
      IMPsz[i] += imps[i].plug[j].sz;
    }
  }
 
  // compare new and old regions to determine which can be copied
  for (int i = 0; i < N_IIF; i++) {
    if (i >= this->N_IIF) {
      log_msg(NULL, 3, 0, "Saved IMP region %d out of range", i);
      continue;
    }
    if (strcmp(imps[i].name, this->IIF[i]->get_type().get_name().c_str()) ) {
      log_msg(NULL, 3, 0, "Saved IMP region %d ionic model does not match that of IMP region %d", i, i);
      continue;
    }
    // format >= 3 stores the de-interleaved per-node size; format 2 the whole SoA block
    const bool per_node_fmt = format >= 3;
    const int  exp_sz = per_node_fmt ? (int) this->IIF[i]->get_sv_per_node_size()
                                     : (int) this->IIF[i]->get_sv_size();
    if (imps[i].sz != exp_sz) {
      log_msg(NULL, 3, 0, "Saved IMP region %d size does not match current IMP size", i);
      continue;
    }
    // a matching size but different SV-layout fingerprint means the model's state
    // variables were reordered/retyped since the checkpoint: skip rather than misread
    if (per_node_fmt && imps[i].fingerprint != this->IIF[i]->sv_fingerprint()) {
      log_msg(NULL, 3, 0, "Saved IMP region %d state-variable layout does not match current model", i);
      continue;
    }
 
    // the IM is good to go
    imps[i].compatible = true;
    for (int j = 0; j < imps[i].nplug; j++)
      for (int k = 0; k < this->IIF[i]->plugins().size(); k++) {
        const int exp_plug_sz = per_node_fmt ? (int) this->IIF[i]->plugins()[k]->get_sv_per_node_size()
                                             : (int) this->IIF[i]->plugins()[k]->get_sv_size();
        const bool plug_fp_ok = !per_node_fmt ||
            imps[i].plug[j].fingerprint == this->IIF[i]->plugins()[k]->sv_fingerprint();
        if (!strcmp(imps[i].plug[j].name, this->IIF[i]->plugins()[k]->get_type().get_name().c_str()) &&
            (imps[i].plug[j].sz == exp_plug_sz) && plug_fp_ok) {
          log_msg(NULL, 3, 0, "Saved IMP region %d plugin %s compatible", i, imps[i].plug[j].name);
          imps[i].plug[j].map        = k;
          imps[i].plug[j].compatible = true;
          break;
        }
      }
  }
 
  // read in region mask
  const size_t global_node_count = static_cast<size_t>(glob_numNode);
  IIF_Mask_t *canMask = static_cast<char *>(malloc(global_node_count*sizeof(this->IIFmask[0]) ));
  f_read_par(canMask, sizeof(IIF_Mask_t), global_node_count, in);
 
  // determine relative IMP data offsets based on canonical ordering
  size_t *offset = static_cast<size_t *>(calloc(global_node_count+1, sizeof(size_t)));
  for (size_t i = 1; i <= global_node_count; i++)
    offset[i] = offset[i-1] + IMPsz[static_cast<int>(canMask[i-1])];
 
  // read in the state variable data
  // TODO: Aurel: Ideally rank0 would read in the data and communicate it, but right now I dont have
  // time to code this up. Thus, at least the ranks will read their data one by one so
  // we dont saturate the file system.
  long SVstart;
  if (my_rank == 0) {
    SVstart = ftell(in->fd);
    f_close(in);
  }
 
  MPI_Bcast(&SVstart, sizeof(long), MPI_BYTE, 0, PETSC_COMM_WORLD);
  global_node_index_t mismatch = 0;
 
  // we need to map from a local algebraic index to a global canonical index
  const sf_mesh & mesh = get_mesh(gid);
  SF::vector<global_node_index_t> loc2canon;
  if(elemwise) {
    const SF::vector<mesh_int_t> & canon_nbr = mesh.get_numbering(SF::NBR_ELEM_SUBMESH);
    loc2canon.resize(canon_nbr.size());
 
    for(size_t i=0; i<canon_nbr.size(); i++) loc2canon[i] = canon_nbr[i];
  }
  else {
    const SF::vector<mesh_int_t> & canon_nbr = mesh.get_numbering(SF::NBR_SUBMESH);
    const SF::vector<mesh_int_t> & alg_idx   = mesh.pl.algebraic_nodes();
    loc2canon.resize(alg_idx.size());
 
    for(size_t i=0; i<alg_idx.size(); i++) loc2canon[i] = canon_nbr[alg_idx[i]];
  }
 
  // loc2canon must hold one canonical index per local node; a shorter array
  // would be indexed out of bounds by the per-node restore below
  assert(loc2canon.size() == static_cast<size_t>(this->numNode));
 
  // the ranks read the file one-by-one
  for (int pid = 0; pid < mpi_size; pid++) {
    if (my_rank == pid) {
      if(in) delete in;
      in = f_open(fname, "r");
 
      for (int i = 0; i < N_IIF; i++) {
        fseek(in->fd, SVstart, SEEK_SET);
        if (format >= 3)
          mismatch += this->IIF[i]->restore_per_node(in, this->N_Nodes[i], this->NodeLists[i],
                                  canMask, offset, imps+i, loc2canon.data());
        else
          mismatch += this->IIF[i]->restore(in, this->N_Nodes[i], this->NodeLists[i],
                                  canMask, offset, imps+i, loc2canon.data());
      }
 
      f_close(in);
    }
    MPI_Barrier(PETSC_COMM_WORLD);
  }
  free(canMask);
  free(IMPsz);
 
  MPI_Reduce(my_rank ? &mismatch : MPI_IN_PLACE, &mismatch, 1,
             opencarp::mpi_datatype<global_node_index_t>(), MPI_SUM, 0, PETSC_COMM_WORLD);
  if ( (my_rank == 0) && mismatch)
    log_msg(NULL, 3, 0, "Number of nonmatching nodes: %jd", printable_int(mismatch));
 
  if ( (my_rank == 0) && (!close) ) {
    if(in) delete in;
    in = f_open(fname, "r");
 
    // position at end of current MIIF data
    fseek(in->fd, SVstart+offset[global_node_count], SEEK_SET);
  }
 
  free(offset);
  for (int i = 0; i < N_IIF; i++)
    free(imps[i].plug);
  free(imps);
 
  double restore_time = timing(t1, t0);
  log_msg(logger, 0, 0, "State restored from file %s in %.3f seconds.\n", fname, restore_time);
 
  return time;
}
 
void MULTI_IF::sv_dump_add_by_name_list( int region, char *imp_name,
                                        char *reg_name, char *sv_lst, char *plg_lst,
                                        char *plg_sv_lst, double t, double dump_dt) {
  char file[8000], svs[1024], plgs[1024], plgsvs[1024];
  char *e, *l, *p, *svnames, *plgnames, *plgsvnames;
 
  strcpy(svs,    sv_lst);
  strcpy(plgs,   plg_lst);
  strcpy(plgsvs, plg_sv_lst);
  svnames    = &svs[0];
  plgnames   = &plgs[0];
  plgsvnames = &plgsvs[0];
 
  // override default setting for dump interval
  this->svd.intv = dump_dt;
 
  // override default start time for dumping (in case of a restart)
  this->svd.t_dump = t;
 
  // parse ionic model sv list
  while ( (e = get_next_list(svnames, ',') ) ) {
    snprintf(file, sizeof file, "%s.%s.bin", reg_name, svnames);
    this->sv_dump_add_by_name(region, imp_name, svnames, reg_name, file);
    svnames = e;
  }
 
  // Information is provided as follows:
  // region[X].plugin        = "PLG_A:PLG_B"
  // region[X].plug_sv_dumps = "m,h,n:x,y"
 
  // parse plugin list
  while ( (p = get_next_list(plgnames, ':') ) ) {
    // get sv list corresponding to current plugin
    l = get_next_list(plgsvnames, ':');
    while ( (e = get_next_list(plgsvnames, ',') ) ) {
      snprintf(file, sizeof file, "%s_%s.%s.bin", reg_name, plgnames, plgsvnames);
      this->sv_dump_add_by_name(region, plgnames, plgsvnames, reg_name, file);
      plgsvnames = e;
    }
    plgnames   = p;
    plgsvnames = l;
  }
}  // sv_dump_add_by_name_list
 
void MULTI_IF::sv_dump_add(int region, const IonType& type, int offset, int size, int dtype,
                           const char *filename, const char *regname) {
  IonIfBase *IF = this->IIF[region];
  int n      = this->svd.n;
 
  // find the proper IMP
  if (IF->get_type() != type) {
    int i;
    for (i = 0; i < IF->plugins().size(); i++)
      if (IF->plugins()[i]->get_type() == type)
        break;
    if (i == IF->plugins().size()) {
      log_msg(logger, 2, 0, "Warning: IMP %s not found in Region %s\n",
              type.get_name().c_str(), regname);
      return;
    } else {
      IF = IF->plugins()[i];
    }
  }
 
  this->svd.active = 1;
  this->svd.n++;
  this->svd.hdls    = static_cast<FILE_SPEC *>(realloc(this->svd.hdls, this->svd.n*sizeof(this->svd.hdls[0])));
  this->svd.fn      = static_cast<char **>(realloc(this->svd.fn, this->svd.n*sizeof(char *)));
  this->svd.reg     = static_cast<int *>(realloc(this->svd.reg,    this->svd.n*sizeof(int)));
  this->svd.svnames = static_cast<char **>(realloc(this->svd.svnames, this->svd.n*sizeof(char *)));
  this->svd.n_dumps = 0;
  this->svd.nwr     = 0;
  this->svd.offset  = static_cast<int *>(realloc(this->svd.offset, this->svd.n*sizeof(int)));
  this->svd.size    = static_cast<int *>(realloc(this->svd.size,   this->svd.n*sizeof(int)));
  this->svd.dlo_vs  = static_cast<int *>(realloc(this->svd.dlo_vs,   this->svd.n*sizeof(int)));
  this->svd.dtype   = static_cast<int *>(realloc(this->svd.dtype,  this->svd.n*sizeof(int)));
  this->svd.svtab   = static_cast<void **>(realloc(this->svd.svtab, this->svd.n*sizeof(void *)));
  this->svd.svsize  = static_cast<size_t *>(realloc(this->svd.svsize, this->svd.n*sizeof(size_t)));
  this->svd.num     = static_cast<node_count_t *>(realloc(this->svd.num, this->svd.n*sizeof(node_count_t)));
  this->svd.reg     = static_cast<int *>(realloc(this->svd.reg,    this->svd.n*sizeof(int)));
 
  if (!get_rank() ) {
#ifdef HDF5
    assert(0);
#else  // ifdef HDF5
    this->svd.hdls[n] = f_open(filename, "w+");
#endif  // ifdef HDF5
  } else {
    this->svd.hdls[n] = NULL;
  }
 
  this->svd.fn[n]      = dupstr(filename);
  this->svd.reg[n]     = region;
  this->svd.svnames[n] = NULL;  // has to be filled in in calling routine
  this->svd.offset[n]  = offset;
  this->svd.size[n]    = size;
  this->svd.dtype[n]   = dtype;
  this->svd.svtab[n]   = IF->get_sv_address();
  this->svd.svsize[n]  = IF->get_sv_size();
  this->svd.num[n]     = IF->get_num_node();
  this->svd.dlo_vs[n] = IF->get_type().dlo_vector_size();
}  // sv_dump_add
 
int MULTI_IF::sv_dump_add_by_name(int region, char *impname,
                                  char *svname, char *regname, char *filename) {
  int offset, size, dtype, added = 0;
  char *svtypename = NULL;
  char *filename_bin;
 
  IonType* type = get_ion_type(std::string(impname));
  filename_bin = strcat(filename, ".bin");
 
  assert(type != NULL);
  if (type->get_sv_offset(svname, &offset, &size) ) {
    type->get_sv_type(svname, &dtype, &svtypename);
    this->sv_dump_add(region, *type, offset, size, dtype, filename_bin, regname);
    this->svd.svnames[this->svd.n-1] = dupstr(svname);
    added                            = 1;
  } else {
    log_msg(NULL, 1, 0, "No state variable added to dump list");
  }
 
  return added;
}
 
void allocate_shared_data(MULTI_IF *miif)
{
  assert(miif->gdata[Vm] != NULL);
 
  for (int i = 0; i < NUM_IMP_DATA_TYPES; i++) {
    for (int n = 0; n < miif->N_IIF; n++)
      if (USED_DAT(miif->IIF[n], imp_data_flag[i]) && (miif->gdata[i] == NULL)) {
        SF::init_vector(&miif->gdata[i], miif->gdata[Vm]);
        break;
      }
  }
}
 
bool MULTI_IF::use_stretch() {
  int num_mech_data = 0;
 
  for (int i = 0; i < NUM_IMP_DATA_TYPES; i++) {
    if ( (Lambda_DATA_FLAG == i) || (delLambda_DATA_FLAG == i) ||
         (Tension_DATA_FLAG == i) || (tension_component_DATA_FLAG == i) ) {
      for (int n = 0; n < this->N_IIF; n++)
        if (USED_DAT(this->IIF[n], imp_data_flag[i]) )
          num_mech_data++;
    }
  }
  return static_cast<bool>(num_mech_data);
}
 
void MULTI_IF::transmem_stim_species(float charge, const char *species,
                                     float beta, int *node, int numnode) {
  float z;
  if (strstr(species, "Ca") != NULL) {
    charge *= 1000;                               // Ca uses micromolar instead of millimolar
    z     = 2;
  } else if (strstr(species, "Na") != NULL) {
    z = 1;
  } else if (strstr(species, "Cl") != NULL) {
    z = -1;
  } else if (strstr(species, "K") != NULL) {
    z = 1;
  } else {
    log_msg(logger, 5, 0, "Unimplemented ion species: %s\n", species);
    exit(1);
  }
  int offset, sz;
  SVgetfcn ion_get;
  SVputfcn ion_put;
  static int warned = 0;
 
  int rank = get_rank();
  SF::vector<int> layout;
  layout_from_count(numnode, layout, PETSC_COMM_WORLD);
  int my_low_idx = layout[rank], my_high_idx = layout[rank+1];
 
  for (int i = 0; i < numnode; i++) {
    if ((node[i] < my_low_idx) || (node[i] >= my_high_idx)) continue; // not on processor
 
    for (int j = 0; j < this->N_IIF; j++) {
      if (this->NodeLists[j] == NULL) continue; // no stim nodes on processor
 
      if ((node[i] < this->NodeLists[j][0]) ||
          (node[i] > this->NodeLists[j][this->N_Nodes[j]-1])) {
        continue;  // not in this IMP
      } else {
        ion_get = this->IIF[j]->get_type().get_sv_offset(species, &offset, &sz);
        if (ion_get == NULL) {
          if (!warned) {
            warned = 1;
            log_msg(logger, 2, 0, "Ion species not present in ionic model: %s\n",
                    species);
          }
          return;
        }
        ion_put = getPutSV(ion_get);
      }
 
      // does this model provide a correct conversion factor?
      double delta_conc;
      cell_geom g = this->IIF[j]->cgeom();
      if (g.sl_i2c != NDEF) {
        delta_conc = charge*g.sl_i2c/z;
      } else {
        // conversion with generic factor
        delta_conc = 10*charge*beta/(FARADAY*z);  // convert to millimolar/L
      }
 
      for (node_index_t k = 0; k < this->N_Nodes[j]; k++)
        if (node[i] == this->NodeLists[j][k])
          ion_put(*this->IIF[j], k, offset, ion_get(*this->IIF[j], k, offset)-delta_conc);
    }
  }
}  // transmem_stim_species
 
void MULTI_IF::MIIF_change_dt( double Dt) {
  this->dt = Dt;
  for (int i = 0; i < this->N_IIF; i++) {
    this->IIF[i]->destroy_luts();
    this->IIF[i]->set_dt((float) Dt);
    this->iontypes[i].get().construct_tables(*this->IIF[i]);
 
    for (int j = 0; j < this->IIF[i]->plugins().size(); j++) {
      auto& plugin = this->IIF[i]->plugins()[j];
      plugin->destroy_luts();
      plugin->set_dt((float) Dt);
      this->plugtypes[i][j].get().construct_tables(*plugin);
    }
  }
}
 
#define MEMFREE(A) free(A)
 
void free_doppel(MULTI_IF *m) {
  assert(m->doppel);
 
  for (int i = 0; i < m->N_IIF; i++) {
    m->IIF[i]->get_type().destroy_ion_if(m->IIF[i]);
  }
 
  for (int j = 0; j < NUM_IMP_DATA_TYPES; j++) {
    if (m->gdata[j])
      delete m->gdata[j];
  }
 
  m->zero_data();
}
 
#undef MEMFREE
 
void doppel_update(MULTI_IF *orig, MULTI_IF *miif_doppel) {
  // removed assertion, need to be able to update bidirectionally
  // assert(doppel->doppel);
 
  for (int i = 0; i < NUM_IMP_DATA_TYPES; i++)
    if (orig->gdata[i])
      SF::init_vector(&miif_doppel->gdata[i], orig->gdata[i]);
 
  for (int i = 0; i < orig->N_IIF; i++) {
    miif_doppel->IIF[i]->copy_SVs_from(*orig->IIF[i], false);
 
    for (int j = 0; j < orig->IIF[i]->plugins().size(); j++)
      miif_doppel->IIF[i]->plugins()[j]->copy_SVs_from(*orig->IIF[i]->plugins()[j], false);
  }
 
  miif_doppel->getRealData();
}
 
void doppel_MIIF(MULTI_IF *orig, MULTI_IF *miif_doppel) {
  *miif_doppel         = *orig;
  miif_doppel->doppel  = true;
 
  // deal with the local portion of the global data
  for (int i = 0; i < NUM_IMP_DATA_TYPES; i++)
    if (orig->gdata[i])
      SF::init_vector(&miif_doppel->gdata[i], orig->gdata[i]);
 
  // copy the state variables
  miif_doppel->IIF = {};
 
  for (int i = 0; i < orig->N_IIF; i++) {
    // Make a new IonIf for this clone and copy state variables and plugins
    miif_doppel->IIF.push_back(orig->IIF[i]->get_type().make_ion_if(orig->IIF[i]->get_target(),
          orig->IIF[i]->get_num_node(), orig->plugtypes[i]));
    miif_doppel->IIF[i]->copy_SVs_from(*orig->IIF[i], true);
    miif_doppel->IIF[i]->copy_plugins_from(*orig->IIF[i]);
  }
}
 
bool isIMPdata(const char *sv) {
  return IMPdataLabel2Index(sv) != -1;
}
 
int IMPdataLabel2Index(const char *sv) {
  int imp_data_id = -1;
 
  for (int i = 0; i < NUM_IMP_DATA_TYPES; i++) {
    if (strcmp(sv, imp_data_names[i]) == 0) {
      imp_data_id = i;
      break;
    }
  }
  return imp_data_id;
}
 
/* copy all of the IMP data from one node to another
 *
 * \param IF   ionic model
 * \param from local IMP number
 * \param to   local IMP number
 */
void dup_IMP_node_state(IonIfBase& IF, node_index_t from, node_index_t to, GlobalData_t **localdata) {
  for (int  i = 0; i < NUM_IMP_DATA_TYPES; i++)
    if (localdata && localdata[i])
      localdata[i][to] = localdata[i][from];
 
  IF.for_each([&](IonIfBase& imp) {
    char **sv_list;
    int sv_list_size = imp.get_type().get_sv_list(&sv_list);
 
    for (int i = 0; i < sv_list_size; i++) {
      char *sv_name = sv_list[i];
      int sv_offset;
      int sv_size;
      SVgetfcn sv_get = imp.get_type().get_sv_offset(sv_name,  &sv_offset, &sv_size);
      if (sv_get == NULL) {
        throw std::runtime_error(std::string(__func__) + " error: " + sv_name + " is not a valid state variable for the " + imp.get_type().get_name() + " model.");
      }
      SVputfcn sv_put = getPutSV(sv_get);
 
      GlobalData_t sv_val = sv_get(imp, from, sv_offset);
      sv_put(imp, to, sv_offset, sv_val);
    }
    free(sv_list);
  });
}
 
}  // namespace limpet