sim_utils.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#include "basics.h"
#include "sim_utils.h"
#include "fem.h"
#include "physics.h"
#include "async_io.h"
#include "SF_init.h"
 
#include <libgen.h>
#include <fstream>
#include <iomanip>
 
 
namespace opencarp {
 
static char  input_dir[1024],        // directory from which to read input
             output_dir[1024],       // directory to which to write results
             postproc_dir[1024],     // postprocessing directory
             current_dir[1024];      // current directory
 
void parse_params_cpy(int argc, char** argv)
{
  // copy the command line parameters that are for carp (i.e. before the first "+")
  int    cpy_argc = argc;
  SF::vector<char*> cpy_argv(cpy_argc, NULL);
 
  for(int i=0; i < cpy_argc; i++) {
    if(strcmp(argv[i], "+") != 0) {
      cpy_argv[i] = dupstr(argv[i]);
    }
    else {
      cpy_argc = i;
      break;
    }
  }
 
  // launch param on the carp command line parameters
  int param_argc = cpy_argc;
  int status;
 
  do {
    status = param(PARAMETERS, &param_argc, cpy_argv.data());
    if ( status==PrMERROR||status==PrMFATAL )
      fprintf( stderr, "\n*** Error reading parameters\n\n");
    else if ( status == PrMQUIT )
      fprintf( stderr, "\n*** Quitting by user's request\n\n");
  } while (status == PrMERROR);
 
  // if --output-setup then loop over PARAMETERS
  //    output value of each PARAMETER
  if (param_globals::output_setup)
    param_output(PARAMETERS, 1);
 
  // free the parameters
  for(int i=0; i<param_argc; i++)
    free(cpy_argv[i]);
 
  // exit on error
  if (status & PrMERROR) exit(status);
 
}
 
 
void register_physics()
{
  // here all the physics can be registered to the physics registry
  // then they should be processed automatically
  if(phys_defined(PHYSREG_LAPLACE))
    user_globals::physics_reg[laplace_phys] = new Laplace();
  if(phys_defined(PHYSREG_EMI))
#if WITH_EMI_MODEL
    user_globals::physics_reg[emi_phys] = new EMI();
#else
    log_msg(NULL, 5, ECHO, "The EMI model was not compiled for this binary.\n");
    exit(EXIT_FAILURE);
#endif
 
  if (phys_defined(PHYSREG_EIKONAL))
    user_globals::physics_reg[eikonal_phys] = new Eikonal();
 
  if(phys_defined(PHYSREG_INTRA_ELEC))
    user_globals::physics_reg[elec_phys] = new Electrics();
}
 
void initialize_physics()
{
  log_msg(0,0,0, "\n    *** Initializing physics ***\n");
 
  //load in the external imp modules
#ifdef HAVE_DLOPEN
  for (int ii = 0; ii < param_globals::num_external_imp; ii++) {
    int loading_succeeded = limpet::load_ionic_module(param_globals::external_imp[ii]);
    assert(loading_succeeded);
  }
#else
  if(param_globals::num_external_imp)
    log_msg(NULL, 4, ECHO,"Loading of external LIMPET modules not enabled.\n"
                          "Recompile with DLOPEN set.\n" );
#endif
 
  // init physics via Basic_physic interface
  for(auto it : user_globals::physics_reg) {
    Basic_physic* p = it.second;
    log_msg(NULL, 0, 0, "Initializing %s ..", p->name);
    p->initialize();
  }
}
 
void destroy_physics()
{
  log_msg(0,0,0, "\n    *** Destroying physics ***\n");
 
  for(auto it : user_globals::physics_reg) {
    Basic_physic* p = it.second;
    log_msg(NULL, 0, 0, "Destroying %s ..", p->name);
    p->destroy();
  }
}
 
// ignore_extracellular stim must be moved to stimulate.cc to be able
// to use all defined set operations instead of defines
 
void ignore_extracellular_stim(Stimulus *st, int ns, int ignore)
{
  // needs to be switch to stim enum types defined in stimulate.h
  for ( int i=0; i<ns; i++ ) {
    int turn_off = 0;
    turn_off += (st[i].stimtype == Extracellular_Ground) && (ignore & NO_EXTRA_GND);
    turn_off += (IsExtraV(st[i])) && (ignore & NO_EXTRA_V);
    turn_off += (st[i].stimtype==Extracellular_I) && (ignore & NO_EXTRA_I);
 
    if (turn_off)  {
      st[i].stimtype = Ignore_Stim;
      log_msg( NULL, 1, 0, "Extracellular stimulus %d ignored for monodomain", i );
    } else if ( st[i].stimtype==Intracellular_I ) {
      st[i].stimtype = Transmembrane_I;
      log_msg( NULL, 1, 0, "Intracellular stimulus %d converted to transmembrane", i );
    }
  }
}
 
int set_ignore_flags( int mode )
{
  if(mode==MONODOMAIN)
    return STM_IGNORE_MONODOMAIN;
  if(mode==BIDOMAIN)
    return STM_IGNORE_BIDOMAIN;
  if(mode==PSEUDO_BIDM)
    return STM_IGNORE_PSEUDO_BIDM;
 
  return IGNORE_NONE;
}
 
 
void check_nullspace_ok()
{
  if(param_globals::floating_ground)
    return;
 
  Stimulus* s = param_globals::stimulus;
 
  for(int i=0; i < param_globals::num_stim; i++) {
    if(s[i].stimtype == Extracellular_Ground ||
       s[i].stimtype == Extracellular_V      ||
       s[i].stimtype == Extracellular_V_OL)
      return;
  }
 
  // for now we only warn, although we should actually stop the run
  log_msg( NULL, 4, 0,"Elliptic system is singular!\n"
                      "Either set floating_ground=1 or use an explicit ground:voltage (stimulus[X].stimtype=3)\n"
                      "Do not trust the elliptic solution of this simulation run!\n");
}
 
void show_build_info()
{
  printf("\n""*** GIT tag:            %s\n",   GIT_COMMIT_TAG);
  printf(    "*** GIT hash:           %s\n",   GIT_COMMIT_HASH);
  printf(    "*** GIT repo:           %s\n",   GIT_PATH);
  printf(    "*** dependency commits: %s\n\n", SUBREPO_COMMITS);
}
 
void check_and_convert_params()
{
  // convert time steps to milliseconds
  param_globals::dt      /= 1000.;
 
  // check parab-solve solution method
  if(param_globals::mass_lumping == 0 && param_globals::parab_solve==0)  {
    log_msg(NULL, 2, ECHO,
            "Warning: explicit solve not possible without mass lumping. \n"
            "Switching to Crank-Nicolson!\n\n");
 
    param_globals::parab_solve = 1;
  }
 
  // check if we have to modify stimuli based on used bidomain setting
  if(!param_globals::extracell_monodomain_stim)
    ignore_extracellular_stim(param_globals::stimulus, param_globals::num_stim,
            set_ignore_flags(param_globals::bidomain));
 
  // check nullspace if necessary
  // if((param_globals::bidomain==BIDOMAIN && !param_globals::ellip_use_pt) ||
  //    (param_globals::bidomain==PSEUDO_BIDM && !param_globals::parab_use_pt))
  //   check_nullspace_ok();
 
  if(param_globals::t_sentinel > 0 && param_globals::sentinel_ID < 0 ) {
    log_msg(0,4,0, "Warning: t_sentinel is set but no sentinel_ID has been specified; check_quiescence() behavior may not be as expected");
  }
 
  if(param_globals::num_external_imp > 0 ) {
    for(int ext_imp_i = 0; ext_imp_i < param_globals::num_external_imp; ext_imp_i++) {
      if(param_globals::external_imp[ext_imp_i][0] != '/') {
        log_msg(0,5,0, "external_imp[%d] error: absolute paths must be used for .so file loading (\'%s\')",
                ext_imp_i, param_globals::external_imp[ext_imp_i]);
        EXIT(1);
      }
    }
  }
 
  if(param_globals::experiment == EXP_LAPLACE && param_globals::bidomain != 1) {
    log_msg(0,4,0, "Warning: Laplace experiment mode requires bidomain = 1. Setting bidomain = 1.");
    param_globals::bidomain = 1;
  }
 
  if(param_globals::num_phys_regions == 0) {
    log_msg(0,4,0, "Warning: No physics region defined! Please set phys_region parameters to correctly define physics.");
 
    if(param_globals::experiment != EXP_LAPLACE) {
      log_msg(0,4,0, "Intra-elec and Extra-elec domains will be derived from fibers.\n");
      param_globals::num_phys_regions = param_globals::bidomain ? 2 : 1;
      param_globals::phys_region = (p_region*) malloc(param_globals::num_phys_regions * sizeof(p_region));
      param_globals::phys_region[0].ptype  = PHYSREG_INTRA_ELEC;
      param_globals::phys_region[0].name   = strdup("Autogenerated intracellular Electrics");
      param_globals::phys_region[0].num_IDs = 0;
 
      if(param_globals::bidomain) {
        param_globals::phys_region[1].ptype   = PHYSREG_EXTRA_ELEC;
        param_globals::phys_region[1].name    = strdup("Autogenerated extracellular Electrics");
        param_globals::phys_region[1].num_IDs = 0;
      }
    } else {
      log_msg(0,4,0, "Laplace domain will be derived from fibers.\n");
      param_globals::num_phys_regions = 1;
      param_globals::phys_region = (p_region*) malloc(param_globals::num_phys_regions * sizeof(p_region));
      param_globals::phys_region[0].ptype  = PHYSREG_LAPLACE;
      param_globals::phys_region[0].name   = strdup("Autogenerated Laplace");
      param_globals::phys_region[0].num_IDs = 0;
    }
  }
 
  if(param_globals::experiment == EXP_LAPLACE && !phys_defined(PHYSREG_LAPLACE)) {
    log_msg(0,4,0, "Warning: Laplace experiment mode requires a laplace physics regions defined.");
 
    int idx = -1;
    if((idx = get_phys_index(PHYSREG_EXTRA_ELEC)) > -1) {
      log_msg(0,4,0, "Converting the defined extracellular-electrics-region to laplace-region.");
      param_globals::phys_region[idx].ptype = PHYSREG_LAPLACE;
    } else if ((idx = get_phys_index(PHYSREG_INTRA_ELEC)) > -1) {
      log_msg(0,4,0, "Converting the defined intracellular-electrics-region to laplace-region.");
      param_globals::phys_region[idx].ptype = PHYSREG_LAPLACE;
    } else {
      param_globals::num_phys_regions += 1;
      param_globals::phys_region = (p_region*) realloc(param_globals::phys_region, param_globals::num_phys_regions * sizeof(p_region));
 
      param_globals::phys_region[param_globals::num_phys_regions - 1].ptype  = PHYSREG_LAPLACE;
      param_globals::phys_region[param_globals::num_phys_regions - 1].name   = strdup("Autogenerated Laplace");
      param_globals::phys_region[param_globals::num_phys_regions - 1].num_IDs = 0;
    }
  }
 
#ifndef WITH_PARMETIS
  if(param_globals::pstrat == 1) {
    log_msg(0,3,0, "openCARP was built without Parmetis support. Swithing to KDtree.");
    param_globals::pstrat = 2;
  }
#endif
 
  // check if we have the legacy stimuli or the new stimuli defined by the user
  bool legacy_stim_set = false, new_stim_set = false;
 
  for(int i=0; i<param_globals::num_stim; i++) {
    Stimulus & legacy_stim = param_globals::stimulus[i];
    Stim & new_stim = param_globals::stim[i];
 
    if(legacy_stim.stimtype || legacy_stim.strength)
      legacy_stim_set = true;
 
    if(new_stim.crct.type || new_stim.pulse.strength)
      new_stim_set = true;
  }
 
  if(legacy_stim_set || new_stim_set) {
    if(legacy_stim_set && new_stim_set) {
      log_msg(0,4,0, "Warning: Legacy stimuli and default stimuli are defined. Only default stimuli will be used!");
    }
    else if (legacy_stim_set) {
      log_msg(0,1,0, "Warning: Legacy stimuli defined. Please consider switching to stimulus definition \"stim[]\"!");
      user_globals::using_legacy_stimuli = true;
    }
  }
  else {
    log_msg(0,4,0, "Warning: No potential or current stimuli found!");
  }
}
 
void set_io_dirs(char *sim_ID, char *pp_ID, IO_t init)
{
  int flg = 0, err = 0, rank = get_rank();
 
  char *ptr = getcwd(current_dir, 1024);
  if (ptr == NULL) err++;
  ptr = getcwd(input_dir, 1024);
  if (ptr == NULL) err++;
  //if (param_globals::experiment == 4 && post_processing_opts == MECHANIC_POSTPROCESS)
  //  {  sim_ID = param_globals::ppID; param_globals::ppID = "POSTPROC_DIR"; }
 
  // output directory
  if (rank == 0)  {
    if (strcmp(sim_ID, "OUTPUT_DIR")) {
      if (mkdir(sim_ID, 0775)) {             // rwxrwxr-x
        if (errno == EEXIST )  {
          log_msg(NULL, 2, 0, "Output directory exists: %s\n", sim_ID);
        } else  {
          log_msg(NULL, 5, 0, "Unable to make output directory\n");
          flg = 1;
        }
      }
    } else if (mkdir(sim_ID, 0775) &&  errno != EEXIST) {
      log_msg(NULL, 5, 0, "Unable to make output directory\n");
      flg = 1;
    }
  }
 
  // terminate?
  if(get_global(flg, MPI_SUM)) { EXIT(-1); }
 
  err += chdir(sim_ID);
  ptr = getcwd(output_dir, 1024);
  if (ptr == NULL) err++;
 
  // terminate?
  if(get_global(err, MPI_SUM)) { EXIT(-1); }
 
  err += chdir(output_dir);
 
  // postprocessing directory
  if (rank == 0 && (param_globals::experiment==EXP_POSTPROCESS))  {
 
    if (strcmp(param_globals::ppID, "POSTPROC_DIR")) {
      if (mkdir(param_globals::ppID, 0775)) {             // rwxrwxr-x
        if (errno == EEXIST )  {
          log_msg(NULL, 2, ECHO, "Postprocessing directory exists: %s\n\n", param_globals::ppID);
        } else  {
          log_msg(NULL, 5, ECHO, "Unable to make postprocessing directory\n\n");
          flg = 1;
        }
      }
    } else if (mkdir(param_globals::ppID, 0775) &&  errno != EEXIST) {
      log_msg(NULL, 5, ECHO, "Unable to make postprocessing directory\n\n");
      flg = 1;
    }
 
  }
 
  if(get_global(flg, MPI_SUM)) { EXIT(-1); }
 
  err += chdir(param_globals::ppID);
  ptr = getcwd(postproc_dir, 1024);
  if (ptr == NULL) err++;
  err = chdir(output_dir);
  if(get_global(err, MPI_SUM)) { EXIT(-1); }
 
  err = set_dir(init);
  if(get_global(err, MPI_SUM)) { EXIT(-1); }
}
 
bool setup_IO(int argc, char **argv)
{
  bool io_node = false;
  int psize = get_size(), prank = get_rank();
 
  if (param_globals::num_io_nodes > 0) {
    // Can't do async IO with only one core
    if (get_size() == 1) {
      log_msg(NULL, 5, 0, "You cannot run with async IO on only one core.\n");
      EXIT(EXIT_FAILURE);
    }
    // Can't do async IO with more IO cores than compute cores
    if (2 * param_globals::num_io_nodes >= psize) {
      log_msg(NULL, 5, 0, "The number of IO cores be less " "than the number of compute cores.");
      EXIT(EXIT_FAILURE);
    }
#if 0
    if (param_globals::num_PS_nodes && param_globals::num_io_nodes > param_globals::num_PS_nodes) {
      LOG_MSG(NULL, 5, 0,
              "The number of IO cores (%d) should not "
              "exceed the number of PS compute cores (%d).\n",
               param_globals::num_io_nodes, param_globals::num_PS_nodes);
      EXIT(-1);
    }
#endif
    // root IO node is global node 0
    io_node = prank < param_globals::num_io_nodes;
 
    MPI_Comm comm;
    MPI_Comm_split(PETSC_COMM_WORLD, io_node, get_rank(), &comm);
    MPI_Comm_set_name(comm, io_node ? "IO" : "compute");
 
    PETSC_COMM_WORLD = comm;    // either the compute world or IO world
 
    prank = get_rank();
 
    MPI_Intercomm_create(comm, 0, MPI_COMM_WORLD, io_node ? param_globals::num_io_nodes : 0,
                         6667, &user_globals::IO_Intercomm);
 
    if(prank != get_rank(user_globals::IO_Intercomm))
      log_msg(NULL, 4, 0, "Global node %d, Comm rank %d != Intercomm rank %d\n",
              get_rank(MPI_COMM_WORLD), get_rank(PETSC_COMM_WORLD),
              get_rank(user_globals::IO_Intercomm));
  } else
    MPI_Comm_set_name(PETSC_COMM_WORLD, "compute");
 
  set_io_dirs(param_globals::simID, param_globals::ppID, OUTPUT);
 
  if((io_node || !param_globals::num_io_nodes) && !prank)
    output_parameter_file("parameters.par", argc, argv);
 
  return io_node;
}
void update_cwd()
{
  getcwd(current_dir, 1024);
}
 
int set_dir(IO_t dest)
{
  int err;
 
  if      (dest==OUTPUT)   err = chdir(output_dir);
  else if (dest==POSTPROC) err = chdir(postproc_dir);
  else if (dest==CURDIR)   err = chdir(current_dir);
  else                     err = chdir(input_dir);
 
  return err;
}
 
void basic_timer_setup()
{
  // if we restart from a checkpoint, the timer_manager will be notified at a later stage
  double start_time = 0.0;
  user_globals::tm_manager = new timer_manager(param_globals::dt, start_time, param_globals::tend);
 
  double end_time    = param_globals::tend;
  timer_manager & tm = *user_globals::tm_manager;
 
  if(param_globals::experiment == EXP_LAPLACE) {
    tm.initialize_singlestep_timer(tm.time, 0, iotm_console,   "IO (console)",  nullptr);
    tm.initialize_singlestep_timer(tm.time, 0, iotm_state_var, "IO (state vars)", nullptr);
    tm.initialize_singlestep_timer(tm.time, 0, iotm_spacedt,   "IO (spacedt)",  nullptr);
  }
  else {
    tm.initialize_eq_timer(tm.time, end_time, 0, param_globals::timedt,  0, iotm_console,   "IO (console)");
    tm.initialize_eq_timer(tm.time, end_time, 0, param_globals::spacedt, 0, iotm_state_var, "IO (state vars)");
    tm.initialize_eq_timer(tm.time, end_time, 0, param_globals::spacedt, 0, iotm_spacedt,   "IO (spacedt)");
  }
 
  if(param_globals::num_tsav) {
    std::vector<double> trig(param_globals::num_tsav);
    for(size_t i=0; i<trig.size(); i++) trig[i] = param_globals::tsav[i];
 
    tm.initialize_neq_timer(trig, 0, iotm_chkpt_list, "instance checkpointing");
  }
 
  if(param_globals::chkpt_intv)
    tm.initialize_eq_timer(param_globals::chkpt_start, param_globals::chkpt_stop, 0,
                           param_globals::chkpt_intv, 0, iotm_chkpt_intv, "interval checkpointing");
 
  if(param_globals::num_trace)
    tm.initialize_eq_timer(tm.time, end_time, 0, param_globals::tracedt, 0, iotm_trace, "IO (node trace)");
}
 
void get_protocol_column_widths(std::vector<int> & col_width, std::vector<int> & used_timer_ids)
{
  char buff[256];
  const short padding = 4;
  Electrics* elec = (Electrics*) get_physics(elec_phys, false);
 
  do {
    snprintf(buff, sizeof buff, "%.3lf", user_globals::tm_manager->time);
    if(col_width[0] < int(strlen(buff)+padding))
      col_width[0] = strlen(buff)+padding;
 
    snprintf(buff, sizeof buff, "%.3d", user_globals::tm_manager->d_time);
    if(col_width[1] < int(strlen(buff)+padding))
      col_width[1] = strlen(buff)+padding;
 
    int col = 2;
    for (size_t tid = 0; tid < used_timer_ids.size(); tid++)
    {
      int timer_id  = used_timer_ids[tid];
      base_timer* t = user_globals::tm_manager->timers[timer_id];
 
      if(t->d_trigger_dur && elec) {
        // figure out value of signal linked to this timer
        double val = 0.;
 
        // determine timer linked to which physics, for now we deal with electrics only
        val = elec->timer_val(timer_id);
 
        snprintf(buff, sizeof buff, "%.3lf", val);
        if(col_width[col] < int(strlen(buff)+padding))
          col_width[col] = strlen(buff)+padding;
      }
      col++;
    }
 
    // advance time
    user_globals::tm_manager->update_timers();
  } while (!user_globals::tm_manager->elapsed());
 
  user_globals::tm_manager->reset_timers();
}
int plot_protocols(const char *fname)
{
  int   err = {0};
  std::ofstream fh;
  const char* smpl_endl = "\n";
 
  if(!get_rank()) {
    fh.open(fname);
 
    // If we couldn't open the output file stream for writing
    if (!fh) {
        // Print an error and exit
        log_msg(0,5,0,"Protocol file %s could not be opened for writing!\n", fname);
        err = -1;
    }
  }
 
  // broadcast and return if err
  if(get_global(err, MPI_SUM))
    return err;
 
  // only rank 0 writes
  if(!get_rank()) {
 
    // collect timer information, label, short label, unit
    std::vector<std::string> col_labels       = {"time", "tick"};
    std::vector<std::string> col_short_labels = {"A",  "B"};
    std::vector<std::string> col_unit_labels  = {"ms", "--" };
    std::vector<int> col_width = {4, 4};
 
    char c_label = {'C'};
    std::string label = {""};
    std::string unit  = {""};
 
    // here we store the IDs of the timers that we care about. currently this are the IO and TS timers
    // and the electricts timers
    std::vector<int> used_timer_ids;
    std::vector<int> used_stim_ids;
 
    Electrics* elec = (Electrics*) get_physics(elec_phys, false);
    if(elec) {
      int sidx = 0;
      for(const stimulus & s : elec->stimuli) {
        if(s.ptcl.timer_id > -1) {
          base_timer* t = user_globals::tm_manager->timers[s.ptcl.timer_id];
          if(t) {
            used_timer_ids.push_back(s.ptcl.timer_id);
            used_stim_ids.push_back(sidx);
          }
        }
 
        sidx++;
      }
    }
 
    // determine longest timer label
    int mx_llen = 0;
    for (size_t tid = 0; tid < used_timer_ids.size(); tid++)
    {
      int timer_id  = used_timer_ids[tid];
      base_timer* t = user_globals::tm_manager->timers[timer_id];
 
      int llen = strlen(t->name);
      mx_llen = llen > mx_llen ? llen : mx_llen;
    }
 
    for (size_t tid = 0; tid < used_timer_ids.size(); tid++)
    {
      int timer_id  = used_timer_ids[tid];
      base_timer* t = user_globals::tm_manager->timers[timer_id];
 
      col_labels.push_back(t->name);
      label = c_label;
      col_short_labels.push_back(label);
 
      if(elec) {
        // search physics for signals linked to timer
        unit = elec->timer_unit(timer_id);
        if(unit.empty()) unit = "--";
        col_unit_labels.push_back(unit);
        col_width.push_back(4);
      }
      c_label++;
    }
 
    get_protocol_column_widths(col_width, used_timer_ids);
 
    // print header + legend first
    fh << "# Protocol header\n#\n" << "# Legend:\n";
    for(size_t i = 0; i<col_short_labels.size(); i++)
    {
      fh << "#" << std::setw(2) << col_short_labels[i] << " = " << std::setw(mx_llen) << col_labels[i];
      fh << " [" << std::setw(10) << col_unit_labels[i] << "]";
 
      if(i >= 2 && used_stim_ids[i-2] > -1) {
        stimulus & s = elec->stimuli[used_stim_ids[i-2]];
 
        if (is_potential(s.phys.type)) {
          if(s.phys.type == GND_ex)
            fh << " ground stim" << smpl_endl;
          else
            fh << " applied: " << std::to_string(s.pulse.strength) << smpl_endl;
        } else {
          fh << smpl_endl;
        }
      } else {
        fh << smpl_endl;
      }
    }
 
    // plot column short labels
    fh << "#";
    for(size_t i = 0; i<col_short_labels.size(); i++)
      fh << std::setw(col_width[i] - 3) <<  col_short_labels[i].c_str() << std::setw(3) << " ";
 
    // plot column units
    fh << smpl_endl << "#";
    for(size_t i = 0; i<col_unit_labels.size(); i++)
      fh << "[" << std::setw(col_width[i]-2) <<  col_unit_labels[i].c_str() << "]";
 
    // step through simulated time period
    fh << smpl_endl << std::fixed;
    do {
      // time and discrete time
      fh << std::setw(col_width[0]) << std::setprecision(3) << user_globals::tm_manager->time;
      fh << std::setw(col_width[1]) << user_globals::tm_manager->d_time;
 
      // iterate over all timers
      int col = 2;
      for (size_t tid = 0; tid < used_timer_ids.size(); tid++)
      {
        int timer_id  = used_timer_ids[tid];
        base_timer* t = user_globals::tm_manager->timers[timer_id];
 
        // type of timer: plain trigger or trigger linked to signal
        if(!t->d_trigger_dur) {
          int On = t->triggered ? 1 : 0;
          fh << std::setw(col_width[col]) << On;
        } else if(elec) {
          // figure out value of signal linked to this timer
          double val = 0.;
 
          // determine timer linked to which physics, for now we deal with electrics only
          val = elec->timer_val(timer_id);
 
          fh << std::setw(col_width[col]) << std::setprecision(3) << val;
        }
        col++;
      }
 
      fh << smpl_endl;
 
      // advance time
      user_globals::tm_manager->update_timers();
    } while (!user_globals::tm_manager->elapsed());
 
    fh.close();
 
    // reset timer to start before actual simulation
    user_globals::tm_manager->reset_timers();
  }
 
  return err;
}
 
void init_console_output(const timer_manager & tm, prog_stats & p)
{
  const char*  h1_prog = "PROG\t----- \t----\t-------\t-------|";
  const char*  h2_prog = "time\t%%comp\ttime\t ctime \t  ETA  |";
  const char*  h1_wc   = "\tELAPS |";
  const char*  h2_wc   = "\twc    |";
 
  p.start = get_time();
  p.last  = p.start;
 
  log_msg(NULL, 0, 0, "%s", h1_prog );
  log_msg(NULL, 0, NONL, "%s", h2_prog );
  log_msg(NULL, 0, 0, "" );
}
 
 
void time_to_string(float time, char* str, short str_size)
{
  int req_hours =  ((int)(time)) / 3600;
  int req_min   = (((int)(time)) % 3600) / 60;
  int req_sec   = (((int)(time)) % 3600) % 60;
 
  snprintf(str, str_size, "%d:%02d:%02d", req_hours, req_min, req_sec);
}
 
void update_console_output(const timer_manager & tm, prog_stats & p)
{
 
  float progress = 100.*(tm.time - tm.start) / (tm.end - tm.start);
  float elapsed_time = timing(p.curr, p.start);
  float req_time = (elapsed_time / progress) * (100.0f - progress);
 
  if(progress == 0.0f)
    req_time = 0.0f;
 
  char elapsed_time_str[256];
  char req_time_str[256];
  time_to_string(elapsed_time, elapsed_time_str, 255);
  time_to_string(req_time,     req_time_str,     255);
 
  log_msg( NULL, 0, NONL, "%.2f\t%.1f\t%.1f\t%s\t%s",
      tm.time,
      progress,
      (float)(p.curr - p.last),
      elapsed_time_str,
      req_time_str);
 
  p.last = p.curr;
 
  // we add an empty string for newline and flush
  log_msg( NULL, 0, ECHO | FLUSH, "");
}
 
void simulate()
{
  // here we want to include all time dependent physics, to check if we have any of those
  bool have_timedependent_phys = (phys_defined(PHYSREG_INTRA_ELEC) || phys_defined(PHYSREG_EIKONAL) || phys_defined(PHYSREG_EMI));
 
  if(!have_timedependent_phys) {
    log_msg(0,0,0, "\n no time-dependent physics region registered, skipping simulate loop..\n");
    return;
  }
 
  log_msg(0,0,0, "\n    *** Launching simulation ***\n");
 
  set_dir(OUTPUT);
 
  if(param_globals::dump_protocol)
    plot_protocols("protocol.trc");
 
  prog_stats prog;
  timer_manager & tm = *user_globals::tm_manager;
  init_console_output(tm, prog);
 
  // main loop
  do {
    // console output
    if(tm.trigger(iotm_console)) {
      // print console
      update_console_output(tm, prog);
    }
 
    // in order to be closer to carpentry we first do output and then compute the solution
    // for the next time slice ..
    if (tm.trigger(iotm_spacedt)) {
      for(const auto & it : user_globals::physics_reg) {
        it.second->output_step();
      }
    }
    // compute step
    for(const auto & it : user_globals::physics_reg) {
      Basic_physic* p = it.second;
      if (tm.trigger(p->timer_idx))
        p->compute_step();
    }
 
    // advance time
    tm.update_timers();
  } while (!tm.elapsed());
 
  log_msg(0,0,0, "\n\nTimings of individual physics:");
  log_msg(0,0,0, "------------------------------\n");
 
  for(const auto & it : user_globals::physics_reg) {
    Basic_physic* p = it.second;
    p->output_timings();
  }
 
}
 
void post_process()
{
  if(param_globals::post_processing_opts & RECOVER_PHIE)  {
    log_msg(NULL,0,ECHO,"\nPOSTPROCESSOR: Recovering Phie ...");
    log_msg(NULL,0,ECHO,  "----------------------------------\n");
 
    // do postprocessing
    int err = postproc_recover_phie();
 
    if(!err)  {
      log_msg(NULL,0,ECHO,"\n-----------------------------------------");
      log_msg(NULL,0,ECHO,  "POSTPROCESSOR: Successfully recoverd Phie.\n");
    }
  }
}
 
Basic_physic* get_physics(physic_t p, bool error_if_missing)
{
  auto it = user_globals::physics_reg.find(p);
 
  if(it != user_globals::physics_reg.end()) {
    return it->second;
  } else {
    if(error_if_missing) {
      log_msg(0,5,0, "%s error: required physic is not active! Usually this is due to an inconsistent experiment configuration. Aborting!", __func__);
      EXIT(EXIT_FAILURE);
    }
 
    return NULL;
  }
}
 
sf_vec* get_data(datavec_t d)
{
  sf_vec* ret = NULL;
 
  if(user_globals::datavec_reg.count(d))
    ret = user_globals::datavec_reg[d];
 
  return ret;
}
 
void register_data(sf_vec* dat, datavec_t d)
{
  if(user_globals::datavec_reg.count(d) == 0) {
    user_globals::datavec_reg[d] = dat;
  }
  else {
    log_msg(0,5,0, "%s warning: trying to register already registered data vector.", __func__);
  }
}
 
void parse_mesh_types()
{
  std::map<mesh_t, sf_mesh> & mesh_registry = user_globals::mesh_reg;
 
  // This is the initial grid we read the hard-disk data into
  mesh_registry[reference_msh] = sf_mesh();
  // we specify the MPI communicator for the reference mesh,
  // all derived meshes will get this comminicator automatically
  mesh_registry[reference_msh].comm = PETSC_COMM_WORLD;
 
  auto register_new_mesh = [&] (mesh_t mt, int pidx) {
    if(!mesh_registry.count(mt)) {
      mesh_registry[mt] = sf_mesh();
      mesh_registry[mt].name = param_globals::phys_region[pidx].name;
    }
    return &mesh_registry[mt];
  };
 
  // based on cli parameters we determine which grids need to be defined
  for(int i=0; i<param_globals::num_phys_regions; i++)
  {
    sf_mesh* curmesh = NULL;
    // register mesh type
    switch(param_globals::phys_region[i].ptype) {
      case PHYSREG_EIKONAL:
      case PHYSREG_INTRA_ELEC:
        curmesh = register_new_mesh(intra_elec_msh, i);
        break;
 
      case PHYSREG_LAPLACE:
      case PHYSREG_EXTRA_ELEC:
        curmesh = register_new_mesh(extra_elec_msh, i);
        break;
#if WITH_EMI_MODEL
      case PHYSREG_EMI:
        curmesh = register_new_mesh(emi_msh, i);
        break;
#endif
 
      default:
        log_msg(0,5,0, "Unsupported mesh type %d! Aborting!", param_globals::phys_region[i].ptype);
        EXIT(EXIT_FAILURE);
    }
 
    if(curmesh) {
      // set mesh unique tags
      for(int j=0; j<param_globals::phys_region[i].num_IDs; j++)
        curmesh->extr_tag.insert(param_globals::phys_region[i].ID[j]);
    }
  }
}
 
void retag_elements(sf_mesh & mesh, TagRegion *tagRegs, int ntr)
{
  if(ntr == 0) return;
  // checkTagRegDefs(ntr, tagRegs);
 
  SF::vector<mesh_int_t> & ref_eidx = mesh.get_numbering(SF::NBR_ELEM_REF);
 
  for (int i=0; i<ntr; i++) {
    tagreg_t type = tagreg_t(tagRegs[i].type);
    SF::vector<mesh_int_t> elem_indices;
 
    if (type == tagreg_list)
      read_indices(elem_indices, tagRegs[i].elemfile, ref_eidx, mesh.comm);
    else {
      geom_shape shape;
      shape.type   = geom_shape::shape_t(tagRegs[i].type);
      shape.p0.x   = tagRegs[i].p0[0];
      shape.p0.y   = tagRegs[i].p0[1];
      shape.p0.z   = tagRegs[i].p0[2];
      shape.p1.x   = tagRegs[i].p1[0];
      shape.p1.y   = tagRegs[i].p1[1];
      shape.p1.z   = tagRegs[i].p1[2];
      shape.radius = tagRegs[i].radius;
 
      bool nodal = false;
      indices_from_geom_shape(elem_indices, mesh, shape, nodal);
    }
 
    if(get_global((long int)elem_indices.size(), MPI_SUM, mesh.comm) == 0)
      log_msg(0,3,0,"Tag region %d is empty", i);
 
    for(size_t j=0; j<elem_indices.size(); j++)
      mesh.tag[elem_indices[j]] = tagRegs[i].tag;
  }
 
  // output the vector?
  if(strlen(param_globals::retagfile))
  {
    update_cwd();
    set_dir(OUTPUT);
 
    int dpn = 1;
    SF::write_data_ascii(mesh.comm, ref_eidx, mesh.tag, param_globals::retagfile, dpn);
 
    // Set dir back to what is was prior to retagfile output
    set_dir(CURDIR);
  }
}
 
size_t renormalise_fibres(SF::vector<mesh_real_t> &fib, size_t l_numelem)
{
  size_t renormalised_count = 0;
 
  // using pragma omp without global OMP control can lead to massive compute stalls,
  // as all cores may be already occupied by MPI, thus they become oversubscribed. Once
  // there is a global OMP control in place, we can activate this parallel for again. -Aurel, 20.01.2022
  // #pragma omp parallel for schedule(static) reduction(+ : renormalised_count)
  for (size_t i = 0; i < l_numelem; i++)
  {
    const mesh_real_t f0 = fib[3*i+0], f1 = fib[3*i+1], f2 = fib[3*i+2];
    mesh_real_t fibre_len = sqrt(f0*f0 + f1*f1 + f2*f2);
 
    if (fibre_len && fabs(fibre_len - 1) > 1e-3) {
      fib[3 * i + 0] /= fibre_len;
      fib[3 * i + 1] /= fibre_len;
      fib[3 * i + 2] /= fibre_len;
      renormalised_count++;
    }
  }
 
  return renormalised_count;
}
 
void setup_meshes(bool require_fibers=true)
{
  log_msg(0,0,0, "\n    *** Processing meshes ***\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(reference_msh) == 1); // There must be a reference mesh
 
  set_dir(INPUT);
 
  // we always read into the reference mesh
  sf_mesh & ref_mesh = mesh_registry[reference_msh];
  MPI_Comm comm = ref_mesh.comm;
 
  int size, rank;
  double t1, t2, s1, s2;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  SF::vector<mesh_real_t>    pts;
  SF::vector<mesh_int_t> ptsidx;
 
  // we add pointers to the meshes that need vertex cooridnates to this list
  std::list< sf_mesh* > ptsread_list;
 
  // read element mesh data
  t1 = MPI_Wtime(); s1 = t1;
  if(verb) log_msg(NULL, 0, 0,"Reading reference mesh: %s.*", basename.c_str());
 
  SF::read_elements(ref_mesh, basename, require_fibers);
  SF::read_points(basename, comm, pts, ptsidx);
 
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  bool check_fibre_normality = true;
  if (check_fibre_normality and ref_mesh.fib.size()>0) {
    t1 = MPI_Wtime();
 
    // make sure that all fibre vectors have unit length
    size_t l_num_fixed_fib = renormalise_fibres(ref_mesh.fib, ref_mesh.l_numelem);
 
    size_t l_num_fixed_she = 0;
    if (ref_mesh.she.size() > 0)
      l_num_fixed_she = renormalise_fibres(ref_mesh.she, ref_mesh.l_numelem);
 
    unsigned long fixed[2] = {(unsigned long) l_num_fixed_fib, (unsigned long) l_num_fixed_she};
    MPI_Allreduce(MPI_IN_PLACE, fixed, 2, MPI_UNSIGNED_LONG, MPI_SUM, comm);
 
    if (fixed[0] + fixed[1] > 0)
      log_msg(NULL, 0, 0, "Renormalised %ld longitudinal and %ld sheet-transverse fibre vectors.", fixed[0], fixed[1]);
 
    t2 = MPI_Wtime();
    if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
  }
 
  if(param_globals::numtagreg > 0) {
    log_msg(0, 0, 0, "Re-tagging reference mesh");
 
    // the retagging requires vertex coordinates, as such we need to read them into
    // the reference mesh
    ptsread_list.push_back(&ref_mesh);
    SF::insert_points(pts, ptsidx, ptsread_list);
 
    retag_elements(ref_mesh, param_globals::tagreg, param_globals::numtagreg);
 
    // we clear the list of meshet to receive vertices
    ptsread_list.clear();
  }
 
  if(verb) log_msg(NULL, 0, 0, "Processing submeshes");
 
  for(auto it = mesh_registry.begin(); it != mesh_registry.end(); ++it) {
    mesh_t grid_type = it->first;
    sf_mesh & submesh = it->second;
 
    if(grid_type != reference_msh) {
      if(verb > 1) log_msg(NULL, 0, 0, "\nSubmesh name: %s", submesh.name.c_str());
      t1 = MPI_Wtime();
 
      if(submesh.extr_tag.size() && grid_type != emi_msh)
        extract_tagbased(ref_mesh, submesh);
      else {
        // all submeshes should be defined on sets of tags, for backwards compatibility
        // we do a fiber based intra_elec_msh extraction if no tags are provided. Also, we
        // could do special treatments of any other physics type here. It would defeat
        // the purpose of the tag-based design, though. -Aurel
        switch(grid_type) {
          case emi_msh:
          case intra_elec_msh: extract_myocardium(ref_mesh, submesh, require_fibers); break;
          default:        extract_tagbased(ref_mesh, submesh);   break;
        }
      }
      t2 = MPI_Wtime();
      if(verb > 1) log_msg(NULL, 0, 0, "Extraction done in %f sec.", float(t2 - t1));
 
      ptsread_list.push_back(&submesh);
    }
  }
 
  // KDtree partitioning requires the coordinates to be present in the mesh data
  if(param_globals::pstrat == 2)
    SF::insert_points(pts, ptsidx, ptsread_list);
 
  for(auto it = mesh_registry.begin(); it != mesh_registry.end(); ++it)
  {
    mesh_t grid_type = it->first;
    sf_mesh & submesh = it->second;
    if(grid_type != reference_msh && grid_type!=emi_msh) {
      if(verb > 2) log_msg(NULL, 0, 0, "\nSubmesh name: %s", submesh.name.c_str());
      SF::vector<mesh_int_t> part;
 
      // generate partitioning vector
      t1 = MPI_Wtime();
      switch(param_globals::pstrat) {
        case 0:
          if(verb > 2) log_msg(NULL, 0, 0, "Using linear partitioning ..");
          break;
 
#ifdef WITH_PARMETIS
        case 1:
        {
          if(verb > 2) log_msg(NULL, 0, 0, "Using Parmetis partitioner ..");
          SF::parmetis_partitioner<mesh_int_t, mesh_real_t> partitioner(param_globals::pstrat_imbalance, 2);
          partitioner(submesh, part);
          break;
        }
#endif
        default:
        case 2: {
          if(verb > 2) log_msg(NULL, 0, 0, "Using KDtree partitioner ..");
          SF::kdtree_partitioner<mesh_int_t, mesh_real_t> partitioner;
          partitioner(submesh, part);
          break;
        }
      }
      t2 = MPI_Wtime();
      if(verb > 2) log_msg(NULL, 0, 0, "Partitioning done in %f sec.", float(t2 - t1));
 
      if(param_globals::pstrat > 0) {
        if(param_globals::gridout_p) {
          std::string out_name = get_basename(param_globals::meshname);
          if(grid_type == intra_elec_msh)      out_name += "_i.part.dat";
          else if(grid_type == extra_elec_msh) out_name += "_e.part.dat";
 
          set_dir(OUTPUT);
          log_msg(0,0,0, "Writing \"%s\" partitioning to: %s", submesh.name.c_str(), out_name.c_str());
          write_data_ascii(submesh.comm, submesh.get_numbering(SF::NBR_ELEM_REF), part, out_name);
        }
 
        t1 = MPI_Wtime();
        SF::redistribute_elements(submesh, part);
        t2 = MPI_Wtime();
        if(verb > 2) log_msg(NULL, 0, 0, "Redistributing done in %f sec.", float(t2 - t1));
      }
 
      t1 = MPI_Wtime();
      SF::submesh_numbering<mesh_int_t, mesh_real_t> sm_numbering;
      sm_numbering(submesh);
      t2 = MPI_Wtime();
      if(verb > 2) log_msg(NULL, 0, 0, "Canonical numbering done in %f sec.", float(t2 - t1));
 
      t1 = MPI_Wtime();
      submesh.generate_par_layout();
      SF::petsc_numbering<mesh_int_t, mesh_real_t> p_numbering(submesh.pl);
      p_numbering(submesh);
      t2 = MPI_Wtime();
      if(verb > 2) log_msg(NULL, 0, 0, "PETSc numbering done in %f sec.", float(t2 - t1));
      if(verb > 2) print_DD_info(submesh);
    }
  }
 
  SF::insert_points(pts, ptsidx, ptsread_list);
  ref_mesh.clear_data();
 
  s2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "All done in %f sec.", float(s2 - s1));
}
 
 
void output_meshes()
{
  bool write_intra_elec = mesh_is_registered(intra_elec_msh) && param_globals::gridout_i;
  bool write_extra_elec = mesh_is_registered(extra_elec_msh) && param_globals::gridout_e;
 
  set_dir(OUTPUT);
  std::string output_base = get_basename(param_globals::meshname);
 
  if(write_intra_elec) {
    sf_mesh & mesh = get_mesh(intra_elec_msh);
 
    if(param_globals::gridout_i & 1) {
      sf_mesh surfmesh;
      if(param_globals::output_level > 1)
        log_msg(0,0,0, "Computing \"%s\" surface ..", mesh.name.c_str());
      compute_surface_mesh(mesh, SF::NBR_SUBMESH, surfmesh);
 
      std::string output_file = output_base + "_i.surf";
      log_msg(0,0,0, "Writing \"%s\" surface: %s", mesh.name.c_str(), output_file.c_str());
      write_surface(surfmesh, output_file);
    }
    if(param_globals::gridout_i & 2) {
      bool write_binary = false;
 
      std::string output_file = output_base + "_i";
      log_msg(0,0,0, "Writing \"%s\" mesh: %s", mesh.name.c_str(), output_file.c_str());
      write_mesh_parallel(mesh, write_binary, output_file.c_str());
    }
  }
 
  if(write_extra_elec) {
    sf_mesh & mesh = get_mesh(extra_elec_msh);
 
    if(param_globals::gridout_e & 1) {
      sf_mesh surfmesh;
      if(param_globals::output_level > 1)
        log_msg(0,0,0, "Computing \"%s\" surface ..", mesh.name.c_str());
 
      compute_surface_mesh(mesh, SF::NBR_SUBMESH, surfmesh);
      std::string output_file = output_base + "_e.surf";
      log_msg(0,0,0, "Writing \"%s\" surface: %s", mesh.name.c_str(), output_file.c_str());
      write_surface(surfmesh, output_file);
    }
    if(param_globals::gridout_e & 2) {
      bool write_binary = false;
      std::string output_file = output_base + "_e";
      log_msg(0,0,0, "Writing \"%s\" mesh: %s", mesh.name.c_str(), output_file.c_str());
      write_mesh_parallel(mesh, write_binary, output_file.c_str());
    }
  }
}
 
[[noreturn]] void cleanup_and_exit()
{
  destroy_physics();
  user_globals::scatter_reg.free_scatterings();
 
  param_clean();
  PetscFinalize();
 
  // close petsc error FD
  if(user_globals::petsc_error_fd)
    fclose(user_globals::petsc_error_fd);
 
  exit(EXIT_SUCCESS);
}
 
char* get_file_dir(const char* file)
{
  char* filecopy = dupstr(file);
  char* dir      = dupstr(dirname(filecopy));
 
  free(filecopy);
  return dir;
}
 
void setup_petsc_err_log()
{
  int rank = get_rank();
  set_dir(OUTPUT);
 
  if(rank == 0) {
    // we open a error log file handle and set it as petsc stderr
    user_globals::petsc_error_fd = fopen("petsc_err_log.txt", "w");
    PETSC_STDERR = user_globals::petsc_error_fd;
  }
  else {
    PetscErrorPrintf = PetscErrorPrintfNone;
  }
}
 
short get_mesh_dim(mesh_t id)
{
  const sf_mesh & mesh = get_mesh(id);
  SF::element_view<mesh_int_t, mesh_real_t> view(mesh, SF::NBR_PETSC);
  short mindim = 3;
 
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    view.set_elem(eidx);
    short cdim = view.dimension();
    if(mindim < cdim) mindim = cdim;
  }
 
