2D code appears to hang without error when running on HPC

[mrhardman]

Using the branch radial-vperp-standard-DKE-with-neutrals-pdf-animations (commit da276a4) for simulations on a HPC cluster, I noted that the simulation appeared to hang after a relatively small number of timesteps, with no error messages. Using the new screen printout that shows the time when the last write to file was completed, the simulation failed to advance to the next nwrite steps after many multiples of the time that the simulation should have taken. The simulation could only be cancelled.

e.g.

 beginning time advance...14:27:41
 finished time step 100  14:32:35
 finished time step 200  14:37:25

But the simulation ran for another hour without further output of any kind.

[mrhardman]

This issue appears to be related to the changes made in #68, as anticipated by the comments on error handling. Running a case which hangs on multiple cores on a single core, we obtain the following error message, which must presumably been suppressed by the try catch block. The underlying error is due to the distribution function going sufficiently negative that the density is negative, and the potential cannot be calculated. @johnomotani FYI.

Activating project at ~/excalibur/moment_kinetics_2D
ERROR: LoadError: UndefVarError: global_size not defined
Stacktrace:
[1] run_moment_kinetics(to::TimerOutputs.TimerOutput, input_dict::Dict{String, Any})
@ moment_kinetics ~/excalibur/moment_kinetics_2D/src/moment_kinetics.jl:104
[2] run_moment_kinetics
@ ~/excalibur/moment_kinetics_2D/src/moment_kinetics.jl:119 [inlined]
[3] run_moment_kinetics(input::String)
@ moment_kinetics ~/excalibur/moment_kinetics_2D/src/moment_kinetics.jl:126
[4] run_moment_kinetics()
@ moment_kinetics ~/excalibur/moment_kinetics_2D/src/moment_kinetics.jl:137
[5] top-level scope
@ ~/excalibur/moment_kinetics_2D/run_moment_kinetics.jl:7
in expression starting at /home/h/hardman/excalibur/moment_kinetics_2D/run_moment_kinetics.jl:7

caused by: DomainError with -0.05520460158933985:
log will only return a complex result if called with a complex argument. Try log(Complex(x)).
Stacktrace:
[1] throw_complex_domainerror(f::Symbol, x::Float64)
@ Base.Math ./math.jl:33
[2] _log(x::Float64, base::Val{:ℯ}, func::Symbol)
@ Base.Math ./special/log.jl:304
[3] log
@ ./special/log.jl:269 [inlined]
[4] macro expansion
@ ~/excalibur/moment_kinetics_2D/src/em_fields.jl:96 [inlined]
[5] macro expansion
@ ~/excalibur/moment_kinetics_2D/src/looping.jl:294 [inlined]
[6] macro expansion
@ ~/excalibur/moment_kinetics_2D/src/em_fields.jl:95 [inlined]
[7] macro expansion
@ ~/excalibur/moment_kinetics_2D/src/looping.jl:294 [inlined]
[8] update_phi!(fields::moment_kinetics.moment_kinetics_structs.em_fields_struct, fvec::moment_kinetics.moment_kinetics_structs.scratch_pdf{5, 3, 6, 3}, z::moment_kinetics.coordinates.coordinate, r::moment_kinetics.coordinates.coordinate, composition::moment_kinetics.input_structs.species_composition, z_spectral::moment_kinetics.chebyshev.chebyshev_info{FFTW.cFFTWPlan{ComplexF64, -1, true, 1, UnitRange{Int64}}, AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, true, 1, UnitRange{Int64}}, Float64}}, r_spectral::moment_kinetics.chebyshev.chebyshev_info{FFTW.cFFTWPlan{ComplexF64, -1, true, 1, UnitRange{Int64}}, AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, true, 1, UnitRange{Int64}}, Float64}})
@ moment_kinetics.em_fields ~/excalibur/moment_kinetics_2D/src/em_fields.jl:62
[9] rk_update!(scratch::Vector{moment_kinetics.moment_kinetics_structs.scratch_pdf{5, 3, 6, 3}}, pdf::moment_kinetics.initial_conditions.pdf_struct, moments::moment_kinetics.initial_conditions.moments_struct, fields::moment_kinetics.moment_kinetics_structs.em_fields_struct, vz::moment_kinetics.coordinates.coordinate, vr::moment_kinetics.coordinates.coordinate, vzeta::moment_kinetics.coordinates.coordinate, vpa::moment_kinetics.coordinates.coordinate, vperp::moment_kinetics.coordinates.coordinate, z::moment_kinetics.coordinates.coordinate, r::moment_kinetics.coordinates.coordinate, rk_coefs::SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, istage::Int64, composition::moment_kinetics.input_structs.species_composition, z_spectral::moment_kinetics.chebyshev.chebyshev_info{FFTW.cFFTWPlan{ComplexF64, -1, true, 1, UnitRange{Int64}}, AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, true, 1, UnitRange{Int64}}, Float64}}, r_spectral::moment_kinetics.chebyshev.chebyshev_info{FFTW.cFFTWPlan{ComplexF64, -1, true, 1, UnitRange{Int64}}, AbstractFFTs.ScaledPlan{ComplexF64, FFTW.cFFTWPlan{ComplexF64, 1, true, 1, UnitRange{Int64}}, Float64}})
@ moment_kinetics.time_advance ~/excalibur/moment_kinetics_2D/src/time_advance.jl:615

