Performance issue in GMGPolar solver

[julianlitz]

"Currently in my tests GMGPolar is slower than our other solver which is unexpected so if you have any ideas to speed things up this could be good. Currently my main hypothesis is that the difference comes from the fact that the other solver allows us to provide an initial guess to reduce the number of iterations whereas GMGPolar doesn't allow this"

@EmilyBourne
Can you run an instance of the solver with

verbose=1 and
solver.printTimings()

so I can take a look at what is causing the issue?

[EmilyBourne]

Here is the output of the first call to GMGPolar:

------------------------------
------- CMake Settings -------
------------------------------
Build: Release
MUMPS: OFF, Likwid: OFF
------------------------------
------ General Settings ------
------------------------------
Maximum number of threads: 8
Across the origin (Interior boundary condition)
A-Take (Stencil Distribution)
Domain geometry mode: Precomputed
Density profile mode: Precomputed
------------------------------
---------- PolarGrid ---------
------------------------------
r ∈ [1.28, 160], θ ∈ [0, 2π]
(nr × nθ) = (33 × 32) → (5 × 4)
Smoother: 5 circles
Splitting radius = 26.08
------------------------------
----- Multigrid settings -----
------------------------------
No Extrapolation
Multigrid Cycle: V(1,1)-Cycle
Number of levels: 4
Residual Norm: Weighted Euclidean (scaled L2 norm)
Tolerances: abs = 1e-08, rel = 0, maxiter = 1000
Full-Multigrid: 2x F(1,1)-Cycle
Preconditioned Conjugate Gradient: Disabled
------------------------------
---- Multigrid Iterations ----
------------------------------
it     ||r_k||      ||r_k||/||rhs||    
0      1.27e-04       4.00e-01           
1      4.87e-07       1.53e-03           
2      3.61e-08       1.14e-04           
3      2.46e-09       7.75e-06           
------------------------------
Total Iterations: 3
Reduction Factor: ρ = 0.019788

------------------
Timing Information
------------------
Setup Time: 0.0138502 seconds
    Create Levels: 0.00314136 seconds
    Smoother: 0.00411497 seconds
    Direct Solver: 0.00640466 seconds
    (Build rhs: 0.00133121 seconds)

Solve Time: 0.0150893 seconds
    Initial Approximation: 0.00768741 seconds
    Multigrid Iterations: 0.00680089 seconds
    Check Convergence: 0.00051269 seconds
    (Check Exact Error: 3.21e-07 seconds)

Average Multigrid Iteration: 0.00226696 seconds
    PreSmoothing: 0.000926368 seconds
    PostSmoothing: 0.000905858 seconds
    Residual: 0.000104623 seconds
    DirectSolve: 0.000120198 seconds
    Other Computations: 0.000209917 seconds
 ---> Save diagnostic 0 at iteration 0 for time 0
 -----> Not distributed quantities 
iter_run = 1 / 10 ; iter_simu = 1; time_simu = 0.5

The average time for a total solver call at this size is 0.0445928 s. The average time for the alternative solver is 0.0243903 s. The total solver call time is calculated with a fenced timer around : https://github.com/gyselax/gyselalibxx/blob/b90aebcad80e82f17cd3e4ed2088ecfa491e1390/src/pde_solvers/gmg_polar_poisson_like_solver.hpp#L247

[julianlitz]

Facts Overview:

• You have a very tiny problem size
• Extrapolation is Off
• PCG isnt used
• 2 FMG + 3 MG Iterations
• No component is sticking out to be particularly slow.
• When compared against FEM solver we are 2x slower.
• Rmin is quite large.

The FEM solver is just a direct solver and thus it is expected that we cannot beat its speed on these small meshes.

If you set maxLevels=1 then you will likely reach similar speeds but then you cannot use extrapolation. If you add Mumps on top you will surely beat the FEM solver here.

GMGPolar scales linearly with the problem size and you will notice the improved performance really only for larger mesh sizes.
I expect that we have a break even point at 65x64 where GMGPolar is the same speed as FEM. For problem sizes larger than that we will be a lot faster than FEM.

[julianlitz]

You could tune the parameters a bit to optimize it for smaller meshes.

E.g. make it a bit faster by reducing the FMG to
1x V-Cycle instead of 2x F-Cycle.

The parameters I gave you are optimized for meshes which are bit larger.

Though I don't think that it will help much. The only reasonable choice to make it faster is setting maxLevels=1.

[julianlitz]

@EmilyBourne Your benchmark is on the left side of the red barrier. I marked the break even point with a red circle :)

[EmilyBourne]

Facts Overview:

• You have a very tiny problem size
• Extrapolation is Off
• PCG isnt used
• 2 FMG + 3 MG Iterations
• No component is sticking out to be particularly slow.
• When compared against FEM solver we are 2x slower.
• Rmin is quite large.

The FEM solver is just a direct solver and thus it is expected that we cannot beat its speed on these small meshes.

If you set maxLevels=1 then you will likely reach similar speeds but then you cannot use extrapolation. If you add Mumps on top you will surely beat the FEM solver here.

GMGPolar scales linearly with the problem size and you will notice the improved performance really only for larger mesh sizes. I expect that we have a break even point at 65x64 where GMGPolar is the same speed as FEM. For problem sizes larger than that we will be a lot faster than FEM.

This is definitely a bit too small. I am running my test case on a small machine so I had to reduce my test case. The problem was originally flagged to me on a 64x64 size mesh on AMD where the direct solver was 4x faster.

