emi.cc

Go to the documentation of this file.

//
// 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
// ----------------------------------------------------------------------------
 
#if WITH_EMI_MODEL
 
#include "petsc_utils.h"
#include "mpi_utils.h"
#include "physics_types.h"
#include "timers.h"
#include "stimulate.h"
#include "electric_integrators.h"
#include "SF_init.h" // for SF::init_xxx()
#include "emi.h"
#include <array>
#include <cmath>
#include <cstdint>
#include <initializer_list>
#include <sys/resource.h>
#include <utility>
 
#ifdef WITH_CALIPER
#include "caliper/cali.h"
#else
#include "caliper_hooks.h"
#endif
 
namespace opencarp {
 
namespace {
 
template<class Assemble>
void assemble_with_exact_preallocation(std::initializer_list<sf_mat*> matrices, Assemble assemble)
{
  bool any_supported = false;
  bool all_supported = true;
  bool saw_matrix = false;
 
  for(sf_mat* mat : matrices) {
    if(mat == nullptr) continue;
 
    saw_matrix = true;
    const bool supported = mat->begin_exact_preallocation();
    any_supported = any_supported || supported;
    all_supported = all_supported && supported;
  }
 
  if(!saw_matrix) return;
 
  // PETSc exact preallocation is a backend-wide mode: either all matrices in
  // this assembly group support it, or none of them do. Mixed support would
  // indicate an inconsistent backend state.
  assert(any_supported == all_supported);
 
  if(all_supported) {
    // First pass collects the sparse graph for PETSc; the second pass below
    // inserts the real numerical values into the exactly preallocated matrix.
    assemble();
    for(sf_mat* mat : matrices) {
      if(mat != nullptr) mat->finalize_exact_preallocation();
    }
    assemble();
    return;
  }
 
  assemble();
}
 
void log_emi_petsc_matrix_preallocation_report(std::initializer_list<std::pair<const char*, sf_mat*>> matrices)
{
#ifdef WITH_PETSC
  PetscBool enabled = PETSC_FALSE;
  PetscOptionsHasName(NULL, NULL, "-mat_view_info", &enabled);
  if(!enabled) return;
 
  // This report estimates PETSc sparse-matrix storage from MatGetInfo().
  // It is not the total process peak memory. Peak RSS includes meshes,
  // vectors, solvers, MPI buffers, temporary assembly data, and allocator
  // overhead, so it should be measured externally, e.g. with
  // mprof run --include-children followed by mprof peak.
  int rank = 0;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  if(rank == 0) {
    log_msg(NULL, 0, 0, "\nEMI PETSc matrix preallocation report");
    log_msg(NULL, 0, 0, "matrix             alloc/used       used_est       allocated_est    overalloc_est");
  }
 
  constexpr double bytes_per_nz = double(sizeof(PetscScalar) + sizeof(PetscInt));
  constexpr double gib = 1024.0 * 1024.0 * 1024.0;
  double total_used_gib = 0.0;
  double total_allocated_gib = 0.0;
 
  for(const auto& item : matrices) {
    auto* petsc_mat = dynamic_cast<SF::petsc_matrix*>(item.second);
    if(petsc_mat == nullptr || petsc_mat->data == nullptr) continue;
 
    MatInfo info;
    MatGetInfo(petsc_mat->data, MAT_GLOBAL_SUM, &info);
 
    const double used = static_cast<double>(info.nz_used);
    const double allocated = static_cast<double>(info.nz_allocated);
    const double ratio = used > 0.0 ? allocated / used : 0.0;
    const double used_gib = used * bytes_per_nz / gib;
    const double allocated_gib = allocated * bytes_per_nz / gib;
    const double overallocated_gib = allocated_gib - used_gib;
    total_used_gib += used_gib;
    total_allocated_gib += allocated_gib;
 
    if(rank == 0) {
      log_msg(NULL, 0, 0,
              "%-18s alloc/used=%5.2f   used_est=%6.2fG allocated_est=%6.2fG overalloc_est=%6.2fG",
              item.first, ratio, used_gib, allocated_gib, overallocated_gib);
    }
  }
 
  if(rank == 0) {
    const double total_ratio = total_used_gib > 0.0 ? total_allocated_gib / total_used_gib : 0.0;
    log_msg(NULL, 0, 0,
            "%-18s alloc/used=%5.2f   used_est=%6.2fG allocated_est=%6.2fG overalloc_est=%6.2fG",
            "TOTAL:", total_ratio, total_used_gib, total_allocated_gib,
            total_allocated_gib - total_used_gib);
  }
 
  struct rusage usage;
  getrusage(RUSAGE_SELF, &usage);
#ifdef __APPLE__
  const double local_peak_rss_gib = double(usage.ru_maxrss) / gib;
#else
  const double local_peak_rss_gib = double(usage.ru_maxrss) * 1024.0 / gib;
#endif
  double summed_peak_rss_gib = 0.0;
  double max_rank_peak_rss_gib = 0.0;
  MPI_Reduce(&local_peak_rss_gib, &summed_peak_rss_gib, 1, MPI_DOUBLE, MPI_SUM, 0, PETSC_COMM_WORLD);
  MPI_Reduce(&local_peak_rss_gib, &max_rank_peak_rss_gib, 1, MPI_DOUBLE, MPI_MAX, 0, PETSC_COMM_WORLD);
 
  if(rank == 0) {
    log_msg(NULL, 0, 0,
            "process peak RSS estimate: summed ranks=%6.2fG max rank=%6.2fG",
            summed_peak_rss_gib, max_rank_peak_rss_gib);
    log_msg(NULL, 0, 0,
            "external peak memory from mprof --include-children is still the recommended whole-run reference.\n");
  }
#else
  (void)matrices;
#endif
}
 
} // namespace
 
void log_mesh_local_element_ranges(const sf_mesh& emi_mesh,
                                   const sf_mesh& emi_surfmesh_w_counter_face,
                                   const sf_mesh& emi_surfmesh_unique_face)
{
  int rank = 0;
  int comm_size = 0;
  MPI_Comm_rank(emi_surfmesh_w_counter_face.comm, &rank);
  MPI_Comm_size(emi_surfmesh_w_counter_face.comm, &comm_size);
 
  const size_t local_emi_elems = emi_mesh.l_numelem;
  const size_t local_both_face_elems = emi_surfmesh_w_counter_face.l_numelem;
  const size_t local_unique_face_elems = emi_surfmesh_unique_face.l_numelem;
 
  std::vector<size_t> all_emi_elems;
  std::vector<size_t> all_both_face_elems;
  std::vector<size_t> all_unique_face_elems;
  if (rank == 0) {
    all_emi_elems.resize(comm_size, 0);
    all_both_face_elems.resize(comm_size, 0);
    all_unique_face_elems.resize(comm_size, 0);
  }
 
  const MPI_Datatype size_mpi_t = mpi_datatype<size_t>();
  MPI_Gather(&local_emi_elems, 1, size_mpi_t,
             rank == 0 ? all_emi_elems.data() : nullptr, 1, size_mpi_t,
             0, emi_surfmesh_w_counter_face.comm);
  MPI_Gather(&local_both_face_elems, 1, size_mpi_t,
             rank == 0 ? all_both_face_elems.data() : nullptr, 1, size_mpi_t,
             0, emi_surfmesh_w_counter_face.comm);
  MPI_Gather(&local_unique_face_elems, 1, size_mpi_t,
             rank == 0 ? all_unique_face_elems.data() : nullptr, 1, size_mpi_t,
             0, emi_surfmesh_w_counter_face.comm);
 
  if (rank != 0) return;
 
  const auto print_min_max = [](const char* label, const std::vector<size_t>& counts) {
    if (counts.empty()) return;
 
    size_t min_val = counts[0];
    size_t max_val = counts[0];
    int min_rank = 0;
    int max_rank = 0;
 
    for (int r = 1; r < static_cast<int>(counts.size()); ++r) {
      if (counts[r] < min_val) {
        min_val = counts[r];
        min_rank = r;
      }
      if (counts[r] > max_val) {
        max_val = counts[r];
        max_rank = r;
      }
    }
 
    log_msg(NULL, 0, 0, "  %s: \n\t\t min=%zu on rank=%d, \n\t\t max=%zu on rank=%d\n",
            label, min_val, min_rank, max_val, max_rank);
  };
  log_msg(NULL, 0, 0, "\n**********************************");
  log_msg(NULL, 0, 0, "min/max number of local-element ranges:");
  print_min_max("emi_mesh", all_emi_elems);
  print_min_max("emi_surfmesh_w_counter_face", all_both_face_elems);
  print_min_max("emi_surfmesh_unique_face", all_unique_face_elems);
  log_msg(NULL, 0, 0, "**********************************");
}
 
#ifdef EMI_DEBUG_MESH
void log_lhs_positive_definite_probe(SF::abstract_matrix<SF_int, SF_real>* mat,
                                     FILE_SPEC logger,
                                     const char* stage,
                                     int num_trials = 5)
{
  (void)logger;
  if (param_globals::flavor != std::string("petsc")) return;
 
  auto* petsc_mat = dynamic_cast<SF::petsc_matrix*>(mat);
  if (petsc_mat == nullptr) return;
 
  Vec x = NULL, y = NULL;
  MatCreateVecs(petsc_mat->data, &x, &y);
 
  PetscRandom rnd = NULL;
  PetscRandomCreate(PETSC_COMM_WORLD, &rnd);
  PetscRandomSetFromOptions(rnd);
 
  int rank = 0;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 
  PetscReal min_q = PETSC_MAX_REAL;
  PetscReal max_q = -PETSC_MAX_REAL;
  PetscInt nonpos_count = 0;
 
  for (int i = 0; i < num_trials; ++i) {
    VecSetRandom(x, rnd);
    MatMult(petsc_mat->data, x, y);
 
    PetscScalar q_scalar = 0.0;
    VecDot(x, y, &q_scalar);
 
    const PetscReal q = PetscRealPart(q_scalar);
    min_q = std::min(min_q, q);
    max_q = std::max(max_q, q);
    if (q <= 0.0) nonpos_count++;
  }
 
  if (rank == 0) {
    PetscPrintf(PETSC_COMM_SELF,
                "%s: SPD probe with %d random vectors: min(x^T A x)=%g, max(x^T A x)=%g, nonpositive=%d\n",
                stage, num_trials, double(min_q), double(max_q), int(nonpos_count));
  }
 
  PetscRandomDestroy(&rnd);
  VecDestroy(&x);
  VecDestroy(&y);
}
 
PetscScalar emi_probe_vector_entry(const PetscInt gid, const int probe_id)
{
  const double x = static_cast<double>(gid + 1);
 
  switch (probe_id) {
    case 0:
      return std::sin(1.0e-3 * x) + 0.5 * std::cos(3.0e-3 * x);
    case 1:
      return std::cos(7.0e-4 * x) - 0.35 * std::sin(2.0e-3 * x);
    default:
      return 0.75 * std::sin(1.3e-3 * x) + 0.25 * std::cos(4.0e-3 * x);
  }
}
 
void log_lhs_operator_probe(SF::abstract_matrix<SF_int, SF_real>* mat,
                            FILE_SPEC logger,
                            const char* stage,
                            int num_probes = 3)
{
  (void)logger;
  if (param_globals::flavor != std::string("petsc")) return;
 
  auto* petsc_mat = dynamic_cast<SF::petsc_matrix*>(mat);
  if (petsc_mat == nullptr) return;
 
  Vec x = NULL, y = NULL;
  MatCreateVecs(petsc_mat->data, &x, &y);
 
  PetscInt i_start = 0, i_end = 0;
  VecGetOwnershipRange(x, &i_start, &i_end);
 
  PetscScalar* x_arr = NULL;
  int rank = 0;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 
  for (int probe_id = 0; probe_id < num_probes; ++probe_id) {
    VecGetArray(x, &x_arr);
    for (PetscInt i = i_start; i < i_end; ++i) {
      x_arr[i - i_start] = emi_probe_vector_entry(i, probe_id);
    }
    VecRestoreArray(x, &x_arr);
 
    MatMult(petsc_mat->data, x, y);
 
    PetscReal y_norm2 = 0.0;
    PetscReal y_norminf = 0.0;
    PetscScalar y_sum = 0.0;
    PetscScalar xAy = 0.0;
    VecNorm(y, NORM_2, &y_norm2);
    VecNorm(y, NORM_INFINITY, &y_norminf);
    VecSum(y, &y_sum);
    VecDot(x, y, &xAy);
 
    PetscScalar weighted_checksum_local = 0.0;
    const PetscScalar* y_arr = NULL;
    VecGetArrayRead(y, &y_arr);
    for (PetscInt i = i_start; i < i_end; ++i) {
      const PetscScalar weight = static_cast<PetscScalar>(i + 1);
      weighted_checksum_local += weight * y_arr[i - i_start];
    }
    VecRestoreArrayRead(y, &y_arr);
 
    PetscScalar weighted_checksum = 0.0;
    MPI_Allreduce(&weighted_checksum_local, &weighted_checksum, 1, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
 
    if (rank == 0) {
      PetscPrintf(PETSC_COMM_SELF,
                  "%s: operator probe %d ||Ax||_2=%g, ||Ax||_inf=%g, sum(Ax)=%g, x^T A x=%g, weighted_checksum=%g\n",
                  stage, probe_id + 1, double(y_norm2), double(y_norminf), double(PetscRealPart(y_sum)),
                  double(PetscRealPart(xAy)), double(PetscRealPart(weighted_checksum)));
    }
  }
 
  VecDestroy(&x);
  VecDestroy(&y);
}
 
void log_rhs_probe(SF::abstract_vector<SF_int, SF_real>* vec,
                   const char* stage)
{
  if (param_globals::flavor != std::string("petsc")) return;
 
  auto* petsc_vec = dynamic_cast<SF::petsc_vector*>(vec);
  if (petsc_vec == nullptr) return;
 
  PetscReal norm2 = 0.0;
  PetscReal norminf = 0.0;
  PetscScalar sum = 0.0;
  VecNorm(petsc_vec->data, NORM_2, &norm2);
  VecNorm(petsc_vec->data, NORM_INFINITY, &norminf);
  VecSum(petsc_vec->data, &sum);
 
  PetscInt i_start = 0, i_end = 0;
  VecGetOwnershipRange(petsc_vec->data, &i_start, &i_end);
 
  PetscScalar weighted_checksum_local = 0.0;
  const PetscScalar* arr = NULL;
  VecGetArrayRead(petsc_vec->data, &arr);
  for (PetscInt i = i_start; i < i_end; ++i) {
    const PetscScalar weight = static_cast<PetscScalar>(i + 1);
    weighted_checksum_local += weight * arr[i - i_start];
  }
  VecRestoreArrayRead(petsc_vec->data, &arr);
 
  PetscScalar weighted_checksum = 0.0;
  MPI_Allreduce(&weighted_checksum_local, &weighted_checksum, 1, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
 
  int rank = 0;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  if (rank == 0) {
    PetscPrintf(PETSC_COMM_SELF,
                "%s: rhs probe ||b||_2=%g, ||b||_inf=%g, sum(b)=%g, weighted_checksum=%g\n",
                stage, double(norm2), double(norminf), double(PetscRealPart(sum)),
                double(PetscRealPart(weighted_checksum)));
  }
}
 
void log_linear_system_probe(SF::abstract_matrix<SF_int, SF_real>* mat,
                             SF::abstract_vector<SF_int, SF_real>* vec,
                             const char* stage,
                             int num_probes = 3)
{
  if (param_globals::flavor != std::string("petsc")) return;
 
  auto* petsc_mat = dynamic_cast<SF::petsc_matrix*>(mat);
  auto* petsc_vec = dynamic_cast<SF::petsc_vector*>(vec);
  if (petsc_mat == nullptr || petsc_vec == nullptr) return;
 
  Vec x = NULL, ax = NULL, residual = NULL;
  MatCreateVecs(petsc_mat->data, &x, &ax);
  VecDuplicate(ax, &residual);
 
  PetscInt i_start = 0, i_end = 0;
  VecGetOwnershipRange(x, &i_start, &i_end);
 
  int rank = 0;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 
  for (int probe_id = 0; probe_id < num_probes; ++probe_id) {
    PetscScalar* x_arr = NULL;
    VecGetArray(x, &x_arr);
    for (PetscInt i = i_start; i < i_end; ++i) {
      x_arr[i - i_start] = emi_probe_vector_entry(i, probe_id);
    }
    VecRestoreArray(x, &x_arr);
 
    MatMult(petsc_mat->data, x, ax);
    VecWAXPY(residual, -1.0, petsc_vec->data, ax);
 
    PetscReal residual_norm2 = 0.0;
    PetscReal residual_norminf = 0.0;
    PetscScalar residual_sum = 0.0;
    PetscScalar xTResidual = 0.0;
    VecNorm(residual, NORM_2, &residual_norm2);
    VecNorm(residual, NORM_INFINITY, &residual_norminf);
    VecSum(residual, &residual_sum);
    VecDot(x, residual, &xTResidual);
 
    PetscScalar weighted_checksum_local = 0.0;
    const PetscScalar* residual_arr = NULL;
    VecGetArrayRead(residual, &residual_arr);
    for (PetscInt i = i_start; i < i_end; ++i) {
      const PetscScalar weight = static_cast<PetscScalar>(i + 1);
      weighted_checksum_local += weight * residual_arr[i - i_start];
    }
    VecRestoreArrayRead(residual, &residual_arr);
 
    PetscScalar weighted_checksum = 0.0;
    MPI_Allreduce(&weighted_checksum_local, &weighted_checksum, 1, MPIU_SCALAR, MPI_SUM, PETSC_COMM_WORLD);
 
    if (rank == 0) {
      PetscPrintf(PETSC_COMM_SELF,
                  "%s: linear-system probe %d ||Ax-b||_2=%g, ||Ax-b||_inf=%g, sum(Ax-b)=%g, x^T(Ax-b)=%g, weighted_checksum=%g\n",
                  stage, probe_id + 1, double(residual_norm2), double(residual_norminf),
                  double(PetscRealPart(residual_sum)), double(PetscRealPart(xTResidual)),
                  double(PetscRealPart(weighted_checksum)));
    }
  }
 
  VecDestroy(&x);
  VecDestroy(&ax);
  VecDestroy(&residual);
}
#endif
 
void set_elec_tissue_properties_emi_volume(MaterialType* mtype, hashmap::unordered_set<int> &extra_tags, hashmap::unordered_set<int> &intra_tags, FILE_SPEC logger)
{
  CALI_CXX_MARK_FUNCTION;
  MaterialType *m = mtype;
  
  // initialize random conductivity fluctuation structure with PrM values
  m->regions.resize(param_globals::num_gregions);
 
  const char* grid_name = "emi_grid_domain";
  log_msg(logger, 0, 0, "Setting up %s tissue poperties for %d regions ..", grid_name,
          param_globals::num_gregions);
 
  char buf[64];
  RegionSpecs* reg = m->regions.data();
 
  // default tags for extra and intra domains
  hashmap::unordered_set<int> extra_tags_default = extra_tags;
  hashmap::unordered_set<int> intra_tags_default = intra_tags;
 
  for (size_t i=0; i<m->regions.size(); i++)
  {
    for (int j=0;j<param_globals::gregion[i].num_IDs;j++)
    {
      int tag = param_globals::gregion[i].ID[j];
 
      // removed tags explicitly defined by user in param_globals::gregion[i>1]
      if(extra_tags_default.find(tag) != extra_tags_default.end())
        extra_tags_default.erase(tag);
 
      if(intra_tags_default.find(tag) != intra_tags_default.end())
        intra_tags_default.erase(tag);
    }
  }
 
  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);
    }
 
    // copy metadata into RegionSpecs
    reg->regname  = strdup(param_globals::gregion[i].name);
    reg->regID    = i;
 
    if(i==0) // default Extracellular region
      reg->nsubregs = extra_tags_default.size();
    if(i==1) // default intracellular region
      reg->nsubregs = intra_tags_default.size();
    if(i>1)  // optional: rest of other param_globals::gregion[i>1] defined by user 
      reg->nsubregs = param_globals::gregion[i].num_IDs;
 
    if(!reg->nsubregs)
      reg->subregtags = NULL;
    else  
    {
      reg->subregtags  = new int[reg->nsubregs];
 
      if(i==0){
        int j = 0;
        for (int tag : extra_tags_default) {
          reg->subregtags[j] = tag;
          j++;
        }
      }
      else if(i==1){
        int j = 0;
        for (int tag : intra_tags_default) {
          reg->subregtags[j] = tag;
          j++;
        }
      }
      else{
        for (int j=0;j<reg->nsubregs;j++)
          reg->subregtags[j] = param_globals::gregion[i].ID[j]; // explicit tags defined by user
      }
    }
 
    // describe material in given region
    elecMaterial *emat  = new elecMaterial();
    emat->material_type = ElecMat;
 
    // Isotropic conductivity is considered in EMI model.
    emat->InVal[0]   = param_globals::gregion[i].g_bath;
    emat->InVal[1]   = param_globals::gregion[i].g_bath;
    emat->InVal[2]   = param_globals::gregion[i].g_bath;
 
    emat->ExVal[0]   = param_globals::gregion[i].g_bath;
    emat->ExVal[1]   = param_globals::gregion[i].g_bath;
    emat->ExVal[2]   = param_globals::gregion[i].g_bath;
 
    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 (strlen(param_globals::gi_scale_vec))
    read_el_scale_vec(param_globals::gi_scale_vec, emi_msh, m->el_scale, m->el_scale_dpn);
}
 
void parabolic_solver_emi::init()
{
  CALI_CXX_MARK_FUNCTION;
  double t0, t1, dur;
  get_time(t0);
  const int log_flag = param_globals::output_level > 1 ? ECHO : 0;
  const auto log_init_timing = [&](const char* label, double start) {
    log_msg(NULL, 0, log_flag, "EMI solver init: %s in %.5f seconds.", label, float(MPI_Wtime() - start));
  };
 
  double phase_t = MPI_Wtime();
  stats.init_logger("par_stats.dat");
 
  // Create/initialise linear solver object: PETSc gets constructed with default settings
  SF::init_solver(&lin_solver);
  log_init_timing("linear solver object", phase_t);
 
  // EMI mesh: After decoupling the interfaces, we obtain a volumetric mesh with new dofs for all vertices.
  sf_mesh & emi_mesh = get_mesh(emi_msh);
  // Surface mesh called emi_surface_counter_msh:
  // This mesh contains both sides of each interface. At this stage its nodes
  // already use the decoupled EMI DOF numbering.
  sf_mesh & emi_surfmesh_w_counter_face = get_mesh(emi_surface_counter_msh);
  sf_mesh & emi_surfmesh_one_side = get_mesh(emi_surface_msh);
  sf_mesh & emi_surfmesh_unique_face = get_mesh(emi_surface_unique_face_msh);
 
  phase_t = MPI_Wtime();
  int max_row_entries_emi = max_nodal_edgecount(emi_mesh);
  log_init_timing("maximum nodal edge counts", phase_t);
 
  int rank;
  MPI_Comm_rank(emi_surfmesh_w_counter_face.comm, &rank);
 
  sf_vec::ltype alg_type = sf_vec::algebraic;
  sf_vec::ltype alg_surface_type = sf_vec::elemwise;
 
  int dpn = 1;
 
  //-----------------------------------------------------------------
  // setup vectors
  //-----------------------------------------------------------------
  phase_t = MPI_Wtime();
  SF::init_vector(&ui,       emi_mesh, dpn, alg_type);
  SF::init_vector(&dui,      emi_mesh, dpn, alg_type);
  SF::init_vector(&ui_pre,   emi_mesh, dpn, alg_type);
  SF::init_vector(&Irhs,     emi_mesh, dpn, alg_type);
  SF::init_vector(&Iij_stim, emi_mesh, dpn, alg_type);
  SF::init_vector(&Iij_temp, emi_mesh, dpn, alg_type);
  SF::init_vector(&vb_both_face, emi_surfmesh_w_counter_face, dpn, alg_surface_type);
  SF::init_vector(&vb_unique_face, emi_surfmesh_unique_face, dpn, alg_surface_type);
  SF::init_vector(&Ib_both_face, emi_surfmesh_w_counter_face, dpn, alg_surface_type);
  SF::init_vector(&Ib_unique_face, emi_surfmesh_unique_face, dpn, alg_surface_type);
  log_init_timing("vectors", phase_t);
 
  // PETSc matrices are reallocated with an exact sparsity graph by replaying
  // the real assembly into MATPREALLOCATOR before the numeric assembly. The
  // scalar hint used here is therefore only a temporary bootstrap value for
  // PETSc; non-PETSc backends keep the conservative scalar bounds.
  const bool use_petsc_exact_preallocation = param_globals::flavor == std::string("petsc");
  const int petsc_initial_prealloc = 1;
  //---------------------------------------------------------------------
  // initialize operator matrices B, Bi, and BsM
  //---------------------------------------------------------------------
  mesh_int_t M = emi_surfmesh_w_counter_face.g_numelem;
  mesh_int_t N = emi_mesh.g_numpts;
  mesh_int_t m = emi_surfmesh_w_counter_face.l_numelem;
  mesh_int_t m_one_side = emi_surfmesh_one_side.l_numelem;
  mesh_int_t M_one_side = emi_surfmesh_one_side.g_numelem; 
  mesh_int_t m_unique_face = emi_surfmesh_unique_face.l_numelem;
  mesh_int_t M_unique_face = emi_surfmesh_unique_face.g_numelem;
  mesh_int_t n = ui->lsize();
 
  if (param_globals::output_level > 1) {
    log_msg(NULL, 0, 0, "\n**********************************");
    log_msg(NULL, 0, 0, "#elements of emi surfmesh unique face: %zu", emi_surfmesh_unique_face.g_numelem);
    log_msg(NULL, 0, 0, "#elements of emi surfmesh one side: %zu", emi_surfmesh_one_side.g_numelem);
    log_msg(NULL, 0, 0, "#elements of emi surfmesh: %zu", emi_surfmesh_w_counter_face.g_numelem);
    log_msg(NULL, 0, 0, "#elements of emi mesh: %zu", emi_mesh.g_numelem);
    log_msg(NULL, 0, 0, "#dofs for emi_mesh: %zu", emi_mesh.g_numpts);
    log_msg(NULL, 0, 0, "#max_row_entries_emi: %zu", max_row_entries_emi);
    log_msg(NULL, 0, 0, "**********************************\n");
    log_mesh_local_element_ranges(emi_mesh, emi_surfmesh_w_counter_face, emi_surfmesh_unique_face);
  }
 
  SF::vector<long int> layout;
  SF::layout_from_count<long int>(emi_surfmesh_w_counter_face.l_numelem, layout, emi_surfmesh_w_counter_face.comm);
  mesh_int_t m_l = layout[rank];
  mesh_int_t n_l = emi_mesh.pl.algebraic_layout()[rank];
  const SF::vector<mesh_int_t> & alg_nod_surface = emi_surfmesh_w_counter_face.pl.algebraic_nodes();
 
  SF::vector<long int> layout_one_side;
  SF::layout_from_count<long int>(emi_surfmesh_one_side.l_numelem, layout_one_side, emi_surfmesh_one_side.comm);
  mesh_int_t m_one_side_l = layout_one_side[rank];
 
  SF::vector<long int> layout_unique_face;
  SF::layout_from_count<long int>(emi_surfmesh_unique_face.l_numelem, layout_unique_face, emi_surfmesh_unique_face.comm);
  mesh_int_t m_unique_face_l = layout_unique_face[rank];
 
  // B, Bi, and BsM have face/volume row spaces whose final PETSc sparsity is
  // collected exactly below by assemble_with_exact_preallocation(). The scalar
  // fallback remains conservative for non-PETSc backends and off-process row
  // insertion through the abstract matrix interface. Resting-potential
  // initialization writes directly into ui and no longer needs a diagonal helper
  // matrix.
  phase_t = MPI_Wtime();
  SF::init_matrix(&B);
  B->init(M, N, m, n, m_l, use_petsc_exact_preallocation ? petsc_initial_prealloc : max_row_entries_emi*2);
  B->zero();
  SF::init_matrix(&Bi);
  Bi->init(M, N, m, n, m_l, use_petsc_exact_preallocation ? petsc_initial_prealloc : max_row_entries_emi*2);
  Bi->zero();
  SF::init_matrix(&BsM);
  BsM->init(N, M, n, m, n_l, use_petsc_exact_preallocation ? petsc_initial_prealloc : max_row_entries_emi*2);
  BsM->zero();
  log_init_timing("EMI coupling matrices", phase_t);
 
  // Build the direct unique-face <-> both-face transfer operators
  phase_t = MPI_Wtime();
  SF::construct_direct_unique_both_operators(operator_unique_to_both_faces,
                                             operator_both_to_unique_face,
                                             map_elem_uniqueFace_to_elem_bothface,
                                             map_elem_uniqueFace_to_elem_oneface,
                                             vec_both_to_one_face,
                                             emi_surfmesh_w_counter_face,
                                             emi_surfmesh_unique_face,
                                             max_row_entries_emi,
                                             dpn,
                                             alg_surface_type);
  log_init_timing("unique/both face transfer operators", phase_t);
 
  // Initialize interpolation operators B/Bi and the scatter operator BsM.
  phase_t = MPI_Wtime();
  assemble_with_exact_preallocation({B, Bi, BsM}, [&]() {
    SF::assemble_restrict_operator(*B, *Bi, *BsM, elemTag_surface_w_counter_mesh, map_vertex_tag_to_dof_petsc, line_face, tri_face, quad_face, emi_surfmesh_w_counter_face, emi_mesh, UM2_to_CM2);
  });
  log_init_timing("restriction operators", phase_t);
 
  //-----------------------------------------------------------------
  // initialize matrices (LHS, K, M_{emi mesh}, M_{surface mesh})
  //-----------------------------------------------------------------
  phase_t = MPI_Wtime();
  SF::init_matrix(&lhs_emi);
  SF::init_matrix(&stiffness_emi);
  SF::init_matrix(&mass_emi);
  SF::init_matrix(&mass_surf_emi);
 
  // These matrices are filled later in rebuild_matrices(). Allocating them
  // after the transfer/restriction operators reduces the EMI init memory peak.
  lhs_emi->init(emi_mesh, dpn, dpn, use_petsc_exact_preallocation ? petsc_initial_prealloc : max_row_entries_emi*3);
  stiffness_emi->init(emi_mesh, dpn, dpn, use_petsc_exact_preallocation ? petsc_initial_prealloc : max_row_entries_emi*2);
  mass_emi->init(emi_mesh, dpn, dpn, use_petsc_exact_preallocation ? petsc_initial_prealloc : max_row_entries_emi*3);
  mass_surf_emi->init(emi_mesh, dpn, dpn, use_petsc_exact_preallocation ? petsc_initial_prealloc : max_row_entries_emi*3);
  log_init_timing("system matrices", phase_t);
 
  // DEBUG: Check mesh sizes and mappings
  #ifdef EMI_DEBUG_MESH
  {
    int local_rank;
    MPI_Comm_rank(emi_mesh.comm, &local_rank);
    
    fprintf(stderr, "RANK %d MESH SIZES: one_side=%zu, counter=%zu, unique=%zu\n",
            local_rank, emi_surfmesh_one_side.l_numelem, 
            emi_surfmesh_w_counter_face.l_numelem, emi_surfmesh_unique_face.l_numelem);
    
    fprintf(stderr, "RANK %d MAP SIZES: uniqueFace_to_oneface=%zu\n",
            local_rank, map_elem_uniqueFace_to_elem_oneface.size());
    fflush(stderr);
  }
  #endif
 
  decltype(map_elem_uniqueFace_to_elem_oneface)().swap(map_elem_uniqueFace_to_elem_oneface);
  decltype(map_elem_uniqueFace_to_elem_bothface)().swap(map_elem_uniqueFace_to_elem_bothface);
  vec_both_to_one_face = SF::vector<mesh_int_t>();
 
  //-----------------------------------------------------------------
  // setup the ionic current
  //-----------------------------------------------------------------
  phase_t = MPI_Wtime();
  // get initial value of the vm and Iion on the face from the ionicOnFace class
  sf_vec* vb_ptr   = get_data(vm_emi_itf);
  sf_vec* Ib_ptr = get_data(iion_emi_itf);
 
  if(!(vb_ptr != NULL && Ib_ptr != NULL)) {
    log_msg(0,5,0, "%s error: global Vb and Ib vectors not properly set up! Ionics seem invalid! Aborting!",
            __func__);
    EXIT(1);
  }
 
  SF::init_vector(&vb, vb_ptr);
  vb->shallow_copy(*vb_ptr);
  SF::init_vector(&Ib, Ib_ptr);
  Ib->shallow_copy(*Ib_ptr);
 
  parab_tech = static_cast<parabolic_solver_emi::parabolic_t>(param_globals::parab_solve_emi);
  log_init_timing("ionic face vectors", phase_t);
 
  dur = timing(t1, t0);
  log_msg(NULL, 0, log_flag, "EMI solver init total in %.5f seconds.", float(dur));
}
 
void parabolic_solver_emi::rebuild_matrices(MaterialType* mtype, limpet::MULTI_IF & miif, SF::vector<stimulus> & stimuli, FILE_SPEC logger)
{
  CALI_CXX_MARK_FUNCTION;
  double start, end, period;
  get_time(start);
  double t0, t1, dur;
  mass_integrator mass_integ;
  mass_integrator mass_integ_emi;
  int dpn = 1;
 
  int            log_flag = param_globals::output_level > 1 ? ECHO : 0;
  MaterialType & mt       = mtype[0];
  const bool have_dbc = have_dbc_stims(stimuli);
  const bool use_petsc_exact_preallocation = param_globals::flavor == std::string("petsc");
  const bool reuse_petsc_fem_preallocation = use_petsc_exact_preallocation && fem_matrices_exact_preallocated;
 
  double Dt = user_globals::tm_manager->time_step;
 
  cond_t    condType = intra_cond;
  sf_mesh & mesh     = get_mesh(emi_msh);
  sf_mesh & emi_surfmesh_w_counter_face    = get_mesh(emi_surface_counter_msh);
 
  // set material and conductivity type
  set_cond_type(mt, condType);
 
  // assemble EMI matrices
  // fill the EMI system
  {
    // emi-> stiffness
    log_msg(NULL, 0, 0, "assemble stiffness matrix");
    get_time(t0);
    elec_stiffness_integrator stfn_integ_emi(mt);
    auto assemble_stiffness_emi = [&]() {
      stiffness_emi->zero();
      SF::assemble_matrix(*stiffness_emi, mesh, stfn_integ_emi);
    };
    if(reuse_petsc_fem_preallocation) {
      assemble_stiffness_emi();
    } else {
      assemble_with_exact_preallocation({stiffness_emi}, assemble_stiffness_emi);
    }
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Computed parabolic stiffness matrix in %.5f seconds.", float(dur));
 
    log_msg(NULL, 0, 0, "assemble mass matrix on the volumetric mesh");
    get_time(t0);
    mass_integrator mass_integ;
    auto assemble_mass_emi = [&]() {
      mass_emi->zero();
      SF::assemble_matrix(*mass_emi, mesh, mass_integ);
    };
    auto* petsc_mass_emi = use_petsc_exact_preallocation ? dynamic_cast<SF::petsc_matrix*>(mass_emi) : nullptr;
    auto* petsc_stiffness_emi = use_petsc_exact_preallocation ? dynamic_cast<SF::petsc_matrix*>(stiffness_emi) : nullptr;
    if(reuse_petsc_fem_preallocation) {
      assemble_mass_emi();
    } else if(petsc_mass_emi != nullptr && petsc_stiffness_emi != nullptr) {
      // mass_emi and stiffness_emi have the same volume FEM sparsity pattern.
      petsc_mass_emi->duplicate_pattern(*petsc_stiffness_emi);
      assemble_mass_emi();
    } else {
      assemble_with_exact_preallocation({mass_emi}, assemble_mass_emi);
    }
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Computed volumetric mass matrix in %.5f seconds.", float(dur));
 
    log_msg(NULL, 0, 0, "assemble LHS matrix and mass matrix on the surface mesh");
    get_time(t0);
    auto assemble_lhs_and_surface_mass = [&]() {
      lhs_emi->zero();
      mass_surf_emi->zero();
      // lhs_emi = (CmM/dt + K)
      SF::assemble_lhs_emi(*lhs_emi, *mass_surf_emi, mesh, emi_surfmesh_w_counter_face, map_vertex_tag_to_dof_petsc, line_face, tri_face, quad_face, stfn_integ_emi, mass_integ_emi, -1., UM2_to_CM2 / Dt);
    };
    if(reuse_petsc_fem_preallocation) {
      assemble_lhs_and_surface_mass();
    } else {
      assemble_with_exact_preallocation({lhs_emi, mass_surf_emi}, assemble_lhs_and_surface_mass);
    }
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Computed parabolic mass matrix in %.5f seconds.", float(dur));
  }
 
 
  bool same_nonzero = false;
 
  // set boundary conditions
  if(have_dbc) {
    log_msg(logger,0,log_flag, "lhs matrix enforcing Dirichlet boundaries.");
    get_time(t0);
 
    if(dbc == nullptr)
      dbc = new dbc_manager(*lhs_emi, stimuli);
    else
      dbc->recompute_dbcs();
 
    dbc->enforce_dbc_lhs();
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "lhs matrix Dirichlet enforcing done in %.5f seconds.", float(dur));
  }
  else {
    log_msg(logger,0,ECHO, "without enforcing Dirichlet boundaries on the lhs matrix!");
    // we are dealing with a singular system
    phie_mat_has_nullspace = true;
  }
 
