bddc.hh

Go to the documentation of this file.

 
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2022-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef BDDC_HH
#define BDDC_HH
 
#include "fem/firstless.hh"
#include "linalg/direct.hh"
#include "linalg/localMatrices.hh"
#include "mg/bddcSpaceTransfer.hh"
#include "utilities/enums.hh"
 
#include "dune/istl/bcrsmatrix.hh"
 
namespace Kaskade::BDDC
{
  using CoarseConstraints = std::vector<std::pair<std::vector<int>,int>>;
 
  struct LocalDof
  {
    int s; // subdomain number
    int i; // local dof index in subdomain
  };
 
  enum InterfaceType { CORNER = 1, EDGE = 2, FACE = 4 };
 
  template <int m=1, class Index=size_t>
  class InterfaceAverages
  {
  public:
    constexpr static int components = m;
 
    InterfaceAverages(std::vector<std::vector<LocalDof>> const& sharedDofs,
                      std::vector<Index> const& subdomainSize,
                      int types = CORNER | FACE);
 
    Interfaces const& interfaces(int id) const;
 
    template <class Scalar=double>
    NumaBCRSMatrix<Dune::FieldMatrix<Scalar,m,m>> coarseConstraint(int id) const;
 
    CoarseConstraints coarseConstraints() const;
 
    std::array<int,3> interfaceCount() const
    {
      return ifCount;
    }
 
  private:
    std::vector<Interfaces> ifs;
    CoarseConstraints ccn;
    std::vector<Index> subdomainSize;
 
    std::vector<std::vector<std::vector<int>>> cifs; // coarse interfaces
    std::array<int,3> ifCount;                       // number of corners, edges, faces
  };
 
  // ----------------------------------------------------------------------------------------------
 
 
  // ----------------------------------------------------------------------------------------------
 
 
  template <int m, class Scalar_=double, class CoarseScalar=Scalar_,
            class Transfer=SpaceTransfer<m,Scalar_> >
  class Subdomain
  {
  public:
    static const int components = m;
 
    using Scalar = Scalar_;
 
    using XEntry = Dune::FieldVector<Scalar,components>;
    using XVector = Dune::BlockVector<XEntry>;
    using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
 
    using TransmissionType = typename Transfer::TransmissionType;
 
 
    using AMatrix = NumaBCRSMatrix<Dune::FieldMatrix<Scalar,components,components> >;
 
    using SMatrix = DynamicMatrix<Dune::FieldMatrix<Scalar,1,1>>;
 
    template <int cComponents>
    Subdomain(int id, AMatrix const& A,
              Interfaces const& interfaces,
              NumaBCRSMatrix<Dune::FieldMatrix<Scalar,cComponents,components>> const& C);
 
    template <class Interface>
    Subdomain(int id, AMatrix const& A, Interface const& ifs);
 
    void imposeDC_inActiveInterface(Vector & Aq, Vector &q);
    void imposeDCq(Vector &q);
 
    void update_rhs(XVector const& f);
 
    void getInterfaceAveragingWeights(int s, XVector& rho) const
    {
      trans.getInterfaceAveragingWeights(s,rho);
    }
 
    void setInterfaceAveragingWeights(int s, XVector const& rho)
    {
      trans.setInterfaceAveragingWeights(s,rho);
    }
 
    std::vector<int> neighbors() const;
    
    void preRestrict();
 
    void getInterfaceResidualUpdate(int s, TransmissionType& dres) const
    {
      trans.getInterfaceResidualUpdate(s,dres);
    }
 
    void setInterfaceResidualUpdate(int s, TransmissionType const& dres)
    {
      trans.setInterfaceResidualUpdate(s,dres);
    }
 
    void restrict();
    
 
    void getS(SMatrix& S) const;
 
    void getSchurResidual(Vector& lambda) const;
 
    void setCorrection(Vector const& c);
    
 
    void preProlongate()
    {
      trans.preProlongate(du);
    }
 
    void getInterfaceValues(int s, TransmissionType& values) const
    {
      trans.getInterfaceValues(s,values);
    }
 
    void setInterfaceValues(int s, TransmissionType const& values)
    {
      trans.setInterfaceValues(s,values);
    }
 
    void setActiveId(bool active){
      activeId = active;
    }
    Dune::FieldVector<double,2> prolongate();
 
    Dune::FieldVector<double,1> updateSearchDirectionNorm();
 
    void updateSearchDirection(double beta);
    
 
    void updateSolution(Scalar alpha);
 
