limexWithoutJens.hh

Go to the documentation of this file.

#ifndef LIMEX_HH
#define LIMEX_HH
 
#include <complex>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <cstdlib>
 
#include <boost/timer/timer.hpp>
 
#include "dune/istl/solvers.hh"
#include "dune/istl/preconditioners.hh"
 
#include "fem/embedded_errorest.hh"
#include "fem/iterate_grid.hh"
#include "timestepping/extrapolation.hh"
#include "timestepping/semieuler.hh"
 
#include "fem/assemble.hh"
#include "fem/istlinterface.hh"
#include "fem/functional_aux.hh"
#include "fem/hierarchicspace.hh"
#include "fem/lagrangespace.hh"
 
#include "linalg/factorization.hh"
#include "linalg/umfpack_solve.hh"
#include "linalg/mumps_solve.hh"
// #include "linalg/superlu_solve.hh"
 
#include "linalg/trivialpreconditioner.hh"
#include "linalg/direct.hh"
#include "linalg/triplet.hh"
#include "linalg/iluprecond.hh"
#include "linalg/hyprecond.hh"
#include "linalg/jacobiPreconditioner.hh"
 
#include "utilities/enums.hh"
 
DirectType xyz = DirectType::UMFPACK;
 
namespace Kaskade
{
  template <class Eq>
  class Limex
  {
  public:
    typedef Eq                                                  EvolutionEquation;
    typedef typename EvolutionEquation::AnsatzVars::VariableSet State;
 
  private:
    typedef SemiLinearizationAt<SemiImplicitEulerStep<EvolutionEquation>> Linearization;
    typedef VariationalFunctionalAssembler<Linearization> Assembler;
 
