PdeErrorCalculator.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/common/PdeErrorCalculator.hpp"
#include "pde/problem/AbstractPdeProblem.hpp"
 
 
PdeErrorCalculator::PdeErrorCalculator(AbstractPdeProblem* pb, CepsReal dt) :
    m_problem    (pb),
    m_solution   (nullptr),
    m_dt         (dt),
    m_refNow     ({0,0,0}),
    m_errNow     ({0,0,0}),
    m_mxTNormsRef({0,0,0}),
    m_mxTNormsErr({0,0,0}),
    m_l1TNormsRef({nullptr,nullptr,nullptr}),
    m_l1TNormsErr({nullptr,nullptr,nullptr}),
    m_ignoreZeroD(false)
{
  CEPS_ABORT_IF(ceps::isNullPtr(pb),
    "Passed ptr to problem is null"
  );
  setAnalyticSolution(pb->getAnalyticSolution());
  m_ignoreZeroD = pb->ignoreZeroDUnknownsForError();
 
  reset();
}
 
PdeErrorCalculator::~PdeErrorCalculator()
{
  for (auto& itg: m_l1TNormsErr)
    ceps::destroyObject(itg);
  for (auto& itg: m_l1TNormsRef)
    ceps::destroyObject(itg);
}
 
void
PdeErrorCalculator::compute(DHVecPtr num, CepsReal time)
{
 
  // Build a vector with the values of the analytic solution
  auto     sd  = m_problem->getSpatialDiscretization();
  DHVecPtr ana = sd->newDofHaloVector();
  sd->evaluateFunctionOnDofs(m_solution,ana,time);
 
  m_refNow[0] = ana->lInfNorm();
  m_refNow[1] = sd->l1Norm(ana);
  m_refNow[2] = sd->l2Norm(ana);
 
  for (CepsUInt i=0;i<3;i++)
  {
    m_mxTNormsRef[i] = std::max(m_mxTNormsRef[i],m_refNow[i]);
    m_l1TNormsRef[i]->update();
  }
 
  ana->addScaled(*num,-1.0);
  if (m_ignoreZeroD)
  {
    ana->getLocalData();
    for (Unknown* u : m_problem->getZeroDUnknowns())
      for (DegreeOfFreedom* dof: sd->getDegreesOfFreedomForUnknown(u))
        if (dof->getOwner()==ceps::getRank())
          (*ana)[dof->getGlobalIndex()] = 0.;
    ana->releaseLocalData();
  }
 
  m_errNow[0] = ana->lInfNorm();
  m_errNow[1] = sd->l1Norm(ana);
  m_errNow[2] = sd->l2Norm(ana);
 
  for (CepsUInt i=0;i<3;i++)
  {
    m_mxTNormsErr[i] = std::max(m_mxTNormsErr[i],m_errNow[i]);
    m_l1TNormsErr[i]->update();
  }
 
  return;
}
 
CepsReal
PdeErrorCalculator::getReferenceNormNow(CepsInt px) const
{
  checkPxPt(px,0);
  return m_refNow[px];
}
 
CepsReal
PdeErrorCalculator::getReferenceNormCumulative(CepsInt px, CepsInt pt) const
{
  checkPxPt(px,pt);
  if (pt==1)
    return m_l1TNormsRef[px]->get();
  return m_mxTNormsRef[px];
}
 
CepsReal
PdeErrorCalculator::getAbsoluteErrorNow(CepsInt px)const
{
  checkPxPt(px,0);
  return m_errNow[px];
}
 
CepsReal
PdeErrorCalculator::getAbsoluteErrorCumulative(CepsInt px, CepsInt pt)const
{
  checkPxPt(px,pt);
  if (pt==1)
    return m_l1TNormsErr[px]->get();
  return m_mxTNormsErr[px];
}
 
CepsReal
PdeErrorCalculator::getRelativeErrorNow(CepsInt px) const
{
  checkPxPt(px,0);
  CepsReal ref = std::max(1E-08,m_refNow[px]);
  return m_errNow[px]/ref;
}
 
CepsReal
PdeErrorCalculator::getRelativeErrorCumulative(CepsInt px, CepsInt pt) const
{
  checkPxPt(px,pt);
  if (pt==1)
    return m_l1TNormsErr[px]->get()/std::max(1E-08,m_l1TNormsRef[px]->get());
  return m_mxTNormsErr[px]/std::max(1E-08,m_mxTNormsRef[px]);
}
 
void
PdeErrorCalculator::setAnalyticSolution(ceps::Function<CepsReal(CepsStandardArgs)>* func)
{
  CEPS_ABORT_IF(ceps::isNullPtr(func),
    "Passed ptr on analytic solution or reference solution is null"
  );
  m_solution = ceps::runtimeCast<ScalarSAFunc*>(func);
}
 
void
PdeErrorCalculator::reset()
{
  m_refNow = m_errNow = m_mxTNormsErr = m_mxTNormsRef = {0.,0.,0.};
  for (CepsUInt i=0;i<3;i++)
  {
    ceps::destroyObject(m_l1TNormsErr[i]);
    ceps::destroyObject(m_l1TNormsRef[i]);
    m_l1TNormsErr[i] = ceps::getNew<TimeIntegrator>(&m_errNow[i],m_dt,4);
    m_l1TNormsRef[i] = ceps::getNew<TimeIntegrator>(&m_refNow[i],m_dt,4);
  }
}
 
void
PdeErrorCalculator::checkPxPt(CepsInt px, CepsInt pt) const
{
  CEPS_ABORT_IF(px>2 or px<0 or pt>1 or pt<0,
    "When requesting error in (l-px,l-pt) norm, px and pt must be either" << std::endl <<
    "     0 for l-infinity norm" << std::endl <<
    "     1 for l-1        norm" << std::endl <<
    "     2 for l-2        norm, but not 2 for pt"
  );
}