ionics.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 "ionics.h"

#include "SF_init.h"

namespace opencarp {

void initialize_sv_dumps(limpet::MULTI_IF *pmiif, IMPregion* reg, int id, double t, double dump_dt);

void Ionics::compute_step()
{
  double t1, t2;
  get_time(t1);

  miif->compute_ionic_current();

  double comp_time = timing(t2, t1);
  this->compute_time += comp_time;

  comp_stats.calls++;
  comp_stats.tot_time += comp_time;

  if(user_globals::tm_manager->trigger(iotm_console))
    comp_stats.log_stats(user_globals::tm_manager->time, false);
}

void Ionics::destroy()
{
  miif->free_MIIF();
}

void Ionics::output_step()
{}

void Ionics::initialize()
{
  double t1, t2;
  get_time(t1);

  set_dir(OUTPUT);

  // initialize generic logger for ODE timings per time_dt
  comp_stats.init_logger("ODE_stats.dat");

  double tstart = 0.;
  sf_mesh & mesh = get_mesh(ion_domain);
  int loc_size = mesh.pl.num_algebraic_idx();

  // create ionic current vector and register it
  sf_vec *IIon;
  sf_vec *Vmv;
  SF::init_vector(&IIon, mesh, 1, sf_vec::algebraic);
  SF::init_vector(&Vmv, mesh, 1, sf_vec::algebraic);
  register_data(IIon, iion_vec);
  register_data(Vmv,  vm_vec);

  // setup miif
  miif = new limpet::MULTI_IF();
  // store IIF_IDs and Plugins in arrays
  miif->name = "myocardium";
  miif->gdata[limpet::Vm]   = Vmv;
  miif->gdata[limpet::Iion] = IIon;

  int rank = get_rank();

  SF::vector<RegionSpecs> rs(param_globals::num_imp_regions);
  for (size_t i = 0; i < rs.size(); i++ ) {
    rs[i].nsubregs   = param_globals::imp_region[i].num_IDs;
    rs[i].subregtags = param_globals::imp_region[i].ID;
    for (int j=0;j<rs[i].nsubregs;j++) {
        if(rs[i].subregtags[j]==-1 && rank==0)
          log_msg(NULL,3,ECHO, "Warning: not all %u IDs provided for imp_region[%u]!\n", rs[i].nsubregs, i);
    }
  }

  SF::vector<int> reg_mask;
  region_mask(ion_domain, rs, reg_mask, false, "imp_region");

  int purkfLen = 1;
  tstart = setup_MIIF(loc_size, param_globals::num_imp_regions, param_globals::imp_region,
                      reg_mask.data(), param_globals::start_statef, param_globals::num_adjustments,
                      param_globals::adjustment, param_globals::dt, purkfLen > 0);

  miif->extUpdateVm = !param_globals::operator_splitting;

  // if we start at a non-zero time (i.e. we have restored a state), we notify the
  // timer manager
  if(tstart > 0.) {
    log_msg(logger, 0, 0, "Changing simulation start time to %.2lf", tstart);
    user_globals::tm_manager->setup(param_globals::dt, tstart, param_globals::tend);
    user_globals::tm_manager->reset_timers();
  }

  set_dir(INPUT);

  this->initialize_time += timing(t2, t1);
}

double Ionics::setup_MIIF(int nnodes, int nreg, IMPregion* impreg, int* mask,
           const char *start_fn, int numadjust, IMPVariableAdjustment *adjust,
           double time_step, bool close)
{
  double tstart = 0;

  miif->N_IIF     = nreg;
  miif->numNode   = nnodes;
  miif->iontypes  = {};
  miif->numplugs  = (int*)calloc( miif->N_IIF, sizeof(int));
  miif->plugtypes = std::vector<limpet::IonTypeList>(miif->N_IIF);
  miif->targets = std::vector<limpet::Target>(miif->N_IIF, limpet::Target::AUTO);

  log_msg(logger,0,ECHO, "\nSetting up ionic models and plugins\n" \
                             "-----------------------------------\n\n" \
                             "Assigning IMPS to tagged regions:" );

  for (int i=0;i<miif->N_IIF;i++) {
    auto pT = limpet::get_ion_type(std::string(impreg[i].im));
    if (pT != NULL)
    {
      miif->iontypes.push_back(*pT);
      log_msg(logger, 0, ECHO|NONL, "\tIonic model: %s to tag region(s)", impreg[i].im);

      if(impreg[i].num_IDs > 0) {
        for(int j = 0; j < impreg[i].num_IDs; j++)
          log_msg(logger,0,ECHO|NONL, " [%d],", impreg[i].ID[j]);
        log_msg(logger,0,ECHO,"\b.");
      }
      else {
        log_msg(logger,0,ECHO, " [0] (implicitely)");
      }
    }
    else {
      log_msg(NULL,5,ECHO, "Illegal IM specified: %s\n", impreg[i].im );
      log_msg(NULL,5,ECHO, "Run bench --list-imps for a list of all available models.\n" );
      EXIT(1);
    }
    if (limpet::get_plug_flag( impreg[i].plugins, &miif->numplugs[i], miif->plugtypes[i]))
    {
      if(impreg[i].plugins[0] != '\0')  {
        log_msg(logger,0, ECHO|NONL, "\tPlug-in(s) : %s to tag region(s)",  impreg[i].plugins);

        for(int j = 0; j < impreg[i].num_IDs; j++)
          log_msg(logger,0,ECHO|NONL, " [%d],", impreg[i].ID[j]);
        log_msg(logger,0,ECHO,"\b.");
      }
    }
    else  {
      log_msg(NULL,5,ECHO,"Illegal plugin specified: %s\n", impreg[i].plugins);
      log_msg(NULL,5,ECHO, "Run bench --list-imps for a list of all available plugins.\n" );
      EXIT(1);
    }
  }

  miif->IIFmask = (limpet::IIF_Mask_t*)calloc(miif->numNode, sizeof(limpet::IIF_Mask_t));

  // The mask is a nodal vector with the region IDs of each node.
  // It is already reduced during the region_mask call to guarantee unique values
  // for overlapping interface nodes
  if (mask) {
    for (int i=0; i<miif->numNode; i++)
      miif->IIFmask[i] = (limpet::IIF_Mask_t) mask[i];
  }

  miif->initialize_MIIF();

  for(int i=0; i<miif->N_IIF; i++)
    if(!miif->IIF[i]->cgeom().SVratio)
      miif->IIF[i]->cgeom().SVratio  = param_globals::imp_region[i].cellSurfVolRatio;

  for (int i=0;i<miif->N_IIF;i++) {
    // the IMP tuning does not handle spaces well, thus we remove them here
    remove_char(impreg[i].im_param, strlen(impreg[i].im_param), ' ');
    miif->IIF[i]->tune(impreg[i].im_param, impreg[i].plugins, impreg[i].plug_param);
  }

  set_dir(INPUT);
  miif->initialize_currents(time_step, param_globals::ode_fac);

  // overriding initial values goes here
  // read in single cell state vector and spread it out over the entire region
  set_dir(INPUT);
  for (int i=0;i<miif->N_IIF;i++) {
    if (impreg[i].im_sv_init && strlen(impreg[i].im_sv_init) > 0)
      if (read_sv(miif, i, impreg[i].im_sv_init)) {
        log_msg(NULL, 5, ECHO|FLUSH, "State vector initialization failed for %s.\n", impreg[i].name);
        EXIT(-1);
      }
  }

  if( !start_fn || strlen(start_fn)>0 )
    tstart = (double) miif->restore_state(start_fn, ion_domain, close);

  for (int i=0; i<numadjust; i++)
  {
    set_dir(INPUT);

    SF::vector<SF_int>  indices;
    SF::vector<SF_real> values;
    bool restrict_to_algebraic = true;

    sf_mesh & imesh = get_mesh(ion_domain);
    std::map<std::string,std::string> metadata;
    read_metadata(adjust[i].file, metadata, PETSC_COMM_WORLD);

    SF::SF_nbr nbr = SF::NBR_REF;
    if(metadata.count("grid") && metadata["grid"].compare("intra") == 0) {
        nbr = SF::NBR_SUBMESH;
    }
    read_indices_with_data(indices, values, adjust[i].file, imesh, nbr, restrict_to_algebraic, 1, PETSC_COMM_WORLD);

    int rank = get_rank();

    for(size_t gi = 0; gi < indices.size(); gi++)
      indices[gi] = SF::local_nodal_to_local_petsc<mesh_int_t,mesh_real_t>(imesh, rank, indices[gi]);

    // debug, output parameters on global intracellular vector
    if(adjust[i].dump) {
      sf_vec* adjPars;
      SF::init_vector(&adjPars, imesh, 1, sf_vec::algebraic);
      adjPars->set(indices, values);

      set_dir(OUTPUT);
      char fname[2085];
      snprintf(fname, sizeof fname, "adj_%s_perm.dat", adjust[i].variable);
      adjPars->write_ascii(fname, false);

      // get the scattering to the canonical permutation
      SF::scattering* sc = get_permutation(ion_domain, PETSC_TO_CANONICAL, 1);
      if(sc == NULL) {
        log_msg(0,3,0, "%s warning: PETSC_TO_CANONICAL permutation needed registering!", __func__);
        sc = register_permutation(ion_domain, PETSC_TO_CANONICAL, 1);
      }

      (*sc)(*adjPars, true);
      snprintf(fname, sizeof fname, "adj_%s_canonical.dat", adjust[i].variable);
      adjPars->write_ascii(fname, false);
    }

    int nc = miif->adjust_MIIF_variables(adjust[i].variable, indices, values);
    log_msg(logger, 0, 0, "Adjusted %d values for %s", nc, adjust[i].variable);
  }

  set_dir(OUTPUT);

  for (int i=0;i<miif->N_IIF;i++)
    initialize_sv_dumps(miif, impreg+i, i, tstart, param_globals::spacedt);

  return tstart;
}

void initialize_sv_dumps(limpet::MULTI_IF *pmiif, IMPregion* reg, int id, double t, double dump_dt)
{
  char  svs[1024], plgs[1024], plgsvs[1024], fname[1024];

  strcpy(svs, reg->im_sv_dumps ? reg->im_sv_dumps : "");
  strcpy(plgs, reg->plugins ? reg->plugins : "");
  strcpy(plgsvs, reg->plug_sv_dumps ? reg->plug_sv_dumps : "");

  if( !(strlen(svs)+strlen(plgsvs) ) )
    return;

  /* The string passed to the "reg_name" argument (#4) of the sv_dump_add
  *     function is supposed to be "region name". It's only purpose is to
  *     provide the base name for the SV dump file, eg: Purkinje.Ca_i.bin.
  *     Thus, we pass:  [vofile].[reg name], Otherwise, dumping SVs in
  *     batched runs would be extremely tedious.
  */
  strcpy(fname, param_globals::vofile);                              // [vofile].igb.gz

  if( !reg->name ) {
    log_msg(NULL, 5, ECHO, "%s: a region name must be specified\n", __func__ );
    exit(0);
  }

  // We want to convert vofile.igb or vofile.igb.gz to vofile.regname
  char* ext_start = NULL;
  ext_start = strstr(fname, "igb.gz");
  if(ext_start == NULL) ext_start = strstr(fname, "igb");
  if(ext_start == NULL) ext_start = fname + strlen(fname);
  strcpy(ext_start, reg->name);

  pmiif->sv_dump_add_by_name_list(id, reg->im, fname, svs, plgs, plgsvs, t, dump_dt);
}

bool check_tags_in_elems(const SF::vector<int> & tags, SF::vector<RegionSpecs> & regspec,
                         const char* gridname, const char* reglist)
{
  bool AllTagsExist = true;

  hashmap::unordered_set<int> tagset;
  tagset.insert(tags.begin(), tags.end());

  // cycle through all user-specified regions
  for (size_t reg=0; reg<regspec.size(); reg++)
    // cycle through all tags which belong to the region
    for (int k=0; k<regspec[reg].nsubregs; k++) {
      // check whether this tag exists in element list
      int n = 0;
      if(tagset.count(regspec[reg].subregtags[k])) n++;
      // globalize n
      int N = get_global(n, MPI_SUM);
      if (N==0) {
        log_msg(NULL, 3, ECHO,
                "Region tag %d in %s[%d] not found in element list for %s grid.\n",
                regspec[reg].subregtags[k], reglist, reg, gridname);
        AllTagsExist = false;
      }
    }

  if (!AllTagsExist) {
    log_msg(NULL, 4, ECHO,"Missing element tags in element file!\n"
                          "Check region ID specs in input files!\n");
  }

  return AllTagsExist;
}

void region_mask(mesh_t meshspec, SF::vector<RegionSpecs> & regspec,
                 SF::vector<int> & regionIDs, bool mask_elem, const char* reglist)
{
  if(regspec.size() == 1) return;

  sf_mesh & mesh = get_mesh(meshspec);
  SF::element_view<mesh_int_t,mesh_real_t> eview(mesh, SF::NBR_PETSC);

  // initialize the list with the default regionID, 0
  size_t rIDsize = mask_elem ? mesh.l_numelem : mesh.l_numpts;

  size_t nelem = mesh.l_numelem;
  const SF::vector<mesh_int_t> & tags = mesh.tag;

  // check whether all specified tags exist in the element list
  check_tags_in_elems(tags, regspec, mesh.name.c_str(), reglist);

  regionIDs.assign(rIDsize, 0);
  hashmap::unordered_map<int,int> tag_to_reg;

  // we generate a map from tags to region IDs. This has many benefits, mainly we can check
  // whether a tag is assigned to multiple regions and simplify our regionIDs filling loop
  for (size_t reg=1; reg < regspec.size(); reg++) {
    int err = 0;
    for (int k=0; k < regspec[reg].nsubregs; k++) {
      int curtag = regspec[reg].subregtags[k];
      if(tag_to_reg.count(curtag)) err++;

      if(get_global(err, MPI_SUM))
        log_msg(0,4,0, "%s warning: Tag idx %d is assigned to multiple regions!\n"
                       "Its final assignment will be to the highest assigned region ID!",
                __func__, curtag);

      tag_to_reg[curtag] = reg;
    }
  }

  SF::vector<int> mx_tag(regionIDs.size(), -1);

  // cycle through the element list
  for(size_t i=0; i<nelem; i++) {
    const mesh_int_t & cur_tag = tags[i];
    // check if current element tag has a custom region ID
    if (tag_to_reg.count(cur_tag))
    {
      int reg = tag_to_reg[cur_tag];

      if (mask_elem) regionIDs[i] = reg;
      else {
        eview.set_elem(i);

        for (int j=0; j < eview.num_nodes(); j++) {
          mesh_int_t n = eview.node(j);
          if (cur_tag > mx_tag[n]) mx_tag[n] = cur_tag;
        }
      }
    }
  }

  if(mask_elem) return;

  for(size_t i=0; i < regionIDs.size(); i++)
    if(mx_tag[i] >= 0) regionIDs[i] = tag_to_reg[mx_tag[i]];

  // for POINT masks we need to do a nodal -> algebraic map
  mesh.pl.reduce(regionIDs, "max");

  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);