  public:
    Limex(GridManager<typename EvolutionEquation::AnsatzVars::Grid>& gridManager_,
          EvolutionEquation& eq_,
          typename EvolutionEquation::AnsatzVars const& ansatzVars_,
          DirectType st=DirectType::MUMPS, PrecondType precondType_ = PrecondType::ILUK,
          int verbosity_ = 1)
    : gridManager(gridManager_), ansatzVars(ansatzVars_), eq(&eq_,0)
    , assembler(ansatzVars.spaces), extrap(0),
          rhsAssemblyTime(0.0), matrixAssemblyTime(0.0), precAssemblyTime(0.0), factorizationTime(0.0), solutionTime(0.0),
          initSolutionTime(0.0), estimateTime(0), directType(st), precondType(precondType_), verbosity(verbosity_)
    {}
 
 
    State const& step(State const& x, double dt, int order,
                      std::vector<std::pair<double,double>> const& tolX)
    {
      boost::timer::cpu_timer timer;
 
      std::vector<double> stepFractions(order+1);
      for (int i=0; i<=order; ++i) 
        stepFractions[i] = 1.0/(i+1); // harmonic sequence for extrapolation
      extrap.clear();
 
      int const nvars = EvolutionEquation::AnsatzVars::noOfVariables;
      int const neq = EvolutionEquation::TestVars::noOfVariables;
      size_t  nnz = assembler.nnz(0,neq,0,nvars,false);
 
      MatrixAsTriplet<double> triplet(nnz);
 
      typedef typename EvolutionEquation::AnsatzVars::template CoefficientVectorRepresentation<0,neq>::type CoefficientVectors;
 
      typedef  AssembledGalerkinOperator<Assembler,0,neq,0,neq> Op;
      State dx(x), dxsum(x), tmp(x);
      double const t = eq.parabolicEquation().time();
 
      eq.parabolicEquation().temporalEvaluationRange(t,t+dt);
 
      bool const iterative = false;      // false: use a direct solver,      true: use an iterative solver
 
 
      // Loop over extrapolation orders. For each order, an Euler solution is computed
      // on an equidistant time grid with shorter and shorter time sub-steps. At the
      // joint final time, the value is extrapolated.
      for (int i=0; i<=order; ++i)
      {
        double const tau = stepFractions[i]*dt;
        eq.setTau(tau);
        eq.parabolicEquation().setTime(t);
 
        // mesh adaptation loop. Just once if i>0.
        gridManager.setVerbosity(verbosity);
        bool accurate = true;
        do {
          int dimension = gridManager.grid().dimension;
 
          // Evaluate and factorize matrix B(t)-tau*J
          dx *= 0;
          timer.start();
          assembler.assemble(Linearization(eq,x,x,dx));
          matrixAssemblyTime += (double)(timer.elapsed().user)/1e9;
          Op A(assembler);
 
          timer.start();
          CoefficientVectors rhs(assembler.rhs());
          CoefficientVectors solution = ansatzVars.zeroCoefficientVector();
 
          if (iterative)
          {
            typedef typename EvolutionEquation::AnsatzVars::template CoefficientVectorRepresentation<0,neq>::type LinearSpaceX;
 
            if ( verbosity>0 )
            {
              std::cout << " eq 0: dof = " << boost::fusion::at_c<0>(dx.data).space().degreesOfFreedom()
                              << " pts in mesh = " << gridManager.grid().size(dimension) << std::endl;
              if ( verbosity>1 )
                std::cout << " dimension of grid = " << dimension << std::endl;
            };
            
            initSolutionTime += (double)(timer.elapsed().user)/1e9;
            elapsedTimeSinceReset = (double)(timer.elapsed().user)/1e9;
 
            // in der Folge finden sich vesrschiedene iterative L�ser
            // und Vorkonditionierer zur Auswahl per Ein-/Auskommentieren.
            // Der Mechanismus muss noch umgestaltet werden mittels
            // passender Parameter aus dem Main()-Programm via itegrate.hh !!!!
 
            int iteSteps = 5000;
            double iteEps = 1.0e-10;
            Dune::InverseOperatorResult res;
 
            switch (precondType)
            {
            case PrecondType::NONE:
            {
              TrivialPreconditioner<Op> trivial;
              precAssemblyTime += (double)(timer.elapsed().user)/1e9-elapsedTimeSinceReset;
              Dune::BiCGSTABSolver<LinearSpaceX> cg(A,trivial,iteEps,iteSteps,verbosity-1);
              cg.apply(solution,rhs,res);
            }
            break;
            case PrecondType::BOOMERAMG:
            {
              int steps = iteSteps;
              int coarsentype = 21;
              int interpoltype = 0;
              int cycleType = 1;
              int relaxType = 3;
              int variant = 0;
              int overlap = 1;
              int syseqn = 1;
              double tol = iteEps;
              double strongThreshold = (dimension==2)?0.25:0.6;
              BoomerAMG<Op> BoomerAMGPrecon(A,steps,coarsentype,interpoltype,tol,cycleType,relaxType, 
                            strongThreshold,variant,overlap,syseqn,verbosity);
              precAssemblyTime += (double)(timer.elapsed().user)/1e9-elapsedTimeSinceReset;
              Dune::BiCGSTABSolver<LinearSpaceX> cg(A,BoomerAMGPrecon,iteEps,iteSteps,verbosity-1);       // 1e-10
              cg.apply(solution,rhs,res);
            }
            break;
            case PrecondType::ILUK:
            {
              int fill_lev = 1;
              ILUKPreconditioner<Op> iluk(A,fill_lev,verbosity);
              precAssemblyTime += (double)(timer.elapsed().user)/1e9-elapsedTimeSinceReset;
              //Dune::CGSolver<LinearSpace> cg(A,iluk,iteEps,iteSteps,verbosity-1);
              Dune::BiCGSTABSolver<LinearSpaceX> cg(A,iluk,iteEps,iteSteps,verbosity-1);
              cg.apply(solution,rhs,res);
            }
            break;
            case PrecondType::JACOBI:
            default:
            {
              JacobiPreconditioner<Op> jacobi(A,1.0);
              precAssemblyTime += (double)(timer.elapsed().user)/1e9-elapsedTimeSinceReset;
              Dune::BiCGSTABSolver<LinearSpaceX> cg(A,jacobi,iteEps,iteSteps,verbosity-1);
              // alternative:  Dune::LoopSolver<LinearSpaceX> cg(A,jacobi,iteEps,iteSteps,verbosity-1);
              cg.apply(solution,rhs,res);
            }
            break;
            }
 
            if ( !(res.converged) || (res.iterations == 2001) ) {
              std::cout << "   no of iterations in cg = " << res.iterations << std::endl;
              std::cout << " convergence status of cg = " << res.converged  << std::endl;
              //std::cout << "A: "  << std::endl;
              //for (size_t ii=0; ii<nnz; ii++)
              //    std::cout   << A.getmat().ridx[ii] << ", "
              //                << A.getmat().cidx[ii] << ", "
              //                << A.getmat().data[ii] << std::endl;
              //std::cout << "A, end: "  << std::endl;
              assert(0);
            }
          }
          else
          {
            directInverseOperator(A,directType).applyscaleadd(1.0,rhs,solution);
            factorizationTime += (double)(timer.elapsed().user)/1e9;
          }
          dx.data = solution.data;
          solutionTime += (double)(timer.elapsed().user)/1e9;
 
 
          timer.start();
          // Mesh adaptation only if requested and only for the full
          // implicit Euler step (first stage).
          if (!tolX.empty() && i==0) 
          {
            // Compute spatial error estimate
            tmp = dx;
            projectHierarchically(tmp);
            tmp -= dx;
 
            // perform mesh adaptation
            accurate = embeddedErrorEstimator(ansatzVars,tmp,x,eq.parabolicEquation().scaling(),tolX,gridManager,verbosity);
            if (!accurate)
              nnz = assembler.nnz(0,neq,0,nvars,false);
            estimateTime += (double)(timer.elapsed().user)/1e9;
          }
        } while (!accurate);
 
        dxsum = dx;
 
        // propagate by linearly implicit Euler
        if (order==0) 
        {
          // insert into extrapolation tableau
          extrap.push_back(dxsum,stepFractions[i]);
 
          // restore initial time
          eq.parabolicEquation().setTime(t);
          if (verbosity>1)
            std::cout << "   end of limex step" << std::endl;
          return extrap.back();
        }
        else
          for (int j=1; j<=i; ++j)
          {
            // Assemble new right hand side tau*f(x_j)
            eq.parabolicEquation().setTime(eq.parabolicEquation().time()+tau);
            tmp = x; tmp += dxsum;
            State zero(x); zero *= 0;
            timer.start();
            assembler.assemble(Linearization(eq,tmp,x,zero),Assembler::RHS|Assembler::MATRIX);
            rhsAssemblyTime += (double)(timer.elapsed().user)/1e9;
 
            Op A(assembler);
 
            timer.start();
            CoefficientVectors rhs(assembler.rhs());
            CoefficientVectors solution(EvolutionEquation::AnsatzVars::template CoefficientVectorRepresentation<0,neq>::init(ansatzVars.spaces));
            solution = 0;
 
            if (iterative)
            {
              typedef typename EvolutionEquation::AnsatzVars::template CoefficientVectorRepresentation<0,neq>::type LinearSpaceX;
 
              initSolutionTime += (double)(timer.elapsed().user)/1e9;
              elapsedTimeSinceReset = (double)(timer.elapsed().user)/1e9;
 
 
              int iteSteps = 5000;
              double iteEps = 1.0e-10;
              Dune::InverseOperatorResult res;
 
              switch (precondType)
              {
                case PrecondType::NONE:
                {
                  TrivialPreconditioner<Op> trivial;
                  precAssemblyTime += (double)(timer.elapsed().user)/1e9-elapsedTimeSinceReset;
                  Dune::BiCGSTABSolver<LinearSpaceX> cg(A,trivial,iteEps,iteSteps,verbosity-1);
                  cg.apply(solution,rhs,res);
                }
                break;
                case PrecondType::BOOMERAMG:
                {
                  int steps = iteSteps;
                  int coarsentype = 21;
                  int interpoltype = 0;
                  int cycleType = 1;
                  int relaxType = 3;
                  int variant = 0;
                  int overlap = 1;
                  int syseqn = 1;
                  double tol = iteEps;
                  int dimension = gridManager.grid().dimension;
                  double strongThreshold = (dimension==2)?0.25:0.6;
                  BoomerAMG<Op> BoomerAMGPrecon(A,steps,coarsentype,interpoltype,tol,cycleType,relaxType,
                                                strongThreshold,variant,overlap,syseqn,verbosity);
                  precAssemblyTime += (double)(timer.elapsed().user)/1e9-elapsedTimeSinceReset;
                  Dune::BiCGSTABSolver<LinearSpaceX> cg(A,BoomerAMGPrecon,iteEps,iteSteps,verbosity-1);       // 1e-10
                  cg.apply(solution,rhs,res);
                }
                break;
                case PrecondType::ILUK:
                {
                  int fill_lev = 1;
                  ILUKPreconditioner<Op> iluk(A,fill_lev,verbosity);
                  precAssemblyTime += (double)(timer.elapsed().user)/1e9-elapsedTimeSinceReset;
                  //Dune::CGSolver<LinearSpace> cg(A,iluk,iteEps,iteSteps,verbosity-1);
                  Dune::BiCGSTABSolver<LinearSpaceX> cg(A,iluk,iteEps,iteSteps,verbosity-1);
                  cg.apply(solution,rhs,res);
                }
                break;
                case PrecondType::JACOBI:
                default:
                {
                  JacobiPreconditioner<Op> jacobi(A,1.0);
                  precAssemblyTime += (double)(timer.elapsed().user)/1e9-elapsedTimeSinceReset;
                  Dune::BiCGSTABSolver<LinearSpaceX> cg(A,jacobi,iteEps,iteSteps,verbosity-1);
                  cg.apply(solution,rhs,res);
                }
                break;
              }
 
              if ( !(res.converged) || (res.iterations == 2001) )
              {
                typedef typename Op::field_type field_type;
                MatrixAsTriplet<field_type> triplet = A.template get<MatrixAsTriplet<field_type> >();
 
                std::cout << "   no of iterations in cg = " << res.iterations << std::endl;
                std::cout << " convergence status of cg = " << res.converged  << std::endl;
                std::cout << "A: "  << std::endl;
                for (size_t ii=0; ii<nnz; ii++)
                  std::cout   << triplet.ridx[ii] << ", " << triplet.cidx[ii] << ", " << triplet.data[ii] << std::endl;
                std::cout << "A, end: "  << std::endl;
                assert(0);
              }
            }
            else
            {
              directInverseOperator(A,directType).applyscaleadd(1.0,rhs,solution);
            }
            dx.data = solution.data;
 
            dxsum += dx;
            solutionTime += (double)(timer.elapsed().user)/1e9;
          }
        // insert into extrapolation tableau
        extrap.push_back(dxsum,stepFractions[i]);
 
        // restore initial time
        eq.parabolicEquation().setTime(t);
      }
 
      return extrap.back();
    }
 
