bench.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
// ----------------------------------------------------------------------------

#define PrMGLOBAL
#include <stdlib.h>
#include <stdio.h>
#include "cmdline.h"
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <limits.h>
#include <math.h>
#include <stdbool.h>
#include <vector>
#include <sstream>
#include <string>
#include <cmath>
#include "basics.h"

#include "MULTI_ION_IF.h"

#include "trace.h"
#include "clamp.h"
#include "ap_analyzer.h"
#include "restitute.h"
#include "stretch.h"
#include "sv_init.h"
#include "bench_utils.h"
#include "target.h"

#include "libgen.h"
#include "build_info.h"

#include "petsc_utils.h" //TODO for EXIT, FPRINT*
#include "SF_init.h" // for SF::init_xxx()


// globals
namespace opencarp {
namespace user_globals {
  SF::scatter_registry scatter_reg;
  std::map<mesh_t, sf_mesh> mesh_reg;
  std::map<SF::quadruple<int>, SF::index_mapping<int> > map_reg;
  // /// the physics
  // std::map<physic_t, Basic_physic*> physics_reg;
  timer_manager* tm_manager;
  // /// important solution vectors from different physics
  // std::map<datavec_t, sf_vec*> datavec_reg;
}  // namespace user_globals
}  // namespace opencarp

using namespace opencarp;
using namespace limpet;

#undef __FUNCT__
#define __FUNCT__ "main"

