FETimedSolver.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/FETimedSolver.hpp"
 
#include "pde/assemblers/AbstractAssembler.hpp"
#include "pde/problem/timed/AbstractTimedPdeProblem.hpp"
 
// #include "ode/solver/RungeKutta.hpp"
 
FETimedSolver::FETimedSolver(AbstractTimedPdeProblem* problem) :
    AbstractTimedPdeSolver(problem),
    m_fe(ceps::runtimeCast<FiniteElements*>(problem->getSpatialDiscretization())),
    m_temp1(m_fe->newDofHaloVector()),
    m_temp2(m_fe->newDofHaloVector())
{}
 
FETimedSolver::~FETimedSolver()
{
}
 
void
FETimedSolver::assembleAndSolve()
{
  CepsReal time   = m_timeStepper->getTime();
  CepsReal dt     = m_timeStepper->getTimeStep();
  CepsUInt iter   = m_timeStepper->getNbTakenTimeSteps();
  CepsUInt nSteps = std::min(iter + 1, m_nbMultiSteps);
 
  // Assemble what is needed
  updateAssemblers();
  // And fill source terms vectors
  fillSourceTermsVector(time);
 
  if (m_type=="SBDF" or m_type=="FBE")
    assembleAndSolveSBDF(nSteps);
  else if (m_type=="CN")
    assembleAndSolveCN(nSteps);
  // else if (m_type == "RK")
  //   updateAndSolveRK();
 
  // Swap pointers for current arrays
  this->swapMultiStepPointers();
 
  // Cheat for convergence tests: if there is an analytic solution,
  // we use the exact solution for first steps of multisteps methods
  if (m_cheatConvergence and m_problem->hasAnalyticSolution() and iter+1<m_nbMultiSteps)
    m_discretization->evaluateFunctionOnDofs(
      m_problem->getAnalyticSolution(),m_uNp1,time+dt);
 
  return;
}
 
