I’m just reaching out to see if anyone has had good experiences with particular solvers when working with om2 output.
Specifically, I need to solve a 2d poisson equation for the scalar potential of a vector field with the aim of separating its divergent/rotational fields. To get my head around how this could work I’ve had a crack at writing a simple Gauss-Seidel style iterative technique with successive over relaxation and some preconditioning on coarser grids. This technique works nicely with xarray/dask as you can update iterate the whole array at once but…. super slow convergence. I can provide the code if that helps get an idea of what I’m trying to do.
Has anyone had success with available direct solvers (e.g. scipy.linalg.lstsq) or multigrid methods (e.g. pyAMG) that work with complex land masks and neumann bcs? or have alternative recommendations?
Back in the day I did this calculation (a Helmholtz decomposition) for vertically-integrated ocean heat transport (see Forget and Ferriera 2019 for the background as to why) in ACCESS-OM2-1 using MatLab’s mldivide function. My code is here (but its messy and old). It seemed to produce reasonable results (see figure). No idea if this is useful or not.
