AbstractStaticPdeSolver.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/static/AbstractStaticPdeSolver.hpp"
#include "pde/problem/static/AbstractStaticPdeProblem.hpp"
#include "pde/assemblers/AbstractAssembler.hpp"
 
AbstractStaticPdeSolver::AbstractStaticPdeSolver(AbstractStaticPdeProblem* pb):
  AbstractPdeSolver(pb),
  m_writeInitialGuess(false)
{
  if (m_doOutput)
    this->initializeWriter();
  m_timeWriter = ceps::getNew<TimeWriter>(m_problem->getOutputFileBase()+"_probes.dat",
                                m_problem->getProbePoints(),
                                m_discretization);
 
  if (pb->getParameters())
    AbstractStaticPdeSolver::setupWithParameters(pb->getParameters());
 
  m_doError = pb->hasAnalyticSolution();
  if (m_doError)
    enableErrorComputation();
 
 
  m_solution = m_problem->getSpatialDiscretization()->newDofHaloVector();
}
 
AbstractStaticPdeSolver::~AbstractStaticPdeSolver()
{
}
 
void 
AbstractStaticPdeSolver::solve()
{
  CEPS_SEPARATOR;
  CEPS_SAYS("Now solving PDE " + m_problem->getProblemName());
  auto lpb = getStaticProblem();
  lpb->getInitialGuess(m_solution);
 
  if (m_writeInitialGuess)
    output(m_solution);
 
  initializeAssemblers();
 
  assembleAndSolve();
 
  if (m_doError)
    m_errors->compute(m_solution);
 
  CEPS_SAYS_NOENDL("  Writing result...");
  output(m_solution);
  CEPS_SAYS_INLINE(" Done");
}
 
void
AbstractStaticPdeSolver::setupWithParameters(InputParameters* params)
{
  m_writeInitialGuess = params->isActiveOption("output initial guess");
}
 
CepsUInt 
AbstractStaticPdeSolver::getExpectedNumberOfOutputs() const
{
  return m_writeInitialGuess ? 2 : 1;
}
 
void
AbstractStaticPdeSolver::output(DHVecPtr solution, CepsBool immediateWriting)
{
  if (m_doOutput)
  {
    for (Unknown* u: m_problem->getUnknowns())
      m_writer->addScalarData(solution,u->getName(),u);
 
    if (m_problem->hasAnalyticSolution())
    {
      auto uNp1 = m_discretization->newDofHaloVector();
      m_discretization->evaluateFunctionOnDofs(
        m_problem->getAnalyticSolution(), uNp1, 0.0);
 
      for (Unknown* u: m_problem->getUnknowns())
        m_writer->addScalarData(uNp1, "REF " + u->getName(), u);
 
    } 
 
    for (auto [key, f] : *m_problem->getSourceTermManager()->getManager())
    {
      auto vec = m_discretization->newDofHaloVector();
      addFieldValuesTo(*f, vec.get(), 0);
 
      for (Unknown* u: m_problem->getSpatialUnknowns())
        m_writer->addScalarData(vec, key + " [" + u->getName() + "]", u);
    }
 
    if (immediateWriting)
      m_writer->write(0.);
    // m_writer->writeSolution(m_problem->getSpatialUnknowns(),solution,0);
  }
  
  if (m_timeWriter->hasSomethingToWrite())
    m_timeWriter->write(0.,solution);
}
 
DHVecPtr
AbstractStaticPdeSolver::getSolution() const
{
  return DHVecPtr(ceps::getNew<DistributedHaloVector>(*(m_solution)));
}
 
void
AbstractStaticPdeSolver::enableErrorComputation()
{
  m_doError = true;
  ceps::destroyObject(m_errors);
  m_errors = ceps::getNew<PdeErrorCalculator>(m_problem);
}
 
CepsArray2<CepsArray3<CepsReal>>
AbstractStaticPdeSolver::getErrors() const
{
  CepsArray2<CepsArray3<CepsReal>> res;
  if (not m_doError)
    return res;
 
  for (CepsUInt norm=0; norm<3; norm++)
  {
    res[0][norm] = m_errors->getAbsoluteErrorNow(norm);
    res[1][norm] = m_errors->getRelativeErrorNow(norm);
  }
  return res;
}
 
// -----------------------------------------------------------------------------------------
 
AbstractStaticPdeProblem*
AbstractStaticPdeSolver::getStaticProblem() const
{
  return ceps::runtimeCast<AbstractStaticPdeProblem*>(m_problem);
}
 
void
AbstractStaticPdeSolver::updateAssemblers()
{
 
  // Assemble operator matrix terms
  m_opMat->zero();
  m_opAsb->setMatrix(m_opMat.get());
  m_opAsb->setScaleFactor(1.0);
  m_opAsb->assemble();
 
  // Assemble only the matrix terms of the Neumann and Robin BCs contributions
  m_bcMat->zero();
  m_bcVec->zero();
  m_bcAsb->setMatrix(m_bcMat.get());
  m_bcAsb->setVector(m_bcVec.get());
  m_bcAsb->setScaleFactor(1.0);
  m_bcAsb->assemble();
 
  // Laplace source terms
  if (m_hasLapSrc)
  {
    m_lapSrcMat->zero();
    m_lapSrcAsb->setMatrix(m_lapSrcMat.get());
    m_lapSrcAsb->assemble();
  }
 
  return;
}