stimulate.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 "stimulate.h"
#include "electrics.h"
 
#include "SF_init.h" // for SF::init_xxx()
#include "petsc_utils.h" // TODO for EXIT
 
namespace opencarp {
 
// stimulus convenience functions
// group stimulus types
namespace stim_info {
  std::set<int> potential_stim {Phi_ex, GND_ex, Phi_ex_ol, Phi_in, Phi_in_ol};
  std::set<int> current_stim   {I_tm, I_ex, I_in};
  std::set<int> dbc_stim       {Phi_ex, GND_ex, Phi_ex_ol, Phi_in, Phi_in_ol};
  std::map<int,std::string> units;
}  // namespace stim_info
 
namespace user_globals {
 
  extern timer_manager* tm_manager;
}  // namespace user_globals
 
void init_stim_info(void)
{
  stim_info::units[I_tm]      = "uA/cm^2";  
  stim_info::units[I_ex]      = "uA/cm^3";  
  stim_info::units[Phi_ex]    = "mV";       
  stim_info::units[GND_ex]    = "mV";       
  stim_info::units[I_in]      = "uA/cm^3";  
  stim_info::units[Phi_ex_ol] = "mV";       
  // stim_info::units[I_tm_grad] = "";          //!<
  // stim_info::units[I_lat]     = "uA/cm^2";   //!< acts as transmembrane current stimulus
  stim_info::units[Vm_clmp]   = "mV";  
  stim_info::units[Phi_in]    = "mV";  
  stim_info::units[Phi_in_ol] = "mV";  
}
 
 
bool is_potential(stim_t type)        
{
  return stim_info::potential_stim.count(type);
}
 
bool is_current(stim_t type)
{
  return stim_info::current_stim.count(type);
}
 
bool is_dbc(stim_t type)
{
  return stim_info::dbc_stim.count(type);
}
 
bool is_extra(stim_t type)
{
  switch(type) {
    default:
    case Illum:
    case I_tm:
    case Vm_clmp:
    case Phi_in:
    case Phi_in_ol:
    case I_in:
      return false;
 
    case I_ex:
    case Phi_ex:
    case Phi_ex_ol:
    case GND_ex:
      return true;
  }
}
 
 
void stimulus::translate(int iidx)
{
  idx = iidx;
 
  // translate protocol
  Stim     & s       = param_globals::stim[idx];
  Stimulus & curstim = param_globals::stimulus[idx];
 
  s.ptcl.bcl      = curstim.bcl;
  s.ptcl.npls     = curstim.npls;
  s.ptcl.start    = curstim.start;
  s.ptcl.duration = curstim.duration;
 
  // translate pulse
  // create a name, old stimulus structure lacks this feature, no naming of pulse waveforms
  std::string pulse_name;
  pulse_name = "pulse_" + std::to_string(idx);
  s.pulse.name = (char*) calloc(pulse_name.length()+1, sizeof(char));
  strcpy(s.pulse.name, pulse_name.c_str());
 
  // old stimulus treats any pulse as truncated exponential shape
  s.pulse.shape = truncExpPulse;
  s.pulse.file  = strdup(curstim.pulse_file);
  s.pulse.bias  = curstim.bias;
  s.pulse.tau_plateau = curstim.tau_plateau;
  s.pulse.tau_edge    = curstim.tau_edge;
  s.pulse.s2          = curstim.s2;
  s.pulse.d1          = curstim.d1;
  s.pulse.strength    = curstim.strength;
 
  // strength is not a scalar, we prescribe a different strength at every node
  s.elec.vtx_fcn      = curstim.vtx_fcn;
 
  // translate circuit
  s.crct.type = curstim.stimtype;
  s.crct.balance = curstim.balance;
  s.crct.total_current = curstim.total_current;
 
  // translate electrode
  s.elec.geomID   = curstim.geometry;
  // s.elec.domain   = curstim.domain;
  s.elec.vtx_file = strdup(curstim.vtx_file);
  s.elec.dump_vtx_file = curstim.dump_vtx_file;
 
  // check whether we have a sound electrode geometry definition
  s.elec.geom_type = 2; // should be region_t block;
  s.elec.p0[0] = curstim.x0 - (curstim.ctr_def?curstim.xd/2.:0.);
  s.elec.p0[1] = curstim.y0 - (curstim.ctr_def?curstim.yd/2.:0.);
  s.elec.p0[2] = curstim.z0 - (curstim.ctr_def?curstim.zd/2.:0.);
  s.elec.p1[0] = s.elec.p0[0] + curstim.xd;
  s.elec.p1[1] = s.elec.p0[1] + curstim.yd;
  s.elec.p1[2] = s.elec.p0[2] + curstim.zd;
 
  // cp overall electrode name
  s.name = strdup(curstim.name);
}
 
void stimulus::setup(int iidx)
{
  idx = iidx;
  const Stim & cur_stim = param_globals::stim[idx];
 
  // label stimulus
  if(cur_stim.name != std::string(""))
    name = cur_stim.name;
  else
    name = "Stimulus_" + std::to_string(idx);
 
  // set up physics of applied stimulus
  phys.setup(idx, associated_intra_mesh, associated_extra_mesh);
 
  // set up stimulus pulse
  pulse.setup(idx);
 
  // set up the spatial region of the stimulus
  electrode.setup(idx, associated_intra_mesh, associated_extra_mesh);
 
  // set up stimulation protocol
  ptcl.setup(idx,name);
}
 
bool stimulus::is_active() const
{
  // closed loop intra-/extracellular potential stimuli are always active
  if(this->phys.type == Phi_ex || this->phys.type == Phi_in)
    return true;
 
  return user_globals::tm_manager->trigger(ptcl.timer_id);
}
 
 
void stim_pulse::setup(int idx)
{
  const Stim & cur_stim = param_globals::stim[idx];
  const Pulse & pulse = cur_stim.pulse;
 
  // assign waveform type and parameters
  waveform_t wvt;
  // set waveform based on circuit type/file/param
  switch (cur_stim.crct.type) {
    case Vm_clmp:
    case GND_ex:
      wvt = constPulse;
      break;
 
    default:
      if (pulse.file && strlen(pulse.file))
        wvt = arbPulse;
      else
        wvt = static_cast<waveform_t>(pulse.shape);
      break;
  }
 
  assign(pulse.strength, cur_stim.ptcl.duration, param_globals::dt, wvt);
  sample_wave_form(*this, idx);
 
  // add unit labels
  std::string sUnit = stim_info::units[(stim_t)cur_stim.crct.type];
  std::string tUnit = "ms";
  this->wave.setUnits(tUnit,sUnit);
}
 
 
void get_stim_list(const char* str_list, std::vector<double> & stlist)
{
  std::string tstim = str_list;
  std::vector<std::string> stims;
  split_string(tstim, ',', stims);
 
  stlist.resize(stims.size());
  size_t widx = 0, err = 0;
 
  for(size_t i=0; i<stims.size(); i++) {
    std::string & stim = stims[i];
 
    if(stim.size()) {
      if(stim[0] == '+') {
        double last = widx == 0 ? 0.0 : stlist[widx-1];
        stim.erase(0, 1); // erase first char
 
        // we use sscanf so that we can actually determine if the string was converted
        // into a double. atof does not allow for error-checking
        double cur = -1.0;
        size_t nread = sscanf(stim.c_str(), "%lf", &cur);
        if(nread)
          stlist[widx++] = cur + last;
      } else {
        double cur = -1.0;
        size_t nread = sscanf(stim.c_str(), "%lf", &cur);
        if(nread)
          stlist[widx++] = cur;
      }
    }
  }
 
  stlist.resize(widx);
 
  if(widx != stims.size())
    log_msg(0, 4, 0, "%s warning: Some values of %s could not be converted in stim times!",
            __func__, str_list);
}
 
 
void stim_protocol::setup(int idx, std::string name)
{
  const Stim & cur_stim = param_globals::stim[idx];
  const Protocol & ptcl = cur_stim.ptcl;
 
  if(strlen(ptcl.stimlist) == 0) {
    start = ptcl.start, npls = ptcl.npls, pcl = ptcl.bcl;
 
    // add timer
    timer_id = user_globals::tm_manager->add_eq_timer(ptcl.start,
               param_globals::tend, npls, pcl, ptcl.duration, name.c_str());
  } else {
    std::vector<double> stims;
    get_stim_list(ptcl.stimlist, stims);
 
    log_msg(0,0,ECHO | NONL, "stim %d: stimulating at times: ", idx);
    for(double t : stims)
      log_msg(0,0, ECHO | NONL, "%.2lf ", t);
    log_msg(0,0,0, "");
    log_msg(0,0,0, "stim %d: .bcl and .npls values will be ignored.", idx);
 
    timer_id = user_globals::tm_manager->add_neq_timer(stims, ptcl.duration, name.c_str());
  }
}
 
 
void stim_physics::setup(int idx, mesh_t intra_mesh, mesh_t extra_mesh)
{
  const Stim & cur_stim = param_globals::stim[idx];
 
  type = stim_t(cur_stim.crct.type);
  unit = stim_info::units[cur_stim.crct.type];
  total_current = cur_stim.crct.total_current;
 
  const short dim = is_extra(stim_t(cur_stim.crct.type)) ?
                    get_mesh_dim(extra_mesh) : get_mesh_dim(intra_mesh);
 
  switch(type) {
    case I_tm:
        scale = param_globals::operator_splitting ? user_globals::tm_manager->time_step : 1.0;
      break;
 
    case I_in:
    case I_ex:
      scale = UM2_to_CM2 * (dim == 2 ? 1.0 : UM_to_CM);
      break;
 
    default: scale = 1.0; break;
  }
}
 
 
void sample_wave_form(stim_pulse& sp, int idx)
{
  const Stim & cur_stim = param_globals::stim[idx];
  const Pulse & pulse = cur_stim.pulse;
 
  switch (sp.wform) {
    case squarePulse: {
      // generate step function using derived class
      sig::Step stepPars;
      stepPars.trig = pulse.trig;
      stepPars.rise = true;
 
      sig::stepFunc step_func(stepPars);
      step_func.sample(sp.wave);
      break;
    }
    case truncExpPulse: {
      // monophasic pulse, first phase lasts for entire pulse duration
      sig::Pulse truncExp;
      truncExp.duration = sp.duration;
      truncExp.d1       = pulse.d1 * sp.duration;  // relative to absolute duration of d1
      truncExp.tau_edge = pulse.tau_edge;
      truncExp.tau_plat = pulse.tau_plateau;
      truncExp.s2r      = pulse.s2;
      truncExp.bias     = pulse.bias;
 
      // from relative duration of d1 we infer mono- or biphasic
      if (pulse.d1 == 1.) {
        // monophasic pulse, first phase lasts for entire pulse duration
        sig::monophasicTruncExpFunc fnc(truncExp);
        fnc.sample(sp.wave);
      } else {
        // biphasic pulse, second phase lasts shorter than pulse duration
        sig::biphasicTruncExpFunc fnc(truncExp);
        fnc.sample(sp.wave);
      }
      break;
    }
    case sinePulse: {
      // generate a sine wave
      sig::sineWave sinePars;
      sinePars.frq   = pulse.freq;
      sinePars.phase = pulse.phase;
      sig::sineFunc sine(sinePars);
      sine.sample(sp.wave);
      break;
    }
    case arbPulse: {
      sig::time_trace file_wave;
      int             _err = 0;
 
      if (!get_rank()) {
        bool unitize = true;
 
        update_cwd();
        set_dir(INPUT);
        _err = file_wave.read_trace(pulse.file, unitize);
        set_dir(CURDIR);
      }
 
      int err = get_global(_err, MPI_MIN);
      if (err) {
        log_msg(user_globals::physics_reg[elec_phys]->logger, MAX_LOG_LEVEL, ECHO,
                "Failed reading pulse for stimulus[%d] from file %s.\n", idx, pulse.file);
      } else {
        get_global(file_wave.f);
        get_global(file_wave.t);
 
        file_wave.resample(sp.wave);
 
        // check length of read in pulse and adjust stimulus settings accordingly
        if (sp.wave.duration() != sp.duration)
          sp.duration = sp.wave.duration();
      }
      break;
    }
    case constPulse: {
      sig::constFunc cnst_fnc(1);
      cnst_fnc.sample(sp.wave);
      break;
    }
    case unsetPulse:
    default:
      assert(0);
  }
 
  /*
  if(sp.wform == APfootPulse)
  {
    // generate AP foot
    APfoot footPars;
    footPars.tau_f = s.duration()/5.;   // use 5 time constants
    APfootFunc APft(footPars);
    s.set_labels("APfoot");
    APft.sample(s);
    s.write_trace();
  }
  */
}
 
bool stimulus::value(double & v) const
{
  if(user_globals::tm_manager->trigger(ptcl.timer_id) == false) {
    // intra-/extracellular potential stimuli are always active, but return 0 when their
    // protocol is not active. Use Phi_ex_ol for extracellular potential stimuli that
    // get removed when they expire.
    if(this->phys.type == Phi_ex || this->phys.type == GND_ex || this->phys.type == Phi_in) {
      v = 0.0;
      return true;
    }
 
    return false;
  }
  else {
    int i = user_globals::tm_manager->trigger_elapse(ptcl.timer_id);
    v = pulse.strength * phys.scale * pulse.wave.f[i];
    return true;
  }
}
 
void stimulus::dump_vtx_file(int idx)
{
  const sf_mesh& mesh = is_extra(phys.type) ? get_mesh(associated_extra_mesh) : get_mesh(associated_intra_mesh);
 
  SF::vector<mesh_int_t> glob_idx(electrode.vertices);
  mesh.pl.globalize(glob_idx);
 
  SF::vector<mesh_int_t> srt_idx;
  SF::sort_parallel(mesh.comm, glob_idx, srt_idx);
 
  size_t num_vtx = get_global(srt_idx.size(), MPI_SUM);
  int    rank    = get_rank();
 
  FILE_SPEC f   = NULL;
  int       err = 0;
 
  if (rank == 0) {
    char dump_name[1024];
 
    if (strlen(name.c_str())) {
      snprintf(dump_name, sizeof dump_name, "%s.vtx", name.c_str());
    } else {
      snprintf(dump_name, sizeof dump_name, "ELECTRODE_%d.vtx", idx);
    }
 
    f = f_open(dump_name, "w");
    if (!f)
      err++;
    else {
      fprintf(f->fd, "%zd\nextra\n", num_vtx);
    }
  }
 
  if (!get_global(err, MPI_SUM)) {
    print_vector(mesh.comm, srt_idx, 1, f ? f->fd : NULL);
  } else {
    log_msg(0, 4, 0, "error: stimulus[%d] cannot be dumped!");
  }
 
  // only root really does that
  f_close(f);
}
 
void stim_electrode::setup(int idx, mesh_t intra_mesh, mesh_t extra_mesh)
{
  const Stim & curstim = param_globals::stim[idx];
  mesh_t grid;
  dump_vtx = curstim.elec.dump_vtx_file;
 
  // based on the stimulus domain type we select grid and physics
  switch(curstim.crct.type) {
    case I_ex:
    case Phi_ex:
    case Phi_ex_ol:
    case GND_ex:
      grid   = extra_mesh; break;
 
    case Illum:
    case Vm_clmp:
    case I_in:
    case I_tm:
    case Phi_in:
    case Phi_in_ol:
      grid   = intra_mesh; break;
 
    case Ignore_stim: return;
 
    default:
      log_msg(0,5,0, "stim_electrode::setup error: Can't determine domain from stim type! Aborting!", __func__);
      EXIT(1);
  }
 
  // get a logger from the physics
  FILE_SPEC logger = NULL;
 
  // this are the criteria for choosing how we extract the stimulus electrode
  bool vertex_file_given = strlen(curstim.elec.vtx_file) > 0;
  bool tag_index_given   = curstim.elec.geomID > -1;
  // the mesh we need for computing the local vertex indices.
  const sf_mesh & mesh   = get_mesh(grid);
 
  if(vertex_file_given && tag_index_given)
    log_msg(0,3,0, "%s warning: More than one stimulus electrode definintions set in electrode %d", __func__, idx);
 
  if(vertex_file_given) {
    definition     = def_t::file_based;
    input_filename = curstim.elec.vtx_file;
 
    log_msg(logger, 0, 0, "Stimulus %d: Selecting vertices from file %s", idx, input_filename.c_str());
 
    set_dir(INPUT);
    warn_when_passing_intra_vtx(input_filename);
 
    // we read the indices. they are being localized w.r.t. the provided numbering. In our
    // case this is always the reference numbering.
    if(curstim.elec.vtx_fcn)
      read_indices_with_data(vertices, scaling, input_filename, mesh, SF::NBR_REF, true, 1, PETSC_COMM_WORLD);
    else
      read_indices(vertices, input_filename, mesh, SF::NBR_REF, true, PETSC_COMM_WORLD);
 
    int gnum_idx = get_global(vertices.size(), MPI_SUM);
    if(gnum_idx == 0) {
      log_msg(0, 5, 0, "Stimulus %d: Specified vertices are not in stimulus domain! Aborting!", idx);
      EXIT(1);
    }
  }
  else if(tag_index_given) {
    definition = def_t::vol_based_tag;
 
    int tag = curstim.elec.geomID;
    log_msg(logger, 0, 0, "Stimulus %d: Selecting vertices from tag %d", idx, tag);
 
    indices_from_region_tag(vertices, mesh, tag);
    // we restrict the indices to the algebraic subset
    SF::restrict_to_set(vertices, mesh.pl.algebraic_nodes());
  }
  else {
    definition = def_t::vol_based_shape;
    log_msg(logger, 0, 0, "Stimulus %d: Selecting vertices from shape.", idx);
 
    geom_shape shape;
    shape.type   = geom_shape::shape_t(curstim.elec.geom_type);
    shape.p0     = curstim.elec.p0;
    shape.p1     = curstim.elec.p1;
    shape.radius = curstim.elec.radius;
 
    bool nodal = true;
    indices_from_geom_shape(vertices, mesh, shape, nodal);
    // we restrict the indices to the algebraic subset
    SF::restrict_to_set(vertices, mesh.pl.algebraic_nodes());
 
    int gsize = vertices.size();
    if(get_global(gsize, MPI_SUM) == 0) {
      log_msg(0,5,0, "error: Empty stimulus[%d] electrode def! Aborting!", idx);
      EXIT(1);
    }
  }
}
 
void dbc_manager::recompute_dbcs()
{
  clear_active_dbc();
 
  for(const stimulus & s : stimuli) {
    if(is_dbc(s.phys.type) && s.is_active()) {
 
      // the local indices of the dbc
      const SF::vector<mesh_int_t> & dbc_idx = s.electrode.vertices;
 
      dbc_data* dbc_buff = new dbc_data();
 
      // the global petsc indices of the dbc
      dbc_buff->nod = new SF::vector<SF_int>(dbc_idx.size());
 
      const SF::vector<mesh_int_t> & petsc_nbr = mat.mesh_ptr()->get_numbering(SF::NBR_PETSC);
 
      size_t widx = 0;
      for(size_t i=0; i<dbc_idx.size(); i++)
        (*dbc_buff->nod)[widx++] = petsc_nbr[dbc_idx[i]];
 
      sf_vec* dbc; SF::init_vector(&dbc, *mat.mesh_ptr(), 1, sf_vec::algebraic);
      SF::init_vector(&dbc_buff->cntr, *mat.mesh_ptr(), 1, sf_vec::algebraic);
 
      if(s.electrode.scaling.size()) {
        // we have spatially heterogenous dirichlet values
        dbc->set(*dbc_buff->nod, s.electrode.scaling);
      } else {
        // we have spatially homogenous dirichlet values
        dbc->set(*dbc_buff->nod, 1.0);
      }
 
      mat.mult(*dbc, *dbc_buff->cntr);
 
      active_dbc[s.idx] = dbc_buff;
    }
  }
}
 
bool dbc_manager::dbc_update()
{
  // check if stim has come offline
  for(const auto & d : active_dbc) {
    if(stimuli[d.first].is_active() == false)
      return true;
  }
 
  // check if stim has come online
  for(const stimulus & s : stimuli) {
    if(is_dbc(s.phys.type) && s.is_active() && active_dbc.count(s.idx) == 0)
      return true;
  }
 
  return false;
}
 
void dbc_manager::enforce_dbc_lhs()
{
  sf_vec* dbc; 
 
  SF::init_vector(&dbc, *mat.mesh_ptr(), mat.dpn_row(), sf_vec::algebraic);
 
  for(auto it = active_dbc.begin(); it != active_dbc.end(); ++it) {
    const SF::vector<SF_int> & dbc_nod = *it->second->nod;
    for (auto dbc_idx : dbc_nod) {
        mat.boundary_idxs.insert(dbc_idx);
    }
    // zero cols and rows in lhs mat
    dbc->set(1.0);
    dbc->set(dbc_nod, 0.0);
    // multiplicate dbc vector from left and right to matrix in order to zero cols and rows
    mat.mult_LR(*dbc, *dbc);
    // now add 1.0 to the diagonal at the dirichlet locations
    dbc->set(0.0);
    dbc->set(dbc_nod, 1.0);
    mat.diag_add(*dbc);
  }
}
 
void dbc_manager::enforce_dbc_rhs(sf_vec &rhs)
{
  for(auto it = active_dbc.begin(); it != active_dbc.end(); ++it) {
    const int & dbc_idx = it->first;
    const dbc_manager::dbc_data & dbc = *(it->second);
    double strength = 0.0;
    bool is_active = stimuli[dbc_idx].value(strength);
 
    // the whole idea of the dbc_manager and its active_dbc is
    // based on only storing the active dbcs. if the stimulus associated
    // to a dbc believed active is actually not active, we want to abort.
    assert(is_active);
 
    rhs.add_scaled(*dbc.cntr, -strength);
 
    if(stimuli[dbc_idx].electrode.scaling.size()) {
      SF::vector<SF_real> scaled_strength(stimuli[dbc_idx].electrode.scaling.size());
 
      size_t widx = 0;
      for(SF_real s : stimuli[dbc_idx].electrode.scaling)
        scaled_strength[widx++] = s * strength;
 
      rhs.set(*dbc.nod, scaled_strength);
    } else {
      rhs.set(*dbc.nod, strength);
    }
  }
}
 
}  // namespace opencarp