uzawa.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) 2002-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef UZAWA_H
#define UZAWA_H
 
 
#include "fem/linearspace.hh"
#include "utilities/timing.hh"
 
#include "dune/istl/preconditioners.hh"
#include "dune/istl/solvers.hh"
 
#include <boost/fusion/sequence.hpp>
 
namespace Kaskade
{
  template <class X, class Y>
  void applyPreconditioner(Dune::Preconditioner<X,Y>& p, X& x, Y& y) {
    p.pre(x,y);
    p.apply(x,y);
    p.post(x);
  }
 
  template<class X, class Y>
  class UzawaSolver : public Dune::InverseOperator<LinearProductSpace<typename X::field_type,
                                                                      boost::fusion::vector<X,Y>>,
                                                   LinearProductSpace<typename X::field_type,
                                                                      boost::fusion::vector<X,Y>>> 
  {
  public:
    using domain_type = LinearProductSpace<typename X::field_type,boost::fusion::vector<X,Y>>;
    using range_type = domain_type;
    typedef typename X::field_type field_type;
    typedef field_type Scalar;
    typedef domain_type Domain;
    typedef range_type Range;
 
    UzawaSolver (Dune::LinearOperator<X,X>& A_, Dune::InverseOperator<X,X>& invA_,
                 Dune::LinearOperator<X,Y>& B_, Dune::LinearOperator<Y,X>& Bt_,
                 Dune::Preconditioner<Y,Y>& precS_,
                 double reduction_, int maxit_, int verbose_)
    : A(A_)
    , invA(invA_)
    , B(B_)
    , Bt(Bt_)
    , precS(precS_)
    , reduction(reduction_)
    , maxit(maxit_)
    , verbose(verbose_)
    {
    }
 
 
    virtual void apply (Domain& x, Range& b, Dune::InverseOperatorResult& res) override
    {
      using namespace boost::fusion;
      using namespace Dune;
 
      ScopedTimingSection timer("Uzawa solve");
 
 
      Timer watch;                // start a timer
 
      InverseOperatorResult resA;
      if (verbose>0) 
        std::cout << "=== UzawaSolver\n";
 
      // Declare variables
      X& u(component<0>(x));        // "velocity" component and multiplier, both as references
      Y& p(component<1>(x));        // such that the provided solution vector is updated in place
      X const& f(component<0>(b));  // right hand side
      Y const& g(component<1>(b));  // and constraint
 
      // We apply preconditioned CG to the Schur complement, with p the computed solution.
      // Block elimination via A yields the Schur system Sp = BA^{-1}f - g with the
      // Schur komplement S = B A^{-1} B'.
      
      X xtmp(u), xtmp2(u);
      Y ytmp(p);
      
      // Compute initial energy p'Sp
      Bt.apply(p,xtmp);                         // xtmp = B'p
      {
        ScopedTimingSection solveA("solve A");
        auto Btp = xtmp;                        // need this copy since rhs arg of invA is non-const
        invA.apply(xtmp2,Btp,resA);             // xtmp2 = A^{-1}B'p
      }
      double energy = spX.dot(xtmp,xtmp2);      // maintain squared energy norm p^T S p
 
 
      // compute initial residual r = B A^{-1}(f-Au-B'p)-(g-Bu)
      Y r(p);
      xtmp = f;
      A.applyscaleadd(-1,u,xtmp);
      Bt.applyscaleadd(-1,p,xtmp);              // xtmp = f - B'p
      {
        ScopedTimingSection solveA("solve A");
        invA.apply(xtmp2,xtmp,resA);            // xtmp2 = A^{-1}(f-Au-B'p)
      }
      B.apply(xtmp2,r);                         // r = B A^{-1}(f-Au-B'p)
      r -= g;                                   // r = B A^{-1}(f-Au-B'p) - g
      B.applyscaleadd(1,u,r);                   // r = B A^{-1}(f-Au-B'p) - (g-Bu)
 
      u += xtmp2;                               // maintain u = A^{-1}(f-B'p) throughout
 
      Y rBar(r);
      {
        ScopedTimingSection preco("Schur preco");
        applyPreconditioner(precS,rBar,r);     // rBar = Shat^{-1} r
      }
      Y q = rBar;
      X z(u);                                  // loop invariant z = A^{-1} B' q
      Y Sq(p);                                 // loop invariant Sq = S q = B z
      
      double sigma = spY.dot(r,rBar);          // approximate Schur error energy r' Shat^{-1} r
      double sigma0 = sigma;                   // i.e. the preconditioned residual
 
 
      // CG iteration
      int i = 0;
      for ( ; i<maxit && sigma>reduction*reduction*energy; i++)
      {
        Bt.apply(q,xtmp);                         // xtmp = B'q
        {
          ScopedTimingSection solveA("solve A");
          invA.apply(z,xtmp,resA);                // z = A^{-1} xtmp = A^{-1} B' q
        }
        B.apply(z,Sq);                            // Sq = B z = S q
        
        double alpha = sigma / spY.dot(q,Sq);     // sigma / (q' S q)
 
        // update solution in background, maintaining u = A^{-1}(f-B'p)
        auto fu = std::async([&]()
        {
          p.axpy(alpha,q);                        // p = p + alpha q
          u.axpy(-alpha,z);                       // u = u - alpha z = u - alpha A^{-1}B'q
        });
 
        energy += alpha*sigma;
        r.axpy(-alpha,Sq);                        // r = r - alpha S q
        {
          ScopedTimingSection preco("Schur preco");
          applyPreconditioner(precS,rBar,r);      // rBar = Shat^{-1} r
        }
 
        fu.wait();                                // join solution update
 
        double sigmaNew = spY.dot(r,rBar);
        double beta = sigmaNew / sigma;
        sigma = sigmaNew;
        q *= beta; q += rBar;                     // q = rBar + beta q
 
        if (verbose>1) 
        {
          if (i%30==0)
            std::printf("%5s %14s %14s\n","Iter","Preco. Res.","Rate");
          std::printf("%5d %14.4e %14.4e\n",i,std::sqrt(sigma),std::sqrt(beta));
        }
      }
      
      // Fill statistics
      res.clear();
      res.iterations = i;
      res.reduction = std::sqrt(sigma/sigma0);
      res.converged = res.reduction<=reduction;
      res.conv_rate = std::pow(res.reduction,1.0/i);
      res.elapsed = watch.elapsed();
 
      if (verbose>0)                 // final print
        printf("=== rate=%g, time=%g, time/it=%g, iter=%d\n",res.conv_rate,res.elapsed,res.elapsed/i,i);
    }
 
    virtual void apply (Domain& x, Range& b, double reduction_, Dune::InverseOperatorResult& res)
    {
      reduction = reduction_;
      this->apply(x,b,res);
    }
 
    void apply(Domain& x, Range& b)
    {
      Dune::InverseOperatorResult tmpResult;
      apply(x,b,tmpResult);
    }
 
    virtual Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }
 