  mindim = get_global(mindim, MPI_MIN, mesh.comm);
 
  return mindim;
}
 
void igb_output_manager::register_output_sync(sf_vec* inp_data,
                                              const mesh_t inp_meshid,
                                              const int dpn,
                                              const char* name,
                                              const char* units,
                                              const SF::vector<mesh_int_t>* idx,
                                              bool elem_data)
{
  sync_io_item IO;
 
  IO.data      = inp_data;
  IO.elem_flag = elem_data;
  IO.restr_idx = idx;
 
  IGBheader regigb;
  const timer_manager & tm = *user_globals::tm_manager;
  const int num_io = tm.timers[iotm_spacedt]->numIOs;
  int       err    = 0;
 
  int gsize = inp_data->gsize();
 
  // if we are restricting, we have to compute the restricted global size
  if(idx != NULL)
    gsize = get_global(int(idx->size()), MPI_SUM);
 
  regigb.x(gsize / dpn);
  regigb.dim_x(regigb.x()-1);
  regigb.inc_x(1);
 
  regigb.y(1); regigb.z(1);
  regigb.t(num_io);
  regigb.dim_t(tm.end-tm.start);
 
  switch(dpn) {
    default:
    case 1: regigb.type(IGB_FLOAT); break;
    case 3: regigb.type(IGB_VEC3_f); break;
    case 4: regigb.type(IGB_VEC4_f); break;
    case 9: regigb.type(IGB_VEC9_f); break;
  }
 
  regigb.unites_x("um"); regigb.unites_y("um"); regigb.unites_z("um");
  regigb.unites_t("ms");
  regigb.unites(units);
 
  regigb.inc_t(param_globals::spacedt);
 
  if(get_rank() == 0) {
    FILE_SPEC file = f_open(name, "w");
    if(file != NULL) {
      regigb.fileptr(file->fd);
      regigb.write();
      delete file;
    }
    else err++;
  }
 
  err = get_global(err, MPI_SUM);
  if(err) {
    log_msg(0,5,0, "%s error: Could not set up data output! Aborting!", __func__);
    EXIT(1);
  }
 
  IO.igb = regigb;
 
  SF::mixed_tuple<mesh_t, int> mesh_spec = {inp_meshid, dpn};
  IO.spec = mesh_spec;
 
  if(elem_data) {
    if(buffmap_elem.find(mesh_spec) == buffmap_elem.end()) {
      sf_vec *inp_copy; SF::init_vector(&inp_copy, inp_data);
      buffmap_elem[mesh_spec] = inp_copy;
    }
  } else {
    if(buffmap.find(mesh_spec) == buffmap.end()) {
      sf_vec *inp_copy; SF::init_vector(&inp_copy, inp_data);
      buffmap[mesh_spec] = inp_copy;
    }
  }
 
  this->sync_IOs.push_back(IO);
}
 