  SF::vector<int> alg_reg(alg_nod.size());
  SF::vector<int> gids   (alg_nod.size());

  for(size_t i=0; i<alg_nod.size(); i++) {
    alg_reg[i] = regionIDs[alg_nod[i]];
    gids[i]    = nbr[alg_nod[i]];
  }

  regionIDs = alg_reg;

  if(param_globals::dump_imp_region) {
    set_dir(OUTPUT);
    write_data_ascii(PETSC_COMM_WORLD, gids, regionIDs, std::string(reglist)+".dat");
  }
}

double Ionics::timer_val(const int timer_id)
{
  double val = std::nan("NaN");
  return val;
}

std::string Ionics::timer_unit(const int timer_id)
{
  std::string s_unit;
  return s_unit;
}

void compute_IIF(limpet::IonIfBase& pIF, limpet::GlobalData_t** impdata, int n)
{
  if (impdata[limpet::Iion] != NULL)
    impdata[limpet::Iion][n] = 0;

  pIF.for_each([&](limpet::IonIfBase& imp) {
    update_ts(&imp.get_tstp());
    imp.compute(n, n + 1, impdata);
  });
}

void* find_SV_in_IMP(limpet::MULTI_IF* miif, const int idx, const char *IMP, const char *SV,
                     int* offset, int* sz, int* plugin_idx)
{
  *plugin_idx = -1;

  if(strcmp(IMP, miif->iontypes[idx].get().get_name().c_str()) == 0) {
    return (void*) miif->iontypes[idx].get().get_sv_offset(SV, offset, sz);
  }
  else {
    for(int k=0; k<miif->numplugs[idx]; k++)
      if(strcmp(IMP, miif->plugtypes[idx][k].get().get_name().c_str()) == 0) {
        *plugin_idx = k;
        return (void*) miif->plugtypes[idx][k].get().get_sv_offset(SV, offset, sz);
      }
  }

  return NULL;
}


