electrics.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
// ----------------------------------------------------------------------------
 
#include "electrics.h"
#include "petsc_utils.h"
#include "timers.h"
#include "stimulate.h"
#include "electric_integrators.h"
 
#include "SF_init.h" // for SF::init_xxx()
 
namespace opencarp {
 
void Electrics::initialize()
{
  double t1, t2;
  get_time(t1);
 
  set_dir(OUTPUT);
 
  // open logger
  logger = f_open("electrics.log", param_globals::experiment != 4 ? "w" : "r");
 
  // setup mappings between extra and intra grids, algebraic and nodal,
  // and between PETSc and canonical orderings
  setup_mappings();
 
  // the ionic physics is currently triggered from inside the Electrics to have tighter
  // control over it
  ion.logger = logger;
  ion.initialize();
 
  // set up Intracellular tissue
  set_elec_tissue_properties(mtype, intra_grid, logger);
  region_mask(intra_elec_msh, mtype[intra_grid].regions, mtype[intra_grid].regionIDs, true, "gregion_i");
 
  if (param_globals::bidomain || param_globals::extracell_monodomain_stim) {
    // set up Extracellular tissue
    set_elec_tissue_properties(mtype, extra_grid, logger);
    region_mask(extra_elec_msh, mtype[extra_grid].regions, mtype[extra_grid].regionIDs, true, "gregion_e");
  }
 
  // add electrics timer for time stepping, add to time stepper tool (TS)
  double global_time = user_globals::tm_manager->time;
  timer_idx = user_globals::tm_manager->add_eq_timer(global_time, param_globals::tend, 0,
                param_globals::dt, 0, "elec::ref_dt", "TS");
 
  // electrics stimuli setup
  setup_stimuli();
 
  // set up the linear equation systems. this needs to happen after the stimuli have been
  // set up, since we need boundary condition info
  setup_solvers();
 
  // the next setup steps require the solvers to be set up, since they use the matrices
  // generated by those
 
  // balance electrodes, we may need the extracellular mass matrix
  balance_electrodes();
  // total current scaling
  constant_total_stimulus_current(stimuli, *parab_solver.mass_i, *ellip_solver.mass_e, ion.miif, logger);
  // initialize the LATs detector
  lat.init(*parab_solver.Vmv, *ellip_solver.phie, 0);
 
  // initialize phie recovery data
  if(strlen(param_globals::phie_rec_ptf) > 0)
    setup_phie_recovery_data(phie_rcv);
 
  // prepare the electrics output. we skip it if we do post-processing
  if(param_globals::experiment != EXP_POSTPROCESS)
    setup_output();
 
  if (param_globals::prepacing_bcl > 0)
    prepace();
  this->initialize_time += timing(t2, t1);
}
 
void set_elec_tissue_properties(MaterialType* mtype, Electrics::grid_t g, FILE_SPEC logger)
{
  MaterialType *m = mtype+g;
 
  // initialize random conductivity fluctuation structure with PrM values
  m->regions.resize(param_globals::num_gregions);
 
  const char* grid_name = g == Electrics::intra_grid ? "intracellular" : "extracellular";
  log_msg(logger, 0, 0, "Setting up %s tissue properties for %d regions ..", grid_name,
          param_globals::num_gregions);
 
  char buf[64];
  RegionSpecs* reg = m->regions.data();
 
  for (size_t i=0; i<m->regions.size(); i++, reg++)  {
    if(!strcmp(param_globals::gregion[i].name, ""))  {
      snprintf(buf, sizeof buf, ", gregion_%d", int(i));
      param_globals::gregion[i].name = dupstr(buf);
    }
 
    reg->regname  = strdup(param_globals::gregion[i].name);
    reg->regID    = i;
    reg->nsubregs = param_globals::gregion[i].num_IDs;
    if(!reg->nsubregs)
      reg->subregtags = NULL;
    else  {
      reg->subregtags  = new int[reg->nsubregs];
      for (int j=0;j<reg->nsubregs;j++) {
        reg->subregtags[j] = param_globals::gregion[i].ID[j];
        if(reg->subregtags[j]==-1)
          log_msg(NULL,3,ECHO, "Warning: not all %u IDs provided for gregion[%u]!\n", reg->nsubregs, i);
      }
    }
 
    // describe material in given region
    elecMaterial *emat  = new elecMaterial();
    emat->material_type = ElecMat;
 
    emat->InVal[0]   = param_globals::gregion[i].g_il;
    emat->InVal[1]   = param_globals::gregion[i].g_it;
    emat->InVal[2]   = param_globals::gregion[i].g_in;
 
    emat->ExVal[0]   = param_globals::gregion[i].g_el;
    emat->ExVal[1]   = param_globals::gregion[i].g_et;
    emat->ExVal[2]   = param_globals::gregion[i].g_en;
 
    emat->BathVal[0] = param_globals::gregion[i].g_bath;
    emat->BathVal[1] = param_globals::gregion[i].g_bath;
    emat->BathVal[2] = param_globals::gregion[i].g_bath;
 
    // convert units from S/m -> mS/um
    for (int j=0; j<3; j++) {
      emat->InVal[j]   *= 1e-3 * param_globals::gregion[i].g_mult;
      emat->ExVal[j]   *= 1e-3 * param_globals::gregion[i].g_mult;
      emat->BathVal[j] *= 1e-3 * param_globals::gregion[i].g_mult;
    }
    reg->material = emat;
  }
 
  if((g == Electrics::intra_grid && strlen(param_globals::gi_scale_vec)) ||
     (g == Electrics::extra_grid && strlen(param_globals::ge_scale_vec)) )
  {
    mesh_t mt = g == Electrics::intra_grid ? intra_elec_msh : extra_elec_msh;
    sf_mesh & mesh = get_mesh(mt);
    const char* file = g == Electrics::intra_grid ? param_globals::gi_scale_vec : param_globals::ge_scale_vec;
 
    size_t num_file_entries = SF::root_count_ascii_lines(file, mesh.comm);
 
    if(num_file_entries != mesh.g_numelem)
      log_msg(0,4,0, "%s warning: number of %s conductivity scaling entries does not match number of elements!",
              __func__, get_mesh_type_name(mt));
 
    // set up parallel element vector and read data
    sf_vec *escale;
    SF::init_vector(&escale, get_mesh(mt), 1, sf_vec::elemwise);
    escale->read_ascii(g == Electrics::intra_grid ? param_globals::gi_scale_vec : param_globals::ge_scale_vec);
 
    if(get_size() > 1) {
      // set up element vector permutation and permute
      if(get_permutation(mt, ELEM_PETSC_TO_CANONICAL, 1) == NULL) {
        register_permutation(mt, ELEM_PETSC_TO_CANONICAL, 1);
      }
      SF::scattering & sc = *get_permutation(mt, ELEM_PETSC_TO_CANONICAL, 1);
      sc(*escale, false);
    }
 
    // copy data into SF::vector
    SF_real* p = escale->ptr();
    m->el_scale.assign(p, p + escale->lsize());
    escale->release_ptr(p);
  }
}
 
void Electrics::setup_mappings()
{
  bool intra_exits = mesh_is_registered(intra_elec_msh), extra_exists = mesh_is_registered(extra_elec_msh);
  assert(intra_exits);
  const int dpn = 1;
 
  // It may be that another physic (e.g. ionic models) has already computed the intracellular mappings,
  // thus we first test their existence
  if(get_scattering(intra_elec_msh, ALG_TO_NODAL, dpn) == NULL) {
    log_msg(logger, 0, 0, "%s: Setting up intracellular algebraic-to-nodal scattering.", __func__);
    register_scattering(intra_elec_msh, ALG_TO_NODAL, dpn);
  }
  if(get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, dpn) == NULL) {
    log_msg(logger, 0, 0, "%s: Setting up intracellular PETSc to canonical permutation.", __func__);
    register_permutation(intra_elec_msh, PETSC_TO_CANONICAL, dpn);
  }
 
  // extracellular mappings
  if(extra_exists) {
    log_msg(logger, 0, 0, "%s: Setting up extracellular algebraic-to-nodal scattering.", __func__);
    register_scattering(extra_elec_msh,  ALG_TO_NODAL, dpn);
    log_msg(logger, 0, 0, "%s: Setting up extracellular PETSc to canonical permutation.", __func__);
    register_permutation(extra_elec_msh, PETSC_TO_CANONICAL, dpn);
    log_msg(logger, 0, 0, "%s: Setting up intra-to-extra scattering.", __func__);
    register_scattering(intra_elec_msh, extra_elec_msh, dpn);
  }
 
  bool check_i2e = false;
  if(check_i2e && extra_exists) {
    sf_mesh & intra_mesh = get_mesh(intra_elec_msh);
    sf_mesh & extra_mesh = get_mesh(extra_elec_msh);
    int rank = get_rank();
 
    SF::scattering & i2e = *get_scattering(intra_elec_msh, extra_elec_msh, SF::NBR_PETSC, 1);
 
    const SF::vector<mesh_int_t> & intra_alg_nod   = intra_mesh.pl.algebraic_nodes();
    const SF::vector<mesh_int_t> & extra_alg_nod   = extra_mesh.pl.algebraic_nodes();
    const SF::vector<mesh_int_t> & extra_petsc_nbr = extra_mesh.get_numbering(SF::NBR_PETSC);
    const SF::vector<mesh_int_t> & intra_ref_nbr   = intra_mesh.get_numbering(SF::NBR_REF);
    const SF::vector<mesh_int_t> & extra_ref_nbr   = extra_mesh.get_numbering(SF::NBR_REF);
 
    // TODO(init) : delete these three at the end of this section?
    sf_vec *intra_testvec; SF::init_vector(&intra_testvec, intra_mesh, 1, sf_vec::algebraic);
    sf_vec *extra_testvec; SF::init_vector(&extra_testvec, extra_mesh, 1, sf_vec::algebraic);
    sf_vec *i2e_testvec; SF::init_vector(&i2e_testvec, extra_mesh, 1, sf_vec::algebraic);
 
    SF_real* id = intra_testvec->ptr();
    for(size_t i=0; i<intra_alg_nod.size(); i++) {
      int lpidx = local_nodal_to_local_petsc(intra_mesh, rank, intra_alg_nod[i]);
      id[lpidx] = intra_ref_nbr[intra_alg_nod[i]];
    }
    intra_testvec->release_ptr(id);
 
    SF_real* ed = extra_testvec->ptr();
    for(size_t i=0; i<extra_alg_nod.size(); i++) {
      int lpidx = local_nodal_to_local_petsc(extra_mesh, rank, extra_alg_nod[i]);
      ed[lpidx] = extra_ref_nbr[extra_alg_nod[i]];
    }
    extra_testvec->release_ptr(ed);
 
    i2e_testvec->set(-1.0);
    i2e.forward(*intra_testvec, *i2e_testvec);
 
    int err = 0;
    for(size_t i=0; i<extra_alg_nod.size(); i++) {
      auto id = i2e_testvec->get(i);
      auto ed = extra_testvec->get(i);
      if(id > -1 && id != ed)
        err++;
    }
 
    if(get_global(err, MPI_SUM))
      log_msg(0,5,0, "Electrics mapping test failed!");
    else
      log_msg(0,5,0, "Electrics mapping test succeeded!");
  }
}
 