On v100 for 64x64 I have:

Solver	Time for 10 steps (s)
Polar FEM	0.375098
GMGPolar	0.659595

• Rmin is quite large.

Is this likely to affect the timings?

Though I don't think that it will help much. The only reasonable choice to make it faster is setting maxLevels=1.

Doesn't this mean that you are only using the solver and not doing any multi-grid?

[EmilyBourne]

What are V1 and V2 in your figure? Should the red line read V1 serial? What is meant by parallel? CPU with N threads? (how many is N?)

[julianlitz]

No Rmin doesn't influence the time it take. Only influence the quality of the solution.
Yes setting maxLevels=1 will result in no multigrid beeing used. Instead we just do a LU Decomp of the original matrix and solve it directly.

[julianlitz]

V1 was the original code you used in the Comparison Paper and V2 is the current GMGPolar code.

I just drew approximately by hand in the values I got on my laptop to get a sense of the relations.
Though I dont know anymore what specific paramters were used as it was just a quick sketch.

Yes the text reads
v1 serial, v1 parallel, v2 serial, v2 parallel

[julianlitz]

Since the code Kokkos code was never tested on GPU before there is probably room for further optimizations.
I made an issue a few months ago about Kokkos thread ranges which I expect to be a bit faster, but likely not much.

[julianlitz]

If you notice a 1:1 balance between the setup time and the solve time, this indicates that the solve phase is already relatively quick and that there really is not much more room to make it faster as half of the time is spend for memory allocations in the setup() phase.

[EmilyBourne]

Turning on PCG with the settings we have currently really ruins the timings for some reason. It took 91.1452 s for the 10 steps!

[EmilyBourne]

Can we avoid setting up so often? At least only as often as we want to update the coefficients?

[julianlitz]

Thats strange. I would expect that with Extrapolation=None it will take the exact same amount of time, as the time per iterations doubles, but you also need 2x less iterations. You can always just paste the timings output here in the chat so I can see whats causing this.

[EmilyBourne]

------------------------------
------- CMake Settings -------
------------------------------
Build: Release
MUMPS: OFF, Likwid: OFF
------------------------------
------ General Settings ------
------------------------------
Maximum number of threads: 8
Across the origin (Interior boundary condition)
A-Take (Stencil Distribution)
Domain geometry mode: Precomputed
Density profile mode: Precomputed
------------------------------
---------- PolarGrid ---------
------------------------------
r ∈ [1.28, 160], θ ∈ [0, 2π]
(nr × nθ) = (33 × 32) → (5 × 4)
Smoother: 5 circles
Splitting radius = 26.08
------------------------------
----- Multigrid settings -----
------------------------------
No Extrapolation
Multigrid Cycle: V(1,1)-Cycle
Number of levels: 4
Residual Norm: Weighted Euclidean (scaled L2 norm)
Tolerances: abs = 1e-08, rel = 0, maxiter = 1000
Full-Multigrid: 2x F(1,1)-Cycle
Preconditioned Conjugate Gradient: Enabled
- PCG Full-Multigrid: 1x V(1,1)-Cycle
- PCG Multigrid Iteration: 1x V(1,1)-Cycle
------------------------------
---- Multigrid Iterations ----
------------------------------
it     ||r_k||      ||r_k||/||rhs||    
0      1.70e-04       4.81e-01           
1      2.85e-05       8.06e-02           
2      2.24e-05       6.35e-02           
3      2.00e-05       5.66e-02           
4      1.09e-05       3.09e-02           
5      3.01e-06       8.53e-03           
6      4.21e-06       1.19e-02           
7      7.28e-06       2.06e-02           
8      7.70e-06       2.18e-02           
9      6.72e-06       1.90e-02           
10     5.23e-06       1.48e-02           
11     3.58e-06       1.01e-02           
12     1.93e-06       5.45e-03           
13     8.12e-07       2.30e-03           
14     1.74e-06       4.92e-03           
15     2.95e-06       8.36e-03           
16     3.88e-06       1.10e-02           
17     4.44e-06       1.26e-02           
18     4.67e-06       1.32e-02           
19     4.64e-06       1.31e-02           
20     4.44e-06       1.26e-02           
21     4.11e-06       1.16e-02           
22     3.70e-06       1.05e-02           
23     3.23e-06       9.15e-03           
24     2.74e-06       7.74e-03           
25     2.22e-06       6.27e-03           
26     1.69e-06       4.77e-03           
27     1.17e-06       3.31e-03           
...      
1000   1.98e-06       5.62e-03           
------------------------------
Total Iterations: 1000
Reduction Factor: ρ = 0.994832

------------------
Timing Information
------------------
Setup Time: 0.00658645 seconds
    Create Levels: 0.00202223 seconds
    Smoother: 0.00345905 seconds
    Direct Solver: 0.000913004 seconds
    (Build rhs: 0.000163373 seconds)

Solve Time: 4.32612 seconds
    Initial Approximation: 0.00558028 seconds
    Preconditioned Conjugate Gradient: 4.28129 seconds
    Check Convergence: 0.0390917 seconds
    (Check Exact Error: 4.9103e-05 seconds)

[julianlitz]

Currently you need to create a new GMGPolar instance each time you want to modify the alpha and beta coefficients and thus .setup() must be called repeatedly.

Currently there is not a way to update the alpha and beta values within the same GMGPolar instance.
It is probably difficult to implement but if it were possible, you would spend quite less time allocating memory.