AbstractPdeProblem.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/problem/AbstractPdeProblem.hpp"
#include "pde/finiteElements/FiniteElements.hpp"
#include "pde/common/functions/FileInterpolatorSAFunc.hpp"
 
AbstractPdeProblem::AbstractPdeProblem(
  Geometry*        geom, 
  InputParameters* params
):
  m_geom              (geom),
  m_parameters        (params),
  m_discr             (nullptr),
  m_ownedDiscr        (false),
  m_analyticSolution  (nullptr),
  m_ownedRefSol       (false),
  m_refSolFiles       (""),
  m_refSolSnapDt      (1.),
  m_functions         (ceps::getNew<FunctionDictionary>()),
  m_boundaryConditions(ceps::getNew<BoundaryConditionManager>(m_functions)),
  m_sourceTerms       (ceps::getNew<SourceTermManager>(m_functions)),
  m_ownedFunctions    (true),
  m_ownedBCs          (true),
  m_ownedSrcs         (true),
  m_outputFileBase    ("./result"),
  m_outputFormat      (CepsOutputFormat::ParaviewSeries),
  m_writeGlobalIDs    (false),
  m_probePoints       ({}),
  m_ignoreZeroDError  (false)
{
  CEPS_ABORT_IF(ceps::isNullPtr(geom),
    "When instantiating pde problem : passed ptr to Geometry is null"
  );
  if (m_parameters)
    AbstractPdeProblem::setupWithParameters(m_parameters);
}
 
AbstractPdeProblem::~AbstractPdeProblem()
{
  for (auto u: m_unknowns)
    ceps::destroyObject(u);
  for (auto ui: m_unknownsInteractions)
    ceps::destroyObject(ui);
  if (m_ownedSrcs)
    ceps::destroyObject(m_sourceTerms);
  if (m_ownedBCs)
    ceps::destroyObject(m_boundaryConditions);
  if (m_ownedFunctions)
    ceps::destroyObject(m_functions);
  if (m_ownedDiscr)
    ceps::destroyObject(m_discr);
  if (m_ownedRefSol)
    ceps::destroyObject(m_analyticSolution);
}
 
InputParameters*
AbstractPdeProblem::getParameters() const
{
  return m_parameters;
}
 
void
AbstractPdeProblem::setupWithParameters(InputParameters* params)
{
  setOutputFileBase (params->getString("output file name",m_outputFileBase));
  if (params->hasKey("output format"))
  {
    CepsString otype = ceps::toUpper(params->getString("output format"));
    if (otype == "MUSICARDIO" or otype=="VTK_LEGACY")
      m_outputFormat = CepsOutputFormat::Music;
    else if (otype == "PARAVIEW_SERIES" or otype == "PSERIES") // default
      m_outputFormat = CepsOutputFormat::ParaviewSeries;
    else if (otype == "CEPS")
      m_outputFormat = CepsOutputFormat::VTKCeps;
    //else if (otype == "MEDIT")
    //  m_outputFormat = CepsOutputFormat::Medit;
    else CEPS_ABORT(
      "Unknown output format : " << otype << ". Valid types are\n" <<
      "     MUSICARDIO, VTK_LEGACY, PARAVIEW_SERIES, PSERIES (default), CEPS"
    );
  }
 
  writeGlobalIndices(params->isActiveOption("output global indices"));
  if (params->hasKey("probe points"))
  {
    auto spoints = ceps::split(params->getString("probe points"),",");
    for (CepsString spoint: spoints)
    {
      std::stringstream ss(spoint);
      CepsReal3D v = ceps::readVertex(ss,"Cannot read probe point with input string " 
                                      + params->getString("probe points"));
      m_probePoints.push_back(v);
    }
  }
 
  m_refSolFiles  = params->getString("reference solution files"        ,m_refSolFiles );
  m_refSolSnapDt = params->getReal  ("reference solution output period",m_refSolSnapDt);
 
  m_ignoreZeroDError = params->isActiveOption("convergence ignore 0D unknowns");
 
}
 
void
AbstractPdeProblem::run()
{
}
 
void 
AbstractPdeProblem::setProblemName(const CepsString& name)
{
  m_name = name;
}
 
CepsString 
AbstractPdeProblem::getProblemName() const
{
  return m_name;
}
 
Geometry*
AbstractPdeProblem::getGeometry() const
{
  return m_geom;
}
 
AbstractDiscretization*
AbstractPdeProblem::getSpatialDiscretization() const
{
  return m_discr;
}
 
void
AbstractPdeProblem::setSpatialDiscretization(AbstractDiscretization* d)
{
  if (m_ownedDiscr)
    ceps::destroyObject(m_discr);
  m_discr      = d;
  m_ownedDiscr = false;
}
 
