preconditioners for complex-valued problems

There is only rudimentary support for block matrices in the Python
interface of DOLFIN. You also have the possibility to use cbc.block. Look
at launchpad, but again only in Python and it requires PyTrilinos.

Johan
On Oct 10, 2012 1:01 PM, "Nico Schlömer" <
<email address hidden>> wrote:

> New question #210832 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/210832
>
> When solving an ordinary real-valued curl-curl problem, some AMG
> preconditioners do an excellent job.
>
> ===================== *snip* =====================
> from dolfin import *
>
> parameters['linear_algebra_backend'] = 'Epetra'
>
> mesh = UnitCube(30, 30, 30)
> W = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
>
> bc = DirichletBC(W, (0, 0, 0), "on_boundary")
>
> rhs = project(Constant([0.0, 0.0, 1.0]), W)
>
> T = Function(W)
> v0 = TestFunction(W)
> u0 = TrialFunction(W)
> a = inner(curl(u0), curl(v0))*dx
> L = -inner(rhs,v0)*dx
> prm = parameters['krylov_solver']
> prm['absolute_tolerance'] = 1E-10
> prm['relative_tolerance'] = 1E-6
> prm['maximum_iterations'] = 10000
> set_log_level(PROGRESS)
> solve(a == L, T, bc,
> solver_parameters={'linear_solver': 'gmres',
> 'preconditioner': 'ml_amg'})
> ===================== *snap* =====================
>
> (Convergence results could be improved even further by setting the default
> parameters to 'Maxwell' instead of 'SA' in
> ./dolfin/la/TrilinosPreconditioner.cpp, but it doesn't seem to be possible
> to modify this at the moment.)
>
> If the vector field A is complex-valued, it is advised to break up A into
> real and imaginary part and do the computation that way.
>
> ===================== *snip* =====================
> from dolfin import *
>
> parameters['linear_algebra_backend'] = 'Epetra'
>
> # Create mesh
> mesh = UnitCube(3, 3, 3)
>
> PN = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
> W = PN * PN
>
> bc = DirichletBC(W, (0, 0, 0, 0, 0, 0), "on_boundary")
>
> fr = project(Constant([0.0, 0.0, 1.0]), PN)
> fi = project(Constant([0.0, 0.0, 1.0]), PN)
>
> T = Function(W)
> vr, vi = TestFunctions(W)
> ur, ui = TrialFunctions(W)
>
> a = inner(curl(ur), curl(vr))*dx + inner(curl(ui), curl(vi))*dx \
> + inner(curl(ur), curl(vi))*dx - inner(curl(ui), curl(vr))*dx
> L = -inner(fr,vr)*dx + inner(fi, vi)*dx \
> -inner(fr,vi)*dx - inner(fi, vr)*dx
> prm = parameters['krylov_solver'] # short form
> prm['absolute_tolerance'] = 1E-10
> prm['relative_tolerance'] = 1E-6
> prm['maximum_iterations'] = 20000
> set_log_level(PROGRESS)
> solve(a == L, T, bc,
> solver_parameters={'linear_solver': 'gmres',
> 'preconditioner': 'ml_amg'})
> #'preconditioner': 'none'})
> ===================== *snap* =====================
>
> However, the preconditioner doesn't seem to work at all.
> A possible reason for this is the fact that the aggregation step in AMG
> mixes up real and imaginary components since at the time of the matrix
> build it has no way of separating those. ML does provide an option for this
> particular case, namely
>
> mlist.set( "PDE equations", 2 );
>
> This one however only works if is the matrix is really a 2x2 block matrix
> of its real and imaginary parts.
> At this point, it becomes important just how 'a' is assembled.
>
> I anyone aware of other possible causes for this problem?
> What is the order in which trial and test functions are called to build
> the matrix?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I wouldn't mind Python-only libraries, but adding just another layer of software is something I'd rather avoid when the actual solution is so much simpler (namely support for complex-valued, sheesh).

I realize that adding support for complex-valued problem would is a major task. One of the principal obstacles I can see is that the workhorse backends (Epetra, PetSC,...) don't support complex.

Anyways. If I was going to think about spending some time here, which of the FEniCS packages would I have to look at?

