Infer dimension of arrays from input to struct. Add global dimension constants.

[mrhardman]

Infer number of field, gyroaveraged field, moment dimensions from call to struct, rather than hardcode size. PR to address spirit of #46, may close #46. Where obvious, we infer moments and field array dimensions from the arrays themselves.

This could be extended by creating global dimension variables, as discussed in issue #46, but this would have to address how to use these variables in code like,

moment_kinetics/moment_kinetics/src/time_advance.jl

Lines 96 to 217 in e99da09

	struct scratch_dummy_arrays
	dummy_s::Array{mk_float,1}
	dummy_sr::Array{mk_float,2}
	dummy_vpavperp::Array{mk_float,2}
	dummy_zrs::MPISharedArray{mk_float,3}
	dummy_zrsn::MPISharedArray{mk_float,3}
	#buffer arrays for MPI
	buffer_z_1::MPISharedArray{mk_float,1}
	buffer_z_2::MPISharedArray{mk_float,1}
	buffer_z_3::MPISharedArray{mk_float,1}
	buffer_z_4::MPISharedArray{mk_float,1}
	buffer_r_1::MPISharedArray{mk_float,1}
	buffer_r_2::MPISharedArray{mk_float,1}
	buffer_r_3::MPISharedArray{mk_float,1}
	buffer_r_4::MPISharedArray{mk_float,1}
	buffer_zs_1::MPISharedArray{mk_float,2}
	buffer_zs_2::MPISharedArray{mk_float,2}
	buffer_zs_3::MPISharedArray{mk_float,2}
	buffer_zs_4::MPISharedArray{mk_float,2}
	buffer_zsn_1::MPISharedArray{mk_float,2}
	buffer_zsn_2::MPISharedArray{mk_float,2}
	buffer_zsn_3::MPISharedArray{mk_float,2}
	buffer_zsn_4::MPISharedArray{mk_float,2}
	buffer_rs_1::MPISharedArray{mk_float,2}
	buffer_rs_2::MPISharedArray{mk_float,2}
	buffer_rs_3::MPISharedArray{mk_float,2}
	buffer_rs_4::MPISharedArray{mk_float,2}
	buffer_rs_5::MPISharedArray{mk_float,2}
	buffer_rs_6::MPISharedArray{mk_float,2}
	buffer_rsn_1::MPISharedArray{mk_float,2}
	buffer_rsn_2::MPISharedArray{mk_float,2}
	buffer_rsn_3::MPISharedArray{mk_float,2}
	buffer_rsn_4::MPISharedArray{mk_float,2}
	buffer_rsn_5::MPISharedArray{mk_float,2}
	buffer_rsn_6::MPISharedArray{mk_float,2}
	buffer_zrs_1::MPISharedArray{mk_float,3}
	buffer_zrs_2::MPISharedArray{mk_float,3}
	buffer_zrs_3::MPISharedArray{mk_float,3}
	buffer_vpavperpzs_1::MPISharedArray{mk_float,4}
	buffer_vpavperpzs_2::MPISharedArray{mk_float,4}
	buffer_vpavperpzs_3::MPISharedArray{mk_float,4}
	buffer_vpavperpzs_4::MPISharedArray{mk_float,4}
	buffer_vpavperpzs_5::MPISharedArray{mk_float,4}
	buffer_vpavperpzs_6::MPISharedArray{mk_float,4}
	buffer_vpavperprs_1::MPISharedArray{mk_float,4}
	buffer_vpavperprs_2::MPISharedArray{mk_float,4}
	buffer_vpavperprs_3::MPISharedArray{mk_float,4}
	buffer_vpavperprs_4::MPISharedArray{mk_float,4}
	buffer_vpavperprs_5::MPISharedArray{mk_float,4}
	buffer_vpavperprs_6::MPISharedArray{mk_float,4}
	# buffer to hold derivative after MPI communicates
	# needs to be shared memory
	buffer_vpavperpzrs_1::MPISharedArray{mk_float,5}
	buffer_vpavperpzrs_2::MPISharedArray{mk_float,5}
	# buffers to hold moment quantities for implicit solves
	implicit_buffer_z_1::MPISharedArray{mk_float,1}
	implicit_buffer_z_2::MPISharedArray{mk_float,1}
	implicit_buffer_z_3::MPISharedArray{mk_float,1}
	implicit_buffer_z_4::MPISharedArray{mk_float,1}
	implicit_buffer_z_5::MPISharedArray{mk_float,1}
	implicit_buffer_z_6::MPISharedArray{mk_float,1}
	# buffers to hold electron for implicit solves
	implicit_buffer_vpavperpz_1::MPISharedArray{mk_float,3}
	implicit_buffer_vpavperpz_2::MPISharedArray{mk_float,3}
	implicit_buffer_vpavperpz_3::MPISharedArray{mk_float,3}
	implicit_buffer_vpavperpz_4::MPISharedArray{mk_float,3}
	implicit_buffer_vpavperpz_5::MPISharedArray{mk_float,3}
	implicit_buffer_vpavperpz_6::MPISharedArray{mk_float,3}
	# buffers to hold ion pdf for implicit solves
	implicit_buffer_vpavperpzrs_1::MPISharedArray{mk_float,5}
	implicit_buffer_vpavperpzrs_2::MPISharedArray{mk_float,5}
	implicit_buffer_vpavperpzrs_3::MPISharedArray{mk_float,5}
	implicit_buffer_vpavperpzrs_4::MPISharedArray{mk_float,5}
	implicit_buffer_vpavperpzrs_5::MPISharedArray{mk_float,5}
	implicit_buffer_vpavperpzrs_6::MPISharedArray{mk_float,5}
	buffer_vzvrvzetazsn_1::MPISharedArray{mk_float,5}
	buffer_vzvrvzetazsn_2::MPISharedArray{mk_float,5}
	buffer_vzvrvzetazsn_3::MPISharedArray{mk_float,5}
	buffer_vzvrvzetazsn_4::MPISharedArray{mk_float,5}
	buffer_vzvrvzetazsn_5::MPISharedArray{mk_float,5}
	buffer_vzvrvzetazsn_6::MPISharedArray{mk_float,5}
	buffer_vzvrvzetarsn_1::MPISharedArray{mk_float,5}
	buffer_vzvrvzetarsn_2::MPISharedArray{mk_float,5}
	buffer_vzvrvzetarsn_3::MPISharedArray{mk_float,5}
	buffer_vzvrvzetarsn_4::MPISharedArray{mk_float,5}
	buffer_vzvrvzetarsn_5::MPISharedArray{mk_float,5}
	buffer_vzvrvzetarsn_6::MPISharedArray{mk_float,5}
	# buffer to hold derivative after MPI communicates
	# needs to be shared memory
	buffer_vzvrvzetazrsn_1::MPISharedArray{mk_float,6}
	buffer_vzvrvzetazrsn_2::MPISharedArray{mk_float,6}
	buffer_vpavperp_1::MPISharedArray{mk_float,2}
	buffer_vpavperp_2::MPISharedArray{mk_float,2}
	buffer_vpavperp_3::MPISharedArray{mk_float,2}
	buffer_vpavperpzr_1::MPISharedArray{mk_float,4}
	buffer_vpavperpzr_2::MPISharedArray{mk_float,4}
	buffer_vpavperpzr_3::MPISharedArray{mk_float,4}
	buffer_vpavperpzr_4::MPISharedArray{mk_float,4}
	buffer_vpavperpzr_5::MPISharedArray{mk_float,4}
	buffer_vpavperpzr_6::MPISharedArray{mk_float,4}
	buffer_vpavperpr_1::MPISharedArray{mk_float,3}
	buffer_vpavperpr_2::MPISharedArray{mk_float,3}
	buffer_vpavperpr_3::MPISharedArray{mk_float,3}
	buffer_vpavperpr_4::MPISharedArray{mk_float,3}
	buffer_vpavperpr_5::MPISharedArray{mk_float,3}
	buffer_vpavperpr_6::MPISharedArray{mk_float,3}
	int_buffer_rs_1::MPISharedArray{mk_int,2}
	int_buffer_rs_2::MPISharedArray{mk_int,2}
	end

