taylorHoodPreconditioner.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) 2023-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef TAYLORHOODPRECONDITIONER_HH
#define TAYLORHOODPRECONDITIONER_HH
 
#include "dune/common/config.h"
 
 
#include "fem/linearspace.hh"
#include "fem/supports.hh"
#include "linalg/symmetricOperators.hh"
#include "linalg/threadedMatrix.hh"
 
#include "dune/grid/common/rangegenerators.hh"
#include "dune/istl/preconditioner.hh"
 
#include <fstream>
 
 
namespace Kaskade
{
  namespace TaylorHoodPreconditionerDetail
  {
    template <int dim, class Real=double>
    class LocalStokesSolver
    {
      using Entry = Dune::FieldMatrix<Real,1,1>;
      using Matrix = DynamicMatrix<Entry>;
 
      template <class R, int n>
      using Vector = Dune::BlockVector<Dune::FieldVector<R,n>>;
      
      using EntryA = Dune::FieldMatrix<Real,dim,dim>;
 
 
    public:
      using MatrixA = DynamicMatrix<EntryA>;
      using MatrixB = DynamicMatrix<Dune::FieldMatrix<Real,1,dim>>;
      
      LocalStokesSolver(MatrixA A, MatrixB B, Real tol);
      
      std::vector<int> const& constraintIndices() const;
 
      template <class R, int nu>
      void apply(Vector<R,dim>& u, Vector<R,1>& p,
                 Vector<R,dim> const& f, Vector<R,1> const& g) const;
      
      Matrix const& getInvK() const { return invK; }
 
    private:
      Matrix invK;            // inverse of local Stokes system
      std::vector<EntryA> D;  // square root of inverse diagonal of A
      std::vector<int> idx;   // row indices of given B that are retained.
    };
  }
 
  
  enum class TaylorHoodPreconditionerFlavour
  {
    SYMMETRIC,
    PRESMOOTH,
    POSTSMOOTH
  };
  