void igb_output_manager::register_output_async(sf_vec* inp_data,
                                               const mesh_t inp_meshid,
                                               const int dpn,
                                               const char* name,
                                               const char* units,
                                               const SF::vector<mesh_int_t>* idx,
                                               bool elem_data)
{
  sf_mesh & mesh = get_mesh(inp_meshid);
  SF::vector<mesh_int_t> ioidx;
  int rank = get_rank();
 
  async_io_item IO;
  IO.data = inp_data;
  IO.restr_idx = idx;
 
  if(elem_data) {
    const SF::vector<mesh_int_t> & nbr = mesh.get_numbering(SF::NBR_ELEM_SUBMESH);
    ioidx.resize(mesh.l_numelem);
    for(size_t i=0; i<mesh.l_numelem; i++)
      ioidx[i] = nbr[i];
  } else {
    const SF::vector<mesh_int_t> & alg_nod = mesh.pl.algebraic_nodes();
    const SF::vector<mesh_int_t> & nbr     = mesh.get_numbering(SF::NBR_SUBMESH);
 
    if(idx == NULL) {
      ioidx.resize(alg_nod.size());
 
      for(size_t i=0; i<alg_nod.size(); i++)
        ioidx[i] = nbr[alg_nod[i]];
    } else {
      ioidx.resize(idx->size());
      IO.restr_petsc_idx.resize(idx->size());
 
      for(size_t i=0; i<idx->size(); i++) {
        mesh_int_t loc_nodal = (*idx)[i];
        ioidx[i] = nbr[loc_nodal];
        IO.restr_petsc_idx[i] = local_nodal_to_local_petsc(mesh, rank, loc_nodal);
      }
    }
  }
 
  int id = async::COMPUTE_register_output(ioidx, dpn, name, units);
  IO.IO_id = id;
 
  this->async_IOs.push_back(IO);
}
 