  set_dir(INPUT);
  log_emi_petsc_matrix_preallocation_report({
      {"B:", B},
      {"Bi:", Bi},
      {"BsM:", BsM},
      {"unique_to_both:", operator_unique_to_both_faces},
      {"both_to_unique:", operator_both_to_unique_face},
      {"stiffness_emi:", stiffness_emi},
      {"mass_emi:", mass_emi},
      {"mass_surf_emi:", mass_surf_emi},
      {"lhs_emi:", lhs_emi},
  });
  if(use_petsc_exact_preallocation) fem_matrices_exact_preallocated = true;
  get_time(t0);
  setup_linear_solver(logger);
  dur = timing(t1,t0);
  log_msg(logger,0,log_flag, "Initializing parabolic solver in %.5f seconds.", float(dur));
  set_dir(OUTPUT);
 
  period = timing(end, start);
}
 
void parabolic_solver_emi::setup_linear_solver(FILE_SPEC logger)
{
  tol    = param_globals::cg_tol_parab;
  max_it = param_globals::cg_maxit_parab;
 
  std::string default_opts;
  std::string solver_file;
  solver_file = param_globals::parab_options_file;
  if (param_globals::flavor == std::string("ginkgo")) {
    default_opts = std::string(
        R"(
{
    "type": "solver::Cg",
    "preconditioner": {
        "type": "solver::Multigrid",
        "min_coarse_rows": 8,
      "max_levels": 16,
        "default_initial_guess": "zero",
        "mg_level": [
            {
                "type": "multigrid::Pgm",
                "deterministic": false
            }
        ],
        "coarsest_solver": {
             "type": "preconditioner::Schwarz",
             "local_solver": {
                 "type": "preconditioner::Jacobi"
             }
        },
        "criteria": [
            {
                "type": "Iteration",
                "max_iters": 1
            }
        ]
    },
    "criteria": [
        {
            "type": "Iteration",
            "max_iters": 100
        },
        {
            "type": "ResidualNorm",
            "reduction_factor": 1e-4
        }
    ]
}
      )");
  } else if (param_globals::flavor == std::string("petsc")) {
    default_opts = std::string("-ksp_type cg -pc_type gamg -options_left");
  }
  lin_solver->setup_solver(*lhs_emi, tol, max_it * 100, param_globals::cg_norm_parab,
                           "parabolic PDE", phie_mat_has_nullspace, logger, solver_file.c_str(),
                           default_opts.c_str());
}
 