void
FETimedSolver::assembleAndSolveSBDF(CepsUInt nSteps)
{
  CepsReal        time = m_timeStepper->getTime();
  CepsReal        dt   = m_timeStepper->getTimeStep();
  CepsUInt        iter = m_timeStepper->getNbTakenTimeSteps();
  const CepsReal* ak   = c_akSBDF[nSteps-1];
  CepsReal        adt  = dt*c_alphaSBDF[nSteps];
 
  CepsBool recomputeLhs = false;
  recomputeLhs          = recomputeLhs or iter+1 <= m_nbMultiSteps;
  recomputeLhs          = recomputeLhs or m_opAsb->isChangingBetweenTimes(time,time+dt);
  recomputeLhs          = recomputeLhs or m_bcAsb->isChangingBetweenTimes(time,time+dt);
 
  // Build LHS = M + adt * OP
  if (recomputeLhs)
  {
    m_lhs->finalize();
    m_dtMat->copy(*m_lhs);
    m_lhs->add(*m_opMat,adt,false);
    m_lhs->add(*m_bcMat,adt,false);
 
    if (not ceps::equals(m_scaleSystem,1.e0))
      *m_lhs *= m_scaleSystem;
  }
 
  // =============================================================================================
  // TODO
  // We must to rewrite this comment for systems of equations with time derivatives only on certain
  // lines (eg bidomain) and remove not needed F^{n-k} for static equations. Complicate ?
  // =============================================================================================
 
  // Build RHS = M*(sum_k a_k*(U^n-k + dt*(k+1)*F^n-k)), 0<k<nSteps with M = derivative
  // matrix
  m_rhs->zero();
 
  // 1. Solutions at previous steps
  m_temp1->zero();
  for (CepsUInt k=0UL; k<nSteps; k++)
    m_temp1->addScaled(*m_uNm[k], ak[k]);
 
  m_temp2->mult(*m_dtMat, *m_temp1);
  m_rhs->add(*m_temp2);
 
  // 2. Regular source terms
  if (m_hasRegSrc)
  {
    // Sources on equations with time derivative
    m_temp1->zero();
    for (CepsUInt k=0UL; k<nSteps; k++)
      m_temp1->addScaled(*m_regSourcesNm[k], ak[k]*dt*(k+1));
 
    m_temp1->pointWiseScale(*m_timeDerEqsVec);
    m_temp2->mult(*m_fe->getMassMatrix(),*m_temp1);
    m_rhs->add(*m_temp2);
 
    // Sources on "static" equations
    m_temp1->zero();
    m_temp1->addScaled(*m_regSourcesNm[0], adt);
    m_temp1->pointWiseScale(*m_noTimeDerEqsVec);
    m_temp2->mult(*m_fe->getMassMatrix(),*m_temp1);
    m_rhs->add(*m_temp2);
  }
 
  // 3. "Laplace" source terms
  if (m_hasLapSrc)
  {
    // Sources on equations with time derivative
    m_temp1->zero();
    for (CepsUInt k=0UL; k<nSteps; k++)
      m_temp1->addScaled(*m_lapSourcesNm[k], ak[k]*dt*(k+1));
 
    m_temp1->pointWiseScale(*m_timeDerEqsVec);
    m_temp2->mult(*m_lapSrcMat,*m_temp1);
    m_rhs->add(*m_temp2);
 
    // Sources on "static" equations
    m_temp1->zero();
    m_temp1->addScaled(*m_lapSourcesNm[0],adt);
    m_temp1->pointWiseScale(*m_noTimeDerEqsVec);
    m_temp2->mult(*m_lapSrcMat,*m_temp1);
    m_rhs->add(*m_temp2);
  }
 
  // 4. The neuman and robin boundary conditions
  m_rhs->addScaled(*m_bcVec,adt);
 
  // Then apply Dirichlet BCs (put 1 on diagonal in the matrix, and the corresponding value into the
  // vector)
  for (auto bc : *(m_problem->getBoundaryConditionManager()->getDirichletBCs()))
    bc.second->applyAsDirichlet(m_lhs, m_rhs, time + dt);
 
  if (not ceps::equals(m_scaleSystem, 1.e0))
    *m_rhs *= m_scaleSystem;
  
  m_rhs->finalize();
  m_lhs->finalize();
 
  // Solve
  m_linearSystem->solve(m_uNp1);
  
  #ifdef CEPS_DEBUG_ENABLED
    m_uNp1->checkNanOrInf("in solution vector ");
  #endif
 
  return;
}
 
