Runtime scaling

[Mikolaj-A-Kowalski]

We need to investigate how the runtime is influenced by:

• I/O
• Number of grid points
• Number of particles

Tasks

• LOBSTER model
  • Scaling with number of grid cells
  • Scaling with number of particles (with CPU and GPU on the same, larger grid)
  • Effect of I/O
• PISCES model
  • Scaling with number of grid cells
  • Effect of I/O
  • Scaling with number of particles

[Mikolaj-A-Kowalski]

The following comments in #262 are really relevant to this issue:
#262 (comment)
#262 (comment)

[Mikolaj-A-Kowalski]

The initial runs for the LOBSTER model has been done. The total runtime (excluding compilation) and total memory allocated on a host are shown in the graphs below.

The raw data in CSV files contains extra details (e.g. total GC time):
cpu_scaling_io.csv
cpu_scaling_noio.csv
gpu_scaling_io.csv
gpu_scaling_noio.csv

The GPU run with I/O did not reach the same grid size as without I/O. This is because I was running them in my home directory and reached the 50 Gb quota 😅 . We should probably rerun the last two points with rds storage.

EDIT: Actually given that the rds is the performant filesystem we should rerun all the I/O runs with rds... it is likely that we will see a performance difference (i.e. rds should be faster).

[Mikolaj-A-Kowalski]

I has some issues with the scaling runs wrt number of particles. The calculations (both CPU and GPU) seemed to crash rather randomly with not very helpful error messages e.g.:

On GPU with CUDA:

ERROR: LoadError: KernelException: exception thrown during kernel execution on device NVIDIA A100-SXM4-80GB
Stacktrace:
  [1] check_exceptions()
    @ CUDA ~/.julia/packages/CUDA/ja0IX/src/compiler/exceptions.jl:39
  [2] synchronize(stream::CuStream; blocking::Bool, spin::Bool)
    @ CUDA ~/.julia/packages/CUDA/ja0IX/lib/cudadrv/synchronization.jl:207
  [3] synchronize
    @ ~/.julia/packages/CUDA/ja0IX/lib/cudadrv/synchronization.jl:194 [inlined]
  [4] synchronize
    @ ~/.julia/packages/CUDA/ja0IX/src/memory.jl:527 [inlined]
  [5] synchronize(x::CuArray{Float64, 3, CUDA.DeviceMemory})
    @ CUDA ~/.julia/packages/CUDA/ja0IX/src/array.jl:493
  [6] (::CUDA.var"#1161#1162"{Float64, Vector{Float64}, Int64, CuArray{Float64, 3, CUDA.DeviceMemory}, Int64, Int64})()
    @ CUDA ~/.julia/packages/CUDA/ja0IX/src/array.jl:613
  [7] #context!#1024
    @ ~/.julia/packages/CUDA/ja0IX/lib/cudadrv/state.jl:168 [inlined]
  [8] context!
    @ ~/.julia/packages/CUDA/ja0IX/lib/cudadrv/state.jl:163 [inlined]
  [9] unsafe_copyto!(dest::Vector{Float64}, doffs::Int64, src::CuArray{Float64, 3, CUDA.DeviceMemory}, soffs::Int64, n::Int64)
    @ CUDA ~/.julia/packages/CUDA/ja0IX/src/array.jl:608
 [10] copyto!
    @ ~/.julia/packages/CUDA/ja0IX/src/array.jl:560 [inlined]
 [11] getindex
    @ ~/.julia/packages/GPUArrays/u6tui/src/host/indexing.jl:52 [inlined]
 [12] getindex
    @ ~/.julia/packages/OffsetArrays/b7qpm/src/OffsetArrays.jl:437 [inlined]
 [13] _getindex
    @ ./abstractarray.jl:1341 [inlined]
 [14] getindex
    @ ./abstractarray.jl:1312 [inlined]
...

On CPU:

ERROR: LoadError: ArgumentError: The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.


Stacktrace:
  [1] assert_bracket
    @ ~/.julia/packages/Roots/ohGnS/src/Bracketing/bracketing.jl:52 [inlined]
  [2] #init_state#53
    @ ~/.julia/packages/Roots/ohGnS/src/Bracketing/bisection.jl:50 [inlined]
  [3] init_state(::Roots.Bisection, F::Roots.Callable_Function{Val{1}, Val{false}, typeof(OceanBioME.Models.CarbonChemistryModel.alkalinity_residual), @NamedTuple{DIC::Float64, Alk::Float64, boron::Float64, sulfate::Float64, fluoride::Float64, silicate::Float64, phosphate::Float64, K1::Float64, K2::Float64, KB::Float64, KW::Float64, KS::Float64, KF::Float64, KP1::Float64, KP2::Float64, KP3::Float64, KSi::Float64}}, x₀::Float64, x₁::Float64, fx₀::Float64, fx₁::Float64)
    @ Roots ~/.julia/packages/Roots/ohGnS/src/Bracketing/bisection.jl:34
  [4] init_state(M::Roots.Bisection, F::Roots.Callable_Function{Val{1}, Val{false}, typeof(OceanBioME.Models.CarbonChemistryModel.alkalinity_residual), @NamedTuple{DIC::Float64, Alk::Float64, boron::Float64, sulfate::Float64, fluoride::Float64, silicate::Float64, phosphate::Float64, K1::Float64, K2::Float64, KB::Float64, KW::Float64, KS::Float64, KF::Float64, KP1::Float64, KP2::Float64, KP3::Float64, KSi::Float64}}, x::Tuple{Float64, Float64})
    @ Roots ~/.julia/packages/Roots/ohGnS/src/Bracketing/bracketing.jl:6
  [5] #init#43
    @ ~/.julia/packages/Roots/ohGnS/src/find_zero.jl:299 [inlined]
  [6] init
    @ ~/.julia/packages/Roots/ohGnS/src/find_zero.jl:289 [inlined]
  [7] solve(𝑭𝑿::Roots.ZeroProblem{typeof(OceanBioME.Models.CarbonChemistryModel.alkalinity_residual), Tuple{Float64, Float64}}, M::Roots.Bisection, p::@NamedTuple{DIC::Float64, Alk::Float64, boron::Float64, sulfate::Float64, fluoride::Float64, silicate::Float64, phosphate::Float64, K1::Float64, K2::Float64, KB::Float64, KW::Float64, KS::Float64, KF::Float64, KP1::Float64, KP2::Float64, KP3::Float64, KSi::Float64}; verbose::Bool, kwargs::@Kwargs{tracks::Roots.NullTracks, atol::Float64})
    @ Roots ~/.julia/packages/Roots/ohGnS/src/find_zero.jl:497
  [8] solve
    @ ~/.julia/packages/Roots/ohGnS/src/find_zero.jl:490 [inlined]
  [9] #find_zero#40
    @ ~/.julia/packages/Roots/ohGnS/src/find_zero.jl:220 [inlined]
 [10] find_zero (repeats 2 times)
    @ ~/.julia/packages/Roots/ohGnS/src/find_zero.jl:210 [inlined]