void igb_output_manager::register_output(sf_vec* inp_data,
                                         const mesh_t inp_meshid,
                                         const int dpn,
                                         const char* name,
                                         const char* units,
                                         const SF::vector<mesh_int_t>* idx,
                                         bool elem_data)
{
  if(param_globals::num_io_nodes == 0)
    register_output_sync(inp_data, inp_meshid, dpn, name, units, idx, elem_data);
  else
    register_output_async(inp_data, inp_meshid, dpn, name, units, idx, elem_data);
}
 
sf_vec* igb_output_manager::fill_output_buffer(const sync_io_item & it)
{
  const SF::mixed_tuple<mesh_t, int> & cspec = it.spec;
  sf_vec* data_vec = it.data;
 
  bool have_perm = it.elem_flag ? have_permutation(cspec.v1, ELEM_PETSC_TO_CANONICAL, cspec.v2):
                                  have_permutation(cspec.v1, PETSC_TO_CANONICAL, cspec.v2);
 
  if(have_perm) {
    sf_vec*   perm_vec = it.elem_flag ? this->buffmap_elem[cspec] : this->buffmap[cspec];
    SF::scattering* sc       = it.elem_flag ? get_permutation(cspec.v1, ELEM_PETSC_TO_CANONICAL, cspec.v2):
                                              get_permutation(cspec.v1, PETSC_TO_CANONICAL, cspec.v2);
    sc->forward(*data_vec, *perm_vec);
    return perm_vec;
  } else {
    return data_vec;
  }
}
 
