Removing the 3rd dimension in the u arrays (or determining they are beneficial)

[jhalpern30]

From conversations with @logan-nc:
The regular u matrix is mpert x mpert x 2, which actually represents a mpert x (2 x mpert) matrix representing the mpert solution vectors, each of size 2 * mpert to contain the mpert displacements and their mpert conjugate momenta. It adds a lot of confusion because the third dimension was just differentiates the two components of the column vector, and the matrix dimensions end up not matching what one would expect based on the paper. This is also used in the asymptotics (ua, ca, vmat, mmat, etc.) and makes things very confusing.

The most obvious example: the entire sing_get_ca function would just be

ca = lu(ua) \ u # where ua is 2mpert x 2mpert and u is mpert x 2mpert

but instead we spend 90% of the function constructing and deconstructing the full 2D matrices from the sliced 3D version

# Build temp1
temp1 = zeros(ComplexF64, 2 * intr.mpert, 2 * intr.mpert)
temp1[1:intr.mpert, :] .= ua[:, :, 1]
temp1[intr.mpert+1:2*intr.mpert, :] .= ua[:, :, 2]

# Built temp2
temp2 = zeros(ComplexF64, 2 * intr.mpert, intr.mpert)
temp2[1:intr.mpert, :] .= odet.u[:, :, 1]
temp2[intr.mpert+1:2*intr.mpert, :] .= odet.u[:, :, 2]

# LU factorization and solve
temp2 .= lu(temp1) \ temp2

# Build ca
ca = zeros(ComplexF64, intr.mpert, intr.mpert, 2)
ca[:, :, 1] .= temp2[1:intr.mpert, :]
ca[:, :, 2] .= temp2[intr.mpert+1:2*intr.mpert, :]

My suspicion is because the most efficient way of computing u' = Lu in sing_der is to first compute the displacement derivative, i.e. du[:, :, 1], and then to use that in the computation of the conjugate momenta derivative since it just reappears in the expression for du[:, :, 2] = Gu[:, :, 1] - adjoint(K) - du{:, :, 1], and this makes the code look more clean than slicing from 1:mpert and mpert+1:2*mpert over and over. But in Julia we do these operations via views already so there's no advantage in this case, we'd just be modifying which slice the view operations are taking at the beginning of this function.

If there are no other obvious reasons why slicing the matrix into a third dimension is beneficial that are still applicable to the Julia code, I think this would be a massive improvement to the code since it would make the fundamental matrices in the code directly line up with what is presented in the DCON papers

[logan-nc]

@parkjk do you have any insight into why this is done in the Fortran?

[parkjk]

I am not sure and perhaps what Jacob mentioned (advantages in du) could be just it. Or it could make it clearer when msol>mpert, in the resistive case, e.g. when the secondary fixup is required, but it is not used any more in DCON.

After all, it was a stylistic choice. It was confusing to me as well but was not a big deal to adjust myself to them along with Fortran extension. Now it is different since you are switching to another language and so it makes sense changing the style in a better way. The only issue is that it may cause more confusion to someone who wants to use both sides, but I can't think of anyone except Nik haha;

[parkjk]

Well, @jmlmir369 may care this change and so is cc'ed here.

[logan-nc]

😆 I'll survive if you change it @jhalpern30 !

@matt-pharr, do you see this being an issue down the road for when we add resistive options? I think Jake's point about using Julia views still applies well, right?

[jhalpern30]

Adding this on here (and for my own sake in case I do this a while from now) - a lot of the complicated index logic that goes on in sing_mmat and sing_vmat that occurs when projecting onto the resonant and nonresonant solution subspaces would be drastically simplified without the third dimension.

More specifically, the r1/r2 and n1/n2 are the indices that project onto each subspace (equivalent to multiplying by R and (I - R) in eq. 23 of Glasser's 2023 PoP paper. In 2D, r1 = r2 = [mn, mn + MN] = the diagonal indices of the identity matrix corresponding to the resonant modes, but in 3D, the 2mpert X 2mpert matrix becomes mpert X 2*mpert X 2 with a single index for r1 = mn that now has to be indexed twice for the third dimension for the mn + MNth mode.

This section of the code is already (in my opinion) one of the most difficult to understand, and this would help a lot, especially to generalize it to 3D

[jhalpern30]

Please take a look at this issue above and the existing repo, and (without writing any code yet) discuss the potential advantages and disadvantages of this PR, focusing on readability of the code and any speedups/optimizations that can be performed using this modification. Provide a determination of whether or not this change should be added to the code, and if so, and implementation plan