void Electrics::compute_step()
{
  double t1, t2;
  get_time(t1);
 
  // if requested, we checkpoint the current state
  checkpointing();
 
  // activation checking
  const double time      = user_globals::tm_manager->time,
               time_step = user_globals::tm_manager->time_step;
  lat.check_acts(time);
  lat.check_quiescence(time, time_step);
 
  // I believe that we need to treat the stimuli in two ways:
  //    - Extracellular potential stimuli (this includes ground) affect the
  //      elliptic solver in a more delicate way, as such, there is a dbc_manager
  //      to take care of that.
  //    - Extracellular and Intracellular current stimuli are applied to the rhs vectors
  //      and can be managed by the stimulate() code directly.
  stimulate_extracellular();
 
  if(param_globals::bidomain == BIDOMAIN)
    ellip_solver.solve(*parab_solver.rhs_parab, *parab_solver.Vmv, *parab_solver.tmp_i1);
 
  clamp_Vm();
 
  // compute ionics update
  ion.compute_step();
 
  stimulate_intracellular();
 
  // store Vm before parabolic step, the full Ic we compute in the output step
  if(param_globals::dump_data & DUMP_IC)
    *parab_solver.Ic = *parab_solver.Vmv;
 
  // solver parabolic system
  parab_solver.solve(*ellip_solver.phie_i);
 
  clamp_Vm();
 
  if(user_globals::tm_manager->trigger(iotm_console)) {
    // output lin solver stats
    parab_solver.stats.log_stats(user_globals::tm_manager->time, false);
    if(param_globals::bidomain == BIDOMAIN)
      ellip_solver.stats.log_stats(user_globals::tm_manager->time, false);
  }
  this->compute_time += timing(t2, t1);
 
  // since the traces have their own timing, we check for trace dumps in the compute step loop
  if(user_globals::tm_manager->trigger(iotm_trace))
    dump_trace( ion.miif, user_globals::tm_manager->time );
 
}
 
void Electrics::output_step()
{
  double t1, t2;
  get_time(t1);
 
  const double time      = user_globals::tm_manager->time,
               time_step = user_globals::tm_manager->time_step;
 
  // for pseudo-bidomain we compute extracellular potential only for output
  if(param_globals::bidomain == PSEUDO_BIDM) {
    ellip_solver.solve(*parab_solver.rhs_parab, *parab_solver.Vmv, *parab_solver.tmp_i1);
    if(user_globals::tm_manager->trigger(iotm_console))
      ellip_solver.stats.log_stats(time, false);
  }
 
  if(param_globals::dump_data & DUMP_IVOL)
    parab_solver.rhs_parab->mult(*parab_solver.Vmv, *parab_solver.Ivol);
 
  if(param_globals::bidomain && (param_globals::dump_data & DUMP_IACT)) {
    parab_solver.rhs_parab->mult(*ellip_solver.phie_i, *parab_solver.Iact);
  }
 
  if(param_globals::dump_data & DUMP_IC) {
    PetscReal *Ic = parab_solver.Ic->ptr(), *Vmv = parab_solver.Vmv->ptr();
 
    for(PetscInt i=0; i < parab_solver.Ic->lsize(); i++)
      Ic[i] = (Ic[i] - Vmv[i]) / (-time_step);
 
    parab_solver.Ic->release_ptr(Ic), parab_solver.Vmv->release_ptr(Vmv);
  }
 
  // recover phie
  if(phie_rcv.pts.size()) {
    recover_phie_std(*parab_solver.Vmv, phie_rcv);
  }
 
  assemble_sv_gvec(gvec, ion.miif);
  output_manager.write_data();
 
  double curtime = timing(t2, t1);
  this->output_time += curtime;
 
  IO_stats.calls++;
  IO_stats.tot_time += curtime;
 
  if(user_globals::tm_manager->trigger(iotm_console))
    IO_stats.log_stats(time, false);
}
 
void Electrics::destroy()
{
  // output LAT data
  lat.output_initial_activations();
 
  // close logger
  f_close(logger);
 
  // close output files
  output_manager.close_files_and_cleanup();
 
  // destroy ionics
  ion.destroy();
}
 
void balance_electrode(elliptic_solver & ellip, SF::vector<stimulus> & stimuli, int balance_from, int balance_to)
{
  log_msg( NULL, 0, 0, "Balancing stimulus %d with %d %s-wise.",balance_from, balance_to,
      is_current(stimuli[balance_from].phys.type) ? "current" : "voltage" );
 
  stimulus & from = stimuli[balance_from];
  stimulus & to   = stimuli[balance_to];
 
  to.pulse     = from.pulse;
  to.ptcl      = from.ptcl;
  to.phys      = from.phys;
  to.pulse.strength *= -1.0;
 
  if (from.phys.type == I_ex)
  {
    // if from is total current, skip volume based adjustment of strength
    // otherwise, calling constant_total_stimulus_current() will undo the balanced scaling of to.pulse.strength
    // constant_total_stimulus_current() will do the scaling based on the volume
    if (!from.phys.total_current) { 
      sf_mat& mass = *ellip.mass_e;
      SF_real vol0 = get_volume_from_nodes(mass, from.electrode.vertices);
      SF_real vol1 = get_volume_from_nodes(mass, to.electrode.vertices);
 
      to.pulse.strength *= fabs(vol0 / vol1);
    }
  }
}
 
void Electrics::balance_electrodes()
{
  for(int i=0; i<param_globals::num_stim; i++) {
    if(param_globals::stim[i].crct.balance != -1) {
      int from = param_globals::stim[i].crct.balance;
      int to   = i;
 
      balance_electrode(this->ellip_solver, stimuli, from, to);
    }
  }
}
 
void Electrics::setup_stimuli()
{
  // initialize basic stim info data (used units, supported types, etc)
  init_stim_info();
  bool dumpTrace = true;
 
  if(dumpTrace) set_dir(OUTPUT);
 
  stimuli.resize(param_globals::num_stim);
  for(int i=0; i<param_globals::num_stim; i++)
  {
    // construct new stimulus
    stimulus & s = stimuli[i];
 
    if(user_globals::using_legacy_stimuli)
      s.translate(i);
 
    s.setup(i);
 
    if(dumpTrace && get_rank() == 0)
      s.pulse.wave.write_trace(s.name+".trc");
  }
}
 
void apply_stim_to_vector(const stimulus & s, sf_vec & vec, bool add)
{
  double val; s.value(val);
  const SF::vector<mesh_int_t> & idx = s.electrode.vertices;
  const int rank = get_rank();
  SF::vector<mesh_int_t> local_idx = idx;
  for (size_t i = 0; i < idx.size(); i++) {
      local_idx[i] = local_nodal_to_local_petsc(*vec.mesh, rank, idx[i]);
  }
  vec.set(local_idx, val, add, true);
}
 
void Electrics::stimulate_intracellular()
{
  parabolic_solver & ps = parab_solver;
 
  // iterate over stimuli
  for(stimulus & s : stimuli) {
    if(s.is_active()) {
      // for active stimuli, deal with the stimuli-type specific stimulus application
      switch(s.phys.type)
      {
        case I_tm: {
          if(param_globals::operator_splitting) {
            apply_stim_to_vector(s, *ps.Vmv, true);
          }
          else {
            SF_real Cm = 1.0;
            timer_manager & tm = *user_globals::tm_manager;
            SF_real sc = tm.time_step / Cm;
 
            ps.Irhs->set(0.0);
            apply_stim_to_vector(s, *ps.Irhs, true);
 
            *ps.tmp_i1 =  *ps.IIon;
            *ps.tmp_i1 -= *ps.Irhs;
            *ps.tmp_i1 *= sc;    // tmp_i1 = sc * (IIon - Irhs)
 
            // add ionic, transmembrane and intracellular currents to rhs
            if(param_globals::parab_solve != parabolic_solver::EXPLICIT)
              ps.mass_i->mult(*ps.tmp_i1, *ps.Irhs);
            else
              *ps.Irhs = *ps.tmp_i1;
          }
          break;
        }
 
        case Illum: {
          sf_vec* illum_vec = ion.miif->gdata[limpet::illum];
 
          if(illum_vec == NULL) {
            log_msg(0,5,0, "Cannot apply illumination stim: global vector not present!");
            EXIT(EXIT_FAILURE);
          } else {
            apply_stim_to_vector(s, *illum_vec, false);
          }
 
          break;
        }
 
        default: break;
      }
    }
  }
}
 
void Electrics::clamp_Vm() {
  for(stimulus & s : stimuli) {
    if(s.phys.type == Vm_clmp && s.is_active())
      apply_stim_to_vector(s, *parab_solver.Vmv, false);
  }
}
 
void Electrics::stimulate_extracellular()
{
  if(param_globals::bidomain) {
    // we check if the DBC layout changed, if so we recompute the matrix and the dbc_manager
    bool dbcs_have_updated = ellip_solver.dbc != nullptr && ellip_solver.dbc->dbc_update();
    bool time_not_final    = user_globals::tm_manager->d_time != user_globals::tm_manager->d_end;
 
    if(dbcs_have_updated && time_not_final)
      ellip_solver.rebuild_stiffness(mtype, stimuli, logger);
 
    ellip_solver.phiesrc->set(0.0);
 
    for(const stimulus & s : stimuli) {
      if(s.is_active() && s.phys.type == I_ex)
        apply_stim_to_vector(s, *ellip_solver.phiesrc, false);
    }
  }
}
 
