AbstractDiscretization.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/AbstractDiscretization.hpp"
#include "pde/problem/AbstractPdeProblem.hpp"
#include "pde/common/functions/SAFunc.hpp"
 
AbstractDiscretization::AbstractDiscretization(AbstractPdeProblem* problem) :
    CepsObject               (),
    m_problem                (problem),
    m_distributedDofs        (new DistributedInfos<DegreeOfFreedom*>()),
    m_dofsFactory            (nullptr),
    m_dofsBuilt              (false),
    m_spatialUnknownsLocation(CepsLocationFlag::ZeroD), // Voluntarily wrong
    m_extraAdjacencyToSend   ({}),
    m_extraAdjacencyToRecv   ({})
{
  CEPS_ABORT_IF(ceps::isNullPtr(m_problem),"Pde problem pointer is nullptr");
  m_geometry = problem->getGeometry();
  CEPS_ABORT_IF(m_problem->getUnknowns().empty(),
    "Problem has no unknowns set ! Make sure to call defineUnknowns() of the problem class"
  );
  m_extraAdjacencyToSend.resize(ceps::getGridSize());
}
 
AbstractDiscretization::~AbstractDiscretization()
{
   // Clear dofs
  for (CepsUInt i=0;i<CepsLocationFlagSize;i++)
    for (auto [gID,dofs] : m_dofsTree[i])
      for (DegreeOfFreedom* dof: dofs)
        ceps::destroyObject(dof);
 
  ceps::destroyObject(m_distributedDofs);
  ceps::destroyObject(m_dofsFactory);
}
 
AbstractPdeProblem*
AbstractDiscretization::getProblem() const
{
  return m_problem;
}
 
DistributedFactory*
AbstractDiscretization::getDofsFactory() const
{
  return m_dofsFactory;
}
 
DistributedInfos<DegreeOfFreedom*>*
AbstractDiscretization::getDistributedDofs() const
{
  return m_distributedDofs;
}
 
const DegreeOfFreedomTree&
AbstractDiscretization::getDegreeOfFreedomTree() const
{
  return m_dofsTree;
}
 
const CepsMap<CepsUnknownIndex,CepsVector<DegreeOfFreedom*>>&
AbstractDiscretization::getDegreesOfFreedomForUnknown()
{
  return m_dofsForUnknown;
}
 
const CepsVector<DegreeOfFreedom*>&
AbstractDiscretization::getDegreesOfFreedomForUnknown(CepsUnknownIndex uId)
{
  return m_dofsForUnknown.at(uId);
}
 
const CepsVector<DegreeOfFreedom*>&
AbstractDiscretization::getDegreesOfFreedomForUnknown(Unknown* u)
{
  return m_dofsForUnknown.at(u->getIdentifier());
}
 
Geometry*
AbstractDiscretization::getGeometry() const
{
  return m_geometry;
}
 
void
AbstractDiscretization::registerSpatialUnkown(Unknown* u)
{
  u->setLocation(m_spatialUnknownsLocation);
}
 
 //=================================================================================================
 
DVecPtr
AbstractDiscretization::newDofVector() const
{
  return m_dofsFactory->getDistributedVector();
}
 
DHVecPtr
AbstractDiscretization::newDofHaloVector() const
{
  return m_dofsFactory->getDistributedHaloVector();
}
 
DMatPtr
AbstractDiscretization::newDofMatrix() const
{
  return m_dofsFactory->getDistributedMatrix();
}
 
DegreeOfFreedom*
AbstractDiscretization::getDof(CepsGlobalIndex globalDofId)
{
  DegreeOfFreedom * dof = m_distributedDofs->getOwnedMapping()->atGlobal(globalDofId,nullptr);
  if (ceps::isNullPtr(dof))
    dof = m_distributedDofs->getHaloMapping()->atGlobal(globalDofId,nullptr);
  return dof;
}
 
const DegreeOfFreedom*
AbstractDiscretization::getDof(CepsGlobalIndex globalDofId) const
{
  const DegreeOfFreedom* dof = m_distributedDofs->getHaloMapping()->atGlobal(globalDofId,nullptr);
  if (ceps::isNullPtr(dof))
    dof = m_distributedDofs->getHaloMapping()->atGlobal(globalDofId,nullptr);
  return dof;
}
 
