AbstractTimedPdeSolver.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 "pde/solver/timed/AbstractTimedPdeSolver.hpp"
#include "pde/problem/timed/AbstractTimedPdeProblem.hpp"
#include "pde/assemblers/AbstractAssembler.hpp"
 
AbstractTimedPdeSolver::AbstractTimedPdeSolver(AbstractTimedPdeProblem* pb):
  AbstractPdeSolver (pb),
  m_uNm             ({}),
  m_type            ("FBE"),
  m_rkType          (""),
  m_nbMultiSteps    (1u),
  m_timeStepper     (nullptr),
  m_regSourcesNm    ({}),
  m_lapSourcesNm    ({}),
  m_rkF             ({}),
  m_rkU             ({}),
  // m_scaleSystem     (1./pb->getTimeStep()),
  m_scaleSystem     (1.),
  m_cheatConvergence(true),
  m_nbIterSnapshot  (1)
{
  m_timeStepper = ceps::getNew<TimeStepper>(pb->getStartTime(),pb->getEndTime(),pb->getTimeStep());
 
  m_nbIterSnapshot = std::max(CepsUInt(1U),CepsUInt(pb->getSnapshotTime()/pb->getTimeStep()));
 
  if (m_doOutput)
    this->initializeWriter();
  m_timeWriter = ceps::getNew<TimeWriter>(m_problem->getOutputFileBase()+"_probes.dat",
                                m_problem->getProbePoints(),
                                m_discretization);
 
 
  if (pb->getParameters())
    AbstractTimedPdeSolver::setupWithParameters(pb->getParameters());
 
  m_doError = pb->hasAnalyticSolution();
  if (m_doError)
    enableErrorComputation();
 
  determineSolverType(m_type);
  AbstractTimedPdeSolver::allocateArrays();
 
  // Matrix of derivative terms 
  m_dtMat = m_discretization->newDofMatrix();
  // Vectors of footprint of diagonal of dtMat
  m_timeDerEqsVec = m_discretization->newDofHaloVector();
  m_noTimeDerEqsVec = m_discretization->newDofHaloVector();
}
 
AbstractTimedPdeSolver::~AbstractTimedPdeSolver()
{
  ceps::destroyObject(m_timeStepper);
}
 
void
AbstractTimedPdeSolver::solve()
{
  CEPS_SAYS("Initialize assemblers for " << std::quoted(m_problem->getProblemName()));
  this->initializeAssemblers();
  getTimedProblem()->getInitialCondition(m_uNm[0]);
 
  // Write initial step
  this->output(m_uNm[0]);
 
  CEPS_SAYS(
    "Now computing PDE solution ("+m_type+(m_nbMultiSteps>1?ceps::toString(m_nbMultiSteps):"")+"):"
  );
  if (ceps::isVerbose())
    displayProgress();
 
  while (not this->finished())
  {
 
    CepsUInt iter = m_timeStepper->getNbTakenTimeSteps();
 
    if (ceps::isProfiling())
      getProfiler()->start("updateLS", "updating system and solving (" + m_problem->getProblemName() + ")");
 
    this->assembleAndSolve();
 
    if (ceps::isProfiling())
    {
      getProfiler()->stop("updateLS");
      getProfiler()->logMemoryUsage("\"PDE solver - after "+ceps::toString(iter)+"th update\"");
    }
 
    // Update time and write solutions
    m_timeStepper->takeOneStep();
 
    if (ceps::isVerbose())
      displayProgress();
    this->output(m_uNp1);
 
    if (getTimedProblem()->canComputeErrorAtTime(m_timeStepper->getTime()))
      m_errors->compute(m_uNp1,m_timeStepper->getTime());
 
    // Swap pointers on solutions at previous steps
    DHVecPtr tmp = m_uNm.back();
    for (CepsUInt k = m_nbMultiSteps-1; k>0;k--)
      m_uNm[k] = m_uNm[k-1];
    m_uNm[0] = m_uNp1;
    m_uNp1   = tmp;
 
  }
 
  // Do a final output, even if it could be a duplicate of previous one
  // (the finished method can be overloaded with something else)
  this->output(m_uNm[0]);
 
  if (ceps::isVerbose())
    std::cout << std::endl;
}
 