void compute_restr_idx(sf_mesh & mesh,
                       SF::vector<mesh_int_t> & inp_idx,
                       SF::vector<mesh_int_t> & idx)
{
  int mpi_rank = get_rank(), mpi_size = get_size();
  const SF::vector<mesh_int_t> & layout = mesh.pl.algebraic_layout();
 
  SF::vector<mesh_int_t> sndbuff;
 
  size_t buffsize = 0;
  idx.resize(0);
 
  for(int pid=0; pid < mpi_size; pid++) {
    if(mpi_rank == pid) {
      sndbuff = inp_idx;
      buffsize = sndbuff.size();
    }
 
    MPI_Bcast(&buffsize, sizeof(size_t), MPI_BYTE, pid, PETSC_COMM_WORLD);
    sndbuff.resize(buffsize);
    MPI_Bcast(sndbuff.data(), buffsize*sizeof(mesh_int_t), MPI_BYTE, pid, PETSC_COMM_WORLD);
 
    mesh_int_t start = layout[mpi_rank], stop = layout[mpi_rank+1];
 
    for(mesh_int_t i : sndbuff) {
      if(i >= start && i < stop)
        idx.push_back(i - start);
    }
  }
 
  binary_sort(idx); unique_resize(idx);
}
 
void compute_restr_idx_async(sf_mesh & mesh,
                             SF::vector<mesh_int_t> & inp_idx,
                             SF::vector<mesh_int_t> & idx)
{
  int mpi_rank = get_rank(), mpi_size = get_size();
  const SF::vector<mesh_int_t> & alg_nod = mesh.pl.algebraic_nodes();
  const SF::vector<mesh_int_t> & nbr     = mesh.get_numbering(SF::NBR_SUBMESH);
 
  hashmap::unordered_map<mesh_int_t,mesh_int_t> amap;
  for(mesh_int_t ii : alg_nod)
    amap[nbr[ii]] = ii;
 
  SF::vector<mesh_int_t> sndbuff;
  size_t buffsize = 0;
  idx.resize(0);
 
  for(int pid=0; pid < mpi_size; pid++) {
    if(mpi_rank == pid) {
      sndbuff = inp_idx;
      buffsize = sndbuff.size();
    }
 
    MPI_Bcast(&buffsize, sizeof(size_t), MPI_BYTE, pid, PETSC_COMM_WORLD);
    sndbuff.resize(buffsize);
    MPI_Bcast(sndbuff.data(), buffsize*sizeof(mesh_int_t), MPI_BYTE, pid, PETSC_COMM_WORLD);
 
    for(mesh_int_t i : sndbuff) {
      if(amap.count(i))
        idx.push_back(amap[i]);
    }
  }
 
  binary_sort(idx); unique_resize(idx);
}
 
void setup_dataout(const int dataout, std::string dataout_vtx, mesh_t grid,
                   SF::vector<mesh_int_t>* & restr, bool async)
{
  sf_mesh & mesh = get_mesh(grid);
 
  switch(dataout) {
 
    case DATAOUT_SURF: {
      sf_mesh surfmesh;
      compute_surface_mesh(mesh, SF::NBR_SUBMESH, surfmesh);
 
      SF::vector<mesh_int_t> idxbuff(surfmesh.con);
      binary_sort(idxbuff); unique_resize(idxbuff);
 
      restr = new SF::vector<mesh_int_t>();
 
      // for sync output, we need restr to hold the local indices in the petsc vectors
      // that have been permuted to canonical numbering. For async, we need the
      // non-overlapping decomposition of indices in NBR_SUBMESH numbering. The petsc indices will be
      // computed at a later stage. The only reason we need to call compute_restr_idx_async,
      // is that surface nodes in NBR_SUBMESH, may
      // reside on partitions where they are not part of the algebraic nodes. thus we need to
      // recommunicate to make sure the data layout is correct. We do not have this problem for
      // DATAOUT_VTX.
      if(!async)
        compute_restr_idx(mesh, idxbuff, *restr);
      else
        compute_restr_idx_async(mesh, idxbuff, *restr);
 
      break;
    }
 
    case DATAOUT_VTX: {
      SF::vector<mesh_int_t> idxbuff;
 
      update_cwd();
 
      set_dir(INPUT);
      read_indices(idxbuff, dataout_vtx, mesh, SF::NBR_REF, true, PETSC_COMM_WORLD);
      set_dir(CURDIR);
 
      restr = new SF::vector<mesh_int_t>();
 
      if(!async) {
        const SF::vector<mesh_int_t> & nbr = mesh.get_numbering(SF::NBR_SUBMESH);
        for(mesh_int_t & i : idxbuff) i = nbr[i];
 
        compute_restr_idx(mesh, idxbuff, *restr);
      } else {
        *restr = idxbuff;
      }
 
      break;
    }
 
    case DATAOUT_NONE:
    case DATAOUT_VOL:
    default: break;
  }
}
 
void Electrics::setup_output()
{
  int rank = get_rank();
  SF::vector<mesh_int_t>* restr_i = NULL;
  SF::vector<mesh_int_t>* restr_e = NULL;
  set_dir(OUTPUT);
 
  setup_dataout(param_globals::dataout_i, param_globals::dataout_i_vtx, intra_elec_msh,
                restr_i, param_globals::num_io_nodes > 0);
 
  if(param_globals::dataout_i)
    output_manager.register_output(parab_solver.Vmv, intra_elec_msh, 1, param_globals::vofile, "mV", restr_i);
 
  if(param_globals::bidomain) {
    setup_dataout(param_globals::dataout_e, param_globals::dataout_e_vtx, extra_elec_msh,
                  restr_e, param_globals::num_io_nodes > 0);
 
    if(param_globals::dataout_i)
      output_manager.register_output(ellip_solver.phie_i, intra_elec_msh, 1, param_globals::phieifile, "mV", restr_i);
    if(param_globals::dataout_e)
      output_manager.register_output(ellip_solver.phie,   extra_elec_msh, 1, param_globals::phiefile,  "mV", restr_e);
  }
 
  if(param_globals::dump_data & DUMP_IC) {
    output_manager.register_output(parab_solver.Ic,   intra_elec_msh, 1, "Ic.igb",  "uA/cm^2", restr_i);
    output_manager.register_output(parab_solver.IIon, intra_elec_msh, 1, "Iion.igb","uA/cm^2", restr_i);
  }
 
  if(param_globals::dump_data & DUMP_IVOL)
    output_manager.register_output(parab_solver.Ivol, intra_elec_msh, 1, "Ivol.igb", "uA", restr_i);
  if(param_globals::dump_data & DUMP_IACT)
    output_manager.register_output(parab_solver.Iact, intra_elec_msh, 1, "Iact.igb", "uA", restr_i);
 
  if(phie_rcv.pts.size())
    output_manager.register_output_sync(phie_rcv.phie_rec, phie_recv_msh, 1, param_globals::phie_recovery_file, "mV");
 
  init_sv_gvec(gvec, ion.miif, *parab_solver.Vmv, output_manager);
 
  if(param_globals::num_trace) {
    sf_mesh & imesh = get_mesh(intra_elec_msh);
    open_trace(ion.miif, param_globals::num_trace, param_globals::trace_node, NULL, &imesh);
  }
 
  // initialize generic logger for IO timings per time_dt
  IO_stats.init_logger("IO_stats.dat");
}
 
void Electrics::dump_matrices()
{
  std::string bsname = param_globals::dump_basename;
  std::string fn;
 
  set_dir(OUTPUT);
 
  // dump monodomain matrices
  if ( param_globals::parab_solve==1 )  {
    // using Crank-Nicolson
    fn = bsname + "_Ki_CN.bin";
    parab_solver.lhs_parab->write(fn.c_str());
  }
  fn = bsname + "_Ki.bin";
  parab_solver.rhs_parab->write(fn.c_str());
 
  fn = bsname + "_Mi.bin";
  parab_solver.mass_i->write(fn.c_str());
 
  if ( param_globals::bidomain )  {
    fn = bsname + "_Kie.bin";
    ellip_solver.phie_mat->write(fn.c_str());
 
    fn = bsname + "_Me.bin";
    ellip_solver.mass_e->write(fn.c_str());
  }
}
 
 
double Electrics::timer_val(const int timer_id)
{
  // determine
  int sidx = stimidx_from_timeridx(stimuli, timer_id);
  double val = 0.0;
  if(sidx != -1) {
    stimuli[sidx].value(val);
  }
  else
    val = std::nan("NaN");
 
  return val;
}
 
std::string Electrics::timer_unit(const int timer_id)
{
  int sidx = stimidx_from_timeridx(stimuli, timer_id);
  std::string s_unit;
 
  if(sidx != -1)
    // found a timer-linked stimulus
    s_unit = stimuli[sidx].pulse.wave.f_unit;
 
  return s_unit;
}
 
int stimidx_from_timeridx(const SF::vector<stimulus> & stimuli, const int timer_id)
{
  // the only electrical quantities linked to a timer are stimuli
  // thus we search for timer links only among stimuli for now
 
  // iterate over stimuli
  for(size_t i = 0; i<stimuli.size(); i++)
  {
    const stimulus & s = stimuli[i];
 
    if(s.ptcl.timer_id == timer_id)
      return s.idx;
  }
 
  // invalid timer index not linked to any stimulus
  return -1;
}
 
void get_kappa(sf_vec & kappa, IMPregion *ir, limpet::MULTI_IF & miif, double k)
{
  double* reg_kappa = new double[miif.N_IIF];
 
  for(int i=0; i<miif.N_IIF; i++)
    reg_kappa[i] = k * miif.IIF[i]->cgeom().SVratio * ir[i].volFrac;
 
  double *kd = kappa.ptr();
 
  for(int i = 0; i < miif.numNode; i++)
    kd[i] = reg_kappa[(int) miif.IIFmask[i]];
 
  kappa.release_ptr(kd);
  delete [] reg_kappa;
}
 
 
void set_cond_type(MaterialType & m, cond_t type)
{
  for(size_t i=0; i < m.regions.size(); i++)  {
    elecMaterial *emat = static_cast<elecMaterial*>(m.regions[i].material);
    emat->g = type;
  }
}
 
void Electrics::setup_solvers()
{
  set_dir(OUTPUT);
  parab_solver.init();
  parab_solver.rebuild_matrices(mtype, *ion.miif, logger);
 
  if (param_globals::bidomain) {
    ellip_solver.init();
    ellip_solver.rebuild_matrices(mtype, stimuli, logger);
  }
 
  if(param_globals::dump2MatLab)
    dump_matrices();
}
 
bool have_dbc_stims(const SF::vector<stimulus> & stimuli)
{
  for(const stimulus & s : stimuli) {
    if(is_dbc(s.phys.type))
      return true;
  }
  return false;
}
 
const char* get_tsav_ext(double time)
{
  int    min_idx  = -1;
  double min_diff = 1e100;
 
  for(int i=0; i<param_globals::num_tsav; i++)
  {
    double diff = fabs(param_globals::tsav[i] - time);
    if(min_diff > diff) {
      min_diff = diff;
      min_idx  = i;
    }
  }
 
  if(min_idx == -1)
    min_idx = 0;
 
  return param_globals::tsav_ext[min_idx];
}
 