void alloc_gvec_data(const int nGVcs, const int nRegs,
                     GVecs *prmGVecs, gvec_data &glob_vecs)
{
  glob_vecs.nRegs = nRegs;

  if (nGVcs) {
    glob_vecs.vecs.resize(nGVcs);
    glob_vecs.plugin_idx.resize(nGVcs);

    for (size_t i = 0; i < glob_vecs.vecs.size(); i++) {
      sv_data &gvec = glob_vecs.vecs[i];

      gvec.name    = dupstr(prmGVecs[i].name);
      gvec.units   = dupstr(prmGVecs[i].units);
      gvec.bogus   = prmGVecs[i].bogus;

      gvec.imps    = (char**) calloc(nRegs, sizeof(char *));
      gvec.svNames = (char**) calloc(nRegs, sizeof(char *));
      gvec.svSizes = (int*)   calloc(nRegs, sizeof(int));
      gvec.svOff   = (int*)   calloc(nRegs, sizeof(int));
      gvec.getsv   = (void**) calloc(nRegs, sizeof(limpet::SVgetfcn));

      for (int j = 0; j < nRegs; j++) {
        if (strlen(prmGVecs[i].imp)) gvec.imps[j] = dupstr(prmGVecs[i].imp);
#ifdef WITH_PURK
        else if (j < param_globals::num_imp_regions)
          gvec.imps[j] = dupstr(param_globals::imp_region[j].im);
        else
          gvec.imps[j] = dupstr(param_globals::PurkIon[j - param_globals::num_imp_regions].im);
#else
        else if (j < param_globals::num_imp_regions)
          gvec.imps[j] = dupstr(param_globals::imp_region[j].im);
#endif
        gvec.svNames[j] = dupstr(prmGVecs[i].ID[j]);
      }
    }
  }
}