void
AbstractTimedPdeSolver::setupWithParameters(InputParameters* p)
{
  m_type = p->getString("pde solver",m_type);
}
 
CepsUInt 
AbstractTimedPdeSolver::getExpectedNumberOfOutputs() const
{
  auto tpb = getTimedProblem();
 
  // Snapshots + initial guess and final output (in case it does not match with snapshot)
  return 2 + CepsUInt(std::floor((tpb->getEndTime()-tpb->getStartTime())/
                                  tpb->getSnapshotTime()));
}
 
void
AbstractTimedPdeSolver::output(DHVecPtr solution, CepsBool immediateWriting)
{
  if (m_doOutput and m_timeStepper->getNbTakenTimeSteps()%m_nbIterSnapshot==0)
  {
    // Unknowns
    for (Unknown* u: m_problem->getUnknowns())
      m_writer->addScalarData(solution,u->getName(),u);
 
    // Reference solution if exists
    if (m_problem->hasAnalyticSolution())
    {
      auto uNp1 = m_discretization->newDofHaloVector();
      m_discretization->evaluateFunctionOnDofs(
        m_problem->getAnalyticSolution(), uNp1, m_timeStepper->getTime());
 
      for (Unknown* u: m_problem->getUnknowns())
        m_writer->addScalarData(uNp1, "Reference " + u->getName(), u);
    } 
 
    // Sources
    for (auto [key, f] : *m_problem->getSourceTermManager()->getManager())
    {
      if (not key.starts_with("__"))
      {
        auto vec = m_discretization->newDofHaloVector();
        addFieldValuesTo(*f, vec.get(), m_timeStepper->getTime());
        for (Unknown* u: m_problem->getSpatialUnknowns())
          m_writer->addScalarData(vec, key + " [" + u->getName() + "]", u);
      }
    }
 
    if (immediateWriting)
      m_writer->write(m_timeStepper->getTime());
    // m_writer->writeSolution(
    //   m_problem->getSpatialUnknowns(),solution,
    //   m_timeStepper->getNbTakenTimeSteps()
    // );
  }
 
  if (m_timeWriter->hasSomethingToWrite())
    m_timeWriter->write(m_timeStepper->getTime(),solution);
}
 
DHVecPtr
AbstractTimedPdeSolver::getSolution() const
{
  return DHVecPtr(ceps::getNew<DistributedHaloVector>(*m_uNp1));
}
 
void
AbstractTimedPdeSolver::enableErrorComputation()
{
  m_doError = true;
  ceps::destroyObject(m_errors);
  CepsReal dt = m_problem->getReferenceSolutionOutputPeriod();
  if (not m_problem->usesReferenceSolution()) 
    dt = getTimedProblem()->getTimeStep();
  m_errors = ceps::getNew<PdeErrorCalculator>(m_problem,dt);
}
 
CepsArray2<CepsArray3<CepsArray3<CepsReal>>>
AbstractTimedPdeSolver::getErrors() const
{
  CepsArray2<CepsArray3<CepsArray3<CepsReal>>> res;
  if (not m_doError)
    return res;
 
  for (CepsUInt normt=0; normt<2; normt++)
  for (CepsUInt norm =0; norm <3; norm ++)
  {
    res[0][normt][norm] = m_errors->getAbsoluteErrorCumulative(norm,normt);
    res[1][normt][norm] = m_errors->getRelativeErrorCumulative(norm,normt);
  }
  for (CepsUInt norm =0; norm <3; norm ++)
  {
    res[0][2][norm] = m_errors->getAbsoluteErrorNow(norm);
    res[1][2][norm] = m_errors->getRelativeErrorNow(norm);
  }
 
  return res;
}
 
// -----------------------------------------------------------------------------------------
 