void Electrics::checkpointing()
{
  const timer_manager & tm = *user_globals::tm_manager;
 
  // regular user selected state save
  if (tm.trigger(iotm_chkpt_list)) {
    char save_fnm[1024];
    const char* tsav_ext = get_tsav_ext(tm.time);
 
    snprintf(save_fnm, sizeof save_fnm, "%s.%s.roe", param_globals::write_statef, tsav_ext);
 
    ion.miif->dump_state(save_fnm, tm.time, intra_elec_msh, false, GIT_COMMIT_COUNT);
  }
 
  // checkpointing based on interval
  if (tm.trigger(iotm_chkpt_intv)) {
    char save_fnm[1024];
    snprintf(save_fnm, sizeof save_fnm, "checkpoint.%.1f.roe", tm.time);
    ion.miif->dump_state(save_fnm, tm.time, intra_elec_msh, false, GIT_COMMIT_COUNT);
  }
}
 
void elliptic_solver::init()
{
  double t0, t1, dur;
  get_time(t0);
  stats.init_logger("ell_stats.dat");
 
  // here we can differentiate the solvers
  SF::init_solver(&lin_solver);
  sf_mesh & extra_mesh = get_mesh(extra_elec_msh);
  sf_vec::ltype alg_type = sf_vec::algebraic;
  const int dpn = 1;
 
  SF::init_vector(&phie, extra_mesh, dpn, alg_type);
  SF::init_vector(&phiesrc, extra_mesh, dpn, alg_type);
  SF::init_vector(&currtmp, extra_mesh, dpn, alg_type);
 
  if(mesh_is_registered(intra_elec_msh)) {
    sf_mesh & intra_mesh = get_mesh(intra_elec_msh);
    SF::init_vector(&phie_i, intra_mesh, dpn, alg_type);
  }
 
  int max_row_entries = max_nodal_edgecount(extra_mesh);
 
  SF::init_matrix(&phie_mat);
  SF::init_matrix(&mass_e);
 
  // alloc stiffness matrix
  phie_mat->init(extra_mesh, dpn, dpn, max_row_entries);
  // alloc mass matrix
  mass_e  ->init(extra_mesh, dpn, dpn, param_globals::mass_lumping ? 1 : max_row_entries);
  dur = timing(t1, t0);
}
 
void elliptic_solver::rebuild_matrices(MaterialType* mtype,
                                       SF::vector<stimulus> & stimuli,
                                       FILE_SPEC logger)
{
  double t0, t1, dur;
  get_time(t0);
  rebuild_stiffness(mtype, stimuli, logger);
  rebuild_mass(logger);
  dur = timing(t1, t0);
}
 
void elliptic_solver::rebuild_stiffness(MaterialType* mtype,
                                        SF::vector<stimulus> & stimuli,
                                        FILE_SPEC logger)
{
  double t0, t1, dur;
  int log_flag = param_globals::output_level > 1 ? ECHO : 0;
 
  MaterialType & mt = mtype[Electrics::extra_grid];
  const bool have_dbc = have_dbc_stims(stimuli);
 
  cond_t condType  = sum_cond;
  set_cond_type(mt, condType);
 
  // get mesh reference
  sf_mesh & mesh = get_mesh(extra_elec_msh);
 
  get_time(t0);
 
  // fill the system
  elec_stiffness_integrator stfn_integ(mt);
 
  phie_mat->zero();
  SF::assemble_matrix(*phie_mat, mesh, stfn_integ);
  phie_mat->scale(-1.0);
 
  dur = timing(t1,t0);
  log_msg(logger,0,log_flag, "Computed ellipitc stiffness matrix in %.3f seconds.", dur);
 
  // set boundary conditions
  if(have_dbc) {
    log_msg(logger,0,log_flag, "Elliptic lhs matrix enforcing Dirichlet boundaries.");
    get_time(t0);
 
    if(dbc == nullptr)
      dbc = new dbc_manager(*phie_mat, stimuli);
    else
      dbc->recompute_dbcs();
 
    dbc->enforce_dbc_lhs();
 
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Elliptic lhs matrix Dirichlet enforcing done in %.3f seconds.", dur);
  }
  else {
    log_msg(logger,1,ECHO, "Elliptic lhs matrix is singular!");
    // we are dealing with a singular system
    phie_mat_has_nullspace = true;
  }
 
  // solver has not been initialized yet
  set_dir(INPUT);
  get_time(t0);
 
  setup_linear_solver(logger);
 
  dur = timing(t1,t0);
  log_msg(logger,0,log_flag, "Initializing elliptic solver in %.5f seconds.", dur);
  set_dir(OUTPUT);
}
 
void elliptic_solver::rebuild_mass(FILE_SPEC logger)
{
  int log_flag = param_globals::output_level > 1 ? ECHO : 0;
  double t0, t1, dur;
  mass_integrator mass_integ;
 
  // get mesh reference
  sf_mesh & mesh = get_mesh(extra_elec_msh);
  get_time(t0);
  mass_e->zero();
 
  if(param_globals::mass_lumping) {
    SF::assemble_lumped_matrix(*mass_e, mesh, mass_integ);
  } else {
    SF::assemble_matrix(*mass_e, mesh, mass_integ);
  }
 
  dur = timing(t1,t0);
  log_msg(logger,0,log_flag, "Computed elliptic mass matrix in %.3f seconds.", dur);
}
 
void elliptic_solver::setup_linear_solver(FILE_SPEC logger)
{
 
  tol    = param_globals::cg_tol_ellip;
  max_it = param_globals::cg_maxit_ellip;
 
  std::string solver_file;
  solver_file = param_globals::ellip_options_file;
  lin_solver->setup_solver(*phie_mat, tol, max_it, param_globals::cg_norm_ellip,
                           "elliptic PDE", phie_mat_has_nullspace,
                           logger, solver_file.c_str(),
                           "-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg "
                           "-pc_hypre_boomeramg_max_iter 1 "
                           "-pc_hypre_boomeramg_strong_threshold 0.0 -options_left");
}
 
void elliptic_solver::solve(sf_mat & Ki, sf_vec & Vmv, sf_vec & tmp_i)
{
  double t0,t1;
  SF::scattering* i2e = get_scattering(intra_elec_msh, extra_elec_msh, 1);
  // assembly of rhs for FE
  if (phiesrc->mag() > 0.0) {
    mass_e->mult(*phiesrc, *currtmp);
    phiesrc->deep_copy(*currtmp);
  }
 
  Ki.mult(Vmv, tmp_i);
 
  bool add = true;
  i2e->forward(tmp_i, *phiesrc, add);
 
  if(dbc != nullptr)
    dbc->enforce_dbc_rhs(*phiesrc);
 
  get_time(t0);
  (*lin_solver)(*phie, *phiesrc);
 
  // treat solver statistics
  auto dur = timing(t1, t0);
  lin_solver->time += dur;
  stats.slvtime += dur;
  stats.update_iter(lin_solver->niter);
  if(lin_solver->reason < 0) {
    log_msg(0, 5, 0,"%s solver diverged. Reason: %s.", lin_solver->name.c_str(),
            petsc_get_converged_reason_str(lin_solver->reason));
    EXIT(1);
  }
 
  add = false;
  i2e->backward(*phie, *phie_i, add);
}
 
void elliptic_solver::solve_laplace()
{
  double t0,t1;
 
  if(dbc != nullptr)
    dbc->enforce_dbc_rhs(*phiesrc);
 
  get_time(t0);
  (*lin_solver)(*phie, *phiesrc);
 
  // treat solver statistics
  auto dur = timing(t1, t0);
  lin_solver->time += dur;
  stats.slvtime += dur;
  stats.update_iter(lin_solver->niter);
 
  if(lin_solver->reason < 0) {
    log_msg(0, 5, 0,"%s solver diverged. Reason: %s.", lin_solver->name.c_str(),
            petsc_get_converged_reason_str(lin_solver->reason));
    EXIT(1);
  }
 
  // phie_i is only set up when we have an IntraMesh registered
  if(is_init(phie_i)) {
    bool add = false;
    SF::scattering* i2e = get_scattering(intra_elec_msh, extra_elec_msh, 1);
    i2e->backward(*phie, *phie_i, add);
  }
}
 
void parabolic_solver::init()
{
  double t0, t1, dur;
  get_time(t0);
  stats.init_logger("par_stats.dat");
 
  // here we can differentiate the solvers
  SF::init_solver(&lin_solver);
 
  sf_vec* vm_ptr   = get_data(vm_vec);
  sf_vec* iion_ptr = get_data(iion_vec);
 
  if(!(vm_ptr != NULL && iion_ptr != NULL)) {
    log_msg(0,5,0, "%s error: global Vm and Iion vectors not properly set up! Ionics seem invalid! Aborting!",
            __func__);
    EXIT(1);
  }
 
  SF::init_vector(&Vmv);
  SF::init_vector(&IIon);
  Vmv-> shallow_copy(*vm_ptr);
  IIon->shallow_copy(*iion_ptr);
 
  if(param_globals::dump_data & DUMP_IC)   SF::init_vector(&Ic  , Vmv);
  if(param_globals::dump_data & DUMP_IVOL) SF::init_vector(&Ivol, Vmv);
  if(param_globals::dump_data & DUMP_IACT) SF::init_vector(&Iact, Vmv);
 
  SF::init_vector(&Diff_term, Vmv);
  SF::init_vector(&old_vm, Vmv);
 
  sf_mesh & intra_mesh = get_mesh(intra_elec_msh);
  sf_vec::ltype alg_type = sf_vec::algebraic;
 
  int dpn = 1;
  SF::init_vector(&kappa_i, intra_mesh, dpn, alg_type);
  SF::init_vector(&tmp_i1, intra_mesh, dpn, alg_type);
  SF::init_vector(&tmp_i2, intra_mesh, dpn, alg_type);
  SF::init_vector(&old_vm, intra_mesh, dpn, alg_type);
 
  if(!param_globals::operator_splitting)
    SF::init_vector(&Irhs, intra_mesh, dpn, alg_type);
 
  // alloc matrices
  int max_row_entries = max_nodal_edgecount(intra_mesh);
 
  SF::init_matrix(&lhs_parab);
  SF::init_matrix(&rhs_parab);
  SF::init_matrix(&mass_i);
 
  rhs_parab->init(intra_mesh, dpn, dpn, max_row_entries);
  mass_i   ->init(intra_mesh, dpn, dpn, param_globals::mass_lumping ? 1 : max_row_entries);
 
  parab_tech = static_cast<parabolic_solver::parabolic_t>(param_globals::parab_solve);
  dur = timing(t1, t0);
}
 
