MPI failure in `anyv` region

[mrhardman]

On commit ab45bf5 (branch https://github.com/mabarnes/moment_kinetics/tree/implicit-fp-collisions-merge-master). I have found an MPI bug that reveals itself when many cores are used. I noticed this because of check-in tests.

E.g., run

$ julia --project -O3 precompile.jl
$ mpirun -n 8 --output-filename mpitest julia --project -O3 -Jmoment_kinetics.so moment_kinetics/test/fokker_planck_tests.jl

and there is output showing errors from an attempt to access an array of different shape to expected in the bounds check https://github.com/mabarnes/moment_kinetics/blob/f4ef4c0942a953f5755bccffd59a9398504d9b0a/moment_kinetics/src/fokker_planck.jl#L434. Commenting these bounds checks just moves the error elsewhere: the key to the problem must be that cores that have nothing to do are trying to access data which doesn't exist for them. I don't believe that we had this particular error before, as I have tried to run with tens of cores in the past, with success. To get the error message below, I commented out the @global_timer macros in fokker_plank.jl.

Error During Test at /home/mrhardman/excalibur-mk/moment_kinetics/moment_kinetics/test/fokker_planck_tests.jl:431
  Got exception outside of a @test
  BoundsError: attempt to access 0×0 Matrix{Float64}
  Stacktrace:
    [1] fokker_planck_collision_operator_weak_form!(ffs_in::Matrix{Float64}, ffsp_in::Matrix{Float64}, ms::Float64, msp::Float64, nussp::Float64, fkpl_arrays::moment_kinetics.fokker_planck_calculus.fokkerplanck_weakform_arrays_struct{SparseArrays.SparseMatrixCSC{Float64, Int64}}, vperp::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, vpa::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, vperp_spectral::moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}, vpa_spectral::moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}; test_assembly_serial::Bool, use_Maxwellian_Rosenbluth_coefficients::Bool, use_Maxwellian_field_particle_distribution::Bool, algebraic_solve_for_d2Gdvperp2::Bool, calculate_GG::Bool, calculate_dGdvperp::Bool, boundary_data_option::moment_kinetics.input_structs.boundary_data_type)
      @ moment_kinetics.fokker_planck ~/excalibur-mk/moment_kinetics/moment_kinetics/src/fokker_planck.jl:434
    [2] fokker_planck_collision_operator_weak_form!
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/src/fokker_planck.jl:425 [inlined]
    [3] macro expansion
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/src/fokker_planck.jl:393 [inlined]
    [4] fokker_planck_self_collision_operator_weak_form!(pdf_in::Matrix{Float64}, ms::Float64, nuss::Float64, fkpl_arrays::moment_kinetics.fokker_planck_calculus.fokkerplanck_weakform_arrays_struct{SparseArrays.SparseMatrixCSC{Float64, Int64}}, vperp::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, vpa::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, vperp_spectral::moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}, vpa_spectral::moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}; boundary_data_option::moment_kinetics.input_structs.boundary_data_type, use_conserving_corrections::Bool)
      @ moment_kinetics.fokker_planck ~/excalibur-mk/moment_kinetics/.julia/packages/TimerOutputs/NRdsv/src/TimerOutput.jl:237
    [5] fokker_planck_self_collision_operator_weak_form!
      @ ~/excalibur-mk/moment_kinetics/.julia/packages/TimerOutputs/NRdsv/src/TimerOutput.jl:230 [inlined]
    [6] (::moment_kinetics.fokker_planck.var"#residual_func!#12"{moment_kinetics.fokker_planck.var"#residual_func!#11#13"{Bool, moment_kinetics.input_structs.boundary_data_type, Matrix{Float64}, Float64, Float64, Float64, moment_kinetics.fokker_planck_calculus.fokkerplanck_weakform_arrays_struct{SparseArrays.SparseMatrixCSC{Float64, Int64}}, moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}, moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}, moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}}})(Fresidual::Matrix{Float64}, Fnew::Matrix{Float64}; krylov::Bool)
      @ moment_kinetics.fokker_planck ~/excalibur-mk/moment_kinetics/moment_kinetics/src/fokker_planck.jl:1012
    [7] residual_func!
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/src/fokker_planck.jl:1011 [inlined]
    [8] newton_solve!(x::Matrix{Float64}, residual_func!::moment_kinetics.fokker_planck.var"#residual_func!#12"{moment_kinetics.fokker_planck.var"#residual_func!#11#13"{Bool, moment_kinetics.input_structs.boundary_data_type, Matrix{Float64}, Float64, Float64, Float64, moment_kinetics.fokker_planck_calculus.fokkerplanck_weakform_arrays_struct{SparseArrays.SparseMatrixCSC{Float64, Int64}}, moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}, moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}, moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}}}, residual::Matrix{Float64}, delta_x::Matrix{Float64}, rhs_delta::Matrix{Float64}, v::Matrix{Float64}, w::Matrix{Float64}, nl_solver_params::moment_kinetics.nonlinear_solvers.nl_solver_info{Matrix{Float64}, Array{Float64, 3}, Vector{Float64}, Vector{Float64}, Nothing, Val{:none}}, solver_type::Val{:vperpvpa}; left_preconditioner::Nothing, right_preconditioner::Function, recalculate_preconditioner::Nothing, coords::@NamedTuple{vperp::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, vpa::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}})
      @ moment_kinetics.nonlinear_solvers ~/excalibur-mk/moment_kinetics/moment_kinetics/src/nonlinear_solvers.jl:451
    [9] macro expansion
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/src/nonlinear_solvers.jl:429 [inlined]
   [10] newton_solve!(x::Matrix{Float64}, residual_func!::Function, residual::Matrix{Float64}, delta_x::Matrix{Float64}, rhs_delta::Matrix{Float64}, v::Matrix{Float64}, w::Matrix{Float64}, nl_solver_params::moment_kinetics.nonlinear_solvers.nl_solver_info{Matrix{Float64}, Array{Float64, 3}, Vector{Float64}, Vector{Float64}, Nothing, Val{:none}}; left_preconditioner::Nothing, right_preconditioner::Function, recalculate_preconditioner::Nothing, coords::@NamedTuple{vperp::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, vpa::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}})
      @ moment_kinetics.nonlinear_solvers ~/excalibur-mk/moment_kinetics/.julia/packages/TimerOutputs/NRdsv/src/TimerOutput.jl:237
   [11] fokker_planck_self_collisions_backward_euler_step!(Fold::Matrix{Float64}, delta_t::Float64, ms::Float64, nuss::Float64, fkpl_arrays::moment_kinetics.fokker_planck_calculus.fokkerplanck_weakform_arrays_struct{SparseArrays.SparseMatrixCSC{Float64, Int64}}, coords::@NamedTuple{vperp::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}, vpa::moment_kinetics.coordinates.coordinate{Vector{Float64}, Vector{Int64}, @NamedTuple{}}}, spectral::@NamedTuple{vperp_spectral::moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}, vpa_spectral::moment_kinetics.gauss_legendre.gausslegendre_info{SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseMatricesCSR.SparseMatrixCSR{1, Float64, Int64}, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, Nothing, Nothing}}, nl_solver_params::moment_kinetics.nonlinear_solvers.nl_solver_info{Matrix{Float64}, Array{Float64, 3}, Vector{Float64}, Vector{Float64}, Nothing, Val{:none}}; test_numerical_conserving_terms::Bool, test_linearised_advance::Bool, test_particle_preconditioner::Bool, use_Maxwellian_Rosenbluth_coefficients_in_preconditioner::Bool, boundary_data_option::moment_kinetics.input_structs.boundary_data_type)
      @ moment_kinetics.fokker_planck ~/excalibur-mk/moment_kinetics/moment_kinetics/src/fokker_planck.jl:1064
   [12] backward_Euler_fokker_planck_self_collisions_test(; vth0::Float64, vperp0::Float64, vpa0::Float64, ngrid::Int64, nelement_vpa::Int64, nelement_vperp::Int64, Lvpa::Float64, Lvperp::Float64, bc_vpa::String, bc_vperp::String, ntime::Int64, delta_t::Float64, test_particle_preconditioner::Bool, test_linearised_advance::Bool, use_Maxwellian_Rosenbluth_coefficients_in_preconditioner::Bool, test_dense_construction::Bool, test_numerical_conserving_terms::Bool, boundary_data_option::moment_kinetics.input_structs.boundary_data_type, print_to_screen::Bool, atol_max::Float64, atol_L2::Float64, atol_dens::Float64, atol_upar::Float64, atol_vth::Float64)
      @ Main.FokkerPlanckTests ~/excalibur-mk/moment_kinetics/moment_kinetics/test/fokker_planck_tests.jl:315
   [13] backward_Euler_fokker_planck_self_collisions_test
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/test/fokker_planck_tests.jl:212 [inlined]
   [14] macro expansion
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/test/fokker_planck_tests.jl:436 [inlined]
   [15] macro expansion
      @ ~/.julia/juliaup/julia-1.10.8+0.x64.linux.gnu/share/julia/stdlib/v1.10/Test/src/Test.jl:1669 [inlined]
   [16] macro expansion
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/test/fokker_planck_tests.jl:431 [inlined]
   [17] macro expansion
      @ ~/.julia/juliaup/julia-1.10.8+0.x64.linux.gnu/share/julia/stdlib/v1.10/Test/src/Test.jl:1577 [inlined]
   [18] macro expansion
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/test/fokker_planck_tests.jl:430 [inlined]
   [19] macro expansion
      @ ~/.julia/juliaup/julia-1.10.8+0.x64.linux.gnu/share/julia/stdlib/v1.10/Test/src/Test.jl:1577 [inlined]
   [20] runtests()
      @ Main.FokkerPlanckTests ~/excalibur-mk/moment_kinetics/moment_kinetics/test/fokker_planck_tests.jl:428
   [21] top-level scope
      @ ~/excalibur-mk/moment_kinetics/moment_kinetics/test/fokker_planck_tests.jl:1155
   [22] include(mod::Module, _path::String)
      @ Base ./Base.jl:495
   [23] exec_options(opts::Base.JLOptions)
      @ Base ./client.jl:323
   [24] _start()
      @ Base ./client.jl:557