void
AbstractPdeProblem::createSpatialDiscretization()
{
  CEPS_ABORT_IF(m_unknowns.empty(), "no unknowns in this problem");
 
  CepsString dtype = "FE 1"; // default value
  if (m_parameters)
    dtype = m_parameters->getString("spatial discretization",dtype);
 
  CepsVector<CepsString> words = ceps::split(dtype);
  CEPS_ABORT_IF(words.size()<2,
    "When creating spatial discretization, not enough parameters to deduce the type of discretization"
  );
 
  if (ceps::toKey(words[0])=="FE")
    m_discr = ceps::getNew<FiniteElements>(this,ceps::toUInt(words[1]));
  else
    CEPS_ABORT("No other discretization than finite elements implemented yet");
 
  // Link dofs to source terms and boundary conditions
  m_sourceTerms       ->setDofsInfos(m_discr->getDistributedDofs());
  m_boundaryConditions->setDofsInfos(m_discr->getDistributedDofs());
 
  m_ownedDiscr = true;
}
 
// -----------------------------------------------------------------------------------------
// UNKNOWNS
    
const CepsVector<Unknown*>&
AbstractPdeProblem::getUnknowns() const
{
  return m_unknowns;
}
 
Unknown*
AbstractPdeProblem::getUnknown(const CepsString& label) const
{
  Unknown* res = nullptr;
  for (Unknown* u:m_unknowns)
    if (ceps::toKey(u->getName())==ceps::toKey(label))
      res = u;
  return res;
}
 
Unknown*
AbstractPdeProblem::getUnknown(CepsUnknownIndex uid) const
{
  Unknown* res = nullptr;
  for (Unknown* u:m_unknowns)
    if (u->getIdentifier()==uid)
      res = u;
  return res;
}
 
const CepsVector<UnknownInteraction*>&
AbstractPdeProblem::getUnknownsInteractions() const
{
  return m_unknownsInteractions;
}
 
CepsVector<Unknown*>
AbstractPdeProblem::getSpatialUnknowns() const
{
  CepsVector<Unknown*> res;
  for (Unknown* u:m_unknowns)
    if (u->isSpatial())
      res.push_back(u);
  return res;
}
 
CepsVector<Unknown*>
AbstractPdeProblem::getZeroDUnknowns() const
{
  CepsVector<Unknown*> res;
 
  for (Unknown* u:m_unknowns)
    if (u->getLocation()==CepsLocationFlag::ZeroD)
      res.push_back(u);
 
  return res;
}
 
CepsBool
AbstractPdeProblem::unknownsInteract(
  Unknown* u1,
  Unknown* u2,
  const CepsSet<CepsAttribute>& attrs
) const
{
  if (u1==u2)
    return true;
  for (UnknownInteraction* ui : m_unknownsInteractions)
  {
    if ((ui->u1==u1 and ui->u2==u2) or (ui->u2==u1 and ui->u1==u2))
      if (attrs.contains(CepsUniversal) or ui->hasOneOfAttributes(attrs)
        or ui->hasAttribute(CepsUniversal) or ui->getNumberOfAttributes()==0)
        return true;
  }
 
  return false;
}
 
CepsBool
AbstractPdeProblem::ignoreZeroDUnknownsForError() const
{
  return m_ignoreZeroDError;
}
 
 
// -----------------------------------------------------------------------------------------
 
CepsBool
AbstractPdeProblem::hasAnalyticSolution() const
{
  return ceps::isValidPtr(m_analyticSolution);
}
 
AbstractPdeProblem::ScalarFunction*
AbstractPdeProblem::getAnalyticSolution() const
{
  return m_analyticSolution;
}
 
void
AbstractPdeProblem::defineAnalyticSolution()
{
  m_analyticSolution = nullptr;
 
  if (not m_refSolFiles.empty())
  {
    setReferenceSolution(m_refSolFiles,m_refSolSnapDt); 
  }
}
 
void
AbstractPdeProblem::setReferenceSolution(const CepsString& baseFiles, CepsReal snapDt)
{
  CEPS_SAYS("Loading reference solution from " + m_refSolFiles + "*");
  m_analyticSolution = ceps::getNew<FileInterpolatorSAFunc>(this,snapDt,m_refSolFiles,3);
  m_ownedRefSol      = true;
}
 
CepsBool
AbstractPdeProblem::usesReferenceSolution() const
{
  return hasAnalyticSolution() and not m_refSolFiles.empty();
}
 
CepsReal
AbstractPdeProblem::getReferenceSolutionOutputPeriod() const
{
  return m_refSolSnapDt;
}
 
FunctionDictionary*
AbstractPdeProblem::getFunctionDictionary() const
{
  return m_functions;
}
 
void
AbstractPdeProblem::setFunctionDictionary(FunctionDictionary* dico)
{
  if (m_ownedFunctions)
    ceps::destroyObject(m_functions);
  m_functions      = dico;
  m_ownedFunctions = false;
}
 
void
AbstractPdeProblem::defineBoundaryConditions()
{
}
 
BoundaryConditionManager*
AbstractPdeProblem::getBoundaryConditionManager() const
{
  return m_boundaryConditions;
}
 
void
AbstractPdeProblem::setBoundaryConditionManager(BoundaryConditionManager* bcm)
{
  if (m_ownedBCs)
    ceps::destroyObject(m_boundaryConditions);
  m_boundaryConditions = bcm;
  m_ownedBCs           = false;
}
 
