Matrix-free-method for an operator depending on input from solver

[Maxi Miller]

I am currently trying to implement the JFNK-method in the matrix-free framework (by following step-37 and step-59) for solving nonlinear equations of the form F(u)=0. The method itself replaces the multiplication of the explicit jacobian J with the krylov vector v during solving with the approximation (F(u+e*v)-F(u))/e, which removes the matrix.

The initial test using a LinearOperator worked, where I modified the vmult-function accordingly (c.f. my earlier questions here), but still I had to form the jacobian at least once (which I wanted to avoid), thus the test of the matrix-free framework.

In theory I just should have to provide the vmult-function, but based on step-37 and step-59 I also have to provide the full operator at least once during initialization, as f.ex. to form the inverse diagonal in step-37. This would be the same as calculating (F(u_i+e*v_i)-F(u_i))/e, but I do not know the value of v_i.

Thus, what would be the best approach of implementing this operator using the MatrixFree-framework in deal.II, without having to form a single matrix, and without knowing the exact value of F? Ideally it can be used for a preconditioner too, such as the GMG in step-37 and step-59.

Thanks!

[Martin Kronbichler]

Dear Maxi,

There is not really a simple answer to your request: It depends on how good a preconditioner you need. Matrix-free methods work best if you can use simple preconditioners in the sense that they only need some matrix entries (if any) and the matrix-vector product (or operator evaluation in the nonlinear case). Note that I also consider multigrid with Chebyshev smoothing as in step-37 and step-59 a simple preconditioner. It is up to you to find out whether your problem admits such a simple form and iterative solvers, such that GMRES, some nonlinear Newton-Krylov or nonlinear Krylov variant, can cope with that and not use too many iterations. We as a mailing list cannot help with the preconditioner because this is a research question that many people would like an answer of. My experience is that matrix-free setups are mostly done for performance reasons (because the mat-vec is cheaper) for well-understood systems, like the Laplacian with multigrid or some linear/nonlinear time-dependent equations with dominating mass matrix terms where simple iterative methods suffice. Sometimes matrix-free is done because one can fit larger problems in memory. But all cases where I know it is used are done when one at least has an expectation of what the matrix and its eigenvalue spectrum look like. If you do not have that, because you have some difficult linear/nonlinear systems, matrix-free setups might help provide you with a black-box view - but they make the preconditioner selection harder, not easier, because you restrict yourself in the choice of preconditioner to the more basic variants.

Now to the technical part: If you need the diagonal of a matrix and rely on an accurate representation, you will need to build the linearization in some way or the other. This is also what step-37 or step-59 do - they compute the local stiffness matrix, throw away everything but the diagonal, and use that inside some solver. Tor the diagonal you would run through all unit vectors with v, i.e., you add (1,0,0,...), (0,1,0,...), and so on. Then you only keep the one entry where v was one. Your approach to do this globally is certainly not going to scale because it is an O(n^2) operation to check every unit vector. In principle, you could do it locally element-by-element similar to step-37 also for nonlinear operators. The question is whether you have the full action inside a single "local_apply" function such that you could start adding unit vectors there. Given a nonlinear field u that you linearize about, you could perturb it with e*v_i and compute the local diagonal according to your formula. From what I understand, you have the operator in an even more abstract way, so collecting everything in one local call with matrix-free facilities would also be non-trivial.

Best,
Martin

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/257bb96a-ba19-48cf-a7c4-49e6f178e417%40googlegroups.com.

[Maxi Miller]

Hei,

based on what I know from literature a GMG-preconditioner with Chebyshev smoothing should work rather well for my case, thus I intended to follow example 37 here. If I can set v_i to 1 for calculating the diagonal, that would simplify the problem quite a lot.

The main reason for testing the matrix-free approach in my case is enormous RAM-consumption when creating the matrix (I have to use a rather dense grid, else there will be errors), and a long build time for the matrix itself.

I hope that makes sense?

To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.