[Dune] Performace of ISTL

[Patrick Leidenberger]

Hi all,

because of the very bad convergence of istl on a very specific FE 
problem I try to solve my linear system with Trilinos. I use a few 
#ifdef's in my Dune code and replace only the call of the istl solver 
with the Trilinos/actecOO solver (the code snippets are below).

For testing I use both different codes on the same problem (~155000 
degrees of freedom), serial:
Dune solver: 57sec, 81 iterations, residual: 1.5402E-08
Trilinos   : 10sec,  6 iterations, residual: 2.8799e-08
             (including copy from dune vector to trilinos vector and back)

Perhaps I use the dune solvers in a wrong way, but I tried different 
settings. Is there a way to get the same performance with dune?

And there is an other question:
Is there no easy way, to get a global, consecutive integer id in a 
parallel DUNE code? I know, there is the global id, which I can compare 
and sort, but if I have such an integer index I can feed my problem to a 
parallel Trilinos solver in a very easy way. An other application for 
such an index will be, to write the solution from different processes in 
the same hdf5 file.
I have a look at the global id and I saw, that it consists of 4 
integers. Is any functionality to send such a global id object via MPI? 
Then can I build my own integer index.

Best regards

Patrick


#ifdef TRILINOS_SOLVER
  #ifdef HAVE_TRILINOS
   
    int vecSize = yi[1].size();
    Epetra_Map map(vecSize, 0, Comm_);

    /** Create trilinos vector. */
    Epetra_Vector yiTrilinos_(map);
    Epetra_Vector xiTrilinos_(map);

    /** Copy from dune to trilinos vector. */
    for (unsigned int vecIter = 0;  vecIter < vecSize; ++vecIter)
    {
      yiTrilinos_[vecIter] = yi[1][vecIter][0];
      xiTrilinos_[vecIter] =  EmpClass::dof1_[1][vecIter][0];
    }
    Epetra_LinearProblem problem(&(*EmpClass::aijTrilinos_), 
&xiTrilinos_, &yiTrilinos_);
    AztecOO solver(problem);
    int options[AZ_OPTIONS_SIZE];
    double params[AZ_PARAMS_SIZE];
    AZ_defaults(options, params);

    /** Solve the linear system. */
    
solver.Iterate(EmpClass::solverMaximumIteration_,EmpClass::solverReduction_);

    /** Copy the solution vector from trilinos to a dune vector. */
    for (unsigned int vecIter = 0;  vecIter < vecSize; ++vecIter)
    {
      EmpClass::dofp1_[1][vecIter][0] = xiTrilinos_[vecIter];
    }
  #endif // HAVE_TRILINOS
#else // TRILINOS_SOLVER

  /**  Initialize defect. */
  EmpVectorTopo defect = yi;

  typedef Dune::SeqJac<EmpMatrixTopo,EmpVectorTopo,
                       EmpVectorTopo,HIERARCHY_LEVEL> Preconditioner;
 
  Preconditioner pr(EmpClass::aij_,1,1);
 
  typedef Dune::MatrixAdapter<EmpMatrixTopo,EmpVectorTopo,EmpVectorTopo> 
Operator;
  Operator op(EmpClass::aij_);
 
  /** \brief Initialize the solver. */
  typedef Dune::GradientSolver<EmpVectorTopo> Solver;
  Solver solver(op, pr,
                EmpClass::solverReduction_,
                EmpClass::solverMaximumIteration_,
                EmpClass::solverVerbose_);
                     
  Dune::InverseOperatorResult r;
 
  /** Apply the solver. */
  solver.apply(EmpClass::dofp1_,yi,r);
#endif // TRILINOS_SOLVER