void parabolic_solver_emi::solve()
{
  switch (parab_tech) {
    case SEMI_IMPLICIT:  solve_semiImplicit();   break;
  }
}
 
void parabolic_solver_emi::solve_semiImplicit()
{
  double t0,t1;
  get_time(t0);
 
  if(dbc != nullptr){  
    CALI_MARK_BEGIN("apply_dbc_rhs");
    dbc->enforce_dbc_rhs(*ui);
    CALI_MARK_END("apply_dbc_rhs");
  }
 
  *ui_pre = *ui;
 
  // -K*u (K is assembled as -K)
  CALI_MARK_BEGIN("stiff_mat");
  stiffness_emi->mult(*ui, *Iij_temp);
  // rhs = Iij + K*u
  *Irhs -= *Iij_temp;
  CALI_MARK_END("stiff_mat");
 
  // add volumetric stimulus currents (I_ex, I_in)
  // rhs = Iij + K*u + M*Iij_stim
  if (Iij_stim->mag() > 0.0) {
    CALI_MARK_BEGIN("stim_application");
    mass_emi->mult(*Iij_stim, *Iij_temp);
    *Irhs -= *Iij_temp;
    CALI_MARK_BEGIN("stim_application");
  }
 
  // rhs = -(Iij + K*u + M*Iij_stim)
  CALI_MARK_BEGIN("rhs_update");
  (*Irhs) *= (-1.0);
  CALI_MARK_END("rhs_update");
 
  // compute step
  CALI_MARK_BEGIN("linear_solve");
  (*lin_solver)(*dui, *Irhs);    
  CALI_MARK_END("linear_solve");
 
  // logfile for solver
  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);
  }
 
  // update solution:: ui_pre = ui_pre + dui
  CALI_MARK_BEGIN("sol_update");
  ui_pre->add_scaled(*dui, 1.0);
  *ui *=0;
  ui->add_scaled(*ui_pre, 1.0);
  CALI_MARK_END("sol_update");
 
  // We need to enforce DBCs again after the solution vector was updated.
  // Otherwise, matrix/solver tolerances allow small nonzero values in dui to accumulate in ui.
  if(dbc != nullptr){  
    CALI_MARK_BEGIN("apply_dbc_rhs");
    dbc->enforce_dbc_rhs(*ui);
    CALI_MARK_END("apply_dbc_rhs");
  }
 
  // treat solver statistics
  stats.slvtime += timing(t1, t0);
  stats.update_iter(lin_solver->niter);
}
 
