PrincetonUniversity/STELLOPT

BEAMS3D: strange behavior when compiling in debug

Closed this issue · 4 comments

kudav commented

When having compiled in debug, e.g with 'make clean_debug', BEAMS3D shows weird behavior. Particularly, GetBCyl in init_vmec does not produce sensible values as it does when compiled with release. One effect is that S_ARR has values of 1.01 inside the Plasma (S<1), but I've also seen it being equal to 4, which is sflxsflxscaleupscaleup=1122.
This did not appear to change when going back on the develop branch a couple of months, so does not seem to be a recent issue.

When running a W7-X VMEC case in debug mode, the following error occurs:

==== backtrace (tid:  87796) ====
 0 0x0000000000016910 __funlockfile()  ???:0
 1 0x00000000009a5cf2 wall_mod_mp_wall_load_txt_()  /raven/u/davku/src/STELLOPT/LIBSTELL/Debug/wall_mod.f90:375
 2 0x0000000000609e05 beams3d_init_()  /raven/u/davku/src/STELLOPT/BEAMS3D/Debug/beams3d_init.f90:424
 3 0x0000000000681efb MAIN__()  /raven/u/davku/src/STELLOPT/BEAMS3D/Debug/../Sources/beams3d_main.f90:44
 4 0x000000000040c922 main()  ???:0
 5 0x000000000003524d __libc_start_main()  ???:0
 6 0x000000000040c83a _start()  /home/abuild/rpmbuild/BUILD/glibc-2.31/csu/../sysdeps/x86_64/start.S:120

This is caused by the following calculation of invDenom:

      DO ik = mystart, myend
         ihit_array(ik) = 0
         DOT00(ik) = V0(ik,1)*V0(ik,1) + V0(ik,2)*V0(ik,2) + V0(ik,3)*V0(ik,3)
         DOT01(ik) = V0(ik,1)*V1(ik,1) + V0(ik,2)*V1(ik,2) + V0(ik,3)*V1(ik,3)
         DOT11(ik) = V1(ik,1)*V1(ik,1) + V1(ik,2)*V1(ik,2) + V1(ik,3)*V1(ik,3)
         d(ik)     = FN(ik,1)*A0(ik,1) + FN(ik,2)*A0(ik,2) + FN(ik,3)*A0(ik,3)
         invDenom(ik) = one / (DOT00(ik)*DOT11(ik) - DOT01(ik)*DOT01(ik))
      END DO

We could bound the denominator eg. by a MIN() statement, skip the calculation for the triangles where this happens or ensure we only include well-defined triangles in the wall model creation. Combining the last two would probably be the most robust and accurate.

What error? Floating point underflow? DOT00 and DOT11 should be positive deffinite. If V1 and V0 are at right angles then DOT01 could underflow but invDenom would still be positive definite.

Sorry, it's a floating-point divide by zero error:
[ravc3321:87756:0:87756] Caught signal 8 (Floating point exception: floating-point divide by zero)

It must be a small triangle issue. DOT00 can equal DOT11 but DOT01 can only equal those if V0=V1 which is no longer a tirangle. So V0 and V1 must be real small.