[mrhardman]

The following code is a Julia script that mocks up how an error should be handled by the current setup in moment_kinetics.jl. Mysteriously, this code below works as intended, calling MPI.Abort and showing an error message for the cause of the problem, even when multiple cores are used.

if abspath(PROGRAM_FILE) == @__FILE__
	using Pkg
	Pkg.activate(".")
	using MPI
	MPI.Init()
	comm = MPI.COMM_WORLD
	println("Hello world, I am $(MPI.Comm_rank(comm)) of $(MPI.Comm_size(comm))")
	try
		
		try
			# set up something that differs across ranks
			if MPI.Comm_rank(comm) > 0
				x = 3
			elseif MPI.Comm_rank(comm) == 0 
				x = -1
			end
			#do a calculation that will fail on some ranks
			y = sqrt(x)
			println( "rank = ",MPI.Comm_rank(comm)," y = ",y)
		finally
			y = 0
		end

	catch e
		
		if MPI.Comm_size(comm) > 1
			println(e)
			MPI.Abort(comm, 1)
		end

		rethrow(e)
	end
end

[mrhardman]

A partial solution to this problem is provided in commits 2cb3835 and 5c087a2. We apply a try-catch statement around offending log() and sqrt() calls that are typically the source of DomainError. This solution appears to be successful in that a call to MPI.Abort() is now made, and the code exits rather than hangs. There may be other errors of this type that are not affected by these changes.

I will leave the issue open.

@johnomotani can you please comment on my fixes? Should the main branch have similar try catch statements inserted?

[johnomotani]

The ERROR: LoadError: UndefVarError: global_size not defined is suspicious. It looks like the problem might be:

• global_size is not exported in communication.jl (this was just an oversight by me, it's missing on master as well).
• As global_size is undefined, the code hits another error in the catch-block at l. 104, before it would actually call MPI.Abort()

  moment_kinetics/src/moment_kinetics.jl

  Lines 101 to 110 in 7975a26

  	catch e
  	# Stop code from hanging when running on multiple processes if only one of them
  	# throws an error
  	if global_size[] > 1
  	println(e)
  	MPI.Abort(comm_world, 1)
  	end
  	rethrow(e)
  	end

• The processes that didn't throw the error carry on to the next MPI call as Abort() was not called, and hang there.

I don't know why no error at all gets printed currently when the code hangs, when it does get printed in serial, but apart from that I think the above makes sense. @mrhardman I'll make a PR for a fix to master that you can cherry-pick.

[johnomotani]

Ah, but I see you already implemented what I hoped was the 'fix' in 2cb3835 and it obviously hasn't fixed everything 😞 I don't know why having the try-catch at a lower level would make any difference though... Seeing your print statements made me realise the lower-level try-catch does help narrow down where the error comes from, which probably just printing the error does not. We can make that better by also printing a stacktrace along with the error (as explained in https://docs.julialang.org/en/v1/manual/stacktraces/#Error-handling). I'll push that to #87, just so you can see what I mean.

[mrhardman]

Besides the error handling issue discussed above, the physics issue here was that of having an instability causing gridscale oscillations that lead to negative densities or temperatures, see https://excalibur-neptune.github.io/Documents/TN-04-3_ConvergenceSolutions2D1VIonDriftKineticEquationWallBoundaryConditionsManufacture.html. There may also have been precompilation issues due to using multiple cores without first precompiling a moment_kinetics.so image. Potentially related issues #137. Closing as stale. Can re-open if similar bugs are reported.