void tags_onFace( hashmap::unordered_map<SF::tuple<mesh_int_t>, 
                      std::pair<SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>,
                      SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>>> & line_face,
                      hashmap::unordered_map<SF::triple<mesh_int_t>, 
                      std::pair<SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>,
                      SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>>> & tri_face,
                      hashmap::unordered_map<SF::quadruple<mesh_int_t>, 
                      std::pair<SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>, 
                      SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>>> & quad_face,
                      hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, mesh_int_t> & map_vertex_tag_to_dof,
                      hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, std::pair<mesh_int_t,mesh_int_t>> & map_vertex_tag_to_dof_petsc,
                      std::vector<std::string> & tags_data)
{
  sf_mesh & emi_surfmesh_w_counter_face = get_mesh(emi_surface_counter_msh);
 
  const SF::vector<mesh_int_t> & tags = emi_surfmesh_w_counter_face.tag;
 
  const SF::vector<mesh_int_t> & rnod = emi_surfmesh_w_counter_face.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<mesh_int_t,mesh_int_t> l2g;
 
  for(size_t i=0; i<rnod.size(); i++){
    l2g[i] = rnod[i];
  }
 
  for(size_t eidx=0; eidx<emi_surfmesh_w_counter_face.l_numelem; eidx++)
  {
    std::vector<mesh_int_t> elem_nodes;
    mesh_int_t tag = emi_surfmesh_w_counter_face.tag[eidx];
    for (int n = emi_surfmesh_w_counter_face.dsp[eidx]; n < emi_surfmesh_w_counter_face.dsp[eidx+1];n++)
    {
      mesh_int_t l_idx = emi_surfmesh_w_counter_face.con[n];
 
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(l2g[l_idx],tags[eidx]);
      mesh_int_t dof = map_vertex_tag_to_dof[Index_tag_old];
      elem_nodes.push_back(dof);
    }
 
    mesh_int_t tag_first = 0;
    mesh_int_t tag_second = 0;
    std::string result_first;
    std::string result_second;
    std::sort(elem_nodes.begin(),elem_nodes.end()); // make the node tuple order-independent
 
    // Extract all surface face pairs separating regions with different material or boundary tags (ionicFaces)
    if(elem_nodes.size()==2){
      SF::tuple<mesh_int_t> key;
 
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      std::pair<SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>,
                SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>> value = line_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
      result_first = std::to_string(tag_first) + ":" + std::to_string(tag_second);
      result_second = std::to_string(tag_second) + ":" + std::to_string(tag_first); 
    }
    else if(elem_nodes.size()==3){
      SF::triple<mesh_int_t> key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      std::pair<SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>,
                SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>>  value = tri_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
      result_first = std::to_string(tag_first) + ":" + std::to_string(tag_second);
      result_second = std::to_string(tag_second) + ":" + std::to_string(tag_first);
    }
    else if(elem_nodes.size()==4){
      SF::quadruple<mesh_int_t> key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      key.v4 = elem_nodes[3];
      std::pair<SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>,
                SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>> value = quad_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
      result_first = std::to_string(tag_first) + ":" + std::to_string(tag_second);
      result_second = std::to_string(tag_second) + ":" + std::to_string(tag_first);
    }
    // check if current element tag has a custom region ID
    tags_data.push_back(result_first);
    tags_data.push_back(result_second);
  }
}
 
