ionics.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#include "ionics.h"
 
#include "SF_init.h"
 
namespace opencarp {
 
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);
  limpet::node_count_t 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();
  // hand down the physics log file, so that what MULTI_IF reports while setting up and
  // restoring state is recorded rather than only echoed to the console
  miif->logger = logger;
  // store IIF_IDs and Plugins in arrays
  miif->name = "myocardium";
  miif->gdata[limpet::Vm]   = Vmv;
  miif->gdata[limpet::Iion] = IIon;
 
  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;
  }
 
  SF::vector<int> reg_mask;
  region_mask(ion_domain, rs, reg_mask, false, "imp_region", logger);
 
  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(limpet::node_count_t 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);
 
      // a single region covers the whole mesh regardless of its ID list
      if(miif->N_IIF == 1) {
        log_msg(logger,0,ECHO, " all");
      }
      else 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, " none");
      }
    }
    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 (limpet::node_index_t 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, false, true);
 
      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);
  }
 
  return tstart;
}
 
bool check_tags_in_elems(const SF::vector<mesh_int_t> & tags, SF::vector<RegionSpecs> & regspec,
                         const char* gridname, const char* reglist, FILE_SPEC logger)
{
  bool AllTagsExist = true;
 
  hashmap::unordered_set<mesh_int_t> 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) {
        if(strcmp(reglist, "gregion_vol"))
          log_msg(logger, 3, ECHO,
                "%s[%d] references tag %d, but no element in the %s grid carries this tag — region will have no elements.\n",
                reglist, reg, regspec[reg].subregtags[k], gridname);
        AllTagsExist = false;
      }
    }
 
  if (!AllTagsExist) {
    log_msg(logger, 4, ECHO,"One or more configured regions are empty. Check that region tag IDs match the tags in your mesh or tagfile.\n");
  }
 
  return AllTagsExist;
}
 
static const char* one_region_per_tag = "Every element tag must be listed in exactly one region.\n";
 
void check_unassigned_tags(const SF::vector<mesh_int_t> & tags,
                           const hashmap::unordered_map<mesh_int_t,int> & tag_to_reg,
                           const char* gridname, const char* reglist, FILE_SPEC logger)
{
  SF::vector<mesh_int_t> unmatched;
  long int n_unassigned = 0;
  for (const mesh_int_t & t : tags) {
    if (!tag_to_reg.count(t)) {
      unmatched.push_back(t);
      n_unassigned++;
    }
  }
 
  n_unassigned = get_global(n_unassigned, MPI_SUM);
  if (!n_unassigned) return;
 
  binary_sort(unmatched);
  unique_resize(unmatched);
  make_global(unmatched, PETSC_COMM_WORLD);
  binary_sort(unmatched);
  unique_resize(unmatched);
 
  std::string taglist;
  for (size_t i = 0; i < unmatched.size(); i++) {
    if (i) taglist += ", ";
    taglist += std::to_string(unmatched[i]);
  }
 
  log_msg(logger, 5, ECHO|FLUSH,
    "%s: %ld element(s) in %s grid carry tags {%s} not assigned to any %s region.\n%s",
    __func__, n_unassigned, gridname, taglist.c_str(), reglist, one_region_per_tag);
  EXIT(EXIT_FAILURE);
}
 
void region_mask(mesh_t meshspec, SF::vector<RegionSpecs> & regspec,
                 SF::vector<int> & regionIDs, bool mask_elem, const char* reglist,
                 FILE_SPEC logger)
{
  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);
 
  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;
 
  // the parameter parser fills IDs missing from a list shorter than num_IDs with -1
  for (size_t reg=0; reg < regspec.size(); reg++)
    for (int k=0; k < regspec[reg].nsubregs; k++)
      if (regspec[reg].subregtags[k] == -1) {
        log_msg(logger, 5, ECHO|FLUSH, "%s: %s[%zu] lists fewer tags than its num_IDs.\n",
                __func__, reglist, reg);
        EXIT(EXIT_FAILURE);
      }
 
  // check whether all specified tags exist in the element list
  check_tags_in_elems(tags, regspec, mesh.name.c_str(), reglist, logger);
 
  // A tag listed in two regions has no well-defined owner, so the region layout is rejected
  // rather than resolved by list order. regspec is identical on all ranks, so every rank
  // reaches the same verdict.
  hashmap::unordered_map<mesh_int_t,int> tag_to_reg;
  for (size_t reg=0; reg < regspec.size(); reg++) {
    for (int k=0; k < regspec[reg].nsubregs; k++) {
      mesh_int_t curtag = regspec[reg].subregtags[k];
      auto it = tag_to_reg.find(curtag);
      if (it != tag_to_reg.end() && it->second != int(reg)) {
        log_msg(logger, 5, ECHO|FLUSH, "%s: tag %ld is listed in both %s[%d] and %s[%zu].\n%s",
                __func__, long(curtag), reglist, it->second, reglist, reg, one_region_per_tag);
        EXIT(EXIT_FAILURE);
      }
      tag_to_reg[curtag] = reg;
    }
  }
 
  check_unassigned_tags(tags, tag_to_reg, mesh.name.c_str(), reglist, logger);
 
  regionIDs.resize(rIDsize);
 
  SF::vector<mesh_int_t> 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];
 
    if (mask_elem) regionIDs[i] = tag_to_reg.at(cur_tag);
    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;
 
  // The tie between the tags meeting at a node must be broken on the tags themselves.
  // Reducing the region IDs instead would let the highest region win at nodes shared
  // between ranks while the highest tag wins everywhere else, making the assignment
  // depend on the partitioning whenever the tags are not listed in ascending order.
  mesh.pl.reduce(mx_tag, "max");
 
  for(size_t i=0; i < regionIDs.size(); i++)
    regionIDs[i] = tag_to_reg.at(mx_tag[i]);
 
  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, limpet::node_index_t 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);
 
    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.plugin_idx = (int*)   calloc(nRegs, sizeof(int));
      gvec.getsv      = (void**) calloc(nRegs, sizeof(limpet::SVgetfcn));
 
      for (int j = 0; j < nRegs; j++) {
        gvec.plugin_idx[j] = -1;
        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, gv.plugin_idx + j);
 
      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, gv.plugin_idx + 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];
 
    // 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];
      int                plugin_idx = gv.plugin_idx[n];
      limpet::IonIfBase* imp     = plugin_idx > -1 ? base_im->plugins()[plugin_idx] : base_im;
 
      for( limpet::node_index_t j=0; j<miif->N_Nodes[n]; j++ ) {
        indices[j] = miif->NodeLists[n][j] + start;
        data[j]    = ((limpet::SVgetfcn)(gv.getsv[n]))( *imp, j, gv.svOff[n]);
      }
 
      bool add = false;
      gv.ordered->set(indices, data, add);
    }
  }
}
 
 
} // namespace opencarp