CardiacSolver.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 "cardiac/solver/CardiacSolver.hpp"
 
CardiacSolver::CardiacSolver(CardiacProblem *problem):
  HeatSolver                    (problem),
  m_activationTracker           (nullptr),
  m_ionSolvers                  ({}),
  m_ionSolvOpts                 ("RL 2"), 
  m_stopAtCompleteActivation    (false),
  m_stopAtCompleteRepolarization(false)
{
  m_hasRegSrc = true; // Ionic models
 
  this->m_scaleSystem = 1.;  // Reverted from 1/dt in HeatSolver
  m_activationTracker = ceps::getNew<ActivationTracker>(this);
 
  // Check parameters
  InputParameters* params = problem->getParameters();
  if (ceps::isValidPtr(params))
  {
    m_stopAtCompleteActivation = params->isActiveOption("stop at complete activation");
    CEPS_WARNS_IF(m_stopAtCompleteActivation and not params->hasKey ("activation time data"),
      "Cannot enable simulation stop at complete activation since no activation time data\n" <<
      "has been provided. Simulation will stop at time " << m_timeStepper->getEndTime()
    );
 
    m_ionSolvOpts = params->getString("ionic model solver",m_ionSolvOpts);
  }
 
  // Initialize ionic solvers
  if (ceps::isProfiling())
    getProfiler()->logMemoryUsage("\"PDE solver - before allocation of ionic data arrays\"");
 
  for (AbstractIonicModel* model: problem->getIonicModels())
    m_ionSolvers.push_back(ceps::getNew<IonicSolver>(m_ionSolvOpts,model,m_fe));
 
 
  if (ceps::isProfiling())
    getProfiler()->logMemoryUsage("\"PDE solver - end of initialization\"");
    
}
 
CardiacSolver::~CardiacSolver()
{
  ceps::destroyObject(m_activationTracker);
  for (auto solver : m_ionSolvers)
    ceps::destroyObject(solver);
}
 
void
CardiacSolver::solve()
{
  HeatSolver::solve();
 
  m_activationTracker->writeActivationMap();
}
 
ActivationTracker*
CardiacSolver::getActivationTracker() const
{
  return m_activationTracker;
}
 
 
CepsBool
CardiacSolver::finished() const
{
  CepsBool ended    = m_timeStepper->atEnd();
  CepsBool allActiv = m_stopAtCompleteActivation and m_activationTracker->allPointsHaveBeenActivated();
 
  return ended or allActiv;
}
 
CardiacProblem*
CardiacSolver::getCardiacProblem() const
{
  return ceps::runtimeCast<CardiacProblem*>(m_problem);
}
 
void 
CardiacSolver::assembleAndSolve()
{
  FETimedSolver::assembleAndSolve();
 
  if (ceps::isProfiling())
    getProfiler()->start("iion","computing ionic currents (CardiacSolver)");
  for (IonicSolver* solver: m_ionSolvers)
  {
    // if (m_type == "RK")
    // solver->takeOneSubStep(tn,dt,m_uNm[0]);
    // else
    solver->takeOneStep(m_timeStepper->getTime(),m_timeStepper->getTimeStep(),m_uNm[0]);
  }
  if (ceps::isProfiling())
    getProfiler()->stop("iion");
}
 
void
CardiacSolver::fillSourceTermsVector(CepsReal tn)
{
  auto cpb = getCardiacProblem();
 
  // Add all non ionic sources (regular and stiffness source terms)
  HeatSolver::fillSourceTermsVector(tn);
 
  // Divide by AmCm
  auto vecUs = cpb->getCardiacUnknowns();
  CepsSet<Unknown*> us(vecUs.begin(),vecUs.end());
  
  scaleByField(m_regSourcesNm[0].get(),*cpb->getCm(),{},us,tn,true);
  scaleByField(m_regSourcesNm[0].get(),*cpb->getAm(),{},us,tn,true);
 
  // Add ionic currents and stim. Ionic models compute directly (iIon+iApp)/cm
  for (IonicSolver* solver: m_ionSolvers)
  {
    solver->addDtv(m_regSourcesNm[0]);
 
    // Advance the ionic model for next step
 
    // if (ceps::isProfiling())
    //   getProfiler()->start("iion","computing ionic currents (CardiacSolver)");
    // if (m_type == "RK")
    // solver->takeOneSubStep(tn,dt,m_uNm[0]);
    // else
      // solver->takeOneStep(tn,dt,m_uNm[0]);
    // if (ceps::isProfiling())
    //   getProfiler()->stop("iion");
  }
 
}
 
void
CardiacSolver::output(DHVecPtr solution, CepsBool immediatWriting)
{
  auto pb = getCardiacProblem();
 
  HeatSolver::output(solution, false);
 
  if (m_doOutput and m_timeStepper->getNbTakenTimeSteps()%m_nbIterSnapshot==0)
  {
    auto vec = m_discretization->newDofHaloVector();
    vec->zero();
    for (IonicSolver* solver: m_ionSolvers)
    {
      solver->addDtv(vec);
    }
 
    for (Unknown* u: m_problem->getSpatialUnknowns())
      if (ceps::contains(pb->getTMVUnknowns(), u))
        m_writer->addScalarData(vec, "IIon (uA) [" + u->getName() + "]", u);
 
    if (immediatWriting)
      m_writer->write(m_timeStepper->getTime());
  }
 
  if(m_timeStepper->getNbTakenTimeSteps()>0)
    m_activationTracker->update(m_timeStepper->getNbTakenTimeSteps(),
                                m_timeStepper->getTime(),m_uNp1,m_uNm[0]);
}