and

moment_kinetics/moment_kinetics/src/time_advance.jl

Lines 219 to 229 in e99da09

	struct advect_object_struct
	vpa_advect::Vector{advection_info{4,5}}
	vperp_advect::Vector{advection_info{4,5}}
	z_advect::Vector{advection_info{4,5}}
	r_advect::Vector{advection_info{4,5}}
	electron_z_advect::Vector{advection_info{4,5}}
	electron_vpa_advect::Vector{advection_info{4,5}}
	neutral_z_advect::Vector{advection_info{5,6}}
	neutral_r_advect::Vector{advection_info{5,6}}
	neutral_vz_advect::Vector{advection_info{5,6}}
	end

where specific structures are created based on the known dimensionality of the code. Making this arbitrary might be more painful than simply upgrading the code everywhere to 3D. Remaining hard-coded dimensional code in moment_kinetics_structs.jl

moment_kinetics/moment_kinetics/src/moment_kinetics_structs.jl

Lines 366 to 373 in e99da09

	struct boundary_distributions_struct
	# knudsen cosine distribution for imposing the neutral wall boundary condition
	knudsen::MPISharedArray{mk_float,3}
	# ion particle r boundary values (vpa,vperp,z,r,s)
	pdf_rboundary_ion::MPISharedArray{mk_float,5}
	# neutral particle r boundary values (vz,vr,vzeta,z,r,s)
	pdf_rboundary_neutral::MPISharedArray{mk_float,6}
	end

and

moment_kinetics/moment_kinetics/src/moment_kinetics_structs.jl