  template <class X>
  class TaylorHoodPreconditioner 
  : public Dune::Preconditioner<X,X>
  {
    using LocalSolver = TaylorHoodPreconditionerDetail::LocalStokesSolver<1>;
 
  public:
    using domain_type = X;
    
    using range_type = domain_type;
    
    using field_type = typename X::field_type;
    using Scalar = field_type;
    
    template <class AGOP,  class VelocitySpace>
    TaylorHoodPreconditioner(AGOP const& K,
                             VelocitySpace const& velocitySpace)
    {
      abort();
    }
    
    template <class MatrixA, class MatrixB, class MatrixBt,  class VelocitySpace>
    TaylorHoodPreconditioner(MatrixA const& A, MatrixB const& B, MatrixBt const& Bt,
                             VelocitySpace const& velocitySpace)
    : dim(VelocitySpace::dim)
    {
      auto const& gridView = velocitySpace.gridView();
      static_assert(MatrixA::block_type::rows==1,"matrix blocks should have one row");
      static_assert(MatrixA::block_type::cols==1,"matrix blocks should have one  col");
      assert(A.N() % dim == 0);
      assert(A.M() % dim == 0);
      assert(A.M() == B.M());
 
      // for each patch around a vertex, store the global velocity dof indices
      patchDofs.resize(gridView.size(dim));
      localInverse.resize(patchDofs.size());
 
      auto patches = computePatches(gridView);
      assert(patches.size()==patchDofs.size());
      
      // for each dof, count the number of incident patches
      velocityPatchCount.resize(A.N()/dim,0);
      pressurePatchCount.resize(B.N(),0);
      
// std::ofstream out("out.m");
// out << "A = [" << A << "];\n";
// out << "B = [" << B << "];\n";
// out.close();
      
      for (size_t i=0; i<patches.size(); ++i) // TODO: parallelize as in domainDecompositionPreconditioner.hh
      {
// Check if this is a boundary patch:
bool boundaryPatch = false;
for (auto const& cell: patches[i])
  for (auto const& face: Dune::intersections(gridView,cell))
    if (face.boundary())
      boundaryPatch = true;
 
        // Find all velocity degrees of freedom with support inside the patch
        auto const& velocityDofs = algebraicPatch(velocitySpace,patches[i]);
        
        // The algebraic patch computation and hence the velocityDofs vector contain indices
        // with respect to the vectorial, i.e. multi-component, velocity space, but A is 
        // flattened out with scalar entries. Thus, we need to compute entries into A.
        std::vector<size_t> vDofs;
        vDofs.reserve(dim*velocityDofs.size());
        for (auto j: velocityDofs)
          for (int k=0; k<dim; ++k)
            vDofs.push_back(dim*j+k);
 
        int const nA = vDofs.size();
        assert(nA==dim*velocityDofs.size());
        
        // Find all pressure constraints which affect any of the velocity dofs on 
        // the patch. This is done by collecting all column indices of B' in rows 
        // given by the velocity dofs.
        std::vector<size_t> pressureDofs;
        for (size_t const i: vDofs)
          for (auto ci=Bt[i].begin(); ci!=Bt[i].end(); ++ci)
            pressureDofs.push_back(ci.index());
        
        std::sort(pressureDofs.begin(),pressureDofs.end());
        pressureDofs.erase(std::unique(pressureDofs.begin(),pressureDofs.end()),
                           pressureDofs.end());
        int nB = pressureDofs.size();
        
        // update the patch count for incident dofs (patch count is w.r.t. vectorial space/indices
        for (size_t uIdx: velocityDofs)
          ++velocityPatchCount[uIdx];
        for (size_t pIdx: pressureDofs)
          ++pressurePatchCount[pIdx];
        
 
        // TODO: implement the method above
 
        // using DenseMatrix = DynamicMatrix<Dune::FieldMatrix<double,1,1>>;
        // TaylorHoodPreconditionerDetail::LocalStokesSolver<1>
        //     invK(submatrix<DenseMatrix>(A,vDofs,vDofs),
        //          submatrix<DenseMatrix>(B,pressureDofs,vDofs),
        //          submatrix<DenseMatrix>(Bt,vDofs,pressureDofs));
        // 
        // // Move over to global data structure
        // patchDofs[i][0] = std::move(velocityDofs);
        // patchDofs[i][1] = std::move(pressureDofs);
        // localInverse[i] = std::move(invK);
      }
    }
    
    virtual void apply(domain_type& x, range_type const& y) override
    {
      apply(x,y,TaylorHoodPreconditionerFlavour::SYMMETRIC);
    }
    
    void apply(domain_type& x, range_type const& y, TaylorHoodPreconditionerFlavour flavour) const
    {
      // Initialize to zero, as we need to accumulate the solution contributions of 
      // all the patch systems.
      x = 0;
      
      // TODO: parallelize
      // for (int k=0; k<patchDofs.size(); ++k)
      for (int k=1; k<2; ++k)
      {
std::cout << "patch " << k << " indices: idxu=[";
for(int i=0; i<patchDofs[k][0].size(); ++i) std::cout << 2*patchDofs[k][0][i]+1 << ' ' << 2*patchDofs[k][0][i]+2 << ' ';
std::cout << "]; idxp=[";
for (int i=0; i<patchDofs[k][1].size(); ++i) std::cout << patchDofs[k][1][i]+1 << ' ';
std::cout << "];\n";
 
        // gather velocity and pressure rhs entries with appropriate weight
        int prolongation = 0;
        std::vector<double> r(dim*patchDofs[k][0].size()+patchDofs[k][1].size());
 
std::cout << "weighted urhs=[";
        // velocity rhs
        for (int i=0; i<patchDofs[k][0].size(); ++i)
        {
          auto uIdx = patchDofs[k][0][i];
          double w = weight(velocityPatchCount[uIdx],prolongation,flavour);
          for (int j=0; j<dim; ++j)
          {
            r[i*dim+j] = component<0>(y)[uIdx][j] * w;
std::cout << r[i*dim+j] << ' ';
          }
        }
std::cout << "];\nweighted prhs=[";
 
 
        // pressure rhs
        for (int i=0; i<patchDofs[k][1].size(); ++i)
        {
          auto pIdx = patchDofs[k][1][i];
          double w = weight(pressurePatchCount[pIdx],prolongation,flavour);
          r[dim*patchDofs[k][0].size()+i] = component<1>(y)[pIdx][0] * w;
std::cout << r[dim*patchDofs[k][0].size()+i] << ' ';
        }
std::cout << "];\n";
 
 
        // solve system
        std::vector<double> z(r.size(),0.0);
        // localInverse[k].mv(r,z);
        
        // write velocity and pressure to solution, with appropriate scaling
        prolongation = 1;
 
        // velocity
std::cout << "u = [";
        for (int i=0; i<patchDofs[k][0].size(); ++i)
        {
          auto uIdx = patchDofs[k][0][i];
          double w = weight(velocityPatchCount[uIdx],prolongation,flavour);
          for (int j=0; j<dim; ++j)
          {
            component<0>(x)[uIdx][j] += z[i*dim+j] * w;
std::cout << z[i*dim+j] * w << ' ';
          }
        }
std::cout << "]; \np = [";
 
        // pressure
        for (int i=0; i<patchDofs[k][1].size(); ++i)
        {
          auto pIdx = patchDofs[k][1][i];
          double w = weight(pressurePatchCount[pIdx],prolongation,flavour);
          component<1>(x)[pIdx][0] += z[dim*patchDofs[k][0].size()+i] * w;
std::cout << z[dim*patchDofs[k][0].size()+i] * w << ' ';
        }
std::cout << "];\n";
 
        // no need to output the constant pressure multiplicator
      }
    }
 
    virtual bool requiresInitializedInput() const
    {
      return true;
    }
    
    virtual void pre (domain_type& x, range_type&) override {}
    virtual void post (domain_type& x) override {}
    virtual Dune::SolverCategory::Category category() const override 
    {
      return Dune::SolverCategory::sequential;
    }
    
  private:    
    // determines restriction (prolongation==0) or prolongation (prolongation==1). See class 
    // documentation for a detailed explanation.
    double weight(int count, int prolongation, TaylorHoodPreconditionerFlavour flavour) const
    {
      assert(count > 0);
      double w = 1.0/count;
      
      if (flavour == TaylorHoodPreconditionerFlavour::SYMMETRIC)
        w = std::sqrt(w);
      else if (flavour==TaylorHoodPreconditionerFlavour::PRESMOOTH && prolongation==0)
        w = 1;
      else if (flavour==TaylorHoodPreconditionerFlavour::POSTSMOOTH && prolongation==1)
        w = 1;
      
      return w;
    }
    
    int dim;
    
    // for each patch around a vertex, store the global velocity and pressure dof indices
    std::vector<std::array<std::vector<size_t>,2>> patchDofs;
 
    // for each patch store the local inverse
    std::vector<LocalSolver> localInverse;
    
    // for each global velocity and pressure dof, store the number of incident patches 
    std::vector<int> velocityPatchCount;
    std::vector<int> pressurePatchCount;
  };
  