void parabolic_solver::rebuild_matrices(MaterialType* mtype, limpet::MULTI_IF & miif, FILE_SPEC logger)
{
  double start, end, period;
  get_time(start);
  double t0, t1, dur;
  mass_integrator mass_integ;
  int dpn = 1;
 
  int            log_flag = param_globals::output_level > 1 ? ECHO : 0;
  MaterialType & mt       = mtype[Electrics::intra_grid];
 
  double Dt = user_globals::tm_manager->time_step;
  get_kappa(*kappa_i, param_globals::imp_region, miif, UM2_to_CM2 / Dt);
 
  cond_t    condType = intra_cond;
  sf_mesh & mesh     = get_mesh(intra_elec_msh);
 
  if( (param_globals::bidomain == MONODOMAIN && param_globals::bidm_eqv_mono) ||
      (param_globals::bidomain == PSEUDO_BIDM) )
    condType = para_cond;
 
  // set material and conductivity type
  set_cond_type(mt, condType);
 
  // fill the system
  {
    get_time(t0);
 
    elec_stiffness_integrator stfn_integ(mt);
    SF::assemble_matrix(*rhs_parab, mesh, stfn_integ);
 
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Computed parabolic stiffness matrix in %.3f seconds.", dur);
    get_time(t0);
 
    if(param_globals::mass_lumping)
      SF::assemble_lumped_matrix(*mass_i, mesh, mass_integ);
    else
      SF::assemble_matrix(*mass_i, mesh, mass_integ);
 
    sf_vec* empty; SF::init_vector(&empty);
    mass_i->mult_LR(*kappa_i, *empty);
 
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Computed parabolic mass matrix in %.3f seconds.", dur);
    delete empty;
  }
 
  // initialize parab lhs
  if(parab_tech != EXPLICIT) {
    lhs_parab->duplicate(*rhs_parab);
    // if we have mass lumping, then the nonzero pattern between Mi and Ki is different
    bool same_nonzero = param_globals::mass_lumping == false;
 
    if (parab_tech==CN) {
      lhs_parab->scale(-param_globals::theta);
      lhs_parab->add_scaled_matrix(*mass_i, 1.0, same_nonzero);
    }
    else if (parab_tech==O2dT) {
      lhs_parab->scale(-0.5);
      mass_i->scale(0.5);
      lhs_parab->add_scaled_matrix(*mass_i, 1.0, same_nonzero);
      lhs_parab->add_scaled_matrix(*mass_i, 1.0, same_nonzero);
      lhs_parab->add_scaled_matrix(*mass_i, 1.0, same_nonzero);
    }
  }
  else {
    SF::init_vector(&inv_mass_diag, Vmv);
    mass_i->get_diagonal(*inv_mass_diag);
 
    SF_real* p = inv_mass_diag->ptr();
 
    for(int i=0; i<inv_mass_diag->lsize(); i++)
      p[i] = 1.0 / p[i];
 
    inv_mass_diag->release_ptr(p);
  }
 
  if(parab_tech == CN || parab_tech == O2dT) {
    set_dir(INPUT);
    get_time(t0);
 
    setup_linear_solver(logger);
 
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Initializing parabolic solver in %.5f seconds.", dur);
    set_dir(OUTPUT);
  }
  period = timing(end, start);
}
 
void parabolic_solver::setup_linear_solver(FILE_SPEC logger)
{
  tol    = param_globals::cg_tol_parab;
  max_it = param_globals::cg_maxit_parab;
 
  std::string solver_file;
  solver_file = param_globals::parab_options_file;
  lin_solver->setup_solver(*lhs_parab, tol, max_it, param_globals::cg_norm_parab,
                           "parabolic PDE", false, logger, solver_file.c_str(),
                           "-pc_type bjacobi -sub_pc_type ilu -ksp_type cg");
}
 
void parabolic_solver::solve(sf_vec & phie_i)
{
  switch (parab_tech) {
    case CN:   solve_CN(phie_i);   break;
    case O2dT: solve_O2dT(phie_i); break;
    default:   solve_EF(phie_i);   break;
  }
}
 
void parabolic_solver::solve_CN(sf_vec & phie_i)
{
  double t0,t1;
  // assembly of rhs for CN
  if (param_globals::bidomain == BIDOMAIN) {
    tmp_i1->deep_copy(phie_i);
    tmp_i1->add_scaled(*Vmv, 1.0 - param_globals::theta);
    rhs_parab->mult(*tmp_i1, *tmp_i2);
  }
  else {
    rhs_parab->mult(*Vmv, *tmp_i2);
    *tmp_i2 *= 1.0 - param_globals::theta;
  }
 
  mass_i->mult(*Vmv, *tmp_i1);
  *tmp_i1 += *tmp_i2;
 
  // add current contributions to rhs
  if(!param_globals::operator_splitting)
    tmp_i1->add_scaled(*Irhs, -1.0);
 
  get_time(t0);
 
  (*lin_solver)(*Vmv, *tmp_i1);
 
  if(lin_solver->reason < 0) {
    log_msg(0, 5, 0,"%s solver diverged. Reason: %s.", lin_solver->name.c_str(),
            petsc_get_converged_reason_str(lin_solver->reason));
    EXIT(1);
  }
 
  // treat solver statistics
  auto dur = timing(t1, t0);
  lin_solver->time += dur;
  stats.slvtime += dur;
  stats.update_iter(lin_solver->niter);
}
 
void parabolic_solver::solve_O2dT(sf_vec & phie_i)
{
  double t0,t1;
  // assembly of rhs for FE
  if (param_globals::bidomain == BIDOMAIN) {
    tmp_i2->deep_copy(phie_i);
    tmp_i2->add_scaled(*Vmv, 0.5);
    rhs_parab->mult(*tmp_i2, *tmp_i1);  // tmp_i1 = K_i(Vm^t * 0.5 + phi_e)
  }
  else {
    rhs_parab->mult(*Vmv, *tmp_i1);
    *tmp_i1 *= 0.5;                   // tmp_i1 = 0.5 * K_i Vm^t
  }
 
  mass_i->mult(*Vmv, *tmp_i2);          // tmp_i2 = M/2 Vm^t
  tmp_i1->add_scaled(*tmp_i2, 4.0);    // tmp_i1 = (2M+K_i/2)Vm^t
  mass_i->mult(*old_vm, *tmp_i2);       // tmp_i2 = M/2 Vm^{t-1}
 
  tmp_i1->add_scaled(*tmp_i2, -1.0);   // tmp_i1 = (2M+K_i/2)Vm^t-M/2 Vm^{t-1}
  *old_vm = *Vmv;
 
  get_time(t0);
 
  // solve
  (*lin_solver)(*Vmv, *tmp_i1);
 
  // treat solver statistics
  stats.slvtime += timing(t1, t0);
  stats.update_iter(lin_solver->niter);
}
 
void parabolic_solver::solve_EF(sf_vec & phie_i)
{
  double t0,t1,t2;
  get_time(t0);
 
  // assembly of rhs for FE
  if (param_globals::bidomain == BIDOMAIN) {
    tmp_i2->deep_copy(phie_i);
    *tmp_i2 += *Vmv;
    rhs_parab->mult(*tmp_i2, *tmp_i1);
  }
  else {
    rhs_parab->mult(*Vmv, *tmp_i1);
  }
 
  *tmp_i1 *= *inv_mass_diag;
  Vmv->add_scaled(*tmp_i1, 1.0);
 
  if(param_globals::operator_splitting == false)
    Vmv->add_scaled(*Irhs, -1.0);
 
  // record rhs timing
  stats.slvtime += timing(t1, t0);
}
 
LAT_detector::LAT_detector()
{
  char* prvSimDir = strlen(param_globals::start_statef) ?
                    get_file_dir(param_globals::start_statef) : NULL;
 
  const char* extn = ".dat";
 
  // if compute_APD we need an extra 2 acts
  int addLATs = param_globals::compute_APD ? 2 : 0;
 
  bool have_sentinel = param_globals::t_sentinel > 0.0;
  bool need_to_add_sentinel = have_sentinel && (param_globals::sentinel_ID < 0);
 
  addLATs += need_to_add_sentinel ? 1 : 0;
  acts.resize(param_globals::num_LATs + addLATs);
 
  int j=0;
  for (int i = 0; i < param_globals::num_LATs; i++ )
  {
    // using Ph only with bidomain runs
    if (param_globals::lats[i].method <= 0 || (param_globals::lats[i].measurand == PHIE && !param_globals::bidomain))  {
      log_msg(NULL, 3, 0, "Phie-based LAT measurement requires bidomain >=1 Ignoring lats[%d].", i);
      continue;
    }
 
    acts[j].method    = (ActMethod)param_globals::lats[i].method;
    acts[j].threshold = param_globals::lats[i].threshold;
    acts[j].mode      = param_globals::lats[i].mode;
    acts[j].all       = param_globals::lats[i].all;
    acts[j].measurand = (PotType)param_globals::lats[i].measurand;
    acts[j].ID        = param_globals::lats[i].ID;
    acts[j].fout      = NULL;
 
    if(param_globals::lats[i].all)  {
      acts[j].fname = (char*) malloc((strlen(param_globals::lats[i].ID)+strlen(extn)+1)*sizeof(char));
      snprintf(acts[j].fname, strlen(param_globals::lats[i].ID)+strlen(extn)+1, "%s%s", param_globals::lats[i].ID, extn);
    }
    else  {
      char prfx[] = "init_acts_";
      int max_len = strlen(prfx) + strlen(param_globals::lats[i].ID) + strlen(extn) + 1;
 
      acts[j].fname = (char*) malloc(max_len*sizeof(char));
      snprintf(acts[j].fname, max_len, "%s%s%s", prfx, param_globals::lats[i].ID, extn);
    }
 
    // restarting
    if(prvSimDir != NULL)  {
      int len_fname = strlen(prvSimDir)+strlen(acts[j].fname)+2;
      acts[j].prv_fname = (char*) malloc(len_fname*sizeof(char));
      snprintf(acts[j].prv_fname, len_fname, "%s/%s", prvSimDir, acts[j].fname);
    }
 
    j++;
  }
 
  if(param_globals::compute_APD) {
    acts[j].method    = ACT_THRESH;    // threshold crossing
    acts[j].threshold = param_globals::actthresh;
    acts[j].mode      =   0;    // upstroke
    acts[j].all       = true;
    acts[j].measurand = VM;    // Vm
    //acts[j].ID      = dupstr("Vm_Activation");
    acts[j].fout      = NULL;
    acts[j].fname     = dupstr("vm_activation.dat");
 
    j++;
    acts[j].method    = ACT_THRESH;       // threshold crossing
    acts[j].threshold = param_globals::recovery_thresh;
    acts[j].mode      = 1;                // repol
    acts[j].all       = true;
    acts[j].measurand = VM;               // Vm
    //(*acts)[j+1].ID        = param_globals::lats[i].ID;
    acts[j].fout      = NULL;
    acts[j].fname     = dupstr("vm_repolarisation.dat");
 
    j++;
  }
 
  // set up sentinel for activity checking
  sntl.activated = have_sentinel;
  sntl.t_start   = param_globals::t_sentinel_start;
  sntl.t_window  = param_globals::t_sentinel;
  sntl.t_quiesc  =-1.;
  sntl.ID        = param_globals::sentinel_ID;
 
  if(need_to_add_sentinel) {
    // add a default LAT detector as sentinel
    acts[j].method    = ACT_THRESH;   // threshold crossing
    acts[j].threshold = param_globals::actthresh;
    acts[j].mode      =   0;          // upstroke
    acts[j].all       = true;
    acts[j].measurand = VM;           // Vm
    //(*acts)[j].ID      = dupstr("Vm_Activation");
    acts[j].fout      = NULL;
    acts[j].fname = dupstr("vm_sentinel.dat");
    // set sentinel index
    sntl.ID = j;
    j++;
  }
 
  if(prvSimDir) free(prvSimDir);
}
 