    void getSolution(XVector& x) const;
    
    XVector getSolution() const;
 
    void set_solver_type(bool cg_solver);
    void updateIteration(int cg_itr);
    void getCorrection(XVector& dx) const;
    void getRawCorrection(XVector& dx) const;
    void getResidual(XVector& rr) const;
    Vector const& getRestrictedResidual() const { return restrictedResidual; }
 
 
    std::vector<bool>  getActiveSubdomains()
    {
      return activeSubdomains_sub;
    }
 
    void setActiveSubdomains_sub(std::vector<bool> & activeSubdomains_)
    {
      activeSubdomains_sub.resize(activeSubdomains_.size());
      // std::cout<< "activeSubdomains_.size() = "<<activeSubdomains_.size() <<std::endl;
      for (int i = 0; i < activeSubdomains_.size(); ++i)
      {
        activeSubdomains_sub[i] = activeSubdomains_[i];
      }
    }
 
 
  private:
    using BMatrix = NumaBCRSMatrix<Dune::FieldMatrix<Scalar,1,1>>;
 
    // These variables describe problem size and connectivity
    int id;                              // my own subdomain number
    int N;
    int nx, nc;                          // number of dofs and of constraints
    Interfaces interfaces;               // list of pairs (s,ids) of adjacent subdomain s and shared dofs
    Transfer trans;                      // manages data conversion and averaging
    
    Vector f;                            // the right hand side
    Vector u;                            // approximate solution [x,lambda]
    Vector r;                            // current residual [f-Au,0]
    Vector rRes;                         // restricted residual
 
    Vector q;                            // search direction for [CG]
    Vector q_pre;                        // search direction for [CG]
    Vector Aq;                           // for [CG]
    Vector r_bar;                        // previous residual for [CG] 
    Scalar alpha;                        // step size for [CG]
    bool cg_solver;                      // cg = false -> steepest descent method with BDDC, cg = true -> CG with BDDC
    int cg_itr; 
 
    Vector du;                           // correction [dx,0]
    Vector Adu;                          // matrix times correction
    
    BMatrix A;
    BMatrix C;
    DirectSolver<Vector,Vector> Binv;    // [ A C^T; C 0 ]^{-1}
    
    bool activeId = true;
 
    std::vector<int> interiorDofs;       // the interior dof indices
    std::vector<int> interfaceDofs;      // the boundary/interface dof indices
    std::vector<bool> activeSubdomains_sub;
    // for debugging purposes
    Vector rawCorrection, residual, restrictedResidual;
 
    std::vector<int> globalIndex_notActive;
 
BMatrix Aii;
DirectSolver<Vector,Vector> AiiInv;
std::vector<int> cIndices;
mutable DynamicMatrix<double> X; // const hack, move computation to constructor
 
    // Solves the homogeneous Dirichlet problem A_ii x_i = r_i for the interior nodes only.
    // Non-zero boundary values x_D can be imposed by exploiting the splitting
    // [Aii AiD] [x_i x_D]^T = r_i <=> Aii x_i = r_i - AiD xD, i.e. modifying the right hand side.
    // The interface nodes of x are not touched, nor is the constraint/multiplier part.
    void solveDirichlet(Vector& x, Vector const& rhs) const;
  };
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
  
  template <class Subdomain>
  class SharedMemoryDomain
  {
  public:
    using Subdom = Subdomain;
    using Scalar = typename Subdom::Scalar;
    using Vector = typename Subdom::Vector;
    using XVector = typename Subdom::XVector;
    static constexpr int m = Subdom::components;
    
    SharedMemoryDomain(std::vector<Subdom>& subdomains);
 
    Dune::FieldVector<double,2> prolongate();
 
    void restrict();
 
    void updateSearchDirection(double beta); // CG
 
    int size() const
    {
      return subdomains->size();
    }
    
    std::vector<DynamicMatrix<Dune::FieldMatrix<Scalar,1,1>>> getS() const;
    std::vector<Vector> getSchurResiduals();
    void setCorrections(std::vector<Vector> const& dcs);
 
    Dune::FieldVector<double,1> updateSearchDirectionNorm(); // CG
    void updateSolution(Scalar alpha);
 
    void set_solver_type(bool cg_solver);
    void updateIteration();
 
    std::vector<bool>  getActiveSubdomains()
    {
      return activeSubdomains;
    }
 