void igb_output_manager::write_data()
{
  SF::vector<float>      restr_buff;
  int rank = get_rank();
  // loop over registered datasets and root-write one by one
  //
  for (sync_io_item & it : sync_IOs) {
    // write to associated file descriptor
    FILE* fd = static_cast<FILE*>(it.igb.fileptr());
 
    // fill the output buffer
    sf_vec* buff = fill_output_buffer(it);
 
    if(it.restr_idx == NULL) {
      buff->write_binary<float>(fd);
    } else {
      const SF::vector<mesh_int_t> & idx = *it.restr_idx;
      SF_real* p = buff->ptr();
 
      restr_buff.resize(idx.size()); restr_buff.resize(0);
 
      for(mesh_int_t ii : idx)
        restr_buff.push_back(p[ii]);
 
      root_write(fd, restr_buff, PETSC_COMM_WORLD);
      buff->release_ptr(p);
    }
  }
 
  // do all A-synchronous output:
  // loop over IDs received from the IO nodes and trigger async output
  //
  for (async_io_item & it : async_IOs) {
    SF_real* p    = it.data->ptr();
    int        ls = it.data->lsize();
    int        id = it.IO_id;
 
    if(it.restr_idx == NULL)
      async::COMPUTE_do_output(p, ls, id);
    else {
      async::COMPUTE_do_output(p, it.restr_petsc_idx, id);
    }
 
    it.data->release_ptr(p);
  }
}
 