CepsVector<DegreeOfFreedom*>
AbstractDiscretization::getOwnedDofs() const
{
  return m_distributedDofs->getOwned();
}
 
CepsVector<DegreeOfFreedom*>
AbstractDiscretization::getHaloDofs() const
{
  return m_distributedDofs->getHalo();
}
 
void
AbstractDiscretization::evaluateFunctionOnDofs(
  ceps::Function<CepsReal(CepsStandardArgs)>* func, 
  DHVecPtr v, 
  CepsReal t
) const
{
  v->zero();
  v->getLocalData();
  for (auto toParse: {m_distributedDofs->getOwned(),
                      m_distributedDofs->getHalo()})
    for (DegreeOfFreedom* dof : toParse)
    {
      CepsStandardArgs args = getStandardArgsFrom(dof);
      args.t = t;
      (*v)[dof->getGlobalIndex()] = func->operator()(args);
    }
  v->releaseLocalData();
}
 
void
AbstractDiscretization::findClosestDof(const CepsReal3D& x, Unknown* u, CepsDofGlobalIndex& dofId, CepsProcId& owner)
{
  auto distanceToX = [&x] (CepsVertex y) {                      
    return std::sqrt(pow(x[0]-y[0],2) + pow(x[1]-y[1],2) + pow(x[2]-y[2],2));
  };
 
  // Find the local minimum on current process
  CepsReal           dMin = std::numeric_limits<CepsReal>::max();
  CepsDofGlobalIndex iMin = -1;
  for(DegreeOfFreedom* dof: m_distributedDofs->getOwned())
    if (dof->getUnknown()==u)
    {
      CepsReal d = distanceToX(*dof->getVertex());
      if (d < dMin)
      {                                                                           
        dMin = d;
        iMin = dof->getGlobalIndex();
      }
    }
 
  CepsVector<CepsReal> myD(ceps::getGridSize(),0.),d(ceps::getGridSize(),0.);
  myD[ceps::getRank()] = dMin;
  MPI_Allreduce(myD.data(),d.data(),ceps::getGridSize(),CEPS_MPI_REAL,MPI_SUM,ceps::getCommunicator());
  owner = ceps::argmin(d);
 
  MPI_Bcast(&iMin,1,CEPS_MPI_GLOBAL_INDEX,owner,ceps::getCommunicator());
  dofId = iMin;
 
  return;
}
 
CepsReal
AbstractDiscretization::getZeroDValue(Unknown* u, DHVecPtr sol)
{
  DegreeOfFreedom* dof   = getDegreesOfFreedomForUnknown(u)[0];
  CepsReal         value = 0.0;
 
  sol->getLocalData();
  if (dof->getOwner() == ceps::getRank())
    value = (*sol)[dof->getGlobalIndex()];
  sol->releaseLocalData();
 
  MPI_Bcast(&value, 1, CEPS_MPI_REAL, dof->getOwner(), ceps::getCommunicator());
  return value;
}
 
CepsReal
AbstractDiscretization::l1Norm(
  DHVecPtr                      v,
  CepsBool                      boundary,
  const CepsSet<CepsAttribute>& attrs,
  const CepsVector<Unknown*>&   unknowns
) 
{
  DHVecPtr av (ceps::getNew<DistributedHaloVector>(*v));
  DHVecPtr one(ceps::getNew<DistributedHaloVector>(*v));
  av ->abs();
  one->fill(1.);
  return this->dotProduct(one,av,boundary,attrs,unknowns);
}
 
CepsReal
AbstractDiscretization::l2Norm(
  DHVecPtr                      v,
  CepsBool                      boundary,
  const CepsSet<CepsAttribute>& attrs,
  const CepsVector<Unknown*>&   unknowns
) 
{
  return std::sqrt(this->dotProduct(v,v,boundary,attrs,unknowns));
}
 
 //=================================================================================================
 
void
AbstractDiscretization::buildDofs()
{
 
  // 1. Build the dofs according to their location and number 
  // (FE P1, FE P2, FV, etc...) from data in the unknown manager and geometry
  this->buildDofsTree();
 
  // 2. Build the local ID to global ID mapping for dofs
  buildDofsDistributedInfo();
 
  // 3. Add dofs for zero D unknowns, if any
  createDofsForZeroDUnknowns();
 
  // 4. Set the interactions between dofs
  setSpatialUnknownsDofsInteractions();
  setZeroDUnknownsDofsInteractions();
 
  // 5. Setup the dofs factory
  setupDofsFactory();
 
  m_dofsBuilt = true;
}
 