void print_act_log(FILE_SPEC logger, const SF::vector<Activation> & acts, int idx)
{
  const Activation & act = acts[idx];
 
  log_msg(logger, 0, 0, "\n");
  log_msg(logger, 0, 0, "LAT detector [%2d]", idx);
  log_msg(logger, 0, 0, "-----------------\n");
 
  log_msg(logger, 0, 0, "Measurand: %s", act.measurand ? "Phie" : "Vm");
  log_msg(logger, 0, 0, "All:       %s", act.all ? "All" : "Only first");
  log_msg(logger, 0, 0, "Method:    %s", act.method==ACT_DT ? "Derivative" : "Threshold crossing");
 
  char buf[64], gt[2], sgn[2];
  snprintf(sgn, sizeof sgn, "%s", act.mode?"-":"+");
  snprintf(gt, sizeof gt, "%s", act.mode?"<":">");
 
  if(act.method==ACT_DT)
    snprintf(buf, sizeof buf, "Maximum %sdf/dt %s %.2f mV", sgn, gt, act.threshold);
  else
    snprintf(buf, sizeof buf, "Intersection %sdf/dt with %.2f", sgn, act.threshold);
 
  log_msg(logger, 0, 0, "Mode:      %s", buf);
  log_msg(logger, 0, 0, "Threshold: %.2f mV\n", act.threshold);
}
 
void LAT_detector::init(sf_vec & vm, sf_vec & phie, int offset, enum physic_t phys_t)
{
  if(!get_physics(phys_t)) {
    log_msg(0,0,5, "There seems to be no EP is defined. LAT detector requires active EP! Aborting LAT setup!");
    return;
  }
 
  // we use the electrics logger for output
  FILE_SPEC logger = get_physics(phys_t)->logger;
 
  // TODO(init): except for the shallow copies, shouldn't these be deleted?
  // When to delete them?
  for(size_t i = 0; i < acts.size(); ++i) {
    acts[i].init   = 1;
    SF::init_vector(&(acts[i].phi));
    acts[i].phi->shallow_copy(!acts[i].measurand ? vm : phie);
    acts[i].offset = offset;
 
    SF::init_vector(&(acts[i].phip), acts[i].phi);
    *acts[i].phip = *acts[i].phi;
 
    // derivative based detector
    if (acts[i].method == ACT_DT) {
      SF::init_vector(&(acts[i].dvp0), acts[i].phi);
      SF::init_vector(&(acts[i].dvp1), acts[i].phi);
      if(acts[i].mode)
        log_msg(NULL,2,0, "Detection of -df/dt|max not implemented, +df/dt|max will be detected.");
    }
 
    // allocate additional local buffers
    acts[i].ibuf   = (int    *)malloc(acts[i].phi->lsize()*sizeof(int));
    acts[i].actbuf = (double *)malloc(acts[i].phi->lsize()*sizeof(double));
 
    if (!acts[i].all)  {
      SF::init_vector(&acts[i].tm, acts[i].phi->gsize(), acts[i].phi->lsize());
      acts[i].tm->set(-1.);
 
      // initialize with previous initial activations
      if(acts[i].prv_fname != NULL)  {
        set_dir(INPUT);
        size_t nread = acts[i].tm->read_ascii(acts[i].prv_fname);
 
        if(nread == 0)
          log_msg(NULL,2,ECHO,"Warning: Initialization of LAT[%2d] failed.", i);
 
        set_dir(OUTPUT);
      }
    }
    else {
      if ( !get_rank() ) {
        // here we should copy over previous file and open in append mode
        if(acts[i].prv_fname!=NULL)  {
          set_dir(INPUT);
          FILE_SPEC in = f_open( acts[i].prv_fname, "r" );
          if(in) {
            log_msg(NULL,2,0, "Copying over of previous activation file not implemented.\n"); f_close(in);
          }
          else
            log_msg(NULL,3,0,"Warning: Initialization in %s - \n"
                             "Failed to read activation file %s.\n", __func__, acts[i].prv_fname);
 
          set_dir(OUTPUT);
        }
        acts[i].fout = f_open( acts[i].fname, acts[i].prv_fname==NULL?"w":"a" );
      }
    }
    print_act_log(logger, acts, i);
  }
 
  sf_mesh & intra_mesh = get_mesh(intra_elec_msh);
  SF::local_petsc_to_nodal_mapping(intra_mesh, this->petsc_to_nodal);
}
 
 
int output_all_activations(FILE_SPEC fp, int *ibuf, double *act_tbuf, int nlacts)
{
  int rank = get_rank(), gacts = 0, numProc = get_size();
 
  if (rank == 0) {
    // rank 0 writes directly to the table
    for (int i=0; i<nlacts; i++)
      fprintf(fp->fd, "%d\t%.6f\n", ibuf[i], act_tbuf[i]);
 
    gacts += nlacts;
 
    SF::vector<int>    buf_inds;
    SF::vector<double> buf_acts;
 
    for (int j=1; j<numProc; j++) {
      int acts = 0;
      MPI_Status status;
      MPI_Recv(&acts, 1, MPI_INT, j, 110, PETSC_COMM_WORLD, &status);
 
      if (acts) {
        buf_inds.resize(acts);
        buf_acts.resize(acts);
 
        MPI_Recv(buf_inds.data(), acts, MPI_INT,    j, 110, PETSC_COMM_WORLD, &status);
        MPI_Recv(buf_acts.data(), acts, MPI_DOUBLE, j, 110, PETSC_COMM_WORLD, &status);
 
        for(int ii=0; ii<acts; ii++)
          fprintf(fp->fd, "%d\t%.6f\n", buf_inds[ii], buf_acts[ii]);
 
        gacts += acts;
      }
    }
    fflush(fp->fd);
  }
  else {
    MPI_Send(&nlacts, 1, MPI_INT, 0, 110, PETSC_COMM_WORLD);
    if (nlacts) {
      MPI_Send(ibuf,     nlacts, MPI_INT,    0, 110, PETSC_COMM_WORLD);
      MPI_Send(act_tbuf, nlacts, MPI_DOUBLE, 0, 110, PETSC_COMM_WORLD);
    }
  }
 
  MPI_Bcast(&gacts, 1, MPI_INT, 0, PETSC_COMM_WORLD);
  return gacts;
}
 
int LAT_detector::check_acts(double tm)
{
  int nacts =  0;
  double *a;
 
  for(Activation* aptr = acts.data(); aptr != acts.end(); aptr++)
  {
    int lacts = 0;
    switch (aptr->method) {
        case ACT_THRESH:
            lacts = check_cross_threshold(*aptr->phi, *aptr->phip, tm,
                    aptr->ibuf, aptr->actbuf, aptr->threshold, aptr->mode);
            break;
 
        case ACT_DT:
            lacts = check_mx_derivative  (*aptr->phi, *aptr->phip, tm,
                    aptr->ibuf, aptr->actbuf, *aptr->dvp0, *aptr->dvp1,
                    aptr->threshold, aptr->mode);
            break;
 
        default:
            break;
    }
 
    if (!aptr->all)
      a = aptr->tm->ptr();
 
    const SF::vector<mesh_int_t> & canon_nbr = get_mesh(intra_elec_msh).get_numbering(SF::NBR_SUBMESH);
 
    for(int j=0; j<lacts; j++) {
      if(aptr->all) {
        int nodal_idx = this->petsc_to_nodal.forward_map(aptr->ibuf[j]);
        aptr->ibuf[j] = canon_nbr[nodal_idx] + aptr->offset;
      }
      else {
        if(a[aptr->ibuf[j]] == -1)
          a[aptr->ibuf[j]] =  aptr->actbuf[j];
      }
    }
 
    if(aptr->all)
      output_all_activations(aptr->fout, aptr->ibuf, aptr->actbuf, lacts);
    else
      aptr->tm->release_ptr(a);
 
    MPI_Allreduce(MPI_IN_PLACE, &lacts, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD);
    nacts += lacts;
 
    aptr->nacts = nacts;
  }
 
  return nacts > 0;
}
 
 
int LAT_detector::check_quiescence(double tm, double dt)
{
  static int savequitFlag =  0;
  int numNodesActivated   = -1;
 
  if(sntl.activated) {
    // initialization
    if(sntl.t_quiesc < 0. && sntl.t_window >= 0.0 ) {
      log_msg(0,0,ECHO | NONL, "================================================================================================\n");
      log_msg(0,0,ECHO | NONL, "%s() WARNING: simulation is configured to savequit() after %.2f ms of quiescence\n", __func__, sntl.t_window);
      log_msg(0,0,ECHO | NONL, "================================================================================================\n");
      sntl.t_quiesc = 0.0;
    }
 
    if(tm >= sntl.t_start && !savequitFlag)
    {
      numNodesActivated = acts[sntl.ID].nacts;
 
      if(numNodesActivated) sntl.t_quiesc = 0.0;
      else                  sntl.t_quiesc += dt;
 
      if(sntl.t_window >= 0.0 && sntl.t_quiesc > sntl.t_window && !savequitFlag) {
        savequitFlag++;
        savequit();
      }
    }
  }
 
  return numNodesActivated;
}
 
 
 
 
int LAT_detector::check_cross_threshold(sf_vec & vm, sf_vec & vmp, double tm,
                          int *ibuf, double *actbuf, float threshold, int mode)
{
  SF_real *c = vm.ptr();
  SF_real *p = vmp.ptr();
  int lsize = vm.lsize();
  int nacts = 0, gnacts = 0;
 
  for (int i=0; i<lsize; i++) {
    int sgn        = 1;
    bool triggered = false;
    if(mode==0)  {// detect +slope crossing
      triggered = p[i] <= threshold && c[i] > threshold; }
    else  {      // detect -slope crossing
      triggered = p[i] >= threshold && c[i] < threshold;
      sgn = -1;
    }
 
    if (triggered)  {
      double tact = tm - param_globals::dt + (threshold-p[i])/(c[i]-p[i])*sgn*param_globals::dt;
      ibuf  [nacts] = i;
      actbuf[nacts] = tact;
      nacts++;
    }
    p[i] = c[i];
  }
 
  vm.release_ptr(c);
  vmp.release_ptr(p);
  return nacts;
}
 