AbstractTimedPdeProblem*
AbstractTimedPdeSolver::getTimedProblem() const
{
  return ceps::runtimeCast<AbstractTimedPdeProblem*>(m_problem);
}
 
void
AbstractTimedPdeSolver::displayProgress() const
{
  if (ceps::isMaster())
  {
    CepsReal t      = m_timeStepper->getTime     ();
    CepsReal tStart = m_timeStepper->getStartTime();
    CepsReal tEnd   = m_timeStepper->getEndTime  ();
    CepsUInt done   = (CepsUInt) std::floor((t-tStart)/(tEnd-tStart)*100);
 
    std::cout << "\r     Simulation time : " << std::setw(10) << t << "/" << std::setw(10) << tEnd << " ms ("
              << std::setw(3) << done << "%)" << std::flush;
  }
}
 
CepsBool
AbstractTimedPdeSolver::finished() const
{
  return m_timeStepper->atEnd();
}
 
void
AbstractTimedPdeSolver::determineSolverType(CepsString opt)
{
  std::stringstream ss(opt);
  ss >> m_type;
 
  if (m_type == "FBE")
  {
    m_nbMultiSteps = 1;
    m_type         = "SBDF";
  }
  else if (m_type == "SBDF")
  {
    m_nbMultiSteps = ceps::readInt(ss,"Unable to read the order of the SBDF method");
    CEPS_ABORT_IF(m_nbMultiSteps > 4,"SBDF solvers of order larger than 4 are not implemented.");
  }
  // else if (m_type == "RK" or m_type == "RUNGE_KUTTA")
  // {
  //   // The solver compatibility with ionic solvers was already checked in CardiacParameters class
  //   ss >> m_rkType;
  //   m_type         = "RK";
  //   m_nbMultiSteps = 1;  // not be confused with RK-substeps
  // }
  else if (m_type == "CN" or m_type=="Crank-Nicolson" or m_type=="Crank-Nicholson")
  {
    m_type         = "CN";
    m_nbMultiSteps = 2;
  }
  else
  {
    CEPS_ABORT ("Timed Pde Solver: Unknown solver type: "
                << m_type << " .Valid types are\n"
                << "   - FBE (for single step Euler method)\n"
                << "   - SBDF <order> where order is an integer between 1 and 4\n"
                // << "   - RK <options> where options are the same as for the line\n"
                << "   - CN (or Crank-Nicolson)");
  }
  return;
 
}
 
