Why does KrylovSolver with amg diverge/converge depending on dof ordering?

Hi, I can provide you with an example how to use AMG with elasticity.
Here the nullspace corresponds to rigid body motions.

class NullspaceMode : public Expression
{
public:
        NullspaceMode(int i) : Expression(3), _i(i) {}

        void eval(Array<double>& values, const Array<double>& x) const
        {
                values[0] = 0.0;
                values[1] = 0.0;
                values[2] = 0.0;

                switch (_i) {
                        case 0: values[0] = 1.0; break; // tx
                        case 1: values[1] = 1.0; break; // ty
                        case 2: values[2] = 1.0; break; // tz
                        case 3: values[1] = -x[2]; values[2] = x[1]; break; // rx
                        case 4: values[0] = -x[2]; values[2] = x[0]; break; // ry
                        case 5: values[0] = -x[1]; values[1] = x[0]; break; // rz
                }
        }

        static PETScVector* generate(int i, const FunctionSpace& V)
        {
                PETScVector* c = new PETScVector(V.dim());
                V.interpolate(*c, NullspaceMode(i));
                return c;
        }

        int _i;
};

...

PETScPreconditioner p("ml_amg");

std::vector<const GenericVector*> nullspace;
for (int i = 0; i < 6; i++) {
 nullspace.push_back(NullspaceMode::generate(i, *(problem->trial_space())));
}
p.set_nullspace(nullspace);
...

PETScKrylovSolver solver("cg", p);
...

Afaik AMG (Hypre and ML) uses by default SSOR for smoothing, which is hidden somewhere in the PETSc code. The results of SSOR are always dependent on the dof ordering, but I'm not sure if it can be that worse.
You can try to adjust the settings (for Hypre AMG) by using this code (or doing some more advanced PETSc coding):

// set some options for Hypre AMG
// see http://www.gnu-darwin.org/www001/ports-1.5a-CURRENT/math/petsc/work/petsc-2.3.3-p0/src/ksp/pc/impls/hypre/hypre.c.html
// for more possibilities
//
// CycleTypes = {"","V","W"};
// CoarsenTypes = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
// MeasureTypes = {"local","global"};
// RelaxTypes = {"Jacobi","sequential-Gauss-Seidel","","SOR/Jacobi","backward-SOR/Jacobi","","symmetric-SOR/Jacobi","","","Gaussian-elimination"};
// InterpTypes = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i", "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};

const char* params[] = {argv[0],
 // "--petsc.pc_hypre_boomeramg_relax_type_up", "SOR/Jacobi",
 // "--petsc.pc_hypre_boomeramg_relax_type_down", "backward-SOR/Jacobi",
 // "--petsc.pc_hypre_boomeramg_max_iter", "1",
 // "--petsc.pc_hypre_boomeramg_rtol", "0.0",
 // "--petsc.pc_hypre_boomeramg_strong_threshold", "0.0",
 // "--petsc.pc_hypre_boomeramg_coarsen_type", "CLJP",
 // "--petsc.pc_hypre_boomeramg_max_levels", "10",
 // "--petsc.pc_hypre_boomeramg_agg_nl", "1",
 // "--petsc.pc_hypre_boomeramg_relax_type_all", "SOR",
 // "--petsc.pc_hypre_boomeramg_print_statistics", "1",
 // "--petsc.pc_hypre_boomeramg_interp_type", "standard",
};
parameters.parse(sizeof(params)/sizeof(char *), const_cast<char**>(params));

Felix

On 23 April 2013 12:36, Marie Rognes
<email address hidden> wrote:
> New question #227292 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/227292
>
>
> I'm trying to set a null space for a PETSc KrylovSolver using "amg" preconditioning, but am running into some issues I don't understand. Sample code below. If dof reordering is turned off, the solver converges quickly in 7 iterations. If dof reordering is turned on, the solver diverges (and in particular fails to converge in 900 iterations)
>

It works by chance with default FFC ordering since it separate the u and p dofs.

> (i) Is this expected, or am I using set_nullspace incorrectly, or ... ?
>
> (ii) Is this related to https://bugs.launchpad.net/dolfin/+bug/1092262 ? If so, is there a way of "attaching the nullspace to the preconditioner"?
>

The null space is required only for smoothed aggregation AMG (e.g. ML
and PETSc GAMG). Hypre does not require the null space.

Making mixed problems work properly with AMG requires support for
block preconditioning. It's being worked on.

Garth

> Any ideas much appreciated.
>
> --
> parameters["reorder_dofs_serial"] = True # Diverges
> #parameters["reorder_dofs_serial"] = False # Converges
>
> # Set-up variational problem
> N = 3
> mesh = UnitCubeMesh(N, N, N)
> V = FunctionSpace(mesh, "CG", 1)
> VU = V*V
>
> (v, u) = TrialFunctions(VU)
> (w, q) = TestFunctions(VU)
> a = (v*w/0.1 + inner(grad(v), grad(w)) + inner(grad(u), grad(w))
> + inner(grad(v), grad(q)) + inner(3.0*grad(u), grad(q)))*dx
> L = 0.1*q*dx
>
> # Assemble matrix/vector
> A = assemble(a)
> b = assemble(L)
>
> # Set-up nullspace vector
> vu = Function(VU)
> nullspace = Vector(vu.vector())
> VU.sub(1).dofmap().set(nullspace, 1.0)
>
> # Set-up solver
> solver = KrylovSolver("gmres", "amg")
> #solver = KrylovSolver("gmres", "jacobi") # Robust, but more iterations...
> solver.set_operator(A)
> solver.set_nullspace([nullspace])
> solver.parameters.monitor_convergence = True
>
> # Solve
> solver.solve(vu.vector(), b)
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

FEniCS no longer uses Launchapd for Questions & Answers. Please
consult the documentation on the FEniCS web page for where and
how to (re)post your question: http://fenicsproject.org/support/

FEniCS no longer uses Launchapd for Questions & Answers. Please
consult the documentation on the FEniCS web page for where and
how to (re)post your question: http://fenicsproject.org/support/