  // ----------------------------------------------------------------------------------------------
  
  template <class AGOP>
  class SymmetricTaylorHoodPreconditioner
  : public SymmetricPreconditioner<typename AGOP::Scalar,
                                   typename AGOP::Domain>
  {
  public:
    using domain_type = typename AGOP::Domain;
    using range_type  = typename AGOP::Range;
    static_assert(std::is_same_v<domain_type,range_type>,"need a symmetric saddle point operator");
    using field_type = typename AGOP::Scalar;
    
    SymmetricTaylorHoodPreconditioner(AGOP const& K_)
    : K(&K_)
    , hatKinv(*K)
    {
    }
    
    void init(domain_type& x, range_type const& y) const 
    {
      hatKinv.apply(x,y,TaylorHoodPreconditionerFlavour::POSTSMOOTH);
    }
    
    virtual void apply(domain_type& x, range_type const& y) override
    {
      // apply presmoother
      auto z = x; 
      z = 0;
      hatKinv.apply(z,y,TaylorHoodPreconditionerFlavour::PRESMOOTH);
      
      // compute residual r = y - Kz
      auto r = y;
      K->applyscaleadd(-1.0,z,r);
      
      // apply postsmoother
      hatKinv.apply(x,r,TaylorHoodPreconditionerFlavour::POSTSMOOTH);
      x += z;
    }
    
  private:
    AGOP const* K;
    
    TaylorHoodPreconditioner<domain_type> hatKinv;
  };
  
}
 
#endif