void
AbstractDiscretization::buildDofsDistributedInfo()
{
  const CepsUInt& rank = ceps::getRank ();
  // ---------------------------------------------------------
  // Count owned and halo
  CepsUInt nbOwned = 0U;
  CepsUInt nbHalo  = 0U;
  for (auto locSlice: m_dofsTree)
    for (auto [gid,setOfDofs]: locSlice)
      for (auto dof: setOfDofs)
        if (dof->getOwner()==rank)
          nbOwned++;
        else
          nbHalo++;
 
  m_distributedDofs->getOwnedMapping()->reserve(nbOwned);
  m_distributedDofs->getHaloMapping ()->reserve(nbHalo);
 
  // ---------------------------------------------------------
  // Create dofs owned and halo, with hash from (u.identifier and dof gid)
  CepsUInt nu = m_problem->getUnknowns().size();
  for (auto locSlice: m_dofsTree)
    for (auto [gid,setOfDofs]: locSlice)
    {
      for (auto dof: setOfDofs)
      {
        CepsHash hash = ceps::hash::get(dof->getUnknown()->getIdentifier()+nu*dof->getDiscrId());
        m_distributedDofs->add(dof,hash,dof->getOwner());
      }
    }
 
  // ---------------------------------------------------------
  // Synchronize each procs
  m_distributedDofs->synchronize(true, &setGlobalIndexFunction<DegreeOfFreedom>);
 
  // ---------------------------------------------------------
  // Construct the "reversed" map (unknown -> dofId)
  m_dofsForUnknown.clear();
 
  for (Unknown* u : m_problem->getUnknowns())
    m_dofsForUnknown[u->getIdentifier()] = {}; 
 
  for (auto toParse : {m_distributedDofs->getOwned(),m_distributedDofs->getHalo()})
    for (DegreeOfFreedom* dof: toParse)
      m_dofsForUnknown[dof->getUnknown()->getIdentifier()].push_back(dof);
 
  return;
}
 
void
AbstractDiscretization::createDofsForZeroDUnknowns()
{
  // The new dofs will be registered after all spatial dofs.
  // So they will be placed in the last rows of the linear system.
  // We give them to the last CPU of the computational grid.
 
  CepsDofGlobalIndex dofId(m_distributedDofs->getTotalNumberOfEntities());
  
  CepsUInt iloc = static_cast<CepsUInt> (CepsLocationFlag::ZeroD);
 
  for (Unknown*u : m_problem->getUnknowns())
    if (u->getLocation()==CepsLocationFlag::ZeroD)
    {
      // The status of the dof is different from spatial dofs.
      // We give it attribute -1, and the unknown's id as geom id.
 
      CepsUnknownIndex uId = u->getIdentifier();
 
      auto dof = ceps::getNew<DegreeOfFreedom>(
        u,CepsSet<CepsAttribute>({CepsUniversal}),
        uId,dofId,ceps::getGridSize()-1
      );
      // Insert dof in tree
      // CepsGlobalIndex           guId  = CepsGlobalIndex(uId);
      // CepsMap<CepsGlobalIndex,CepsSet<DegreeOfFreedom*>> tmp;
      // tmp[guId] = {};
      // tmp[guId].insert(dof);
      m_dofsTree[iloc][CepsGlobalIndex(uId)] = {dof};
      // Insert dof in distributed info
      m_distributedDofs->add(dof,dofId,dofId,ceps::getGridSize()-1);
      
      // Insert dof in reverse map
      m_dofsForUnknown[u->getIdentifier()].push_back(dof);
      dofId++; // next zeroD unknown 
 
    }
 
  m_distributedDofs->synchronize(false,nullptr);
}
 
