Add W_v output eigenvalues fix

[bford26]

Potential Bug find and Fix.

Comparisons between W_t and reconstructed matrix from eigen values and vectors was not working.

Matrix reconstruction tests were done with this script. (both before and after changes are included in comments of script)
dcon_W_reconstruction_test.py

The changes to dcon/free.f and dcon/dcon_netcdf.f now reflect that 'ev' is defined as explicitly real eigenvalues.

This can be seen from the call to zheev line305, ev needs to be real zheev doc

[jhalpern30]

Question on this - I believe Wv should be Hermitian, but can accumulate numerical errors (like Wp). In vacuum, it is forced Hermitian with lsymz (line 526 of vacuum_vac.f), which is set to true for most of the gpec example inputs. From my understanding, it seems like this fix assumes lsymz = t, correct? But this is technically still a quantity exposed to the user. So what happens with this fix if lsymz = f?

If this breaks, perhaps another option is to use zgeev for the Wv matrix in free.f as well to avoid making any Hermitian assumptions, since this isn't a numerically costly section of the code so zheev vs. zgeev isn't going to make much of a difference in terms of runtime. But open to thoughts on this

[matt-pharr]

I agree with @jhalpern30

[bford26]

So not sure if this fixes the lsymz = false case. (i tested it with lsymz false and got similar output, so that means it's fine? The simple fix might just be using a new variable ev_real ?

The versions i tested reinstate ev as COMPLEX(r8). From the outputs i test, the ev_real addition seems to work best.

  CALL zheev('V','U',mpert,wv,mpert,ev_real,work,lwork,rwork,info)
  ev = CMPLX(ev_real, 0.0_r8, KIND=r8)

vs

  lwork=2*mpert+1
  CALL zgeev('V','V',mpert,wv,mpert,ev,
 $     vl,mpert,vr,mpert,work2,lwork,rwork2,info)
  eindex(1:mpert)=(/(ipert,ipert=1,mpert)/)
  CALL bubble(REAL(ev),eindex,1,mpert)
  tt=ev
  DO ipert=1,mpert
     wv(:,ipert)=vr(:,eindex(mpert+1-ipert))
     ev(ipert)=tt(eindex(mpert+1-ipert))
  ENDDO

Test Script:
test_Wreconstruct_script.py

Test outputs:
w_reconstruct_outputs.txt

[logan-nc]

@jhalpern30 please respond. Id like to merge or close this before releasing v1.6

[jhalpern30]

Sorry all, this got lost in other Github notifications.

I think this is fine - the lsymz = f case is quite niche even if there's a case we're missing here, it would only occur if the Hermitian errors in Wv are large because the vacuum grids were under resolved anyway.

[jhalpern30]

Looks like it just needs to pass the Fortran build checks

[StuartBenjamin]

My agent (extremely powerful) was going through the code-base and made the following remark on the issue this PR is dealing with:

"
Problem: zheev returns real eigenvalues, but ev is declared COMPLEX(r8). The eigenvalues get packed pairwise into the real and imaginary parts of the first half of ev, with stale values in the rest. The garbled values are what's written as W_v_eigenvalue. Also, lwork still holds 2*mpert+1 from the preceding zgeev, but work has 2*mpert-1 elements. #256 fixes DCON's ev, but not lwork, and not RDCON, which has the same code in rdcon/free.f.

Fix (DCON is the same, plus lwork):

!     rdcon/free.f, declarations
      REAL(r8), DIMENSION(mpert) :: evr
!     before the zheev call (replaces it)
      lwork=2*mpert-1
      CALL zheev('V','U',mpert,wv,mpert,evr,work,lwork,rwork,info)
      ev=evr

All other outputs (Δ′, W_t, eigenvectors) are bit-identical on the DIII-D resistive example.
"

I haven't fact checked it but thought I'd share regardless - seems like a few more call-sites need the fix

[matt-pharr]

@StuartBenjamin That does appear to be the case. I will go through all our zheev calls and make sure that this is fixed on this branch.

[matt-pharr]

OK @logan-nc @jhalpern30 @StuartBenjamin this problem is much larger than I thought, actually. I have found a couple other instances of problems with ZHEEV calls besides the one Stuart's agent found. That is going to be a much larger bugfix PR. I recommend merging this. I will open the bugfix PR once I make progress on that, which I will work on tonight.

This also makes me wonder if there are any other instances of LAPACK/LSODE/ZVODE misuse. Since there is no interface for them the compiler doesn't check, unfortunately. I will look into if there is any way to check these things automatically.