void igb_output_manager::close_files_and_cleanup()
{
  if(get_rank() == 0) {
    // loop over registered datasets and close fd
    for(sync_io_item & it : sync_IOs) {
      FILE* fd = static_cast<FILE*>(it.igb.fileptr());
      fclose(fd);
    }
  }
 
  for(auto it = buffmap.begin(); it != buffmap.end(); ++it)
    delete it->second;
 
  for(auto it = buffmap_elem.begin(); it != buffmap_elem.end(); ++it)
    delete it->second;
 
  // we resize the arrays and clear the maps so that we are safe when calling
  // close_files_and_cleanup multiple times.
  sync_IOs.resize(0);
  async_IOs.resize(0);
  buffmap.clear();
  buffmap_elem.clear();
}
 
IGBheader* igb_output_manager::get_igb_header(const sf_vec* vec)
{
  for(sync_io_item & it : sync_IOs) {
    if(it.data == vec)
      return &it.igb;
  }
 
  return NULL;
}
 
void output_parameter_file(const char *fname, int argc, char **argv)
{
  const int   max_line_len = 128;
  const char* file_sep = "#=======================================================";
 
  // make sure only root executes this function
  if(get_rank() != 0)
    return;
 
  FILE_SPEC out = f_open(fname, "w");
  fprintf(out->fd, "# CARP GIT commit hash: %s\n", GIT_COMMIT_HASH);
  fprintf(out->fd, "# dependency hashes:    %s\n", SUBREPO_COMMITS);
  fprintf(out->fd, "\n");
 
  // output the command line
  char line[8196] = "# ";
 
  for (int j=0; j<argc; j++) {
    strcat(line, argv[j]);
    if(strlen(line) > max_line_len) {
      fprintf(out->fd, "%s\n", line);
      strcpy(line, "# ");
    } else
      strcat(line, " ");
  }
 
  fprintf(out->fd, "%s\n\n", line);
  set_dir(INPUT);
 
  // convert command line to a par file
  for (int i=1; i<argc; i++) {
    std::string argument(argv[i]);
    if (argument == "+F" || argument.find("_options_file")!= std::string::npos) {
 
      std::string init = "";
      if (argument.find("_options_file")!= std::string::npos) {
        fprintf(out->fd, "%s = %s\n", argument.substr(1).c_str(), argv[i+1]);
        init = "#";
      }
      fprintf(out->fd, "%s>>\n", file_sep);
      // import par files
      i++;
      fprintf(out->fd, "## %s ##\n", argv[i]);
      FILE *in = fopen(argv[i], "r");
      while (fgets(line, 8196, in))
        fprintf(out->fd, "%s%s", init.c_str(), line);
      fclose(in);
      fprintf(out->fd, "\n##END of %s\n", argv[i]);
      fprintf(out->fd, "%s<<\n\n", file_sep);
    }
    else if(argv[i][0] == '-')
    {
      bool prelast       = (i==argc-1);
      bool paramFollows  = !prelast && ((argv[i+1][0] != '-') ||
                           ((argv[i+1][1] >= '0') && (argv[i+1][1] <= '9')));
 
      // strip leading hyphens from command line opts
      // assume options do not start with numbers
      if(paramFollows) {
        // nonflag option
        char *optcpy = strdup(argv[i+1]);
        char *front  = optcpy;
        // strip {} if present for arrays of values
        while(*front==' ' && *front) front++;
        if(*front=='{') {
          while(*++front == ' ');
          char *back = optcpy+strlen(optcpy)-1;
          while(*back==' ' && back>front) back--;
          if(*back == '}')
            *back = '\0';
        }
        if (strstr(front, "=") != nullptr)  // if "=" is find then we need ""
          fprintf(out->fd, "%-40s= \"%s\"\n", argv[i]+1, front);
        else
          fprintf(out->fd, "%-40s= %s\n", argv[i]+1, front);
        free(optcpy);
        i++;
      }
      else // a flag was specified
        fprintf(out->fd, "%-40s= 1\n", argv[i]);
    }
  }
  f_close(out);
}
 
void savequit()
{
  if(!get_physics(elec_phys)) return;
 
  set_dir(OUTPUT);
 
  double time = user_globals::tm_manager->time;
  char save_fn[512];
 
  snprintf(save_fn, sizeof save_fn, "exit.save.%.3f.roe", time);
  log_msg(NULL, 0, 0, "savequit called at time %g\n", time);
 
  Electrics* elec = (Electrics*) get_physics(elec_phys);
  elec->ion.miif->dump_state(save_fn, time, intra_elec_msh, false, GIT_COMMIT_COUNT);
 
  cleanup_and_exit();
}
 
}  // namespace opencarp