void
AbstractDiscretization::setZeroDUnknownsDofsInteractions()
{
 
  for (UnknownInteraction* ui: m_problem->getUnknownsInteractions())
  {
    Unknown* u1 = ui->u1;
    Unknown* u2 = ui->u2;
    if (not (u1->isSpatial() and u2->isSpatial()))
    {
      // At least one of those two has only one element
      for (DegreeOfFreedom* dof1: m_dofsForUnknown[u1->getIdentifier()])
      for (DegreeOfFreedom* dof2: m_dofsForUnknown[u2->getIdentifier()])
      {
        CepsBool okAttributeI = ui->hasAttribute(CepsUniversal) 
                             or ui->getNumberOfAttributes()==0;
        CepsBool okAttribute1 = dof1->hasAttribute(CepsUniversal) or (not okAttributeI and 
                                dof1->hasOneOfAttributes(ui->getAttributes()));
        CepsBool okAttribute2 = dof2->hasAttribute(CepsUniversal) or (not okAttributeI and 
                                dof2->hasOneOfAttributes(ui->getAttributes()));
 
        if (okAttributeI or (okAttribute1 and okAttribute2))
        {
          dof1->addNeighbor(dof2);
          dof2->addNeighbor(dof1);
        }
      }
 
    }
  }
 
}
 
void
AbstractDiscretization::setupDofsFactory()
{
 
  ceps::barrier();
  // Get information on number of dofs
  CepsUInt nDofs  = m_distributedDofs->getTotalNumberOfEntities();
  CepsUInt nLDofs = m_distributedDofs->getNumberOfOwned();
 
  // Prepare the non zero template
  CepsVector<CepsInt> dNnz,oNnz;
 
  dNnz.resize(nLDofs,0U);
  oNnz.resize(nLDofs,0U);
 
  // Neighbours can be built asymetrically, so the oNnz is not correct if deduced directly
  CepsVector<CepsMap<CepsDofGlobalIndex,CepsUInt>> othersONnz(ceps::getGridSize());
 
  for (CepsUInt j=0U; j<nLDofs; ++j)
  {
    // Fetch dof adjacency seen by other cpus
    DegreeOfFreedom*            dof   = m_distributedDofs->getOwned(j);
    CepsDofGlobalIndex          dofId = dof->getGlobalIndex();
    CepsSet<CepsDofGlobalIndex> extraNgb(m_extraAdjacencyToRecv[dofId]);
 
    for (DegreeOfFreedom *ndof: m_distributedDofs->getOwned(j)->getNeighbors())
    {
      CepsProcId         owner  = ndof->getOwner();
      CepsDofGlobalIndex ndofId = ndof->getGlobalIndex();
 
      extraNgb.erase(ndofId);
      if (owner!=ceps::getRank())
      {
        oNnz[j]++;
        if (not othersONnz[owner].contains(ndofId))
          othersONnz[owner][ndofId] = 1;
        else
          othersONnz[owner][ndofId]++;
      }
      else
        dNnz[j]++;
    }
 
    // +1 in diagonal because we need to count the dof itself
    // + the dofs relationships that were seen by other processors only
    dNnz[j] += 1+extraNgb.size();
 
  }
 
  // Some interactions can be missing in the halo, when other dofs are
  // distributed with two other processes
  for (DegreeOfFreedom* hdof: m_distributedDofs->getHalo())
  {
    CepsProcId hOwner = hdof->getOwner();
    for (DegreeOfFreedom* ndof: hdof->getNeighbors())
    {
      CepsProcId         nOwner = ndof->getOwner();
      CepsDofGlobalIndex nDofId = ndof->getGlobalIndex();
      if (hOwner!=nOwner)
      {
        if (not othersONnz[nOwner].contains(nDofId))
          othersONnz[nOwner][nDofId] = 1;
        else
          othersONnz[nOwner][nDofId]++;
      }
    }
  }
 
 
  // Send othersONnz data. Not optimal but
  // From cpu1 to cpu2
  ceps::barrier();
  CepsProcId rank     = ceps::getRank();
  CepsUInt   gridSize = ceps::getGridSize();
  for (CepsProcId cpu1 = 0; cpu1<gridSize; cpu1++)
  for (CepsProcId cpu2 = 0; cpu2<gridSize; cpu2++)
  {
 
  if (cpu1!=cpu2 and (rank==cpu1 or rank==cpu2))
  {
    CepsMap<CepsDofGlobalIndex,CepsUInt> tmpMap; 
    if (rank==cpu1)
    {
      tmpMap = othersONnz[cpu2];
    }
 
    // For cpu2, otherONnz[cpu2] is empty before this call (NO!!!!!!!)
    ceps::sendMap(cpu1,cpu2,tmpMap,gridSize*gridSize,cpu1*gridSize+cpu2);
    if (rank==cpu2)
    {
      for (auto [row,nNz] : tmpMap)
      {
        // get local index on cpu2. The global index may not have been registered though...
        // we have to do it manually
        oNnz[row-m_distributedDofs->getGlobalIndexOffset(cpu2)] += nNz;
      }
    }
  }
  ceps::barrier();
  }
 
 
  for (CepsUInt j=0U; j<nLDofs; ++j)
  {
    // Sometimes the patch the othersONnz leads to unused allocated space, that can, in some cases,
    // even exceed the total size of the matrix (e.g. a lagrangian coefficient)
    oNnz[j] = std::min(oNnz[j],CepsInt(nDofs-nLDofs));
 
    // CEPS_ABORT_IF(CepsUInt(dNnz[j]+oNnz[j])>nDofs,
    //   "Total number of nnz on line " << j << " exceeds size of matrix.\n"
    //   "   Global size: " << nDofs << " x " << nDofs << "\n" <<
    //   "   Local size: " << nLDofs << " x " << nLDofs << "\n" <<
    //   "   In diagonal block: " << dNnz[j] << "\n" <<
    //   "   In off-diagonal blocks: " << oNnz[j] << "\n" <<
    //   "   Total number of coefficients on this row: " << dNnz[j]+oNnz[j]  << "\n" <<
    //   "   othersONnz " << othersONnz[0] << " " << othersONnz[1] << " " << othersONnz[2] << "\n"
    // );
  }
 
  CepsVector<CepsIndex> haloIndices;
  haloIndices.reserve(m_distributedDofs->getHalo().size());
  for (DegreeOfFreedom* dof: m_distributedDofs->getHalo())
  {
    haloIndices.push_back(dof->getGlobalIndex());
  }
 
  // Get the factory
  m_dofsFactory = ceps::getNew<DistributedFactory>(nDofs,nDofs,nLDofs,nLDofs,dNnz,oNnz);
  m_dofsFactory->setHaloIndices(haloIndices.data(),haloIndices.size(),haloIndices.data(),haloIndices.size());
 
}
 