On Wed, Oct 10, 2012 at 7:11 PM, Nico Schlömer
<email address hidden> wrote:
> Question #210832 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210832
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
> I wouldn't mind Python-only libraries, but adding just another layer of
> software is something I'd rather avoid when the actual solution is so
> much simpler (namely support for complex-valued, sheesh).
>
> I realize that adding support for complex-valued problem would is a
> major task. One of the principal obstacles I can see is that the
> workhorse backends (Epetra, PetSC,...) don't support complex.
>
> Anyways. If I was going to think about spending some time here, which of
> the FEniCS packages would I have to look at?
>

We're working on support for block preconditioning via the C++
interface, and contributions in this direction would be welcome. I
would suggest looking at PETSc FieldSplit and how this could be
integrated nicely into the DOLFIN abstractions.

Garth

> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Adding complex number support would require to look at all FEniCS
packages. Not sure this is the path you want to take...

Block systems with block preconditioners are pretty common and to me it
seems like this is really what you need.

Johan

On 10/10/2012 08:11 PM, Nico Schlömer wrote:
> Question #210832 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210832
>
> Status: Answered => Open
>
> Nico Schlömer is still having a problem:
> I wouldn't mind Python-only libraries, but adding just another layer of
> software is something I'd rather avoid when the actual solution is so
> much simpler (namely support for complex-valued, sheesh).
>
> I realize that adding support for complex-valued problem would is a
> major task. One of the principal obstacles I can see is that the
> workhorse backends (Epetra, PetSC,...) don't support complex.
>
> Anyways. If I was going to think about spending some time here, which of
> the FEniCS packages would I have to look at?
>

@Garth
I'm coming more from the Trilinos-side of things, and I'm aware of http://trilinos.sandia.gov/packages/docs/dev/packages/epetra/doc/html/classEpetra__FEVbrMatrix.html as a block matrix implemenation. It's only efficient for a block size larger than five or ten or so, so the intention is really for larger multiphysics problems.
Otherwise one can use ordinary Epetra_(FE)CrsMatrices, having in mind that some routines might need additional hints about the matrix structure (see my first post).

@Johan
Complex-valued systems are also pretty common, at least as far as someone's concerned who's had the pleasure to work in quantum dynamics (Schrödinger!). Also, everything that happens in the frequency domain is complex-valued (Maxwell!).
Anyways, I've never seen a consistent implementation of anything that incorporates complex values, and although FEniCS is the best software I've seen so far in many respects, I'm not going to take it out on her.
I'm sure block systems are the lower fruit, especially since there are some people interested in this.

I'll look at thing with my block matrix hat on.

So I would say there are two different notions of "block matrices" here:

   * Instead of hosting scalars, the entries of a matrix are more or less small k*k blocks. All of the data is held in one underlying data structure, an ordinary nk*nk matrix if you want. This format doesn't expose much any of its structure, and in particular it is hardly possible to treat the blocks separately.

   * Creating a new data structure that can be fed with, e.g., four linear operators A, B, C, D to be composed to something like
     [ A, B ]
     [ C, D ]
    The entities A, B, C, D could still be treated separately which would make it possible, for exampe, to calculate the Schur complement w.r.t. A or D. The Trilinos package Teko has gone some lengths here and may be worth looking at (design-wise, too).

I believe what you are talking about is no. 2 when you mention things like block preconditioners.
What I had in mind for the faked complex-valued problem was no. 1 though. One could certainly formulate a problem with such tightly coupled parameters such as Re(u) and Im(u) in the form [A B; C D], but only think about ILU and fill ins for a k*k block-diagonal matrix.

On Thu, Oct 11, 2012 at 12:15 AM, Nico Schlömer
<email address hidden> wrote:
> Question #210832 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/210832
>
> Nico Schlömer posted a new comment:
> So I would say there are two different notions of "block matrices" here:
>
> * Instead of hosting scalars, the entries of a matrix are more or
> less small k*k blocks. All of the data is held in one underlying data
> structure, an ordinary nk*nk matrix if you want. This format doesn't
> expose much any of its structure, and in particular it is hardly
> possible to treat the blocks separately.
>

I experimented with this for elasticity:

    https://code.launchpad.net/~garth-wells/dolfin/block-dofmaps

but came to the conclusion (for the application I considered) that it
was not worth the effort and increase in complexity. It's simple for a
library that handles just one problem, but when generating code for
arbitrary PDEs with arbitrary combinations of function spaces it's not
straightforward to identify the local 'node' blocks for compressed row
storage. Complex numbers would be easier to handle because they always
come in pairs.

