I'm developing an FEM solver for a coupled system. I have diffusion and potential equations which result in positive definite matrices for each equation, but the coupling makes the overall system matrix non-symmetric, non-positive-definite.
The solver is planned to run entirely on a GPU. Mesh and data representation is highly optimized for that and standard libraries are therefore not very useful, also because one would need to transfer between GPU and CPU, or had to use CPU-optimized data formats that would not perform well on the GPU.
I plan to use multigrid later on, so I currently use jacobi iterations to solve the equation. That works, let's say, semi-well. It converges not very fast and when errors get small, there are some signs of other problems (which I do not fully understand currently).
I'm now thinking about my next step. Should I try some CG-based algorithms? Or what would you try or investigate next?