void init_sv_gvec(gvec_data& GVs, limpet::MULTI_IF* miif, sf_vec & tmpl,
                  igb_output_manager & output_manager)
{
  GVs.inclPS = false;
  // int num_purk_regions = GVs->inclPS ? purk->ion.N_IIF : 0;
  int num_purk_regions = 0;
  int nRegs            = param_globals::num_imp_regions + num_purk_regions;

  alloc_gvec_data(param_globals::num_gvecs, nRegs, param_globals::gvec, GVs);

#ifdef WITH_PURK
  if (GVs->inclPS) sample_PS_ionSVs(purk);
#endif

  for (unsigned int i = 0; i < GVs.vecs.size(); i++) {
    sv_data & gv = GVs.vecs[i];
    int noSV    = 0;

    for (int j = 0; j < miif->N_IIF; j++) {
      gv.getsv[j] = find_SV_in_IMP(miif, j, gv.imps[j], gv.svNames[j], gv.svOff + j, gv.svSizes + j, &GVs.plugin_idx[i]);

      if (gv.getsv[j] == NULL) {
        log_msg(NULL, 3, ECHO, "\tWarning: SV(%s) not found in region %d\n", gv.svNames[j], j);
        noSV++;
      }
    }

    SF::init_vector(&gv.ordered, &tmpl);
    output_manager.register_output(gv.ordered, intra_elec_msh, 1, gv.name, gv.units);

#ifdef WITH_PURK
    // same procedure for Purkinje
    if (GVs->inclPS) {
      IF_PURK_PROC(purk) {
        MULTI_IF* pmiif = &purk->ion;
        for (int j = miif->N_IIF; j < nRegs; j++) {
          gv.getsv[j] = find_SV_in_IMP(pmiif, j - miif->N_IIF, gv.imps[j],
                                       gv.svNames[j], gv.svOff + j, gv.svSizes + j);

          if (gv.getsv[j] == NULL) {
            LOG_MSG(NULL, 3, ECHO, "\tWarning: state variable \"%s\" not found in region %d\n", gv.svNames[j], j);
            noSV++;
          }
        }
        RVector_dup(purk->vm_pt, &gv.orderedPS_PS);
        MYO_COMM(purk);
      }
      RVector_dup(purk->vm_pt_over, &gv.orderedPS);
      initialize_grid_output(grid, NULL, tmo, intra_elec_msh, GRID_WRITE, 0., 1., gv.units, gv.orderedPS,
                             1, gv.GVcName, -purk->npt, param_globals::output_level);
    }
#endif

    if (noSV == nRegs) {
      log_msg(NULL, 5, ECHO, "\tError: no state variables found for global vector %d\n", i);
      log_msg(NULL, 5, ECHO, "Run bench --imp=YourModel --imp-info to get a list of all parameters.\\n");
      exit(1);
    }
  }
}