void EMI::initialize()
{
  double t1, t2;
  get_time(t1);
 
  set_dir(OUTPUT);
 
  // open logger
  logger = f_open("emi.log", param_globals::experiment != 4 ? "w" : "r");
  const int verb = param_globals::output_level;
  const auto log_init_timing = [&](const char* label, double start) {
    if(verb)
      log_msg(logger, 0, ECHO, "EMI init: %s in %.5f seconds.", label, float(MPI_Wtime() - start));
  };
 
  double phase_t = MPI_Wtime();
  // Mesh processing step: convert the input mesh into 
  // - an EMI mesh (with discontinuities at gap junctions and membranes) 
  // - and generate the corresponding surface mesh
  CALI_MARK_BEGIN("mesh_setup");
  setup_EMI_mesh();
  log_init_timing("mesh setup", phase_t);
 
  // setup mappings between extra and intra grids, algebraic and nodal,
  // and between PETSc and canonical orderings
  phase_t = MPI_Wtime();
  setup_mappings();
  log_init_timing("mesh mappings", phase_t);
 
  // the ionicOnFace physics is currently triggered from inside the emi to have tighter
  // control over it
  phase_t = MPI_Wtime();
  ion.logger = logger;
 
  ion.set_surface_mesh_data(parab_solver.line_face,
                            parab_solver.tri_face,
                            parab_solver.quad_face,
                            parab_solver.map_vertex_tag_to_dof);
 
  // builds per-face adjacency tags
  std::vector<std::string> tags_data;                    
  tags_onFace(parab_solver.line_face,
              parab_solver.tri_face,
              parab_solver.quad_face,
              parab_solver.map_vertex_tag_to_dof,
              parab_solver.map_vertex_tag_to_dof_petsc,
              tags_data);
  CALI_MARK_END("mesh_setup");
  log_init_timing("ionic face metadata", phase_t);
 
  ion.set_tags_onFace(tags_data);
 
  ion.set_face_region_data(parab_solver.intra_tags, parab_solver.map_elem_uniqueFace_to_tags);
  phase_t = MPI_Wtime();
  ion.initialize();
  log_init_timing("ionic model initialization", phase_t);
 
  phase_t = MPI_Wtime();
  // set up tissue properties on the extra and intracellular domains
  set_elec_tissue_properties_emi_volume(mtype_vol, parab_solver.extra_tags, parab_solver.intra_tags, logger);
  // In EMI, the default extra/intra compartment tags are intentionally handled
  // as the implicit regions 0 and 1, so the generic "unassigned" warning is
  // not useful here.
  region_mask(emi_msh, mtype_vol[0].regions, mtype_vol[0].regionIDs, true, "gregion_vol", false);
  
  // 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");
  log_init_timing("tissue properties and timers", phase_t);
 
  // EMI stimuli setup
  CALI_MARK_BEGIN("stimulus_setup");
  phase_t = MPI_Wtime();
  param_globals::operator_splitting = 0; // EMI does not use operator splitting, so keep stimulus scaling at the monodomain/default setting.
  setup_stimuli();
  log_init_timing("stimuli", phase_t);
  CALI_MARK_END("stimulus_setup");
 
  // set up the linear equation systems. this needs to happen after the stimuli have been
  // set up, since we need boundary condition info
  CALI_MARK_BEGIN("solver_setup");
  phase_t = MPI_Wtime();
  setup_solvers();
  log_init_timing("solver setup", phase_t);
  CALI_MARK_END("solver_setup");
  
  phase_t = MPI_Wtime();
  // Balance paired electrodes before total-current scaling.
  balance_electrodes();
  // total current scaling
  scale_total_stimulus_current(stimuli, *parab_solver.mass_emi, *parab_solver.mass_surf_emi, logger);
  log_init_timing("stimulus current scaling", phase_t);
 
  sf_mesh & emi_mesh = get_mesh(emi_msh);
  sf_mesh & emi_surfmesh_w_counter_face = get_mesh(emi_surface_counter_msh);
 
  phase_t = MPI_Wtime();
  // Initialize ui from the ionic resting membrane voltage:
  // - vb is defined on the unique-face ionic layout.
  // - The direct unique -> both operator expands it to the both-face layout.
  // - assign_resting_potential_from_ionic_models_on_myocyte writes the matching
  //   membrane voltage into intracellular volume DOFs and zero into extracellular DOFs.
  parab_solver.operator_unique_to_both_faces->mult(*parab_solver.vb, *parab_solver.vb_both_face);
  SF::assign_resting_potential_from_ionic_models_on_myocyte(*parab_solver.ui,
                    parab_solver.vb_both_face,
                    parab_solver.elemTag_emi_mesh,
                    parab_solver.map_vertex_tag_to_dof_petsc, 
                    parab_solver.line_face, parab_solver.tri_face, parab_solver.quad_face, 
                    emi_surfmesh_w_counter_face, emi_mesh);  
 
  *parab_solver.vb_unique_face = *parab_solver.vb;
  log_init_timing("initial membrane state projection", phase_t);
 
  CALI_MARK_BEGIN("output_setup");
  phase_t = MPI_Wtime();
  // prepare the electrics output. we skip it if we do post-processing
  if(param_globals::experiment != EXP_POSTPROCESS)
    setup_output();
  log_init_timing("output setup", phase_t);
  CALI_MARK_END("output_setup");
 
  const double init_dur = timing(t2, t1);
  this->initialize_time += init_dur;
  if(verb)
    log_msg(logger, 0, ECHO, "EMI init total in %.5f seconds.", float(init_dur));
}
 
void EMI::setup_mappings()
{
  bool emi_exits = mesh_is_registered(emi_msh);
  assert(emi_exits);
  const int dpn = 1;
 
  if(get_scattering(emi_msh, ALG_TO_NODAL, dpn) == NULL) 
  {
    log_msg(logger, 0, 0, "%s: Setting up intracellular algebraic-to-nodal scattering.", __func__);
    register_scattering(emi_msh, ALG_TO_NODAL, dpn);
  }
  if(get_permutation(emi_msh, PETSC_TO_CANONICAL, dpn) == NULL) 
  {
    log_msg(logger, 0, 0, "%s: Setting up intracellular PETSc to canonical permutation.", __func__);
    register_permutation(emi_msh, PETSC_TO_CANONICAL, dpn);
  }
}
 
void EMI::compute_step()
{
  double t1, t2;
  get_time(t1);
 
  // activation checking
  const double time      = user_globals::tm_manager->time,
               time_step = user_globals::tm_manager->time_step;
 
  const int verb = param_globals::output_level;
  // We treat stimuli by type:
  //    - Potential stimuli (Phi_ex, GND_ex, Phi_ex_ol, Phi_in, Phi_in_ol) are
  //      managed by a dbc_manager and applied to the left- and right-hand side.
  //    - Current stimuli (I_ex, I_in) are applied to the right hand side,
  //      while (I_tm) is applied to the vector Ib directly.
  CALI_MARK_BEGIN("apply_dbc_lhs");
  apply_dbc_stimulus();
  CALI_MARK_END("apply_dbc_lhs");
 
  // compute ionics update
  CALI_MARK_BEGIN("ion_compute");
  ion.compute_step();
  CALI_MARK_END("ion_compute");
 
  CALI_MARK_BEGIN("apply_stim");
  apply_current_stimulus();
  CALI_MARK_END("apply_stim");
 
  // convert Ib -> Irhs
  parab_solver.operator_unique_to_both_faces->mult(*parab_solver.Ib, *parab_solver.Ib_both_face);
  parab_solver.BsM->mult(*parab_solver.Ib_both_face, *parab_solver.Irhs);
  
  // solver parabolic system
  CALI_MARK_BEGIN("parab_solve");
  parab_solver.solve();
  {
    // v_b = B_i * u
    parab_solver.B->mult(*parab_solver.ui, *parab_solver.vb_both_face);
    // Direct both -> unique mapping
    parab_solver.operator_both_to_unique_face->mult(*parab_solver.vb_both_face, *parab_solver.vb_unique_face);
    *parab_solver.vb = *parab_solver.vb_unique_face;
  }
  CALI_MARK_END("parab_solve");
 
  if(user_globals::tm_manager->trigger(iotm_console)) {
    // output lin solver stats
    parab_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 EMI::output_step()
{
  double t1, t2;
  get_time(t1);
  
  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(user_globals::tm_manager->time, false);
}
 
void EMI::destroy()
{
  CALI_CXX_MARK_FUNCTION;
  // close logger
  f_close(logger);
 
  // close output files
  output_manager.close_files_and_cleanup();
 
  // destroy ionics
  ion.destroy();
}
 
void EMI::setup_stimuli()
{
  // initialize basic stim info data (used units, supported types, etc)
  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);
 
    // we associate to the EMI mesh. this is needed for the stim_phys and stim_electrode setups.
    s.associated_intra_mesh = emi_msh, s.associated_extra_mesh = emi_msh;
 
    s.setup(i);
 
    if (s.phys.type == Illum) {
      log_msg(0, MAX_LOG_LEVEL, ECHO, "Stimulus of type Illum (=6) is not implemented in EMI. Abort.");
      EXIT(EXIT_FAILURE);
    }
 
    // Depending on the stimulus type, we make sure to only stimulate the correct regions of the mesh:
    // Extracellular stimuli restrict to all DOFs of the extracellular region (including the extracellular side of the membrane).
    // Equivalent for intracellular stimuli.
    // Stimuli that act directly on the membrane restrict to DOFs with ptsData > 0, i.e. membrane, gap junctions, and complex junctions.
    if (is_extra(s.phys.type)) {
      SF::vector<mesh_int_t> extra_vertices;
      const sf_mesh&         mesh = get_mesh(s.associated_extra_mesh);
 
      // Gather vertices from emi_msh using extra_tags
      indices_from_region_tags(extra_vertices, mesh, parab_solver.extra_tags);
 
      // Restrict electrode vertices to extra region
      restrict_to_set(s.electrode.vertices, extra_vertices);
    } else if (s.phys.type == I_tm) {
      const sf_mesh&         mesh = get_mesh(emi_msh);
      SF::restrict_to_membrane(s.electrode.vertices, dof2ptsData, mesh);
    } else {
      SF::vector<mesh_int_t> intra_vertices;
      const sf_mesh&         mesh = get_mesh(s.associated_intra_mesh);
 
      // Gather vertices from emi_msh using intra_tags
      indices_from_region_tags(intra_vertices, mesh, parab_solver.intra_tags);
 
      // Restrict electrode vertices to intra region
      restrict_to_set(s.electrode.vertices, intra_vertices);
    }
 
    if (s.electrode.dump_vtx) {
      set_dir(OUTPUT);
      s.dump_vtx_file(i);
    }
 
    if(param_globals::stim[i].pulse.dumpTrace && get_rank() == 0) {
      set_dir(OUTPUT);
      s.pulse.wave.write_trace(s.name+".trc");
    }
 
  }
}
 
void EMI::apply_dbc_stimulus()
{
  parabolic_solver_emi& ps = parab_solver;
 
  // Rebuild only if the active DBC set changed. Time-dependent DBC values with
  // the same constrained DOFs are applied through enforce_dbc_rhs().
  bool dbcs_have_updated = ps.dbc != nullptr && ps.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) {
    parab_solver.rebuild_matrices(mtype_vol, *ion.miif, stimuli, logger);
  }
}
 
void EMI::apply_current_stimulus()
{
  parabolic_solver_emi& ps = parab_solver;
  ps.Iij_stim->set(0.0);
 
  // iterate over stimuli
  for(stimulus & s : stimuli) {
    if(s.is_active()) {
      switch (s.phys.type) {
        case I_tm: {
          apply_stim_to_vector(s, *ps.Iij_temp, true);
          ps.Bi->mult(*ps.Iij_temp, *ps.Ib_both_face);
          ps.operator_both_to_unique_face->mult(*ps.Ib_both_face, *ps.Ib_unique_face);
          ps.Ib->add_scaled(*ps.Ib_unique_face, -0.5); // compensate the two-sided contribution produced by Bi.
        } break;
 
        case I_ex:
        case I_in: {
          apply_stim_to_vector(s, *ps.Iij_stim, true);
        } break;
 
        default: break;
      }
    }
  }
}
 
void EMI::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;
 
      log_msg(NULL, 0, 0, "Balancing stimulus %d with %d %s-wise.", from, to,
              is_current(stimuli[from].phys.type) ? "current" : "voltage");
 
      stimulus& s_from = stimuli[from];
      stimulus& s_to   = stimuli[to];
 
      s_to.pulse = s_from.pulse;
      s_to.ptcl  = s_from.ptcl;
      s_to.phys  = s_from.phys;
      s_to.pulse.strength *= -1.0;
 
      if (s_from.phys.type == I_ex || s_from.phys.type == I_in) {
        // if from is total current, skip volume based adjustment of strength
        // otherwise, scale_total_stimulus_current() will undo the balanced scaling of to.pulse.strength
        // scale_total_stimulus_current() will do the scaling based on the volume
        if (!s_from.phys.total_current) {
          sf_mat& mass = *parab_solver.mass_emi;
          SF_real vol0 = get_volume_from_nodes(mass, s_from.electrode.vertices);
          SF_real vol1 = get_volume_from_nodes(mass, s_to.electrode.vertices);
 
          s_to.pulse.strength *= fabs(vol0 / vol1);
        }
      }
    }
  }
}
 
void EMI::scale_total_stimulus_current(SF::vector<stimulus>& stimuli,
                                       sf_mat&               mass_vol,
                                       sf_mat&               mass_surf,
                                       FILE_SPEC             logger)
{
  for (stimulus & s : stimuli){
    if(is_current(s.phys.type) && s.phys.total_current){
      switch (s.phys.type) {
        case I_in:
        case I_ex: {
          // compute affected volume in um^3
          SF_real vol = get_volume_from_nodes(mass_vol, 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/I_in there is an additional um^3 to cm^3 scaling in phys.scale,
          // since I_ex/I_in is expected to be in uA/cm^3. Therefore, we need to compensate for that to arrive at uA later.
          assert(vol > 0);
          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);
        } break;
 
        case I_tm: {
          // In the EMI model, I_tm only affects the membrane. Therefore, we compute
          // the affected membrane surface in um^2 using the membrane mass matrix.
          // The electrode vertices are resticted to the membrane during setup, hence this function returns a surface already.
          SF_real surf = get_volume_from_nodes(mass_surf, s.electrode.vertices);
 
          // convert to cm^2
          assert(surf > 0);
          surf /= 1.e8;
 
          // 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);
        } break;
 
        default: break;
      }
    }
  }
}
 
