[dune-pdelab] Dune::PDELab::ISTLBackend_NOVLP_CG_AMG_SSOR with ALUGrid

[Eike Mueller]

Hi Rebecca,

Rebecca Neumann wrote:
> Hey Eike,
> 
> with the P0GhostConstraints I used ISTLBackend_BCGS_AMG_SSOR and 
> ISTLBackend_OVLP_BCGS_SSORk, this works for YaspGrid as well as ALU.
> Why should the normal solver not work for FV?
> 

This is what I do to solve with the FV scheme:

/* *** Construct finite element space P0 *** */
     typedef Dune::PDELab::P0LocalFiniteElementMap<Coord, Real, dim> FEM;
     FEM fem(Dune::GeometryType(Dune::GeometryType::cube,dim));

/* *** Construct grid function space *** */
     typedef Dune::PDELab::OverlappingConformingDirichletConstraints CON;
     typedef Dune::PDELab::ISTLVectorBackend<1> VBE;
     typedef Dune::PDELab::GridFunctionSpace<GV, FEM, CON, VBE> GFS;
     GFS gfs(gv,fem);

/* *** Construct grid operator space *** */
     typedef HelmholtzLocalOperator<aReference,bReference,fRHS> LOP;
     LOP lop(aref,
             bref,
             frhs,
             model_parameters);
     typedef Dune::PDELab::ISTLBCRSMatrixBackend<1,1> MBE;
     typedef Dune::PDELab::EmptyTransformation CC;
     typedef Dune::PDELab::GridOperatorSpace<GFS, GFS, LOP, CC, CC, MBE> 
GOS;
     GOS gos(gfs,gfs,lop);

/* *** Initialise solution vector *** */
     typedef typename GFS::template VectorContainer<Real>::Type U;
     U u(gfs, 0.0);

/* *** Construct constraints *** */
     typedef BCType<GV> B;
     B b(gv,model_parameters);
     typedef typename GFS::template ConstraintsContainer<Real>::Type C;
     C c;
     Dune::PDELab::constraints(b,gfs,c);

/* *** Construct linear solver *** */
     typedef Dune::PDELab::ISTLBackend_OVLP_CG_SSORk<GFS,C> LS;
     LS ls(gfs,c,maxiter,smoothsteps, verbose);

/* *** Solve linear problem *** */
     Dune::PDELab::StationaryLinearProblemSolver<GOS,LS,U>
       slp(gos,u,ls, model_parameters.solver.resreduction);
     slp.apply();
     result = ls.result();

it works on one core, but not on four, the error is largest at the point 
where all four domains meet, in the centre of the cube. Maybe I'm not 
setting up the constraints correctly? I should say that this is for DUNE 
version 2.1, not for the trunk.

> With the normal OverlappingConformingDirichletConstraints the code runs 
> as well, but the solution at the dirichlet points is wrong.
> Does that really give you the right solutions?
> 
I've just tested it again, and the numerical solution for a test problem 
in a unit cube, implemented with ALUGrid, agrees very well with the 
analytical solution, in particular I can't see large errors at the 
processor domain boundaries.

Cheers,

Eike

> Best Regards,
> Rebecca
>