void
FETimedSolver::assembleAndSolveCN(CepsUInt nSteps)
{
  CepsReal time = m_timeStepper->getTime();
  CepsReal dt   = m_timeStepper->getTimeStep();
  CepsUInt iter = m_timeStepper->getNbTakenTimeSteps();
  CepsReal adt  = dt * c_alphaCN[nSteps];
 
  CepsBool recomputeLhs = false;
  recomputeLhs          = recomputeLhs or iter + 1 <= m_nbMultiSteps;
  recomputeLhs          = recomputeLhs or m_opAsb->isChangingBetweenTimes(time, time + dt);
  recomputeLhs          = recomputeLhs or m_bcAsb->isChangingBetweenTimes(time, time + dt);
  
  // Build LHS = M + adt * OP
  if (recomputeLhs)
  {
    m_lhs->finalize();
    m_dtMat->copy(*m_lhs);
    m_lhs->add(*m_opMat, adt, false);
    m_lhs->add(*m_bcMat, adt, false);
 
    if (not ceps::equals(m_scaleSystem, 1.e0))
      *m_lhs *= m_scaleSystem;
  }
  
  // =============================================================================================
  // TODO
  // We must to rewrite this comment for systems of equations with time derivatives only on certain
  // lines (eg bidomain) and remove not needed F^{n-k} for static equations. Complicate ?
  // =============================================================================================
 
  // Build RHS = -0.5*dt*K*V_n + M*(V_n + 1.5*dt*F_n - 0.5* dt * F_{n-1})
  m_rhs->zero();
 
  // 1. Solutions at previous steps
  m_temp1->zero();
  m_temp1->add(*m_uNm[0]);
  m_temp2->mult(*m_dtMat, *m_temp1);
  m_rhs->add(*m_temp2);
 
  // 2. Regular source terms
  if (m_hasRegSrc)
  {
    // Sources on equations with time derivative
    m_temp1->zero();
    m_temp1->addScaled(*m_regSourcesNm[0], 1.5 * dt);
    m_temp1->addScaled(*m_regSourcesNm[1], -0.5 * dt);
 
    m_temp1->pointWiseScale(*m_timeDerEqsVec);
    m_temp2->mult(*m_fe->getMassMatrix(), *m_temp1);
    m_rhs->add(*m_temp2);
 
    // Sources on "static" equations
    m_temp1->zero();
    m_temp1->addScaled(*m_regSourcesNm[0], adt);
    m_temp1->pointWiseScale(*m_noTimeDerEqsVec);
    m_temp2->mult(*m_fe->getMassMatrix(), *m_temp1);
    m_rhs->add(*m_temp2);
  }
 
  // 3. "Laplace" source terms
  if (m_hasLapSrc)
  {
    // Sources on equations with time derivative
    m_temp1->zero();
    m_temp1->addScaled(*m_lapSourcesNm[0], 1.5 * dt);
    m_temp1->addScaled(*m_lapSourcesNm[1], -0.5 * dt);
 
    m_temp1->pointWiseScale(*m_timeDerEqsVec);
    m_temp2->mult(*m_lapSrcMat, *m_temp1);
    m_rhs->add(*m_temp2);
 
    // Sources on "static" equations
    m_temp1->zero();
    m_temp1->addScaled(*m_lapSourcesNm[0], adt);
    m_temp1->pointWiseScale(*m_noTimeDerEqsVec);
    m_temp2->mult(*m_lapSrcMat, *m_temp1);
    m_rhs->add(*m_temp2);
  }
  
  // 4. The neuman and robin boundary conditions
  m_rhs->addScaled(*m_bcVec, adt);
 
  // 5. half Op term
  m_temp1->zero();
  m_temp1->addScaled(*m_uNm[0], -0.5 * dt);
  m_temp2->mult(*m_opMat, *m_temp1);
  m_rhs->add(*m_temp2);
 
  // Then apply Dirichlet BCs (put 1 on diagonal in matrix, and the corresponding value into the
  // vector)
  for (auto bc : *(m_problem->getBoundaryConditionManager()->getDirichletBCs()))
    bc.second->applyAsDirichlet(m_lhs, m_rhs, time + dt, m_scaleSystem);
 
  if (not ceps::equals(m_scaleSystem, 1.e0))
    *m_rhs *= m_scaleSystem;
  
  m_rhs->finalize();
  m_lhs->finalize();
 
  // Solve
  m_linearSystem->solve(m_uNp1);
 
  #ifdef CEPS_DEBUG_ENABLED
    m_uNp1->checkNanOrInf("in solution vector ");
  #endif
 
  return;
}
 
/*void
FETimedSolver::setupLinearSystemRK()
{
  // LHS matrix is the mass matrix + boundary conditions in this case
 
  // Copy mass matrix from finite elements
  ceps::destroyObject(m_lhs.get());
  m_lhs = DMatPtr(ceps::getNew<DistributedMatrix>(*(m_fe->getMassMatrix()),true));
  m_linearSystem->setLhsMatrix(m_lhs);
 
  // Neumann and Robin boundary conditions
  FEDivKGradBCAssembler basb(m_fe);
  basb.setMatrix(m_lhs)
  basb.assemble();
 
  // Then apply Dirichlet BCs, for the matrix only
  m_lhs->finalize();
  for (auto bc : m_problem->getBoundaryConditionsManager()->getDirichletBCs())
    bc->applyAsDirichlet(m_lhs,nullptr);
 
  m_lhs->finalize();
  return;
}*/
 
/* void
// FETimedSolver::buildRHSVectorRK ()
// {
//   // RHS = -SV^n_k + M(V^n+dt*sum_j b_jkF(time^n_k))
//   m_rhs->mult (*m_lhs, *m_regSourcesNm[0]);
//   m_temp->zero ();
//   m_temp->mult (*m_stiffMat, *m_VnM[0]);
//   m_rhs->addScaled (*m_temp, -1.);
// }*/