[jagoosw]

Was this with the kelp particles or just normal particles?

[Mikolaj-A-Kowalski]

Was this with the kelp particles or just normal particles?

That would be SugerKelpParticles (so kelp particles).
Set up here in the script.

Originally there were only 5 initialised at different depths. Now the number is arbitrary with, but I kept the range of depths the same. I was thinking the problem may stem from that?

[jagoosw]

How many particles are you going up to? That error is from the carbon chemistry solve (and we so far haven't found a good way to catch errors which occur inside GPU kernels). If the number of particles is very large then it might be making the DIC leave the "normal" range which the chemistry model can cope with, or possibly its taking all the nutrients and making the phytoplankton go negative (although I'd expect an error from the light attenuation model if it is that).

[Mikolaj-A-Kowalski]

How many particles are you going up to?

The aim is of course as much as possible ;-)
But the errors appeared for:

• 1024 particles on CPU
• 128 particles on GPU

If I understand what you are saying the problems may arise from the concentration of the particles, right?
In that case perhaps a simple solution would be to distribute them in x and y axis as well?
The other suspicious thing in the modification I made are the scalefactors I think. If I understand correctly they set a 'weight' of each particle... in that case I should have probably decreased it when increasing the number of particles to keep the total quantity of the kelp constant?

[Mikolaj-A-Kowalski]

It seems to have been the total weight indeed!

-particles = SugarKelpParticles(n; grid, scalefactors = fill(2000, n))
+particles = SugarKelpParticles(n; grid, scalefactors = fill(2000/n, n))

Makes the calculation not crash (on a GPU) will check a CPU shortly (EDIT: CPU is a go 🎆 )

[jagoosw]

Changing the scale factor is a good solution!

[Mikolaj-A-Kowalski]

Here is the result for the scaling with the number of kelp particles:

But following the discussion during today’s meeting we should rerun it with CPU and GPU on the same grid size.

[Mikolaj-A-Kowalski]

But following the discussion during today’s meeting we should rerun it with CPU and GPU on the same grid size.

I would say that the conclusion is that for the LOBSTER model there is no difference in behaviour between CPU and GPU calculation (other than total runtime).

[jagoosw]

I have a question, I've been using the sugar kelp particles and it makes simple runs a lot slower. When I @profile nothing stands out, do you have any suggestions what to look for?

[Mikolaj-A-Kowalski]

I have a question, I've been using the sugar kelp particles and it makes simple runs a lot slower. When I @Profile nothing stands out, do you have any suggestions what to look for?

It is very difficult to say anything without seeing the case or the profile data. Just couple of things that first come to mind:

• Did you try to side-by-side compere @profile results from a run with and without particles? If all functions would scale by the same amount it would be very weird... something must increase more and it should contain clues about what is going on.
• I assume you are talking about the Profiling.@profile and a CPU run? Maybe trying a GPU run and using CUDA.@profile external=false could be helpful by checking if it is only CPU problem or it effects GPU as well. In the latter case CUDA's macro will give results on a higher granularity (kernels basically), which for me is always easier to interpret.
• Could it be that particles introduce extra I/O? I am thinking that with the Julia's sampling profiler no samples should collected when inside the system calls and perhaps no samples are collected when inside I/O operations... this is just a speculation though. I haven't verified that it indeed behaves like that.

[jagoosw]

• Did you try to side-by-side compere @profile results from a run with and without particles? If all functions would scale by the same amount it would be very weird... something must increase more and it should contain clues about what is going on.

Thank you for this suggestion, it turns out that it spends a significant extra time computing the tendencies for the kelp timestepping because they involve a root finding (which I forgot happened) so its quite a long operation so I guess this is to be expected

[johnryantaylor]

@jagoosw, where is the root finding happening? I wonder if there is a way that we can avoid that (even if it introduces additional assumptions).

[jagoosw]

@jagoosw, where is the root finding happening? I wonder if there is a way that we can avoid that (even if it introduces additional assumptions).

They have a very complicated formulation for photosynthesis which requires solving for a parameter, I have tried replacing the find_zero with an inhouse newton raphson and it seems to be a lot faster (since it converges in ~10 steps). I'm working on it in #308

[AdelekeBankole]

Runtime scaling for PISCES