void assemble_sv_gvec(gvec_data & gvecs, limpet::MULTI_IF *miif)
{
  for(size_t i=0; i<gvecs.vecs.size(); i++ ) {
    sv_data & gv         = gvecs.vecs[i];
    int       plugin_idx = gvecs.plugin_idx[i];

    // set to the defalt value
    gv.ordered->set(gv.bogus);
    SF_int start, stop;
    gv.ordered->get_ownership_range(start, stop);

    for( int n = 0; n<miif->N_IIF; n++ ) {
      if( !gv.getsv[n] ) continue;

      SF::vector<SF_real> data   (miif->N_Nodes[n]);
      SF::vector<SF_int>  indices(miif->N_Nodes[n]);

      limpet::IonIfBase* base_im = miif->IIF[n];
      limpet::IonIfBase* imp     = plugin_idx > -1 ? base_im->plugins()[plugin_idx] : base_im;

      for( int j=0; j<miif->N_Nodes[n]; j++ ) {
        indices[j] = miif->NodeLists[n][j] + start;
        data[j]    = ((limpet::SVgetfcn)(gv.getsv[n]))( *miif->IIF[n], j, gv.svOff[n]);
      }

      bool add = false;
      gv.ordered->set(indices, data, add);
    }
  }
}


} // namespace opencarp