IonicSolver.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   CEPS is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "ode/solver/IonicSolver.hpp"
#include "ode/solver/AllOdeSolvers.hpp"
 
IonicSolver::IonicSolver(
  CepsString              solverOpts, 
  AbstractIonicModel*     m,
  AbstractDiscretization* sDiscr
):
  m_model        (m),
  m_fLin         (nullptr),
  m_fNLin        (nullptr),
  m_sDiscr       (sDiscr),
  m_solver       (nullptr),
  m_y            (nullptr),
  m_dtyLin       (nullptr),
  m_dtyNLin      (nullptr),
  m_v            (nullptr),
  m_dtv          (nullptr),
  m_tauOk        (false),
  m_localToGlobal({}),
  m_size         (0)
{
  CEPS_ABORT_IF(ceps::isNullPtr(m),
    "Pointer to ionic model is null"
  );
  CEPS_ABORT_IF(ceps::isNullPtr(sDiscr),
    "Pointer to spatial discretization is null"
  );
 
  m_model->setLinkedSolver(this);
 
  // Convert the members of the ionic model into ceps::MemberFunction callable objects
  m_fLin  = ceps::newFunction(this,&IonicSolver::computeRatesLin   );
  m_fNLin = ceps::newFunction(this,&IonicSolver::computeRatesNonLin);
 
  computeIndexMapping();
  createSolvers(solverOpts);
  initializeSolver();
}
 
IonicSolver::~IonicSolver()
{
  ceps::destroyObject(m_solver);
 
  ceps::destroyObject(m_fLin );
  ceps::destroyObject(m_fNLin);
 
  ceps::destroyTabular(m_v      );
  ceps::destroyTabular(m_y      );
  ceps::destroyTabular(m_dtyLin );
  ceps::destroyTabular(m_dtyNLin);
  ceps::destroyTabular(m_dtv);
}
 
void
IonicSolver::setActionPotential(DHVecPtr vPde)
{
  vPde->getLocalData();
  for (CepsUInt i=0;i<m_size;i++)
    *(m_v+i) = m_model->convertVoltageFromCepsUnit((*vPde)[m_localToGlobal[i]]);
  vPde->releaseLocalData(); 
 
}
 
void
IonicSolver::takeOneStep(CepsReal t, CepsReal dt, DHVecPtr vPde)
{
  setActionPotential(vPde);
  m_solver->takeOneStep(t,dt);
  m_tauOk = false;
}
 
void
IonicSolver::takeOneSubStep(CepsReal t, CepsReal dt, DHVecPtr vPde)
{
  setActionPotential(vPde);
  m_solver->takeOneSubStep(t,dt);
  m_tauOk = false;
}
 
void
IonicSolver::addDtv(DHVecPtr dtv, CepsReal) const
{
 
  dtv->getLocalData();
  for (CepsUInt i=0; i<m_size; i++)
  {
    CepsUInt gId = m_localToGlobal[i];
    CepsReal val = m_model->convertDtvToCepsUnit(m_dtv[i]);
    (*dtv)[gId] += val;
  }
  dtv->releaseLocalData();
}
 
void
IonicSolver::fillPdeVm(DHVecPtr vPde) const
{
  vPde->getLocalData();
  for (CepsUInt i=0;i<m_size;i++)
    (*vPde)[m_localToGlobal[i]] = m_model->convertVoltageToCepsUnit(*(m_v+i));
  vPde->releaseLocalData(); 
}
 
CepsReal
IonicSolver::getVm(CepsBool convert) const
{
  if (convert)
    return m_model->convertVoltageToCepsUnit(*m_v);
  return *m_v;
}
 
CepsReal*
IonicSolver::getStates() const
{
  return m_y;
}
 
CepsReal*
IonicSolver::getDtv() const
{
  return m_dtv;
}
 
CepsUInt
IonicSolver::getNbMultiSteps() const
{
  return m_solver->getNbMultiSteps();
}
 
// Protected ---------------------------------------------------------------------------------------
 
void
IonicSolver::computeIndexMapping()
{
  auto dofs   = m_sDiscr->getDegreesOfFreedomForUnknown(m_model->getUnknown());
  auto mattrs = m_model->getAttributes();
 
  for (DegreeOfFreedom* dof : dofs)
  {
    if (dof->getOwner()==ceps::getRank() and 
       (mattrs.empty() or dof->hasOneOfAttributes(mattrs) or mattrs.contains(CepsUniversal)))
      m_localToGlobal.push_back(dof->getGlobalIndex());
  }
}
 