void
AbstractPdeProblem::defineSourceTerms()
{
}
 
SourceTermManager*
AbstractPdeProblem::getSourceTermManager() const
{
  return m_sourceTerms;
}
 
void
AbstractPdeProblem::setSourceTermManager(SourceTermManager* stm)
{
  if (m_ownedSrcs)
    ceps::destroyObject(m_sourceTerms);
  m_sourceTerms = stm;
  m_ownedSrcs   = false;
}
 
void
AbstractPdeProblem::initializeEquation()
{
  this->defineUnknowns();
  this->createSpatialDiscretization();
  CEPS_SEPARATOR;
  CEPS_SAYS("Setting up elements of the problem...");
  this->defineSourceTerms();
  this->defineBoundaryConditions();
  this->defineAnalyticSolution();
}
 
// -----------------------------------------------------------------------------------------
 
CepsString
AbstractPdeProblem::getOutputFileBase() const
{
  return m_outputFileBase;
}
 
void
AbstractPdeProblem::setOutputFileBase(CepsString s)
{
  m_outputFileBase = ceps::getDir(s) + ceps::getBaseName(s);
}
 
CepsOutputFormat
AbstractPdeProblem::getOutputFormat() const
{
  return m_outputFormat;
}
 
void
AbstractPdeProblem::setOutputFormat(CepsOutputFormat opt)
{
  m_outputFormat = opt;
}
 
CepsBool
AbstractPdeProblem::writesGlobalIndices() const
{
  return m_writeGlobalIDs;
}
 
void
AbstractPdeProblem::writeGlobalIndices(CepsBool opt)
{
  m_writeGlobalIDs = opt;
}
 
CepsVector<CepsReal3D>
AbstractPdeProblem::getProbePoints() const
{
  return m_probePoints;
}
 
CepsReal
AbstractPdeProblem::getStartTime() const
{
  return 0;
}
 
CepsReal
AbstractPdeProblem::getEndTime() const
{
  return 0;
}
 
 
CepsReal
AbstractPdeProblem::getSnapshotTime() const
{
  return 0;
}
 
// -----------------------------------------------------------------------------------------
 
 
void
AbstractPdeProblem::addUnknown(const CepsString& label, CepsSet<CepsAttribute> attrs,
                               CepsLocationFlag loc, const CepsString& unit)
{
  Unknown* u = ceps::getNew<Unknown>(m_unknowns.size(),label);
  if (attrs.empty())
    u->addAttribute(CepsUniversal);
  else
    u->setAttributes(attrs.begin(),attrs.end());
  u->setUnit(unit);
  u->setLocation(loc);
  m_unknowns.push_back(u);
}
 
void
AbstractPdeProblem::addZeroDUnknown(CepsString label, const CepsString& unit)
{
  Unknown* u = ceps::getNew<Unknown>(m_unknowns.size(),label);
  u->addAttribute(CepsUniversal);
  u->setUnit(unit);
  u->setLocation(CepsLocationFlag::ZeroD);
  m_unknowns.push_back(u);
}
 
void 
AbstractPdeProblem::addUnknownInteraction(CepsString label1, CepsString label2, CepsSet<CepsAttribute> attrs)
{
  Unknown* u1 = nullptr;
  Unknown* u2 = nullptr;
 
  for (Unknown* u: m_unknowns)
  {
    if (ceps::toKey(u->getName())==ceps::toKey(label1))
      u1 = u;
    if (ceps::toKey(u->getName())==ceps::toKey(label2))
      u2 = u;
  }
 
  CEPS_ABORT_IF(ceps::isNullPtr(u1),
    "When setting up interactions between unknowns : cannot find unknown " << std::quoted(label1)
  );
  CEPS_ABORT_IF(ceps::isNullPtr(u2),
    "When setting up interactions between unknowns : cannot find unknown " << std::quoted(label2)
  );
  CEPS_ABORT_IF(u1==u2,
    "When setting up interactions between unknowns : unknwons are the same!" 
  );
 
  // Check if given attributes concern the unknowns...
  if (not(attrs.empty() or attrs.contains(CepsUniversal)))
  {
    CEPS_ABORT_IF(not u1->hasAttribute(CepsUniversal) and not u1->hasOneOfAttributes(attrs),
      "When setting up interactions between unknowns : unknown " << std::quoted(label1) << 
      "\n   is not defined on any of the given attributes"
    );
    CEPS_ABORT_IF(not u2->hasAttribute(CepsUniversal) and not u2->hasOneOfAttributes(attrs),
      "When setting up interactions between unknowns : unknown " << std::quoted(label2) << 
      "\n   is not defined on any of the given attributes"
    );
  }
 
  UnknownInteraction* ui = ceps::getNew<UnknownInteraction>();
  ui->u1 = u1;
  ui->u2 = u2;
  ui->setAttributes(attrs.begin(),attrs.end());
  m_unknownsInteractions.push_back(ui);
 
  u1->addUnknownInteraction(u2);
  u2->addUnknownInteraction(u1);
 
}