int main(int argc, char *argv[]) {

  double dt = 10.e-3;   // time step needed for the setup of the LUT's (ms)

  // specify the number of regions and their names
  const int num_region = 1;
  enum { R1 }
  Region;

  // specifiy the ionic model for each region
  IonType* RegionDef;

  // specify the plugins for each region
  int num_plugins;
  IonTypeList RegionPlug = {};

  // specify the region for each node
  char *Regions;


  struct gengetopt_args_info params;

  // we allocate the vectors on the heap since they will be freed by the MIIF
  sf_vec* Vmv;
  sf_vec* I_ion;
  
  FILE *fhdls[NUM_IMP_DATA_TYPES+1];
  GVEC_DUMP gvd;
  IOCtrl io;

  char *PetscDBfile = COMPAT_PETSC_NULLPTR;  // could be used to pass in ~/.petscrc
  char *help_msg    = COMPAT_PETSC_NULLPTR;
  initialize_PETSc(&argc, argv, PetscDBfile, help_msg);
  COMPAT_PetscOptionsInsertString("-options_left no");

  if (cmdline_parser(argc, argv, &params) != 0)
    exit(1);

  // performance measurements
  double t1, t2, t3, t4, t5, t6, t7, t8;
  struct timeval tt1;

  event_timing timings[N_TIMINGS];
  initialize_timings(timings);

  // start measurementstrace_IA
  get_time(t1);

  // assign command line params and depending variables
  dt = params.dt_arg;
  double i_stim      = params.stim_curr_arg*dt;
  double bcl         = params.bcl_arg;
  double stim_start  = params.stim_start_arg;
  double stim_dur    = params.stim_dur_arg;
  int    nstim       = params.numstim_arg ? params.numstim_arg : -1;
  double dt_out      = params.dt_out_arg;
  double start_out   = params.start_out_arg;
  int    numNode     = params.num_arg;
  int    stim_assign = params.stim_assign_flag;
  bool   extVmUpdate = params.ext_vm_update_given;
  bool   analyze_AP  = params.APstatistics_given || params.restitute_given;
  bool   restitute   = params.restitute_given;
  bool   validate    = params.validate_given;
  bool   APclamp     = params.AP_clamp_file_given;
  bool   do_trace    = restitute ? params.res_trace_flag : (!params.no_trace_flag);
  TrgList stim_lst;
  float   duration = determine_duration(&params, &stim_lst);

  io.w2file = validate ? 1 : params.fout_given;
  io.wbin   = validate ? 1 : params.bin_flag;

  if(get_size() > numNode) {
    numNode = get_size();

    log_msg(NULL, 3, 0, "Warning: Number of nodes is less than number of processes used!");
    log_msg(NULL, 3, 0, "Setting number of nodes to %d.", numNode);
  }

  // stimulate with trace from file
  trace stim_trace;
  bool  tr_stim = false;
  if (params.stim_file_given) {
    char *fname = dupstr(params.stim_file_arg);
    read_trace(&stim_trace, fname);
    resample_trace(&stim_trace, dt);
    tr_stim  = true;
    stim_dur = stim_trace.dur;
    if (bcl < stim_trace.dur)
      log_msg(NULL, 4, 0, "BCL less than stimulus duration!");
  }

  // illuminate with trace from file
  trace light_trace;
  bool tr_light = false;
  if( params.light_file_given ) {
    char* fname = dupstr( params.light_file_arg );
    read_trace( &light_trace, fname );
    resample_trace( &light_trace, dt );
    tr_light = true;

    int n_clipped_to_0 = 0;
    for(int i=0; i<light_trace.N; i++) {
      if( light_trace.s[i] < 0.0 ) {
        light_trace.s[i] = 0.0;
        n_clipped_to_0++;
      }
    }
    if( n_clipped_to_0 ) {
      FPRINTF( WORLD stderr, "WARNING: %d irradiance values in %s were negative (clipped to 0.0)\n", n_clipped_to_0, fname );
    }
  }

  // state variable clamp
  Clamp *sv_cl;
  int SVclamp = process_sv_clamps(params.clamp_SVs_arg, params.SV_clamp_files_arg, &sv_cl, dt);

  // action potential clamp
  Clamp ap_cl;
  if (APclamp){
    const char *Vm_clamp = "Vm";
    initialize_sv_clamp(&ap_cl, Vm_clamp, params.AP_clamp_file_arg, dt);
  }

  // voltage clamp experiment
  Clamp cl;
  bool  IsVclamp = initialize_clamp(&cl, params.clamp_arg, params.clamp_ini_arg,
                                    params.clamp_start_arg, params.clamp_dur_arg,
                                    params.clamp_file_arg,
                                    !params.clamp_ini_given ? CLAMP_TM_IDX : 0, &duration);
  if (IsVclamp) {
    if (APclamp) {
      log_msg(NULL, 0, 0,  "info: cannot clamp voltage and have AP clamp at same time");
      exit(1);
    }
    if (cl.transient) {
      log_msg(NULL, 0, 0,  "info: Vm will be clamped to %f mV over the interval [%f, %f]",
              params.clamp_arg, params.clamp_start_arg, params.clamp_start_arg + params.clamp_dur_arg);
    } else {
      log_msg(NULL, 0, 0,  "info: Vm will be clamped from %f to %f mV over the interval [%f, %f]",
              params.clamp_ini_arg, params.clamp_arg, params.clamp_start_arg,
              params.clamp_start_arg + params.clamp_dur_arg);
    }
  }

  action_potential AP;
  if (analyze_AP) initialize_AP_analysis(&AP);

  Regions = static_cast<char *>(calloc(numNode, sizeof(char) ));

  SF::init_vector(&Vmv, numNode, -1);
  SF::init_vector(&I_ion, Vmv);

  int *stim_list = static_cast<int *>(malloc(numNode*sizeof(int) ));
  for (int i = 0; i < numNode; i++) {
    Regions[i]   = R1;
    stim_list[i] = i;
  }

  for (unsigned int i = 0; i < params.load_module_given; i++) {
#ifdef HAVE_DLOPEN
    load_ionic_module(params.load_module_arg[i]);
    if (!params.imp_given) {

      char *imp = strdup(basename(params.load_module_arg[i]) );
      char *dot = strrchr(imp, '.');
      if (dot)
        *dot = 0;
      params.imp_arg = imp;
    }
#else // ifdef HAVE_DLOPEN
    log_msg(NULL, 5, 0, "Compile with USE_DLOPEN to support dynamic module loading.");
    EXIT(-1);
#endif // ifdef HAVE_DLOPEN
  }

  if (get_ion_type(std::string(params.imp_arg)) == NULL) {
    fprintf(stderr, "Illegal IMP specified: %s\n", params.imp_arg);
    exit(1);
  } else {
    RegionDef = get_ion_type(std::string(params.imp_arg));
  }
  if (!get_plug_flag(params.plug_in_arg, &num_plugins, RegionPlug)) {
    fprintf(stderr, "Illegal plugin specified: %s\n", params.plug_in_arg);
    exit(1);
  }
  // Fetch target enum from command line string
  Target target = get_target_from_string(std::string(params.target_arg));
  // Automatic target selection
  if (target == Target::AUTO) {
    target = RegionDef->select_target(target);
  }
  // Check whether target exists and was generated for the given model and plugin
  if (target == Target::UNKNOWN) {
    fprintf(stderr, "Unkown target: %s\n", params.target_arg);
    fprintf(stderr, "Available targets are: %s\n", get_target_list_string().c_str());
    exit(1);
  }
  // Check if the main imp was generated for the requested target
  if (RegionDef->select_target(target) == Target::UNKNOWN) {
    fprintf(stderr, "Model %s was not generated for target %s\n",
        params.imp_arg, get_string_from_target(target).c_str());
    exit(1);
  }
  // Check the same thing if a plugin is specified
  if (!RegionPlug.empty()) {
    if (RegionPlug[0].get().select_target(target) == Target::UNKNOWN) {
      fprintf(stderr, "Plugin %s was not generated for target %s\n",
          params.plug_in_arg, get_string_from_target(target).c_str());
      exit(1);
    }
  }

  // data structure to manage multiple IIF's
  MULTI_IF MIIF;
  MULTI_IF doppel, *cMIIF = &MIIF;

  MIIF.N_IIF       = num_region;
  MIIF.iontypes    = {*RegionDef};
  MIIF.targets     = {target};
  MIIF.IIFmask     = Regions;
  MIIF.numplugs    = &num_plugins;
  MIIF.plugtypes   = {RegionPlug};
  MIIF.numNode     = Vmv->lsize();
  MIIF.gdata[Vm]   = Vmv;
  MIIF.gdata[Iion] = I_ion;

  if (extVmUpdate)
    MIIF.extUpdateVm = true;

  double setup_time = timing(t2, t1);
  update_timing(timings+SETUP_IDX, setup_time);

  MIIF.initialize_MIIF();

  if (params.list_imps_flag || params.plugin_outputs_flag) {
    print_models(params.plugin_outputs_flag);
    exit(0);
  }

  if (params.imp_info_flag) {
    print_param_help(RegionDef, RegionPlug);
    exit(0);
  }

  fprintf(stderr, "\n*** GIT tag:              %s\n", GIT_COMMIT_TAG);
  fprintf(stderr, "*** GIT hash:             %s\n", GIT_COMMIT_HASH);
  fprintf(stderr, "*** GIT repo:             %s\n", GIT_PATH);
  fprintf(stderr, "*** dependency commits:   %s\n\n", SUBREPO_COMMITS);

  if (params.buildinfo_flag) {
    fprintf(stderr, "\n*** -buildinfo flag is deprecated and will be removed.\n");
    fprintf( stderr, "\n*** The build information is printed by default. \n\n");
    exit(0);
  }

  MIIF.IIF[0]->tune(params.imp_par_arg, params.plug_in_arg, params.plug_par_arg);//FIXME

  // now compute the tables
  MIIF.initialize_currents(dt, 1);

  if (analyze_AP) {
    // initialize AP analyzer
    double Vm0 = getCellVal(MIIF.gdata[Vm], 0);
    for (int i = 0; i < VM_HIST_LEN; i++)
      AP.vm_trc[i] = Vm0;

    char ap_stats_fname[1024];
    snprintf(ap_stats_fname, sizeof ap_stats_fname, "%s_AP_stats.dat", params.fout_arg);
    AP.stats = fopen(ap_stats_fname, "wt");
    print_AP_stats_header(&AP, AP.stats);
  }

  restitution r;
  int n_dopple;
  double *t_dopple;
  if (restitute) {
    char *res_file = dupstr(params.res_file_arg);
    restitution_trigger_list(res_file, &r, params.restitute_arg, &n_dopple, &t_dopple);
    if (r.dur > duration) duration = r.dur;

    // output restitution data
    char ap_rstats_fname[1024];
    snprintf(ap_rstats_fname, sizeof ap_rstats_fname, "%s_APD_restitution.dat", params.fout_arg);
    AP.rstats = fopen(ap_rstats_fname, "wt");
    print_AP_stats_header(&AP, AP.rstats);
  }

  // read initial state vector from file and/or assignment
  double t = 0.;
  // if (params.restore_given)
  //   t = MIIF.restore_state( params.restore_arg, &pl, true);

  // illumination stuff
  double light_irrad = params.light_irrad_arg;
  double light_start = params.light_start_arg;
  double light_end = params.duration_arg;
  int    light_nstim = params.light_numstim_arg;
  double light_bcl = params.light_bcl_arg;
  double light_dur = params.light_dur_arg;

  if( light_irrad < 0 ) {
    FPRINTF( WORLD stderr, "WARNING: irradiance value %.3f makes no sense (clipped to 0.0)\n", light_irrad );
    light_irrad = 0.0;
  }
  TrgList light_lst;
  if(params.light_times_given)  //user specified list
    determine_stim_list(params.light_times_arg, &light_lst, false);

  if (params.SV_init_given)
    initial_SVs(&MIIF, params.SV_init_arg, params.imp_arg, params.plug_in_arg, numNode);

  timer_manager tmo(dt, t, duration);
  tmo.timers.resize(N_TIMERS);

  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0., CON_TM_IDX, "CON_TIMER");
  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0., SVD_TM_IDX, "SVD_TIMER");

  if (params.save_time_arg)
    tmo.initialize_singlestep_timer(params.save_time_arg, 0., STA_TM_IDX, "STA_TIMER");

  if (params.save_ini_time_given)
    tmo.initialize_singlestep_timer(params.save_ini_time_arg, 0., SSV_TM_IDX, "SSV_TIMER");

  // use a trigger for stimulation
  if (params.restitute_given) {
    std::vector<double> trg;
    trg.assign(r.trigs.lst, r.trigs.lst + r.trigs.n);
    tmo.initialize_neq_timer(trg, stim_dur, STM_TM_IDX, "STIM_TIMER");

    if (params.res_state_vector_given) {
      trg.assign(r.saveState.lst, r.saveState.lst + r.saveState.n);
      tmo.initialize_neq_timer(trg, 0, RES_SAVE_TM_IDX, "RES_SAVE_TIMER");
    }

    if (!strcmp(params.restitute_arg, "S1S2f") ) {
      trg.assign(t_dopple, t_dopple + n_dopple);
      tmo.initialize_neq_timer(trg, 1, DOPPLE_TM_IDX, "DOPPLE_TIMER");
    }
  }
  else if (params.stim_times_given) {
    std::vector<double> trg;
    trg.assign(stim_lst.lst, stim_lst.lst + stim_lst.n);
    tmo.initialize_neq_timer(trg, stim_dur, STM_TM_IDX, "STIM_TIMER");
  }
  else {
    tmo.initialize_eq_timer(stim_start, stim_start+duration, nstim, bcl, stim_dur, STM_TM_IDX, "STIM_TIMER");
  }

  // illumination timers
  if ( MIIF.gdata[illum] != NULL && !tr_light ) {
    if ( params.light_times_given ) {
      std::vector<double> trg;
      trg.assign(light_lst.lst, light_lst.lst + light_lst.n);
      tmo.initialize_neq_timer(trg, light_dur, LIGHT_TM_IDX, "LIGHT_TIMER", NULL);
    } else {
      tmo.initialize_eq_timer(light_start, light_end, light_nstim,
                              light_bcl, light_dur, LIGHT_TM_IDX, "LIGHT_TIMER", NULL);
    }
  }

  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0, TRACE_TM_IDX, "TRACE_TIMER");

  if (IsVclamp && cl.transient)
    tmo.initialize_singlestep_timer(params.clamp_start_arg, params.clamp_dur_arg, CLAMP_TM_IDX, "CLAMP_TIMER");
  if (!restitute && params.doppel_on_given) {
    tmo.initialize_singlestep_timer(params.doppel_on_arg, params.doppel_dur_arg, DOPPLE_TM_IDX, "DOPPLE_TIMER");
  }

  // read initial state from single cell file and spread it out to all cells
  if (params.read_ini_file_given) {
    if (read_sv(&MIIF, R1, params.read_ini_file_arg) != 0) {
      return 1;  // exit with error code.
    }
  }

  // stretch experiment
  stretch s;
  memset(&s, 0, sizeof(stretch));
  if (params.strain_given && MIIF.use_stretch()) {
    initializePulseStretch(params.strain_arg, params.strain_time_arg,
                           params.strain_dur_arg, params.strain_rate_arg,
                           params.strain_rate_arg, &s);

    sf_vec* lambdavec; SF::init_vector(&lambdavec, MIIF.gdata[Vm]);
    MIIF.gdata[Lambda] = lambdavec;
    SF_real *l = lambdavec->ptr();

    for (int i = 0; i < lambdavec->lsize(); i++)
      l[i] = s.pulse.sr;

    lambdavec->release_ptr(l);
  }

  // set up assigning stimulus current to ion concentrations
  std::stringstream stim_species(params.stim_species_arg);
  std::stringstream stim_ratios(params.stim_ratios_arg);
  std::string specie, ratio;
  std::vector<std::string> species_list;
  std::vector<float> ratios_list;
  float ratio_sum = 0;      
  if (stim_assign) {
    // split stim_species by ":" and then iterate over them
    while(std::getline(stim_species, specie, ':'))
    {
       species_list.push_back(specie);
    }
    while(std::getline(stim_ratios, ratio, ':'))
    {
       ratios_list.push_back(std::stof(ratio));
    }
    if (species_list.size() !=  ratios_list.size()) {
      fprintf(stderr, "Number of ion species for which stimulus should be assigned does not match number of ratio entries: %s | %s\n", params.stim_species_arg, params.stim_ratios_arg);
      exit(1);    
    }
    for (auto i : ratios_list ) 
      ratio_sum += i;
    if (std::abs(ratio_sum - 1) >= 0.001)
      log_msg(NULL, 3, 0, "Warning: Sum of ion species ratios for stimulus assignement != 1. Normalizing to sum %.1f.", ratio_sum);
  }

  // dumping of state variables goes here
  char sv_dump_names[1024];
  if (!params.fout_given && !params.validate_flag)
    snprintf(sv_dump_names, sizeof sv_dump_names, "%s_%s", params.fout_arg, params.imp_arg);
  else
    strcpy(sv_dump_names, params.fout_arg);

  if (params.validate_flag) {
    dump_all(&MIIF, R1, params.imp_arg, params.plug_in_arg, t, dt, sv_dump_names);
  } else {
    MIIF.sv_dump_add_by_name_list(R1, params.imp_arg, sv_dump_names,
                             params.imp_sv_dump_arg, params.plug_in_arg,
                             params.plug_sv_dump_arg, t, dt_out);
  }

  // dumping of LUTs
  if (params.dump_lut_flag)
    MIIF.dump_luts_MIIF(true);

  // initialize tracing
  int trace_nodes[] = {0};
  if (do_trace)
    open_trace(&MIIF, sizeof(trace_nodes)/sizeof(trace_nodes[0]), trace_nodes, &params.trace_no_arg, NULL);

  // dumping global vectors to file?
  open_globalvec_dump(fhdls, &gvd, &MIIF, sv_dump_names, &io);

  // initialization time
  double ini_time = timing(t3, t2);
  update_timing(timings + INIT_IDX, ini_time);
  doppel_MIIF(&MIIF, &doppel);

  int gpu_Vm_update = 1;
  double counter = 0.0;
  double epsilon = 0.0000001f;
  /***********************************************/
  /***************** MAIN LOOP *******************/
  /***********************************************/

  // Print simulation target
  fprintf(stderr, "Running simulation on target %s\n", get_string_from_target(target).c_str());
  while (!tmo.elapsed())
  {
    get_time(t4);
    t = tmo.time;

    if (tmo.triggered_now(DOPPLE_TM_IDX) ) {
      printf("[%f] Switching to new ion structure\n", tmo.time);
      doppel_update(&MIIF, &doppel);
      cMIIF = &doppel;
    }

    if (IsVclamp && cl.transient && tmo.trigger(CLAMP_TM_IDX) )
      cMIIF->extUpdateVm = true;
    else
      cMIIF->extUpdateVm = extVmUpdate;

    if (IsVclamp)
      clamp_signal(cMIIF, &cl, &tmo);

    if (APclamp)
      AP_clamp(&ap_cl, &tmo, cMIIF->gdata[Vm], tmo.triggered_now(STM_TM_IDX) );

    for (int i = 0; i < SVclamp; i++) {
      bool trigger;
      if (params.SV_I_trigger_flag)
        trigger =  tmo.triggered_now(STM_TM_IDX);
      else
        trigger = !tmo.d_time;

      sv_clamp(sv_cl+i, &tmo, cMIIF, trigger);
    }

    // output selected state variables
    cMIIF->dump_svs(tmo.timers[SVD_TM_IDX]);
    globalvec_dump(fhdls, &gvd, cMIIF, &tmo, &io, numNode);

    // save state of all cells to binary file
    if (tmo.trigger(STA_TM_IDX) )
      cMIIF->dump_state(params.save_file_arg, tmo.time, mesh_t::unset_msh, false, 0);

    // save state of cell 0 to text file
    if (tmo.trigger(SSV_TM_IDX) )
      save_sv(cMIIF, R1, params.save_ini_file_arg);

    // save state within an restitution protocol
    if (tmo.trigger(RES_SAVE_TM_IDX) )
      restitution_save_sv(cMIIF, R1, &r, &AP);

    // add stim current
    if (tmo.trigger(STM_TM_IDX) ) {
      gpu_Vm_update = 0;
      if (tr_stim) {
        i_stim = stim_trace.s[tmo.trigger_elapse(STM_TM_IDX)]*dt;
      }
      else if (params.stim_volt_given) {
        SF_real *d = cMIIF->gdata[Vm]->ptr();
        // peak_vm = last_vm = d[0];
        i_stim = (params.stim_volt_arg-d[0])/params.resistance_arg;
        cMIIF->gdata[Vm]->release_ptr(d);
      }

#ifndef ALGEBRAIC
      *(cMIIF->gdata[Vm]) += i_stim;
#endif // ifndef ALGEBRAIC

      if (stim_assign) {
        int iIon = 0;
        for (auto ion : species_list)
        {
          cMIIF->transmem_stim_species( -i_stim*(ratios_list[iIon]/ratio_sum), ion.c_str(), params.surface_to_volume_arg, stim_list, numNode);
          FPRINTF( WORLD stderr, "[stim_assign] @ t = %.2f ms; magnitude: %fpA/pF; type: %s\n", tmo.time, i_stim, ion.c_str() );
          iIon++;
        }
      }
    }

    if ( MIIF.gdata[illum] != NULL ) {
      MIIF.gdata[illum]->set(tmo.trigger(LIGHT_TM_IDX) ? light_irrad : 0.0);

      if(tmo.triggered_now(LIGHT_TM_IDX) && light_irrad) {
        FPRINTF( WORLD stderr, "[stim] illum ON  @ t = %.2f ms; E_e = %.3f mW/mm^2\n", tmo.time, light_irrad );
      }
      if(tmo.trigger_end(LIGHT_TM_IDX) && light_irrad) {
        FPRINTF( WORLD stderr, "[stim] illum OFF @ t = %.2f ms\n", tmo.time );
      }
    }

    // strain
    if (params.strain_given && cMIIF->use_stretch())
      apply_stretch(cMIIF, &s, &tmo);

    if (analyze_AP) {
      if (restitute) AP.pmat = r.trigs.pmat[AP.beat];
      check_events(getCellVal(cMIIF->gdata[Vm], 0), &AP, &tmo);
    }

    get_time(t5);


    if (is_gpu(target)) {
      if (gpu_Vm_update == 0) {
        cMIIF->compute_ionic_current(1, 1);
      } else if (gpu_Vm_update == 1) {
        cMIIF->compute_ionic_current(fabs(counter - 0.0) < epsilon || fabs(counter - 1.0) < epsilon  , fabs((counter+dt) - 1.0) < epsilon);
      }
    }
    else {
      cMIIF->compute_ionic_current();
    }

    // update Vm from Iion
    if (extVmUpdate)
      cMIIF->gdata[Vm]->add_scaled(*cMIIF->gdata[Iion], -cMIIF->dt);

    double ode_time = timing(t6, t5);
    update_timing(timings+ODE_IDX, ode_time);

    if (do_trace && tmo.trigger(TRACE_TM_IDX)) {
      dump_trace(cMIIF, t);
    }

    if (!restitute && tmo.trigger_end(DOPPLE_TM_IDX) ) {
      log_msg(NULL, 0, 0, "[%f] Switching back to original ion structure\n", tmo.time);
      cMIIF = &MIIF;
    }

    tmo.update_timers();

    double loop_time = timing(t7, t4);
    update_timing(timings+LOOP_IDX, loop_time);

    gpu_Vm_update = 1;
    counter = counter + dt;
    if ( fabs(counter - 1.0) < epsilon ) {
      counter = 0.0;
    }
  }//loop end

  // loop_current time
  double main_time = timing(t8, t3);

  // write header so that we know which files belong together
  const char *ExpID = params.fout_arg;
  write_dump_header(&gvd, &MIIF.svd, ExpID);

  // clean up everything
  if (tr_stim) free_trace(&stim_trace);
  if (analyze_AP) cleanup_AP_analysis(&AP);
  if (do_trace) close_trace(&MIIF);
  MIIF.close_svs_dumps();
  free_doppel(&doppel);
  MIIF.free_MIIF();

  close_globalvec_dump(fhdls, &gvd, &io);

  log_msg(NULL, 0, 0, "\n\n\nAll done!\n\n");

  // print timing data
  event_timing *t_setup = timings+SETUP_IDX;
  log_msg(NULL, 0, 0, "setup time          %.6f s",   (double)t_setup->tot);
  event_timing *t_init = timings+INIT_IDX;
  log_msg(NULL, 0, 0, "initialization time %.6f s",   (double)t_init->tot);
  event_timing *t_loop = timings+LOOP_IDX;
  log_msg(NULL, 0, 0, "main loop time      %.6f s",   (double)t_loop->tot);
  event_timing *t_ode = timings+ODE_IDX;
  log_msg(NULL, 0, 0, "total ode time      %.6f s\n", (double)t_ode->tot);

  log_msg(NULL, 0, 0, "mn/avg/mx loop time %.6f %.6f %.6f ms",
          (double)t_loop->mn*1e3, (double)t_loop->avg*1e3, (double)t_loop->mx*1e3);
  log_msg(NULL, 0, 0, "mn/avg/mx ODE time  %.6f %.6f %.6f ms",
          (double)t_ode->mn*1e3, (double)t_ode->avg*1e3, (double)t_ode->mx*1e3);

  log_msg(NULL, 0, 0, "real time factor    %.6f\n\n", duration/(t_ode->tot*1e3) );


  CHKERRQ(PetscFinalize());
  return 0;
} // main