This problem is a blocker for me to proceed with development.

[mrhardman]

@johnomotani Suggestions appreciated!

[johnomotani]

This seems odd.

• Does it happen on every MPI process?
• It seems like there's a zero-size array. I don't see any reason zero-size arrays should be allocated or passed around. Suggest tracking back which array it is, and what its size is at the top level (in test/fokker_planck_tests.jl).

[mrhardman]

This seems odd.

• Does it happen on every MPI process?
• It seems like there's a zero-size array. I don't see any reason zero-size arrays should be allocated or passed around. Suggest tracking back which array it is, and what its size is at the top level (in test/fokker_planck_tests.jl).

In the example above, only on process 5, 6 and 7. In that example, it was ffsp_in, but I have also seen it affect the Rosenbluth potentials. I will look again for more detail.

[johnomotani]

The algorithm that sets up the anyv parallelism is allowed to not include some processes. Sometimes this can be more efficient, e.g. if a grid is 3x4 it might parallelise better with 6 processes than 8... Not sure it's done the best thing for this test, but the algorithm has decided to use 5 processes instead of all 8. anyv operations are supposed to be inside a loop over {s,r,z} (or over {r,z} once #361 is merged) so that the unused processes just don't do anything, and will never try to access anyv-shared-memory arrays. So the fix is to just put the function call inside an @loop_s_r_z is ir iz begin ... end.

Fixed in #382 by 2a1191e.