// Assign a deterministic global element numbering for a surface mesh based on
// element type, tag, and sorted node ids. This only affects output ordering.
static void assign_deterministic_elem_numbering(sf_mesh & mesh)
{
  const int KEY_SIZE = 6; // type, tag, n1, n2, n3, n4
  int rank = 0, size = 0;
  MPI_Comm_rank(mesh.comm, &rank);
  MPI_Comm_size(mesh.comm, &size);
 
  auto make_key = [&](size_t i) {
    std::array<mesh_int_t, KEY_SIZE> k;
    k.fill(-1);
    k[0] = static_cast<mesh_int_t>(mesh.type[i]);
    k[1] = mesh.tag[i];
 
    int nn = 0;
    if (mesh.type[i] == SF::Line) nn = 2;
    else if (mesh.type[i] == SF::Tri) nn = 3;
    else if (mesh.type[i] == SF::Quad) nn = 4;
    else nn = 0;
 
    std::vector<mesh_int_t> nodes;
    nodes.reserve(nn);
    size_t off = mesh.dsp[i];
    for (int j = 0; j < nn; j++) {
      nodes.push_back(mesh.con[off + j]);
    }
    std::sort(nodes.begin(), nodes.end());
    for (int j = 0; j < (int)nodes.size(); j++) {
      k[2 + j] = nodes[j];
    }
    return k;
  };
 
  // Pack local keys
  std::vector<mesh_int_t> local_keys(mesh.l_numelem * KEY_SIZE, -1);
  for (size_t i = 0; i < mesh.l_numelem; i++) {
    auto k = make_key(i);
    for (int j = 0; j < KEY_SIZE; j++) local_keys[i * KEY_SIZE + j] = k[j];
  }
 
  // Gather sizes
  std::vector<int> counts(size, 0), displs(size, 0);
  int local_count = (int)local_keys.size();
  MPI_Allgather(&local_count, 1, MPI_INT, counts.data(), 1, MPI_INT, mesh.comm);
  int total = 0;
  for (int r = 0; r < size; r++) {
    displs[r] = total;
    total += counts[r];
  }
 
  std::vector<mesh_int_t> all_keys;
  if (rank == 0) all_keys.resize(total, -1);
  const MPI_Datatype key_mpi_t = mpi_datatype<mesh_int_t>();
  MPI_Gatherv(local_keys.data(), local_count, key_mpi_t,
              rank == 0 ? all_keys.data() : nullptr, counts.data(), displs.data(), key_mpi_t,
              0, mesh.comm);
 
  // Build sorted global keys on rank 0
  std::vector<std::array<mesh_int_t, KEY_SIZE>> sorted_keys;
  if (rank == 0) {
    const int nkeys = total / KEY_SIZE;
    sorted_keys.resize(nkeys);
    for (int i = 0; i < nkeys; i++) {
      std::array<mesh_int_t, KEY_SIZE> k;
      for (int j = 0; j < KEY_SIZE; j++) k[j] = all_keys[i * KEY_SIZE + j];
      sorted_keys[i] = k;
    }
    std::sort(sorted_keys.begin(), sorted_keys.end());
  }
 
  // Broadcast sorted keys
  int nkeys = 0;
  if (rank == 0) nkeys = (int)sorted_keys.size();
  MPI_Bcast(&nkeys, 1, MPI_INT, 0, mesh.comm);
  std::vector<mesh_int_t> flat_sorted(nkeys * KEY_SIZE, -1);
  if (rank == 0) {
    for (int i = 0; i < nkeys; i++) {
      for (int j = 0; j < KEY_SIZE; j++) flat_sorted[i * KEY_SIZE + j] = sorted_keys[i][j];
    }
  }
  MPI_Bcast(flat_sorted.data(), (int)flat_sorted.size(), key_mpi_t, 0, mesh.comm);
 
  // Reconstruct sorted_keys on all ranks
  if (rank != 0) {
    sorted_keys.resize(nkeys);
    for (int i = 0; i < nkeys; i++) {
      std::array<mesh_int_t, KEY_SIZE> k;
      for (int j = 0; j < KEY_SIZE; j++) k[j] = flat_sorted[i * KEY_SIZE + j];
      sorted_keys[i] = k;
    }
  }
 
  // Assign deterministic element numbering (both REF and SUBMESH)
  SF::vector<mesh_int_t> & nbr_ref = mesh.register_numbering(SF::NBR_ELEM_REF);
  SF::vector<mesh_int_t> & nbr_sub = mesh.register_numbering(SF::NBR_ELEM_SUBMESH);
  nbr_ref.resize(mesh.l_numelem);
  nbr_sub.resize(mesh.l_numelem);
  for (size_t i = 0; i < mesh.l_numelem; i++) {
    auto k = make_key(i);
    auto it = std::lower_bound(sorted_keys.begin(), sorted_keys.end(), k);
    if (it == sorted_keys.end() || *it != k) {
      log_msg(0, 5, 0, "deterministic numbering failed to find key (rank %d, elem %zu)", rank, i);
      EXIT(1);
    }
    mesh_int_t gid = (mesh_int_t)(it - sorted_keys.begin());
    nbr_ref[i] = gid;
    nbr_sub[i] = gid;
  }
}
 