> * Creating a new data structure that can be fed with, e.g., four linear operators A, B, C, D to be composed to something like
> [ A, B ]
> [ C, D ]
> The entities A, B, C, D could still be treated separately which would make it possible, for exampe, to calculate the Schur complement w.r.t. A or D. The Trilinos package Teko has gone some lengths here and may be worth looking at (design-wise, too).
>

This is a current priority. My initial focus will be on PETSc and
later Trilinos since in my experience PETSc is better for very large
problems in parallel, is better documented and the PETSc devs are more
responsive. Plus PETSc provides a nice interface to the useful parts
of Trilinos (i.e. ML).

Garth

> I believe what you are talking about is no. 2 when you mention things like block preconditioners.
> What I had in mind for the faked complex-valued problem was no. 1 though. One could certainly formulate a problem with such tightly coupled parameters such as Re(u) and Im(u) in the form [A B; C D], but only think about ILU and fill ins for a k*k block-diagonal matrix.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

> So I would say there are two different notions of "block matrices" here:
>
> * Creating a new data structure that can be fed with, e.g., four linear operators A, B, C, D to be composed to something like
> [ A, B ]
> [ C, D ]
> The entities A, B, C, D could still be treated separately which would make it possible, for exampe, to calculate the Schur complement w.r.t. A or D. The Trilinos package Teko has gone some lengths here and may be worth looking at (design-wise, too).
>
> I believe what you are talking about is no. 2 when you mention things like block preconditioners.
> What I had in mind for the faked complex-valued problem was no. 1 though. One could certainly formulate a problem with such tightly coupled parameters such as Re(u) and Im(u) in the form [A B; C D], but only think about ILU and fill ins for a k*k block-diagonal matrix.

I always appear when I hear cbc.block mentioned ;)

Here's an example that follows your #2 above. It looks like it
converges fine in just a few iterations (8 on 30x30x30). That said, I
think of cbc.block mainly as a tool for experimentation, not for
production runs.

===

from dolfin import *
from block import block_mat, block_vec
from block.algebraic.trilinos import ML
from block.iterative import LGMRES

mesh = UnitCube(30, 30, 30)
PN = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
bc = DirichletBC(PN, (0, 0, 0), "on_boundary")

fr = project(Constant([0.0, 0.0, 1.0]), PN)
fi = project(Constant([0.0, 0.0, 1.0]), PN)

u, v = TrialFunction(PN), TestFunction(PN)

a = inner(curl(u), curl(v))*dx
A = assemble(a, bcs=bc)

Lr = -inner(fr,v)*dx + inner(fi, v)*dx
Li = -inner(fr,v)*dx - inner(fi, v)*dx

br = assemble(Lr, bcs=bc)
bi = assemble(Li, bcs=bc)

AA = block_mat([[A,-A], [A, A]])
bb = block_vec([br, bi])

mlA = ML(A)
AAp = block_mat([[mlA, 0], [0, mlA]])
# or: AAp = AA.scheme('jacobi', inverse=ML)

AAinv = LGMRES(AA, precond=AAp, tolerance=1e-6, show=1)

set_log_level(PROGRESS)
xx = AAinv * bb

print 'Residual norm reduced by', (AA*xx-bb).norm('l2') / bb.norm('l2')
===

@Garth
I agree that a matrix of blocks would only make sense if the interactions between the components are pointwise, such as the case for real-imaginary.
I wouldn't know what a reasonable API for this would look like, though, so I guess I'm going to have to look at your elasticity code.

@Joachim
Thanks for your example code!
My actual operator is `curlcurlA + iA`, so I guess the block matrix would be looking more complex for me, but anyhow this will be a good makeshift solution.

> @Joachim
> Thanks for your example code!
> My actual operator is `curlcurlA + iA`, so I guess the block matrix would be looking more complex for me, but anyhow this will be a good makeshift solution.

In case you actually run it, there's a small bug in my example.
Dirichlet rows should be zero in off-diagonal blocks. This should do
the trick:

A0 = A.copy()
bc.zero(A0)
AA = block_mat([[A, -A0], [A0, A]])

-j.

@joachim
The solution xx is a pair of Epetra-wrapped objects. Can they be inserted back into a Dolfin-friendly format?

> @joachim
> The solution xx is a pair of Epetra-wrapped objects. Can they be inserted back into a Dolfin-friendly format?

They are Dolfin vectors (from the Epetra backend in this case). You can do f.x.

Tr = Function(PN, xx[0])
Ti = Function(PN, xx[1])

-j.