    void setActiveSubdomains(std::vector<int> & activeIds)
    {
      std::cout << "activeIds.size() = " << activeIds.size() <<std::endl;
      activeSubdomains.resize(size());
      for (int i = 0; i < activeIds.size(); ++i)
      {
        activeSubdomains[activeIds[i]] = true;
      }
    }
    std::array<int,6> traffic() const
    {
      return std::array<int,6>{exchangeCount,transfers,initTraffic,restrictionTraffic,
                               prolongationTraffic,coarseGridTraffic};
    }
    
  private:
    std::vector<Subdom>* subdomains;
    std::vector<std::pair<int,int>> neighbors;
    std::vector<bool> activeSubdomains;
 
    // Support for monitoring the amount of data exchanged between subdomains.
    std::mutex mutex;
    int exchangeCount = 0;        // number of iterations
    int transfers = 0;            // number of bidirectional data transfers
    int initTraffic = 0;          // amount of data (in bytes) sent during initialization
    int restrictionTraffic = 0;   // amount of data (in bytes) sent during restriction
    int prolongationTraffic = 0;  // amount of data (in bytes) sent during prolongation
    int coarseGridTraffic = 0;    // amount of data (in bytes) sent to coarse grid solver
  };
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  template <class Domain>
  class KKTSolver
  {
  public:
    using Scalar = typename Domain::Scalar;
    using Constraint = std::tuple<int,int,int>;
 
  private:
    using Vector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
 
  public:
    KKTSolver(Domain& domain_,
              CoarseConstraints const& constr)
    : domain(&domain_)
    , constraints(domain->size())
    , localConstraintsSize(domain->size())
    , sOffsets(constr.size()+1)
    {
      // Compute the offsets of individual slack vectors s_ij into the global slack vector
      // and the constraints in which each subdomain is involved.
      sOffsets[0] = 0;
      for (int k=0; k<constr.size(); ++k)
      {
        auto const& [s,m] = constr[k];
        for (auto i: s)
        {
          assert(0<=i && i<domain->size());
          constraints[i].push_back(k);
          localConstraintsSize[i] += m;
        }
        sOffsets[k+1] = m;
      }
      std::partial_sum(begin(sOffsets),end(sOffsets),begin(sOffsets));
 
      // Create Schur complement matrix B. This is sparse, since each constraint, i.e.
      // slack variable, is only affected by two subdomains.
      int const ns = sOffsets.back();                                 // Total number of slacks.
      if (ns == 0)                                                    // If there are no coarse constraints,
        return;                                                       // we're done - there's nothing to do.
      NumaCRSPatternCreator<> sCreator(ns,ns,false,2);
      for (auto const& ss: constraints)                               // For each subdoman find all pairs
        for (int sRow: ss)                                            // of incident slacks and add the
          for (int sCol: ss)                                          // corresponding dense block into
            sCreator.addDenseBlock(sOffsets[sRow],sOffsets[sRow+1],   // the sparsity pattern of S.
                                   sOffsets[sCol],sOffsets[sCol+1]);
 
      using Entry = Dune::FieldMatrix<Scalar,1,1>;
      NumaBCRSMatrix<Entry> S(sCreator);                              // And then create the matrix.
 
      // Assemble the Schur complement matrix.
      using SortedIdx = std::vector<std::pair<int,int>>;
      LocalMatrices<Entry,false,SortedIdx,SortedIdx> lm(S);
      auto localSchurComplements = domain->getS();
      for (int i=0; i<domain->size(); ++i)
      {
        SortedIdx rows;
        int localIndex = 0 ;
        for (int k: constraints[i])
          for (int l=sOffsets[k]; l<sOffsets[k+1]; ++l)
            rows.push_back({l,localIndex++});
        assert(std::is_sorted(rows.begin(),rows.end(),FirstLess()));
        lm.push_back(rows,rows);
        lm.back() = localSchurComplements[i];                         // add local stiffness matrix Si
      }
      lm.scatter();
      S *= -1.0;                                                      // take care of correct Schur complement sign
 
      Sinv = DirectSolver<Vector,Vector>(S,DirectType::UMFPACK,
                                         MatrixProperties::SYMMETRICSTRUCTURE);
    }
 