  private:
    Dune::SeqScalarProduct<X> spX;
    Dune::SeqScalarProduct<Y> spY;
 
    Dune::LinearOperator<X,X>& A;
    Dune::InverseOperator<X,X>& invA;
 
    Dune::LinearOperator<X,Y>& B;
    Dune::LinearOperator<Y,X>& Bt;
    Dune::Preconditioner<Y,Y>& precS;
 
    double reduction;
    int    maxit;
    int    verbose;
  };
  
  // ----------------------------------------------------------------------------------------------
  
  template <class MatA, class SolA, class MatB, class MatBt, class PrecS>
  auto uzawa(MatA& A, SolA& invA, MatB& B, MatBt& Bt, PrecS& precS, 
             double reduction=1e-3, int maxit=1000, int verbose=0)
  {
    using X = typename MatA::domain_type;
    using Y = typename MatB::range_type;
    return UzawaSolver<X,Y>(A,invA,B,Bt,precS,reduction,maxit,verbose);
  }
  
  template <class Scalar, class Sequence>
  auto nestUzawa(LinearProductSpace<Scalar,Sequence> const& z)
  {
    using namespace boost::fusion;
    using Z = LinearProductSpace<Scalar,Sequence>;
    using X = LinearProductSpace<Scalar,vector<std::remove_const_t<typename Z::Component<0>>>>;
    using Y = LinearProductSpace<Scalar,vector<std::remove_const_t<typename Z::Component<1>>>>;
    using Seq = vector<X,Y>;
    return LinearProductSpace<Scalar,Seq>(Seq(X(component<0>(z)),
                                              Y(component<1>(z))));
  }
  
  template <class Scalar, class Sequence>
  auto unnestUzawa(LinearProductSpace<Scalar,Sequence> const& z)
  {
    using namespace boost::fusion;
    using Z = LinearProductSpace<Scalar,Sequence>;
    using X = std::remove_const_t<typename std::remove_const_t<typename Z::Component<0>>::Component<0>>;
    using Y = std::remove_const_t<typename std::remove_const_t<typename Z::Component<1>>::Component<0>>;
    using Seq = boost::fusion::vector<X,Y>;
    return LinearProductSpace<Scalar,Seq>(Seq(component<0>(component<0>(z)),
                                              component<0>(component<1>(z))));
  }
  
}
 
#endif