int LAT_detector::check_mx_derivative(sf_vec & vm, sf_vec & vmp, double tm,
                        int *ibuf, double *actbuf, sf_vec & dvp0, sf_vec & dvp1,
                        float threshold, int mode)
{
  int  nacts = 0, gnacts = 0;
  double   tact, dt2 = 2 * param_globals::dt;
  int  lsize = vm.lsize();
  SF_real  ddv0, ddv1, dv, dvdt;
  SF_real *c, *p, *pd0, *pd1;
 
  c   = vm.ptr();
  p   = vmp.ptr();
  pd0 = dvp0.ptr();
  pd1 = dvp1.ptr();
 
  for (int i=0; i<lsize; i++ )  {
    dv   = (c[i]-p[i]);
    dvdt = dv/param_globals::dt;
    ddv0 = pd1[i]-pd0[i];
    ddv1 = dv    -pd1[i];
 
    if (dvdt>=threshold && ddv0>0 && ddv1<0)  {
      tact = tm-dt2+(ddv0/(ddv0-ddv1))*param_globals::dt;
      ibuf  [nacts] = i;
      actbuf[nacts] = tact;
      nacts++;
    }
    p  [i] = c[i];
    pd0[i] = pd1[i];
    pd1[i] = dv;
  }
 
  vm  .release_ptr(c);
  vmp .release_ptr(p);
  dvp0.release_ptr(pd0);
  dvp1.release_ptr(pd1);
 
  return nacts;
}
 
void LAT_detector::output_initial_activations()
{
  SF::scattering * sc = get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, 1);
  assert(sc != NULL);
 
  bool forward = true;
 
  for (size_t i = 0; i < acts.size(); i++)  {
    if (is_init(acts[i].tm)) {
      (*sc)(*acts[i].tm, forward);
      acts[i].tm->write_ascii(acts[i].fname, false);
    }
  }
}
 
void Electrics::prepace() {
  log_msg(NULL, 0, 0, "Using activation times from file %s to distribute prepacing states\n",
          param_globals::prepacing_lats);
  log_msg(NULL, 0, 0, "Assuming stimulus strength %f uA/uF with duration %f ms for prepacing\n",
          param_globals::prepacing_stimstr, param_globals::prepacing_stimdur);
 
  limpet::MULTI_IF* miif = this->ion.miif;
 
  const sf_mesh &   mesh = get_mesh(intra_elec_msh);
  sf_vec* read_lats; SF::init_vector(&read_lats, mesh, 1, sf_vec::algebraic);
 
  // read in the global distributed vector of all activation times
  set_dir(INPUT);
  size_t numread = read_lats->read_ascii(param_globals::prepacing_lats);
  if (numread == 0) {
    log_msg(NULL, 5, 0, "Failed reading required LATs! Skipping prepacing!");
    return;
  }
  set_dir(OUTPUT);
 
  SF::scattering* sc = get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, 1);
  assert(sc != NULL);
 
  // permute in-place to petsc permutation
  bool forward = false;
  (*sc)(*read_lats, forward);
 
  // take care of negative LAT values
  {
    PetscReal* lp = read_lats->ptr();
    for(int i=0; i<read_lats->lsize(); i++)
      if(lp[i] < 0.0) lp[i] = param_globals::tend + 10.0;
 
    read_lats->release_ptr(lp);
  }
 
  // make LATs relative and figure out the first LAT
  // so we know when to save state of each point
  SF_real LATmin = read_lats->min();
 
  if(LATmin < 0.0) {
    log_msg(0,3,0, "LAT data is not complete. Skipping prepacing.");
    return;
  }
 
  SF_real offset  = floor(LATmin / param_globals::prepacing_bcl) * param_globals::prepacing_bcl;
  SF_real last_tm = param_globals::prepacing_bcl * param_globals::prepacing_beats;
 
  // compute read_lats[i] = last_tm - (read_lats[i] - offset)
  *read_lats += -offset;
  *read_lats *= -1.;
  *read_lats += last_tm;
 
  miif->getRealData();
  SF_real *save_tm = read_lats->ptr();
  SF_real *vm      = miif->gdata[limpet::Vm]->ptr();
 
  for (int ii = 0; ii < miif->N_IIF; ii++) {
    if (!miif->N_Nodes[ii]) continue;
 
    // create sorted array of save times.
    SF::vector<SF::mixed_tuple<double,int>> sorted_save(miif->N_Nodes[ii]); // v1 = time, v2 = index
    for (int kk = 0; kk < miif->N_Nodes[ii]; kk++) {
      sorted_save[kk].v1 = save_tm[miif->NodeLists[ii][kk]];
      sorted_save[kk].v2 = kk;
    }
    std::sort(sorted_save.begin(), sorted_save.end());
 
    size_t lastidx = sorted_save.size() - 1;
    int    paced   = sorted_save[lastidx].v2;  // IMP index of latest node
    int    csav    = 0;
 
    for (double t = 0; t < sorted_save[lastidx].v1; t += param_globals::dt) {
      if (fmod(t, param_globals::prepacing_bcl) < param_globals::prepacing_stimdur &&
          t < param_globals::prepacing_bcl * param_globals::prepacing_beats - 1)
        miif->ldata[ii][limpet::Vm][paced] += param_globals::prepacing_stimstr * param_globals::dt;
 
      compute_IIF(*miif->IIF[ii], miif->ldata[ii], paced);
 
      // Vm update always happens now outside of the imp
      miif->ldata[ii][limpet::Vm][paced]    -= miif->ldata[ii][limpet::Iion][paced] * param_globals::dt;
      vm[miif->NodeLists[ii][paced]] = miif->ldata[ii][limpet::Vm][paced];
 
      while (csav < miif->N_Nodes[ii] - 1 && t >= sorted_save[csav].v1)
        dup_IMP_node_state(*miif->IIF[ii], paced, sorted_save[csav++].v2, miif->ldata[ii]);
    }
 
    // get nodes which may be tied for last
    while (csav < miif->N_Nodes[ii] - 1)
      dup_IMP_node_state(*miif->IIF[ii], paced, sorted_save[csav++].v2, miif->ldata[ii]);
    // ipdate global Vm vector
    for (int k = 0; k < miif->N_Nodes[ii]; k++) vm[miif->NodeLists[ii][k]] = miif->ldata[ii][limpet::Vm][k];
  }
 
  read_lats->release_ptr(save_tm);
  miif->gdata[limpet::Vm]->release_ptr(vm);
  miif->releaseRealData();
}
 
 
void recover_phie_std(sf_vec & vm, phie_recovery_data & rcv)
{
  if (!rcv.pts.size())
    return;
 
  int rank = get_rank();
 
  if(!get_physics(elec_phys)) {
    log_msg(0,0,5, "There seems to be no EP is defined. Phie recovery requires active EP! Aborting!");
    return;
  }
 
  Electrics* elec = static_cast<Electrics*>(get_physics(elec_phys));
  sf_mat & Ki = *elec->parab_solver.rhs_parab;
 
  const sf_mesh &                imesh   = get_mesh(intra_elec_msh);
  const SF::vector<mesh_int_t> & alg_nod = imesh.pl.algebraic_nodes();
 
  SF_int start, end;
  vm.get_ownership_range(start, end);
 
  if(!rcv.Im) {
    SF::init_vector(&rcv.Im, &vm);
    SF::init_vector(&rcv.dphi, &vm);
  }
 
  SF_int r_start, r_end;
  rcv.phie_rec->get_ownership_range(r_start, r_end);
 
  SF_real *ph_r = rcv.phie_rec->ptr();
 
  // use minimum distance to ensure r>0
  // consistent with the line source approximation, the "cable radius"
  // is used as a lower limit for the source-field point distance
  float minDist = 2. / param_globals::imp_region[0].cellSurfVolRatio;  // radius in um
 
  Ki.mult(vm, *rcv.Im);
  int numpts = rcv.pts.size() / 3;
  Point fpt, cpt;
 
  for (int j=0; j<numpts; j++) {
    fpt = rcv.pts.data() + j*3;
 
    *rcv.dphi = *rcv.Im;
    SF_real* dp = rcv.dphi->ptr();
 
    for (size_t i = 0; i<alg_nod.size(); i++)
    {
      mesh_int_t loc_nodal_idx = alg_nod[i];
      mesh_int_t loc_petsc_idx = local_nodal_to_local_petsc(imesh, rank, loc_nodal_idx);
      cpt = imesh.xyz.data()+loc_nodal_idx*3;
 
      double r = dist(fpt, cpt) + minDist;
      dp[loc_petsc_idx] /= r;
    }
 
    rcv.dphi->release_ptr(dp);
 
    SF_real phi = rcv.dphi->sum() / 4. / M_PI / rcv.gBath;
    if ( (j>=r_start) && (j<r_end) )
      ph_r[j-r_start] = phi;
  }
 
  rcv.phie_rec->release_ptr(ph_r);
}
 