    void solve()
    {
      int const ns = sOffsets.back();                 // Total number of slacks.
 
      if (ns==0)
        return;
      // std::cout << "ns---> " << ns <<std::endl;
      // Obtain right hand side of the Schur complement system. This is the sum of right
      // hand sides from the subdomains, scattered into appropriate positions in the
      // right hand side vector f.
      Timings& timer = Timings::instance();
      timer.start("get Schur rhs");
      Vector f(ns);
      getRightHandSide(f);
      timer.stop("get Schur rhs");
 
  
      std::vector<int> d(ns,1);
      for (int i=0; i<domain->size(); ++i)
      {
        if(not domain->getActiveSubdomains()[i])   // vector d is a {0,1} where zero is for the interfaces 
        {                                          // belongs to inactive subdomains
          for (auto const k: constraints[i])          
          {                                            
            int const off = sOffsets[k];               
            int const nsk = sOffsets[k+1]-off;
            for (int l=0; l<nsk; ++l)
              d[l+off] = 0;
          }
        }
      }
 
      for (int i=0; i<ns; ++i)                        // assign zero on the rhs for those belongs to  
        f[i] *= d[i];                                 // inactive subdomains via vector d
 
      timer.start("solve Schur system");
      Vector s(ns);                                   // Create zero initial slack vector
      Sinv.apply(s,f);                                // and solve the global Schur system.
      timer.stop("solve Schur system");
      
 
      for (int i=0; i<ns; ++i)                        // assign zero on the slacks corresponding to  
        s[i] *= d[i];                                 // inactive subdomains via vector d
 
 
      timer.start("set corrections");
      std::vector<Vector> dcs(domain->size());
      for (int i=0; i<domain->size(); ++i)
      {
        if(domain->getActiveSubdomains()[i])
        {
          Vector si(localConstraintsSize[i]);           // gather the slacks corresponding to the
          int pos = 0;                                  // constraints affecting active subdomain i
          for (auto const k: constraints[i])            // by concatenating all the slacks
          {                                             // corresponding to constraints in which
            int const off = sOffsets[k];                // subdomain i is involved.
            int const nsk = sOffsets[k+1]-off;
            for (int l=0; l<nsk; ++l)
              si[pos++] = s[l+off];
          }
          dcs[i] = si;
        }
      }
      domain->setCorrections(dcs);
      timer.stop("set corrections");
    }
 
 
  private:
 
    void getRightHandSide(Vector& f) const
    {
      f = 0;
      auto floc = domain->getSchurResiduals();
      
      for (int i=0; i<domain->size(); ++i)
      {
        if(domain->getActiveSubdomains()[i])
        {
          int pos = 0;                                  // scatter the lambda solution component
          for (auto const k: constraints[i])            // into the right hand side vector
          {
            int const off = sOffsets[k];
            int const nsk = sOffsets[k+1]-off;
            for (int l=0; l<nsk; ++l)
              f[l+off] += floc[i][pos++];
          }
        }
      }
    }
 
 
    Domain* domain;
    std::vector<std::vector<int>> constraints;        // involved constraints by subdomain id
    std::vector<int> localConstraintsSize;            // total number of constraints for subdomain id
    std::vector<int> sOffsets;                        // offset of slack vector in total slacks
    DirectSolver<Vector,Vector> Sinv;                 // inverse of Schur complement
  };
 
 
  template <class Subdomain>
  class BDDCSolver // TODO: implement in terms of SymmetricPreconditioner
  {
  public:
 
  
    using Subdom = Subdomain;
    using Scalar = typename Subdom::Scalar;
    using Vector = typename Subdom::Vector;
    static constexpr int m = Subdom::components;
 
    using XEntry = Dune::FieldVector<Scalar,m>;
    using XVector = Dune::BlockVector<XEntry>;
 
    BDDCSolver(std::vector<Subdomain>& subdomains,
               CoarseConstraints const& constraints,
               std::vector<int> & activeIds,
               bool cg_solver=false,
               bool BDDC_verbose=false);
 
    void update_rhs(std::vector<XVector> & Fs)
    {
      parallelFor(0,subdomains->size(),[&](int i)
      {
        (*subdomains)[i].update_rhs(Fs[i]);
      });
    }
    double solve();
 
    std::array<int,6> traffic() const
    {
      return domain.traffic();
    }
 
  private:
    using Domain = SharedMemoryDomain<Subdomain>;
    std::vector<Subdomain>* subdomains;
    Domain domain;
    CoarseConstraints constraints;
    KKTSolver<Domain> coarseSolver;
    std::vector<XVector> const Fs;
    bool cg_solver;
    double sigma_pre = 0;
    double BDDC_verbose = false;
  };
 
}
 
#endif