void EMI::setup_output()
{
  std::string output_base = get_basename(param_globals::meshname);
 
  set_dir(INPUT);
  const bool write_binary =
      SF::fileExists(std::string(param_globals::meshname) + ".belem") ||
      SF::fileExists(std::string(param_globals::meshname) + ".bpts");
 
  set_dir(OUTPUT);
 
  const int gridout_emi = param_globals::gridout_emi;
 
  // write entire mesh
  sf_mesh & mesh = get_mesh(emi_msh);
  if(gridout_emi & 2) {
    std::string output_file = output_base + "_e";
    log_msg(0, 0, 0, "Writing \"%s\" mesh: %s (%s)", mesh.name.c_str(), output_file.c_str(), write_binary ? "binary" : "text");
    const double t0 = MPI_Wtime();
    write_mesh_parallel(mesh, write_binary, output_file.c_str());
    log_msg(0, 0, 0, "Wrote \"%s\" mesh in %.5f seconds.", mesh.name.c_str(), float(MPI_Wtime() - t0));
  }
  else if(param_globals::output_level > 1) {
    log_msg(0, 0, 0, "Skipping \"%s\" mesh output.", mesh.name.c_str());
  }
  // register output for overall phi on the entire mesh
  output_manager.register_output(parab_solver.ui, emi_msh, 1, param_globals::phiefile, "mV", NULL);
 
  sf_mesh & mesh_m = get_mesh(emi_surface_unique_face_msh);
  mesh_m.name = "Membrane";
  if(gridout_emi & 1) {
    std::string output_file = output_base + "_m";
    log_msg(0, 0, 0, "Writing \"%s\" mesh: %s (%s)", mesh_m.name.c_str(), output_file.c_str(), write_binary ? "binary" : "text");
    const double t0 = MPI_Wtime();
    write_mesh_parallel(mesh_m, write_binary, output_file.c_str());
    log_msg(0, 0, 0, "Wrote \"%s\" mesh in %.5f seconds.", mesh_m.name.c_str(), float(MPI_Wtime() - t0));
  }
  else if(param_globals::output_level > 1) {
    log_msg(0, 0, 0, "Skipping \"%s\" mesh output.", mesh_m.name.c_str());
  }
  // register output for Vm on membrane interface
  // ensure deterministic element ordering in output
  if(get_permutation(emi_surface_unique_face_msh, ELEM_PETSC_TO_CANONICAL, 1) == NULL) {
    register_permutation(emi_surface_unique_face_msh, ELEM_PETSC_TO_CANONICAL, 1);
  }
  output_manager.register_output(parab_solver.vb_unique_face, emi_surface_unique_face_msh, 1, param_globals::vofile, "mV", NULL, true);
 
  if(param_globals::num_trace) {
    sf_mesh & imesh = get_mesh(emi_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 EMI::dump_matrices()
{
  std::string bsname = param_globals::dump_basename;
  std::string fn;
 
  set_dir(OUTPUT);
 
  fn = bsname + "_lhs.bin";
  parab_solver.lhs_emi->write(fn.c_str());
 
  fn = bsname + "_K.bin";
  parab_solver.stiffness_emi->write(fn.c_str());
 
  fn = bsname + "_B.bin";
  parab_solver.B->write(fn.c_str());
 
  fn = bsname + "_Bi.bin";
  parab_solver.Bi->write(fn.c_str());
 
  fn = bsname + "_BsM.bin";
  parab_solver.BsM->write(fn.c_str());
 
  fn = bsname + "_M.bin";
  parab_solver.mass_emi->write(fn.c_str());
 
  fn = bsname + "_Ms.bin";
  parab_solver.mass_surf_emi->write(fn.c_str());
 
}
 
double EMI::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 EMI::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;
}
 
void EMI::setup_solvers()
{
  set_dir(OUTPUT);
  const int log_flag = param_globals::output_level > 1 ? ECHO : 0;
  double t0 = MPI_Wtime();
  parab_solver.init();
  log_msg(logger, 0, log_flag, "EMI setup_solvers: parabolic solver init in %.5f seconds.", float(MPI_Wtime() - t0));
  t0 = MPI_Wtime();
  parab_solver.rebuild_matrices(mtype_vol, *ion.miif, stimuli, logger);
  log_msg(logger, 0, log_flag, "EMI setup_solvers: matrix assembly and linear solver setup in %.5f seconds.", float(MPI_Wtime() - t0));
 
  if(param_globals::dump2MatLab)
    dump_matrices();
}
 
void extract_unique_tag(SF::vector<mesh_int_t>& unique_tags)
{
  MPI_Comm comm = SF_COMM;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  binary_sort(unique_tags);
  unique_resize(unique_tags); // unique_resize is done locally(for each rank)
  make_global(unique_tags, comm);
  binary_sort(unique_tags);
  unique_resize(unique_tags);
}
 
void compute_tags_per_rank(int num_tags, SF::vector<mesh_int_t>& num_tags_per_rank)
{
  MPI_Comm comm = SF_COMM;
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
  divide(num_tags, size, num_tags_per_rank);
}
 
void EMI::setup_EMI_mesh() 
{
  log_msg(0,0,0, "\n    *** Processing EMI mesh ***\n");
 
  const std::string basename = param_globals::meshname;
  const int verb = param_globals::output_level;
  std::map<mesh_t, sf_mesh> & mesh_registry = user_globals::mesh_reg;
  assert(mesh_registry.count(emi_msh) == 1); 
 
  set_dir(INPUT);
 
  sf_mesh & emi_mesh = mesh_registry[emi_msh];
  sf_mesh & emi_surfmesh_one_side = mesh_registry[emi_surface_msh];
  sf_mesh & emi_surfmesh_w_counter_face = mesh_registry[emi_surface_counter_msh];
  sf_mesh & emi_surfmesh_unique_face = mesh_registry[emi_surface_unique_face_msh];
 
  MPI_Comm comm = emi_mesh.comm;
 
  int size, rank;
  double t1, t2, s1, s2;
  const double total_setup_t0 = MPI_Wtime();
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  //-----------------------------------------------------------------
  // Validate mesh dimension: EMI only supports 3D volumetric meshes
  //-----------------------------------------------------------------
  if (emi_mesh.l_numelem > 0) {
    // Check first local element type (all elements should be same type in CARP)
    SF::elem_t first_elem = emi_mesh.type[0];
    
    if (first_elem == SF::Line || first_elem == SF::Tri || first_elem == SF::Quad) {
      const char* type_name = (first_elem == SF::Line) ? "1D (Line)" :
                             (first_elem == SF::Tri)  ? "2D (Tri)" :
                                                         "2D (Quad)";
      if (rank == 0) {
        log_msg(0, 5, 0, "\n*** ERROR: EMI model requires a 3D volumetric mesh!");
        log_msg(0, 5, 0, "*** Current mesh element type: %s", type_name);
        log_msg(0, 5, 0, "*** EMI only supports 3D element types: Tetra, Pyramid, Prism, Hexa");
        log_msg(0, 5, 0, "*** Please provide a 3D mesh with volume elements.\n");
      }
      EXIT(EXIT_FAILURE);
    }
  }
 
  //-----------------------------------------------------------------
  // Step 1: READ *.intra and *.extra 
  //-----------------------------------------------------------------
  int total_num_tags = 0;
  if(verb) log_msg(NULL, 0, 0,"\nReading tags for extra and intra regions from input files");
  t1 = MPI_Wtime();
  {
    SF::vector<mesh_int_t> unique_extra_tags;
    SF::vector<mesh_int_t> unique_intra_tags;
 
    if(verb) log_msg(NULL, 0, 0,"Read extracellular tags");
    read_indices_global(unique_extra_tags,basename+".extra", comm);
    for(mesh_int_t tag:unique_extra_tags){
      parab_solver.extra_tags.insert(tag);
    }
 
    if(verb) log_msg(NULL, 0, 0,"Read intracellular tags");
    read_indices_global(unique_intra_tags,basename+".intra", comm);
    for(mesh_int_t tag:unique_intra_tags){
      parab_solver.intra_tags.insert(tag);
    }
 
    total_num_tags = parab_solver.extra_tags.size() + parab_solver.intra_tags.size();
    if(total_num_tags < size){
      log_msg(0,5,0, "\nThe number of unique tags on EMI mesh is smaller than number of processors!");
      EXIT(EXIT_FAILURE);
    }
    if(verb) log_msg(NULL, 0, 0,"\nextra_tags=%lu, intra_tags=%lu", parab_solver.extra_tags.size(), parab_solver.intra_tags.size());
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 2: READ points from *.pts
  //-----------------------------------------------------------------
  SF::vector<mesh_real_t> pts;
  SF::vector<mesh_int_t>  ptsidx;
  SF::vector<mesh_int_t>  ptsData;
  if(verb) log_msg(NULL, 0, 0,"\nReading points with data on each vertex");
  t1 = MPI_Wtime();
  SF::read_points(basename, comm, pts, ptsidx);
  ptsData.resize(ptsidx.size());
  assert(ptsidx.size()==ptsData.size());
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  std::list< sf_mesh* > meshlist;
  meshlist.push_back(&emi_mesh);
  
  //-----------------------------------------------------------------
  // Step 3: Distribute mesh based on tag or *.part
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0,"\nDistribute mesh based on tags");
  // should be replaced by scotch/pt-scotch for efficiency
  t1 = MPI_Wtime(); 
  distribute_elements_based_tags(emi_mesh, total_num_tags);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 4: insert coordinates to the mesh
  //-----------------------------------------------------------------
  // insert points
  if(verb) log_msg(NULL, 0, 0, "\nInserting points");
  t1 = MPI_Wtime();
  SF::insert_points(pts, ptsidx, meshlist);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 5: extract ptsData based on the location of each vertex
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "\nCompute location of all DOFs on EMI mesh (inner DOFs, membrane, gap junction) saved into ptsData from original mesh");
  t1 = MPI_Wtime();
  compute_ptsdata_from_original_mesh( emi_mesh, 
                                      SF::NBR_REF,
                                      vertex2ptsdata, 
                                      parab_solver.extra_tags,
                                      parab_solver.intra_tags);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 6: extract faces and counterfaces on the interfaces (mem/gap) 
  // and mark the unique faces on surface mesh
  // generate required mapping between surface mesh and unique faces
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "\nExtract EMI surface mesh");
  hashmap::unordered_map<mesh_int_t, SF::emi_index_rank<mesh_int_t>> unused_map_elem_oneface_to_elem_uniqueFace;
  extract_face_based_tags(emi_mesh, SF::NBR_REF, vertex2ptsdata,
                          parab_solver.line_face, 
                          parab_solver.tri_face, 
                          parab_solver.quad_face, 
                          parab_solver.extra_tags, 
                          parab_solver.intra_tags, 
                          emi_surfmesh_one_side, emi_surfmesh_w_counter_face, emi_surfmesh_unique_face,
                          parab_solver.map_elem_uniqueFace_to_elem_oneface,
                          unused_map_elem_oneface_to_elem_uniqueFace);
  meshlist.push_back(&emi_surfmesh_one_side);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
  //-----------------------------------------------------------------
  // Step 7: register surface meshes based on the {typeof}_face with counter face
  //-----------------------------------------------------------------
  compute_surface_mesh_with_counter_face(emi_surfmesh_w_counter_face, SF::NBR_REF,
                                         parab_solver.line_face, 
                                         parab_solver.tri_face, 
                                         parab_solver.quad_face);
 
  compute_surface_mesh_with_unique_face(emi_surfmesh_unique_face, SF::NBR_REF,
                                         parab_solver.line_face, 
                                         parab_solver.tri_face, 
                                         parab_solver.quad_face,
                                         parab_solver.map_elem_uniqueFace_to_tags);
 
  //-----------------------------------------------------------------
  // Step 8: create a map between emi_surfmesh_w_counter_face and emi_surfmesh (both -> one)
  //-----------------------------------------------------------------
  SF::create_reverse_elem_mapping_between_surface_meshes(parab_solver.line_face,
                          parab_solver.tri_face,
                          parab_solver.quad_face,
                          parab_solver.vec_both_to_one_face,
                          comm);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 9: global number of interfaces only on emi_surfmesh without counter face
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "\ncompute global number of interface");
  size_t global_count_surf = 0;
  size_t numelem_surface = emi_surfmesh_one_side.l_numelem;
  size_t local_count_surf = numelem_surface;
  MPI_Reduce(&local_count_surf, &global_count_surf, 1, mpi_datatype<size_t>(), MPI_SUM, 0, MPI_COMM_WORLD);
  if(verb && rank==0) fprintf(stdout, "global number of interfaces = %zu\n", global_count_surf);
 
  //-----------------------------------------------------------------
  // Step 10: submesh_numbering on the emi_mesh and generate parallel layout 
  //-----------------------------------------------------------------
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_mesh);
    emi_mesh.generate_par_layout();
  }
 
  SF::meshdata<mesh_int_t, mesh_real_t>     tmesh_backup_old = emi_mesh;
 
  //-----------------------------------------------------------------
  // Step 11: make a map between key(vertex,tag)-> value(dof) after decoupling interfaces defined on emi mesh
  //-----------------------------------------------------------------
  // During interface decoupling and the introduction of new degrees of freedom (DoFs),
  // we assume that the mesh partitioning is performed at least on a per-tag basis.
  // This means that all elements sharing the same tag number belong to the same rank.
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "\ndecouple emi interfaces");
  if(verb) log_msg(NULL, 0, 0, "\tcompute map oldIdx to dof");
  compute_map_vertex_to_dof(emi_mesh, SF::NBR_REF, vertex2ptsdata, parab_solver.extra_tags, parab_solver.map_vertex_tag_to_dof);
  
  //-----------------------------------------------------------------
  // Step 12: complete the map between key(vertex,tag)-> value(dof) for counter faces defined on emi mesh
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "\tcomplete map oldIdx to dof with counter interface");
  // add to the map the counter part of the interface
  complete_map_vertex_to_dof_with_counter_face(parab_solver.line_face, parab_solver.tri_face, parab_solver.quad_face, parab_solver.map_vertex_tag_to_dof);
  
  //-----------------------------------------------------------------
  // Step 13: update EMI mesh with new DOFs, so interface decoupling is applied on the EMI mesh
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "\tupdate mesh with dof");
  update_emi_mesh_with_dofs(emi_mesh, SF::NBR_REF, parab_solver.map_vertex_tag_to_dof, parab_solver.dof2vertex);  
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 14: initialize map_vertex_tag_to_dof_petsc where the key is (vertex,tag) and the value is (dof,petsc)
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "\nInitialize petsc =0 for  map<oldIdx,tag> -><dof, petsc>");
  // Iterate over map to assign the (oldIndx, tag) -> (dof, petsc) atm petsc = 0
  for(const auto & key_value : parab_solver.map_vertex_tag_to_dof) 
  {
    mesh_int_t gIndex_old  = key_value.first.first;
    mesh_int_t tag_old  = key_value.first.second;
    mesh_int_t dof = key_value.second;
 
    std::pair <mesh_int_t,mesh_int_t> dof_petsc = std::make_pair(dof,-1); // petsc numbering ...
    parab_solver.map_vertex_tag_to_dof_petsc.insert({key_value.first,dof_petsc});
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 15: insert the coordinates on new DOFs to the EMI mesh
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Inserting points and ptsData of dofs to emi_mesh");
  insert_points_ptsData_to_dof(tmesh_backup_old, emi_mesh, SF::NBR_REF, parab_solver.dof2vertex, vertex2ptsdata, dof2ptsData);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 16: submesh_numbering on the emi_mesh and generate parallel layout after decoupling the interfaces on the EMI mesh
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Generating unique PETSc numberings");
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_mesh);
    emi_mesh.generate_par_layout(); 
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
  
  //-----------------------------------------------------------------
  // Step 17: register petsc_numbering on EMI mesh
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Generating unique PETSc numberings");
  {
    SF::petsc_numbering<mesh_int_t,mesh_real_t> petsc_numbering(emi_mesh.pl);
    petsc_numbering(emi_mesh);
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 18: complete map_vertex_tag_to_dof_petsc where the key is (vertex,tag) and the value is (dof,petsc) for counter faces
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "Updating the map between indices to PETSc numberings");
  t1 = MPI_Wtime();
  update_map_indices_to_petsc(emi_mesh, SF::NBR_REF, SF::NBR_PETSC, parab_solver.extra_tags, parab_solver.map_vertex_tag_to_dof_petsc, parab_solver.dof2vertex, parab_solver.elemTag_emi_mesh);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 19: update values of map {type_of}_face with new DOFs after decoupling
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Updating surface mesh with dof");
  update_faces_on_surface_mesh_after_decoupling_with_dofs(emi_surfmesh_one_side, parab_solver.map_vertex_tag_to_dof, parab_solver.line_face, parab_solver.tri_face, parab_solver.quad_face);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 20: register emi_surfmesh for SF::NBR_ELEM_REF, NBR_REF, SF::NBR_SUBMESH
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Layout for element of EMI surfmesh");
  SF::vector<mesh_int_t> & emi_surfmesh_elem =  emi_surfmesh_one_side.register_numbering(SF::NBR_ELEM_REF);
  SF::vector<long int> layout;
  SF::layout_from_count<long int>(emi_surfmesh_one_side.l_numelem, layout, emi_surfmesh_one_side.comm);
  size_t count = layout[rank+1] - layout[rank];
  emi_surfmesh_elem.resize(count);
  for (int i = 0; i < count; ++i){
    emi_surfmesh_elem[i] = layout[rank]+i;
  }
 
  emi_surfmesh_one_side.localize(SF::NBR_REF);
  emi_surfmesh_one_side.register_numbering(SF::NBR_SUBMESH);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 21: register emi_surfmesh_w_counter_face for SF::NBR_ELEM_REF, NBR_REF, SF::NBR_SUBMESH
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Layout for element of EMI surfmesh w counter face");
  SF::vector<mesh_int_t> & emi_surfmesh_counter_elem =  emi_surfmesh_w_counter_face.register_numbering(SF::NBR_ELEM_REF);
  SF::vector<long int> layout_counter;
  SF::layout_from_count<long int>(emi_surfmesh_w_counter_face.l_numelem, layout_counter, emi_surfmesh_w_counter_face.comm);
  size_t count_counter = layout_counter[rank+1] - layout_counter[rank];
  emi_surfmesh_counter_elem.resize(count_counter);
  for (int i = 0; i < count_counter; ++i){
    emi_surfmesh_counter_elem[i] = layout_counter[rank]+i;
  }
  emi_surfmesh_w_counter_face.localize(SF::NBR_REF);
  emi_surfmesh_w_counter_face.register_numbering(SF::NBR_SUBMESH);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
 
  //-----------------------------------------------------------------
  // Step 22: register emi_surfmesh_unique_face for SF::NBR_ELEM_REF, NBR_REF, SF::NBR_SUBMESH
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Layout for element of EMI unique-face surfmesh");
  SF::vector<mesh_int_t> & emi_surfmesh_unique_elem =  emi_surfmesh_unique_face.register_numbering(SF::NBR_ELEM_REF);
  SF::vector<long int> layout_unique;
  SF::layout_from_count<long int>(emi_surfmesh_unique_face.l_numelem, layout_unique, emi_surfmesh_unique_face.comm);
  size_t count_unique = layout_unique[rank+1] - layout_unique[rank];
  emi_surfmesh_unique_elem.resize(count_unique);
  for (int i = 0; i < count_unique; ++i){
    emi_surfmesh_unique_elem[i] = layout_unique[rank]+i;
  }
  emi_surfmesh_unique_face.localize(SF::NBR_REF);
  emi_surfmesh_unique_face.register_numbering(SF::NBR_SUBMESH);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 23: insert coordinates to one-sided, both-face, and unique-face surface meshes
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Inserting points to EMI surfmesh");
  insert_points_to_surface_mesh(tmesh_backup_old, emi_surfmesh_one_side, SF::NBR_REF, parab_solver.dof2vertex, parab_solver.extra_tags, parab_solver.elemTag_surface_mesh);
  insert_points_to_surface_mesh(tmesh_backup_old, emi_surfmesh_w_counter_face, SF::NBR_REF, parab_solver.dof2vertex, parab_solver.extra_tags, parab_solver.elemTag_surface_w_counter_mesh);
  insert_points_to_surface_mesh(tmesh_backup_old, emi_surfmesh_unique_face, SF::NBR_REF, parab_solver.dof2vertex);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 24: submesh_numbering and PETSc numbering for all EMI surface meshes
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Generating submesh_numbering and PETSc numberings for surface mesh");
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_surfmesh_one_side);
    emi_surfmesh_one_side.generate_par_layout(); 
 
    SF::petsc_numbering<mesh_int_t,mesh_real_t> petsc_numbering(emi_surfmesh_one_side.pl);
    petsc_numbering(emi_surfmesh_one_side);
  }
 
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_surfmesh_w_counter_face);
    emi_surfmesh_w_counter_face.generate_par_layout();
 
    SF::petsc_numbering<mesh_int_t,mesh_real_t> petsc_numbering(emi_surfmesh_w_counter_face.pl);
    petsc_numbering(emi_surfmesh_w_counter_face);
  }
 
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_surfmesh_unique_face);
    emi_surfmesh_unique_face.generate_par_layout();
 
    SF::petsc_numbering<mesh_int_t,mesh_real_t> petsc_numbering(emi_surfmesh_unique_face.pl);
    petsc_numbering(emi_surfmesh_unique_face);
    // Ensure deterministic element ordering across MPI for output
    assign_deterministic_elem_numbering(emi_surfmesh_unique_face);
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 25: assign PETSc numbering to map_vertex_tag_to_dof_petsc where the key is (vertex,tag) and value is (dof,petsc)
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "assign PETSc numbering for new faces");
  SF::assign_petsc_on_counter_face(parab_solver.map_vertex_tag_to_dof_petsc,comm);
  {
    hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, std::pair<mesh_int_t,mesh_int_t>>::iterator it;
    for (it = parab_solver.map_vertex_tag_to_dof_petsc.begin(); it != parab_solver.map_vertex_tag_to_dof_petsc.end(); it++)
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old = it->first;
      std::pair <mesh_int_t,mesh_int_t> dof_petsc = it->second;
      parab_solver.dof2petsc[dof_petsc.first] = dof_petsc.second;
      parab_solver.petsc2dof[dof_petsc.second] = dof_petsc.first;
    }
  }
  
  // DEBUG: Check for invalid PETSc indices after exchange
  #ifdef EMI_DEBUG_MESH
  {
    int invalid_count = 0;
    for (const auto& [key, val] : parab_solver.map_vertex_tag_to_dof_petsc) {
      if (val.second < 0) {
        invalid_count++;
        if (invalid_count <= 3) {
          fprintf(stderr, "RANK %d INVALID: vertex=%ld tag=%ld dof=%ld petsc=%ld\n", 
                  rank, (long)key.first, (long)key.second, (long)val.first, (long)val.second);
        }
      }
    }
    fprintf(stderr, "RANK %d: After Step 26: %d invalid PETSc indices out of %zu total\n", 
            rank, invalid_count, parab_solver.map_vertex_tag_to_dof_petsc.size());
    fflush(stderr);
  }
  #endif
  
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
  
  //-----------------------------------------------------------------
  // Step 26: update counter faces in map{type_of}_face
  //-----------------------------------------------------------------
  added_counter_faces_to_map(parab_solver.line_face, parab_solver.tri_face, parab_solver.quad_face);
  const double total_setup = MPI_Wtime() - total_setup_t0;
  log_msg(0,0,0, "Total setup_EMI_mesh processing time: %.5f sec.", float(total_setup));
  
  log_msg(0,0,0, "\n    *** EMI mesh processing Done ***\n");
}
 
