Should we make Conjugate Gradient + FFTBased-preconditioner the default for HydrostaticFreeSurfaceModel?

[navidcy]

Seems that #2412 suggests that this is the fastest way for grids that don't have regular horizontal spacing or have an immersed boundary.

Yet, the free surface solver constructor falls back to the MatrixBasedSolver in those cases.

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl

Line 101 in 95206ae

default_method = is_horizontally_regular(grid) ? :FastFourierTransform : :HeptadiagonalIterativeSolver

Should we change that?

[simone-silvestri]

Sounds good! we have to leave it only for LatitudeLongitudeGrid

[simone-silvestri]

I wonder about CPU execution though... because the MatrixSolver with ILU was not tested and that in my experience seemed always to be the fastest method

[glwagner]

These benchmarks:

#2412

show a factor of 2 speed up on the CPU. In principle, the FFT is faster than a matrix multiply which may explain it? Or, perhaps the benchmarks test the wrong thing so we should re-run them.

As explained on that PR, some performance is left on the table, because the fastest combination would use an FFT-based preconditioner with a matrix multiply to compute the LHS.

[simone-silvestri]

Yeah, but the ILU preconditioner doesn't seem to be tested, which I remember was quite good to reduce the number of iterations. So in case of a large time step that would be quite efficient even if we have to solver two triangular linear systems each iteration.

but anyways, maybe it's not faster, it was just I was wondering how it would compare

[glwagner]

Hmm, well I guess we should add that combination to

https://github.com/CliMA/Oceananigans.jl/blob/main/benchmark/benchmark_hydrostatic_model.jl

and re-run the benchmarks.

[simone-silvestri]

I think we can close this as we have the SplitExplicitFreeSurface as a default now