void
IonicSolver::createSolvers(CepsString opts)
{
 
  std::istringstream iss(opts);
 
  CEPS_ABORT_IF(opts.empty(),
    "Could not read options of ionic model solver"
  );
 
  CepsString type;
  CepsUInt   order;
  iss >> type;
  type = ceps::toUpper(type);
 
  if (type == "EULER" or type == "EXPLICIT_EULER")
  {
    m_solver = ceps::getNew<ExplicitEulerOdeSolver>();
  }
 
  else if (type == "IMPLICIT_EULER" or type == "FBE")
  {
    m_solver = ceps::getNew<FBEOdeSolver>();
  }
 
  else if (type == "RUSH_LARSEN" or type == "RL")
  {
    order = ceps::readInt(iss,"Ionic solver: missing order for Rush-Larsen ODE solver");
    m_solver = ceps::getNew<RushLarsenOdeSolver>(order);
  }
 
  else if (type == "EXP_ADAMS_BASHFORTH" or type == "EAB")
  {
    order = ceps::readInt(iss,"Ionic solver: missing order for Adams-Bashforth ODE solver");
    m_solver = ceps::getNew<ExponentialAdamsBashforthOdeSolver>(order);
  }
 
  else if (type == "RUNGE_KUTTA" or type == "RK")
  {
    CepsString opt;
    iss >> opt;
    CEPS_ABORT_IF(opt.empty(),
      "Ionic solver : missing type of Runge-Kutta ODE solver"
    );
    m_solver = ceps::getNew<RungeKuttaOdeSolver>(opt);
  }
 
  else if (type == "SBDF" or type == "BACKWARD_DIFFERENTIATION")
  {
    order = ceps::readInt(iss,"Ionic solver: missing order for SBDF ODE solver");
    m_solver = ceps::getNew<SBDFOdeSolver>(order);
  }
 
  else
    CEPS_ABORT("Unknown solver type. Valid types are  \n"
               << " - EXPLICIT_EULER (or just EULER)\n"
               << " - IMPLICIT_EULER (or FBE)\n"
               << " - RUSH_LARSEN <1-4> (or RL <1-4>)\n"
               << " - EXP_ADAMS_BASHFORTH <1-4> (or EAB <1-4>)\n"
               << " - BACKWARD_DIFFERENTIATION <1-4> (or SBDF <1-4>)\n"
               << " - RUNGE_KUTTA <opt> (or RK <opt>)\n"
               << "     where opt can be an integer between 1 and 5 or\n"
               << "                      Heun2 Heun3\n"
               << "                      Ralston3 SSP3 38 Ralston4\n"
               << "                      Nystrom5 Fehlberg5.");
 
}
  
void
IonicSolver::initializeSolver()
{
 
  // Allocate arrays  
  m_size = m_localToGlobal.size();
  CepsUInt nVars = m_model->getNbStateVariables();
 
  m_y       = ceps::newArray<CepsReal>(m_size*nVars);
  m_dtyLin  = ceps::newArray<CepsReal>(m_size*nVars);
  m_dtyNLin = ceps::newArray<CepsReal>(m_size*nVars);
  m_v       = ceps::newArray<CepsReal>(m_size);
  m_dtv     = ceps::newArray<CepsReal>(m_size);
 
  // Build the initial condition vector, also sets the PDE vector
  for (CepsUInt i=0;i<m_size;i++)
  {
    m_model->getInitialCondition(m_v+i,m_y+i*nVars);
    //m_dtv[i] = 0.; 
  }
  // Computes m_dtv with provided initial conditions (not necessarily steady state)
  this->computeModelRates(0.,m_y);
  m_tauOk = false;
 
  // Initialize solver with initial condition data and link to linear/non linear evaluation calls
  m_solver->initialize(m_y,m_size*nVars,m_fNLin,m_fLin,false);
 
}
 
void
IonicSolver::computeModelRates(CepsReal t, CepsReal* Y)
{
 
  CepsUInt nVars = m_model->getNbStateVariables();
 
  // Loop on all dofs
  for (CepsUInt i=0; i<m_size; i++)
  {
    
    CepsReal* y       = Y         + i*nVars;
    CepsReal* v       = m_v       + i;
    CepsReal* dtyLin  = m_dtyLin  + i*nVars;
    CepsReal* dtyNLin = m_dtyNLin + i*nVars;
    CepsReal* dtv     = m_dtv     + i;
 
    DegreeOfFreedom* dof = m_sDiscr->getDof(m_localToGlobal[i]);
 
    m_model->computeRates(t,y,v,dtyLin,dtyNLin,dtv,dof);
  }
 
  m_tauOk = true;
}
 
void
IonicSolver::computeRatesLin(CepsReal t, CepsReal* y, CepsReal* res, CepsUInt n)
{
  if (not m_tauOk)
    computeModelRates(t,y);
  std::copy(m_dtyLin,m_dtyLin+n,res);
}
 
void
IonicSolver::computeRatesNonLin(CepsReal t, CepsReal* y, CepsReal* res, CepsUInt n)
{
  if (not m_tauOk)
    computeModelRates(t,y);
  std::copy(m_dtyNLin,m_dtyNLin+n,res);
}