[dune-pdelab] Dune::PDELab::ISTLBackend_NOVLP_CG_AMG_SSOR with ALUGrid
Eike Mueller
E.Mueller at bath.ac.uk
Mon Mar 26 17:23:31 CEST 2012
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
>
More information about the dune-pdelab
mailing list