void
AbstractTimedPdeSolver::allocateArrays()
{
 
  m_uNp1 = m_discretization->newDofHaloVector();
  for (CepsUInt i=0;i<m_nbMultiSteps;i++)
    m_uNm.push_back(m_discretization->newDofHaloVector());
 
  for (CepsUInt i=0;i<m_nbMultiSteps;i++)
    m_regSourcesNm.push_back(m_discretization->newDofHaloVector());
 
  for (CepsUInt i=0;i<m_nbMultiSteps;i++)
    m_lapSourcesNm.push_back(m_discretization->newDofHaloVector());
 
      
  // Allocate arrays for the RK method
  // if (m_type=="RK")
  // {
  //   RungeKuttaOdeSolver rks (m_RKType);
  //   for (CepsUInt k = 0U; k < rks.nSubSteps (); k++)
  //   {
  //     m_rkF.push_back(m_discretization->newDofHaloVector());
  //     m_rkU.push_back(m_discretization->newDofHaloVector());
  //   }
  // }
  
}
 
 
void
AbstractTimedPdeSolver::updateAssemblers()
{
  CepsUInt iter = m_timeStepper->getNbTakenTimeSteps();
  CepsReal time = m_timeStepper->getTime();
  CepsReal dt   = m_timeStepper->getTimeStep();
 
  // Assemble derivative matrix terms
  if (iter == 0 or m_opAsb->isChangingBetweenTimes(time, time + dt))
  {
    m_dtMat = m_discretization->newDofMatrix();
    m_opAsb->setMatrix(m_dtMat.get());
    m_opAsb->setAssemblingFlag(CepsAssemblingFlag::TimeDerivative);
    m_opAsb->setScaleFactor(1.0);
    m_opAsb->assemble(time + dt);
 
    // I couldn't think of an easier way
    m_dtMat->getDiagonalFootPrintAsDistributedVector(*m_timeDerEqsVec);
    m_noTimeDerEqsVec->zero();
    m_noTimeDerEqsVec->addScaled(*m_timeDerEqsVec, -1.0);
    *m_noTimeDerEqsVec += 1.0;
  }
 
  // Assemble operator matrix terms
  if (iter == 0 or m_opAsb->isChangingBetweenTimes(time, time + dt))
  {
    m_opMat = m_discretization->newDofMatrix();
    m_opAsb->setMatrix(m_opMat.get());
    m_opAsb->setAssemblingFlag(CepsAssemblingFlag::SpatialOperator);
    m_opAsb->setScaleFactor(1.0);
    m_opAsb->assemble(time + dt);
  }
 
  // Assemble only the matrix terms of the Neumann and Robin BCs contributions
  if (iter == 0 or m_bcAsb->isChangingBetweenTimes(time, time+dt))
  {
    m_bcMat = m_discretization->newDofMatrix();
    m_bcAsb->setMatrix(m_bcMat.get());
  }
  else {
    m_bcAsb->setMatrix(nullptr);
  }
 
  // if (ceps::isNullPtr(m_bcVec.get()))
  m_bcVec = m_discretization->newDofHaloVector();
  m_bcAsb->setVector(m_bcVec.get());
  m_bcAsb->setScaleFactor(1.0);
  m_bcAsb->assemble(time + dt);
 
  // Laplace source terms
  if (m_hasLapSrc)
    if (iter == 0 or m_bcAsb->isChangingBetweenTimes(time, time + dt))
    {
      m_lapSrcMat = m_discretization->newDofMatrix();
      m_lapSrcAsb->setMatrix(m_lapSrcMat.get());
      m_lapSrcAsb->assemble(time);
    }
 
  return;
}
 
void
AbstractTimedPdeSolver::fillSourceTermsVector(CepsReal tn)
{
  m_problem->getSourceTermManager()->actualizeAll(tn);
 
  // Regular source terms
  if (m_hasRegSrc)
  {
    m_regSourcesNm[0]->zero();
    // Update source terms and add them to the sources vectors at t^n
    for (ScalarSourceTerm *src : m_problem->getSourceTermManager()->asVector())
      if (not src->hasOption(CepsSourceTermFlag::Laplace)
      and not src->hasOption(CepsSourceTermFlag::Stimulation))
        addFieldValuesTo(*src,m_regSourcesNm[0].get(),tn);
    m_regSourcesNm[0]->finalize();
 
    CEPS_CHECK_IS_NAN(m_regSourcesNm[0]->dot(*m_regSourcesNm[0]), "source term contains NaN values");
  }
 
  // Laplace source terms
  if (m_hasLapSrc)
  {
    m_lapSourcesNm[0]->zero();
    // Update source terms and add them to the sources vectors at t^n
    for (ScalarSourceTerm *src : m_problem->getSourceTermManager()->asVector())
      if (src->hasOption(CepsSourceTermFlag::Laplace))
        addFieldValuesTo(*src,m_lapSourcesNm[0].get(),tn);
    m_lapSourcesNm[0]->finalize();
 
    CEPS_CHECK_IS_NAN(m_lapSourcesNm[0]->dot(*m_lapSourcesNm[0]), "laplace source term contains NaN values");
  }
}
 
void
AbstractTimedPdeSolver::swapMultiStepPointers()
{
  if (m_hasRegSrc)
    std::rotate(m_regSourcesNm.rbegin(),m_regSourcesNm.rbegin()+1,m_regSourcesNm.rend());
 
  if (m_hasLapSrc)
    std::rotate(m_lapSourcesNm.rbegin(),m_lapSourcesNm.rbegin()+1,m_lapSourcesNm.rend());
}