void
AbstractDiscretization::extractValuesForUnknownInternal(
  const CepsUnknownIndex&                    uId,
  DVecPtr                                    vec,
  CepsVector<CepsIndex>&                     indices,
  CepsVector<CepsReal>&                      values,
  CepsBool                                   sendToMaster,
  CepsBool                                   returnGeomIndices,
  std::function<CepsBool(DegreeOfFreedom*)>& selector
)
{
  // Firstly, test if the u exists
  CepsBool hasUnknown = false;
  for (Unknown* u: m_problem->getUnknowns())
    hasUnknown = hasUnknown or u->getIdentifier()==uId;
 
  if (not hasUnknown or not m_dofsForUnknown.contains(uId))
    return;
 
  CepsVector<CepsIndex> resI;
  CepsVector<CepsReal>  resV;
 
  // Be careful of sizes ! If regroup on master: remove halo nodes
  CepsVector<DegreeOfFreedom*> dofs  = m_dofsForUnknown.at(static_cast<CepsSize>(uId));
 
  CepsUInt resSize = 0u;
  for (DegreeOfFreedom* dof:dofs)
    if (selector(dof))
      resSize++;
 
  // Resize to correct size the IndicesValues struct
  resI.reserve(resSize);
  for (DegreeOfFreedom* dof:dofs)
    if (selector(dof))
      resI.push_back(dof->getGlobalIndex());
 
  // Get the correct entries of the vector
  resV.resize(resSize,0.);
  vec->getLocalData();
  vec->getValues(resV.data(),resSize,resI.data());
  vec->releaseLocalData();
 
  // Overwrite with geom indices if needed
  if (returnGeomIndices)
  {
    resI.clear();
    resI.reserve(resSize);
    for (DegreeOfFreedom* dof:dofs)
      if (selector(dof))
        resI.push_back(this->getDofSpatialId(dof));
  }
 
  if (sendToMaster)
  {
    ceps::allGatherv(resI,resI.size(),indices,true);
    ceps::allGatherv(resV,resV.size(),values ,true);
  }
  else
  {
    values  = resV;
    indices = resI;
  }
 
  return;
}