int postproc_recover_phie()
{
  int err = 0, rank = get_rank();
 
  if(!get_physics(elec_phys) || !mesh_is_registered(intra_elec_msh)) {
    log_msg(0,0,5, "There seems to be no EP is defined. Phie recovery requires active EP! Aborting!");
    return 1;
  }
 
  sf_mesh & imesh = get_mesh(intra_elec_msh);
  Electrics* elec = static_cast<Electrics*>(get_physics(elec_phys));
  phie_recovery_data & phie_rcv = elec->phie_rcv;
 
  // we close the files of the default electrics if there are any open
  elec->output_manager.close_files_and_cleanup();
 
  // register output
  set_dir(POSTPROC);
  igb_output_manager phie_rec_out;
  phie_rec_out.register_output(phie_rcv.phie_rec, phie_recv_msh, 1,
         param_globals::phie_recovery_file, "mV");
 
  // Buffer for Vm data
  sf_vec* vm = get_data(vm_vec); assert(vm);
 
  // set up igb header and point fd to start of Vm file
  set_dir(OUTPUT);
  IGBheader vm_igb;
  if(rank == 0) {
    FILE_SPEC file = f_open(param_globals::vofile, "r");
    if(file != NULL) {
      vm_igb.fileptr(file->fd);
      vm_igb.read();
 
      if(vm_igb.x() != vm->gsize()) {
        log_msg(0,4,0, "%s error: Vm dimension does not fit to %s file. Aborting recovery! \n",
                __func__, param_globals::vofile);
        err++;
      }
 
      delete file;
    }
    else err++;
  }
 
  err = get_global(err, MPI_MAX);
 
  if(err == 0) {
    FILE* fd = static_cast<FILE*>(vm_igb.fileptr());
 
    // number of data slices
    const int num_io = user_globals::tm_manager->timers[iotm_spacedt]->numIOs;
 
    // scatterers
    SF::scattering* petsc_to_canonical = get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, 1);
    assert(petsc_to_canonical != NULL);
 
    // loop over vm slices and recover phie
    for(int i=0; i<num_io; i++) {
      log_msg(0,0,0, "Step %d / %d", i+1, num_io);
      size_t nread = vm->read_binary<float>(fd);
 
      if(nread != size_t(vm->gsize())) {
        log_msg(0,3,0, "%s warning: read incomplete data slice! Aborting!", __func__);
        err++;
        break;
      }
 
      // permute vm_buff
      bool forward = false;
      (*petsc_to_canonical)(*vm, forward);
 
      // do phie computation
      recover_phie_std(*vm, phie_rcv);
 
      phie_rec_out.write_data();
    }
 
    phie_rec_out.close_files_and_cleanup();
  }
  return err;
}
 
void setup_phie_recovery_data(phie_recovery_data & data)
{
  if(!get_physics(elec_phys) ) {
    log_msg(0,0,5, "There seems to be no EP is defined. Phie recovery requires active EP! Aborting!");
    return;
  }
 
  int rank = get_rank(), size = get_size();
  Electrics* elec = static_cast<Electrics*>(get_physics(elec_phys));
 
  sf_mesh & imesh = get_mesh(intra_elec_msh);
  const std::string basename = param_globals::phie_rec_ptf;
  SF::vector<mesh_int_t>  ptsidx;
 
  set_dir(INPUT);
  SF::read_points(basename, imesh.comm, data.pts, ptsidx);
  make_global(data.pts, imesh.comm); // we want all ranks to have all points
 
  // set up parallel layout of recovery points
  SF::vector<mesh_int_t> layout;
  layout_from_count(mesh_int_t(ptsidx.size()), layout, imesh.comm);
 
  // set up petsc_vector for recovered potentials
  SF::init_vector(&data.phie_rec, layout[size], layout[rank+1]-layout[rank], 1, sf_vec::algebraic);
 
  // get conductivty
  SF::vector<RegionSpecs> & intra_regions = elec->mtype[Electrics::intra_grid].regions;
  data.gBath = static_cast<elecMaterial*>(intra_regions[0].material)->BathVal[0];
}
 
void Laplace::initialize()
{
  int rank = get_rank();
 
  assert(param_globals::bidomain == BIDOMAIN);
  double t1, t2;
  get_time(t1);
 
  // set up Extracellular tissue
  set_elec_tissue_properties(mtype, Electrics::extra_grid, NULL);
  region_mask(extra_elec_msh, mtype[Electrics::extra_grid].regions,
              mtype[Electrics::extra_grid].regionIDs, true, "gregion_e");
 
  // set up a subset of the complete electrical mappings
  int dpn = 1;
  register_permutation(extra_elec_msh, PETSC_TO_CANONICAL, dpn);
 
  if(mesh_is_registered(intra_elec_msh)) {
    // set up Intracellular tissue
    set_elec_tissue_properties(mtype, Electrics::intra_grid, NULL);
    region_mask(intra_elec_msh, mtype[Electrics::intra_grid].regions,
                mtype[Electrics::intra_grid].regionIDs, true, "gregion_i");
 
    register_permutation(intra_elec_msh, PETSC_TO_CANONICAL, dpn);
    register_scattering(intra_elec_msh, extra_elec_msh, dpn);
  }
 
  // set up stimuli
  init_stim_info();
  stimuli.resize(param_globals::num_stim);
 
  for(int i=0; i<param_globals::num_stim; i++) {
    // construct new stimulus
    stimulus & s = stimuli[i];
 
    if(user_globals::using_legacy_stimuli)
      s.translate(i);
 
    s.setup(i);
 
    if(s.phys.type == Phi_ex) {
      s.pulse.wform = constPulse;
      sample_wave_form(s.pulse, i);
    }
  }
 
  set_dir(OUTPUT);
 
  ellip_solver.init();
  ellip_solver.rebuild_stiffness(mtype, stimuli, logger);
 
  if(param_globals::dump2MatLab) {
    std::string bsname = param_globals::dump_basename;
    std::string fn;
 
    set_dir(OUTPUT);
    fn = bsname + "_Kie.bin";
    ellip_solver.phie_mat->write(fn.c_str());
  }
 
  // the laplace solver executes only once, thus we need a singlestep timer
  timer_idx = user_globals::tm_manager->add_singlestep_timer(0.0, 0.0, "laplace trigger", nullptr);
 
  SF::vector<mesh_int_t>* restr_i = NULL;
  SF::vector<mesh_int_t>* restr_e = NULL;
 
  setup_dataout(param_globals::dataout_e, param_globals::dataout_e_vtx, extra_elec_msh,
                restr_e, param_globals::num_io_nodes > 0);
  if(param_globals::dataout_e)
    output_manager.register_output(ellip_solver.phie, extra_elec_msh, 1, param_globals::phiefile, "mV", restr_e);
 
  if(mesh_is_registered(intra_elec_msh)) {
    setup_dataout(param_globals::dataout_i, param_globals::dataout_i_vtx, intra_elec_msh,
                  restr_i, param_globals::num_io_nodes > 0);
    if(param_globals::dataout_i)
      output_manager.register_output(ellip_solver.phie_i, intra_elec_msh, 1, param_globals::phieifile, "mV", restr_i);
  }
 
  this->initialize_time += timing(t2, t1);
 
  this->compute_step();
}
 
void Laplace::destroy()
{}
 
void Laplace::compute_step()
{
  // Laplace compute might be called multiple times, we want to run only once..
  if(!ellip_solver.lin_solver) return;
 
  double t0, t1, dur;
  log_msg(0,0,0, "Solving Laplace problem ..");
 
  get_time(t0);
  ellip_solver.solve_laplace();
  dur = timing(t1,t0);
 
  log_msg(0,0,0, "Done in %.5f seconds.", dur);
 
  ellip_solver.stats.log_stats(user_globals::tm_manager->time, false);
  this->compute_time += timing(t1, t0);
  set_dir(OUTPUT);
  output_manager.write_data();
  output_manager.close_files_and_cleanup();
 
  // we clear the elliptic matrices and solver to save some memory when computing
  // the laplace solution on-the-fly
  delete ellip_solver.mass_e;     ellip_solver.mass_e     = NULL;
  delete ellip_solver.phie_mat;   ellip_solver.phie_mat   = NULL;
  delete ellip_solver.lin_solver; ellip_solver.lin_solver = NULL;
}
 
void Laplace::output_step()
{}
 
double Laplace::timer_val(const int timer_id)
{
  int sidx = stimidx_from_timeridx(stimuli, timer_id);
  double val = 0.0;
 
  if(sidx != -1) stimuli[sidx].value(val);
  else           val = std::nan("NaN");
  return val;
}
 
std::string Laplace::timer_unit(const int timer_id)
{
  int sidx = stimidx_from_timeridx(stimuli, timer_id);
  std::string s_unit;
  if(sidx != -1) s_unit = stimuli[sidx].pulse.wave.f_unit;
  return s_unit;
}
 
void constant_total_stimulus_current(SF::vector<stimulus> & stimuli,
                                     sf_mat & mass_i,
                                     sf_mat & mass_e,
                                     limpet::MULTI_IF *miif,
                                     FILE_SPEC logger)
{
  hashmap::unordered_map<mesh_int_t,mesh_int_t> alg_idx_map;
 
  for(stimulus & s : stimuli) {
    if(is_current(s.phys.type) && s.phys.total_current) {
      // extracellular current injection
      if (s.phys.type == I_ex) {
        // compute affected volume in um^3
        SF_real vol = get_volume_from_nodes(mass_e, s.electrode.vertices);
 
        // s->strength holds the total current in uA, compute current density in uA/cm^3
        // Theoretically, we don't need to scale the volume to cm^3 here since we later 
        // multiply with the mass matrix and we get um^3 * uA/um^3 = uA.
        // However, for I_ex there is an additional um^3 to cm^3 scaling in phys.scale,
        // since I_e is expected to be in uA/cm^3. Therefore, we need to compensate for that to arrive at uA later.
        float scale = 1.e12/vol;
 
        s.pulse.strength *= scale;
 
        log_msg(logger,0,ECHO,
                "%s [Stimulus %d]: current density scaled to %.4g uA/cm^3\n",
                s.name.c_str(), s.idx, s.pulse.strength);
      }
      else if (s.phys.type == I_tm) {
        // compute affected volume in um^3
        SF_real vol = get_volume_from_nodes(mass_i, s.electrode.vertices);
        const sf_mesh & imesh = get_mesh(intra_elec_msh);
        const SF::vector<mesh_int_t> & alg_nod = imesh.pl.algebraic_nodes();
 
        if(alg_idx_map.size() == 0) {
          mesh_int_t lidx = 0;
          for(mesh_int_t n : alg_nod) {
            alg_idx_map[n] = lidx;
            lidx++;
          }
        }
 
        SF_real surf = 0.0;
        for(mesh_int_t n : s.electrode.vertices) {
          if(alg_idx_map.count(n)) {
            mesh_int_t lidx = alg_idx_map[n];
            int r = miif->IIFmask[lidx];
            // surf = vol*beta [1/um], surf is in [um^2]
            surf = vol * miif->IIF[r]->cgeom().SVratio * param_globals::imp_region[r].volFrac;
            //convert to cm^2
            surf /= 1.e8;
            break;
          }
        }
        surf = get_global(surf, MPI_MAX, PETSC_COMM_WORLD);
 
        // scale surface density now to result in correct total current
        s.pulse.strength /= surf;
        log_msg(logger, 0, ECHO,
                "%s [Stimulus %d]: current density scaled to %.4g uA/cm^2\n",
                s.name.c_str(), s.idx, s.pulse.strength);
      }
    }
  }
}
 
 
 
}  // namespace opencarp