Try using Krylov to solve the Poisson equation

[glwagner]

This PR adds support for Krylov.jl and aims to test it out for solving the Poisson equation.

cc @xkykai @amontoison

Closes #3803

[glwagner]

To do:

• Implement mul! for the Laplacian operator
• Implement a wrapper type for Field that subtypes AbstractVector?
• Implement ldiv! for the FFT-based preconditioner (after the first two are done), planning to use ldiv=true with cg

This example illustrates implementing a custom operator:

https://jso.dev/Krylov.jl/dev/matrix_free/

Here is the documentation for implementing preconditioners in Krylov:

https://jso.dev/Krylov.jl/dev/preconditioners/

In Krylov, we want to remove broadcasting.

To make Field look like a vector we just need

struct QuasiVector{FT, F} <: AbstractVector{FT}
    field :: F
    function QuasiVector(field)
        FT = eltype(field)
        F = typeof(field)
        return new{FT, F}(field)
    end
end

[amontoison]

@glwagner I remember now why the right-hand side b is an AbstractVector{FC}: I use the type FC for type stability with the keyword arguments (tolerances, etc.) as well as default values (sqrt(eps(real(FC)))).

[glwagner]

@glwagner I remember now why the right-hand side b is an AbstractVector{FC}: I use the type FC for type stability with the keyword arguments (tolerances, etc.) as well as default values (sqrt(eps(real(FC)))).

True but then you could use AbstractArray{FT} similarly right?

[amontoison]

Yes, I can use AbstractArray{FT}, but it should represent a vector (whatever format he has).
For common array types, it always leads to an AbstractVector, so I never found the need to remove it.
I was always hesitant to relax this constraint because the type of b is used for allocating the "vectors" in the workspace of a Krylov method.
If a user makes a mistake and provides a matrix with one column instead of a vector, it will fail internally with an unfriendly error message.

[glwagner]

If a user makes a mistake and provides a matrix with one column instead of a vector, it will fail internally with an unfriendly error message.

Ah for sure. I was just remarking that if you need the eltype, you can use AbstractArray{FT} (or you can simply call eltype, which is preferred by YASGuide --- type parameters should be used for dispatch only).

I think it's ok if you keep AbstractVector. We can create a wrapper for Krylov.

[amontoison]

@glwagner @xkykai
I added a example in the documentation for you:
https://jso.dev/Krylov.jl/previews/PR927/custom_workspaces

[glwagner]

@glwagner @xkykai I added a example in the documentation for you: https://jso.dev/Krylov.jl/previews/PR927/custom_workspaces

Just a comment --- the Poisson operator is usually denoted $\nabla^2$ (rather than just $\nabla$ which denotes a gradient and yields a vector):

otherwise this is great!

[amontoison]

Thanks!
Good catch for typo with \nabla, I will fix it.