    std::vector<std::pair<double,double> > estimateError(State const& x,int i, int j) const
      {
      assert(extrap.size()>1);
 
      std::vector<std::pair<double,double> > e(ansatzVars.noOfVariables);
 
      relativeError(typename EvolutionEquation::AnsatzVars::Variables(),extrap[i].data,
          extrap[j].data,x.data,
          ansatzVars.spaces,eq.parabolicEquation().scaling(),e.begin());
 
      return e;
      }
 
    template <class OutStream>
    void reportTime(OutStream& out) const {
      out << "Limex time: " << matrixAssemblyTime << "s matrix assembly\n"
          << "            " << rhsAssemblyTime << "s rhs assembly\n"
          << "            " << precAssemblyTime << "s preconditioner assembly\n"
          << "            " << factorizationTime << "s factorization\n"
          << "            " << initSolutionTime << "s init solution\n"
          << "            " << solutionTime << "s solution\n"
          << "            " << estimateTime << "s estimate\n"
          << std::flush;
    }
 
    void advanceTime(double dt)
    {
      eq.parabolicEquation().setTime(eq.time()+dt);
    }
 
 
  private:
    GridManager<typename EvolutionEquation::AnsatzVars::Grid>& gridManager;
    typename EvolutionEquation::AnsatzVars const& ansatzVars;
    SemiImplicitEulerStep<EvolutionEquation>      eq;
    Assembler                                     assembler;
 
  public:
    ExtrapolationTableau<State>  extrap;
    double rhsAssemblyTime, matrixAssemblyTime, precAssemblyTime, factorizationTime, solutionTime;
    double elapsedTimeSinceReset, initSolutionTime, estimateTime;
    DirectType directType;
    PrecondType precondType;
    int verbosity;
  };
} // namespace Kaskade
#endif