void distribute_elements_based_tags(SF::meshdata<mesh_int_t, mesh_real_t>& mesh,
                                    int total_num_tags)
{
  MPI_Comm comm = mesh.comm;
  int size, rank;
  double t1, t2, s1, s2;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
  const int verb = param_globals::output_level;
 
  if(total_num_tags < size)
  {
    PetscPrintf(PETSC_COMM_WORLD,"\nThe number of processors should be less than the number of tags, size = %d & ntags = %d !!!\n",
                size, total_num_tags);
    cleanup_and_exit();
  }
 
  if(verb==10) log_msg(NULL, 0, 0,"\ncompute the number of tags which belongs to one rank");
  // compute the number of tags per rank
  SF::vector<mesh_int_t> ntags_per_rank;
  t1 = MPI_Wtime();
  compute_tags_per_rank(total_num_tags, ntags_per_rank); 
  t2 = MPI_Wtime();
  if(verb==10) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  // generate global destination w.r.t. tags
  // redistribute elements based on new partition 
  t1 = MPI_Wtime(); 
  SF::vector<mesh_int_t> part_based_Tags(mesh.l_numelem);
  partition_based_tags(total_num_tags, mesh.tag, ntags_per_rank, part_based_Tags);
  SF::redistribute_elements(mesh,part_based_Tags);
  // permute elements locally first based on the tag then element index
  permute_mesh_locally_based_on_tag_elemIdx(mesh);
  t2 = MPI_Wtime(); 
  if(verb==10) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
}
 
// In the EMI model we cannot have partition boundaries crossing cells.
// Therefore this function assigns a partition to each tag. All elements with this tag will go to the same partition.
void partition_based_tags(int num_tags,
                          SF::vector<mesh_int_t> tag,
                          SF::vector<mesh_int_t> num_tags_per_rank,
                          SF::vector<mesh_int_t>& part)
{
  const int verb = param_globals::output_level;
  MPI_Comm comm = SF_COMM;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  // we need to have tags_to_rank as a map globally.
  hashmap::unordered_map<mesh_int_t, mesh_int_t> tags_to_rank_map; 
  // Try to load the mapping from a .part file and fall back on the internal mapping if there is no such file.
  // If the mapping is read from file, the num_tags_per_rank argument is not used.
  if (!load_partitions_from_file(tags_to_rank_map, num_tags, comm)) {
    if(verb==10) log_msg(NULL, 0, 0,"\ncompute the number of unique tags");
    SF::vector<mesh_int_t> unique_tags = tag;
    double t1 = MPI_Wtime();
    extract_unique_tag(unique_tags);
    double t2 = MPI_Wtime();
    if(verb==10) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
    if (unique_tags.size() != static_cast<size_t>(num_tags)) {
      log_msg(0,5,0, "\nerror: the number of tags in the EMI mesh does not match the total tags from *.extra and *.intra");
      EXIT(1);
    }
    map_tags_to_rank(size, unique_tags, num_tags_per_rank, tags_to_rank_map);
  }
 
  for (size_t i = 0; i < part.size(); ++i) {
    if(tags_to_rank_map.count(tag[i]))
      part[i] = tags_to_rank_map[tag[i]];
  }
}
 
void map_tags_to_rank(int size, const SF::vector<mesh_int_t> & unique_tags, const SF::vector<mesh_int_t> & num_tags_per_rank, hashmap::unordered_map<mesh_int_t, mesh_int_t> &tags_to_rank_map)
{
  // the goal is to map the tag as a key -> rank, should be global one
  mesh_int_t t = 0;  // tag index
  for (size_t r = 0; r < size; ++r) {
    for (size_t count = 0; count < num_tags_per_rank[r];) {
      mesh_int_t tag =  unique_tags[t];
 
      // // here we consider each tag belongs to only one rank.
      tags_to_rank_map.insert({tag, r});  
      count +=1; 
      t+=1;
    }
  }
}
 
// This function creates a tag-to-rank mapping based only on the numerical order of the tags.
bool load_partitions_from_file(hashmap::unordered_map<mesh_int_t, mesh_int_t>& tags_to_rank_map,
                               int expected_num_tags,
                               MPI_Comm comm)
{
  int size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  SF::vector<int> parts;
  const std::string basename = param_globals::meshname;
  FILE* fd = fopen((basename + ".part").c_str(), "r");
  if (fd != NULL) { // check if the file exists otherwise read_indices_global() crash
    read_indices_global(parts, basename + ".part", comm);
    fclose(fd);
  } 
  
  if (parts.size() == 0) return false;  // file does not exist or is empty
 
  if (parts.size() % 2 != 0) {
    log_msg(0,5,0, "\nThe part file should contain 2 lines per tag, one for the tag number and the next for its associated partition number.!");
    EXIT(1);
  }
 
  int min_part = std::numeric_limits<int>::max();
  int max_part = std::numeric_limits<int>::min();
  for(int i = 0; i < parts.size(); i+=2) {
    min_part = std::min(min_part, parts[i + 1]);
    max_part = std::max(max_part, parts[i + 1]);
    tags_to_rank_map.insert({parts[i], parts[i + 1]});
  }
 
  if (tags_to_rank_map.size() != static_cast<size_t>(expected_num_tags)) {
    log_msg(0,5,0, "\nerror: the number of tags in the .part file does not match the total tags from *.extra and *.intra");
    EXIT(1);
  }
 
  if (min_part < 0 || max_part >= size) {
    log_msg(0,5,0,
            "\nerror: EMI partition file %s.part is incompatible with this run.\n"
            "The file contains partition IDs in [%d, %d], but the current MPI communicator has %d rank(s).\n"
            "Remove/regenerate the .part file or run with a matching number of MPI tasks.",
            basename.c_str(), min_part, max_part, size);
    EXIT(1);
  }
 
  if (rank == 0 && max_part + 1 != size) {
    log_msg(0,3,0,
            "Warning: EMI partition file %s.part uses %d partition ID(s), but the current run uses %d MPI rank(s).",
            basename.c_str(), max_part + 1, size);
  }
 
  return true;
}
 
void permute_mesh_locally_based_on_tag_elemIdx(SF::meshdata<mesh_int_t, mesh_real_t>& mesh)
{
  mesh.globalize(SF::NBR_REF);
  SF::meshdata<mesh_int_t, mesh_real_t>     tmesh = mesh;
  SF::vector<mesh_int_t> perm;
  interval(perm, 0, mesh.tag.size());
 
  SF::vector<mesh_int_t>  tags    = tmesh.tag;
  SF::vector<mesh_int_t>& elemIdx = tmesh.get_numbering(SF::NBR_ELEM_REF);
  binary_sort_sort_copy(tags, elemIdx, perm);
  permute_mesh(tmesh, mesh, perm);
 
  mesh.localize(SF::NBR_REF);
}
 
}  // namespace opencarp
 
#endif