Lines 335 to 342 in e99da09

	struct pdf_struct
	#ion particles: s + r + z + vperp + vpa
	ion::pdf_substruct{5}
	# electron particles: r + z + vperp + vpa
	electron::Union{electron_pdf_substruct{4},Nothing}
	#neutral particles: s + r + z + vzeta + vr + vz
	neutral::pdf_substruct{6}
	end

can be more easily generalised, but the value of doing so might be limited if the rest of the code does not benefit from similar clarifications of meaning of number of dimensions.

[johnomotani]

Why not create the global constants? I get it might be a lot of work to consistently use them everywhere right away, but if they exist we can build them in gradually where it makes sense. Given the current design of the code, arrays have to have a hard-coded number of dimensions - and I think that would be very hard to change - so we might as well be explicit. For example, we could then also use the values in function signatures to verify that we are doing the right thing, e.g.

function foo(ion_moment::AbstractArray{mk_float,ion_moment_ndim})
    ...
end

Just on the variable names, I'd prefer to use something like ndim rather than n to avoid ambiguity. n_field sounds like it should be the number of fields, etc. I'd suggest something like field_ndim, ion_moment_ndim, etc.

[mrhardman]

Why not create the global constants?

I agree this could be nice for the instances which represent the full pdf, but I noticed that many instances (e.g., the buffer arrays) we are defining some slice of the pdf. If we slice 1, 2, 3, 4 ... dimensions then we just end up with ndim_pdf - 1, etc, which seems just as bad. What would you do about this?

Just on the variable names, I'd prefer to use something like ndim rather than n to avoid ambiguity. n_field sounds like it should be the number of fields, etc. I'd suggest something like field_ndim, ion_moment_ndim, etc.

Here I followed the existing convention for dimension variables in moment_kinetics_structs.jl, we can change this, but this means changing more than what my commits above touch. Do you agree with that?

[johnomotani]

Why not create the global constants?

I agree this could be nice for the instances which represent the full pdf, but I noticed that many instances (e.g., the buffer arrays) we are defining some slice of the pdf. If we slice 1, 2, 3, 4 ... dimensions then we just end up with ndim_pdf - 1, etc, which seems just as bad. What would you do about this?

How many is 'many'? I think it depends on the case, and there may not be a perfect answer. My gut reaction is to say define just the full dimension sizes ion_pdf_ndim, etc. as global constants. Where we need slices, do whatever is convenient in that place. The buffer arrays I would probably not bother to touch for now - as they are explicitly named for the dimensions included, they might as well have hard-coded dimension sizes. If they are refactored to things like 'slice of a moment at a z-boundary', there I might just leave the dimension sizes as type parameters because you can't do arithmetic in type parameters - e.g. buffer::MPISharedArray{mk_float, ion_moment_ndim-1} is not valid - and it doesn't seem worth defining global constants for each slice. For the advection_structs I'd suggest a similar thing - where the number of dimensions is one of the 'standard' ones, use the global constant, but where it is ion_moment_ndim+1, etc., just leave as an arbitrary type parameter.

For function arguments, the workarounds for the restrictions on type parameters are less awful, e.g.

function foo(ion_moment_slice::MPISharedArray{mk_float,N}) where {N}
    @boundscheck N == ion_moment_ndim - 1 || throw DimensionMismatch("wrong array size")
    ...
end

Just on the variable names, I'd prefer to use something like ndim rather than n to avoid ambiguity. n_field sounds like it should be the number of fields, etc. I'd suggest something like field_ndim, ion_moment_ndim, etc.

Here I followed the existing convention for dimension variables in moment_kinetics_structs.jl, we can change this, but this means changing more than what my commits above touch. Do you agree with that?

Yes, I'd be in favour of more descriptive names - I don't think I chose the original ones well. The originals are not too bad when they're just type parameters within a single struct, where it's fairly clear from the context what they mean, but as global constants they would be pretty awful names.

[mrhardman]

OK, I agree with your comments. I have made an initial set of name choices (they are a cyclic permutation of yours -- I made the choices before reading your last comment, happy to change to another cyclic permutation), but I haven't started using them outside of moment_kinetics_structs.jl. Can you take a look and comment if you want different choices?

[johnomotani]

Names are fine 👍
Need to remove some redundant template parameters - see comments below.

[mrhardman]

Names are fine 👍 Need to remove some redundant template parameters - see comments below.

Right -- so you want to move away from inferred types entirely where possible and just use the global variables? If so, that's fine by me.

[johnomotani]

Right -- so you want to move away from inferred types entirely where possible and just use the global variables? If so, that's fine by me.

Yes. I think it's nicest to minimise the template parameters where possible, so where the global parameters exist we should use them.

[mrhardman]

Right -- so you want to move away from inferred types entirely where possible and just use the global variables? If so, that's fine by me.

Yes. I think it's nicest to minimise the template parameters where possible, so where the global parameters exist we should use them.

The latest version should implement all your comments.