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

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)

(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"?

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)
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)

## Question information

Language:
English Edit question
Status:
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Felix Ospald (felix-ospald) said on 2013-04-30: #1

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

 Revision history for this message Garth Wells (garth-wells) said on 2013-04-30: #2

On 23 April 2013 12:36, Marie Rognes
> New question #227292 on DOLFIN:
>
>
> 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)
> 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.

 Revision history for this message Anders Logg (logg) said on 2013-05-09: #3

consult the documentation on the FEniCS web page for where and
how to (re)post your question: http://fenicsproject.org/support/

 Revision history for this message Anders Logg (logg) said on 2013-05-09: #4