Bugfix/debug fixes

BEAMS3D/Sources/beams3d_diagnostics.f90

@@ -37,12 +37,10 @@ SUBROUTINE beams3d_diagnostics
       LOGICAL, ALLOCATABLE     :: partmask(:), partmask2(:,:), partmask2t(:,:)
       INTEGER, ALLOCATABLE  :: int_mask(:), int_mask2(:,:)
       INTEGER, ALLOCATABLE  :: dist_func(:,:,:)
-      REAL, ALLOCATABLE     :: real_mask(:),vllaxis(:),vperpaxis(:), nlost(:), norbit(:), tlow(:), thigh(:)
+      REAL, ALLOCATABLE     :: vllaxis(:),vperpaxis(:), nlost(:), norbit(:), tlow(:), thigh(:)
       REAL, ALLOCATABLE     :: help3d(:,:,:)
-#if defined(MPI_OPT)
-      INTEGER :: mystart, mypace
-      REAL(rprec), ALLOCATABLE :: buffer_mast(:,:), buffer_slav(:,:)
-#endif
+      REAL, ALLOCATABLE     :: dbl_help(:)
+      INTEGER :: mystart
 !-----------------------------------------------------------------------
 !     Begin Subroutine
 !-----------------------------------------------------------------------
@@ -67,12 +65,16 @@ SUBROUTINE beams3d_diagnostics
       IF (ALLOCATED(norbit)) DEALLOCATE(norbit)
       IF (ALLOCATED(tlow)) DEALLOCATE(tlow)
       IF (ALLOCATED(thigh)) DEALLOCATE(thigh)
+      IF (ALLOCATED(dbl_help)) DEALLOCATE(dbl_help)
+      IF (ALLOCATED(int_mask)) DEALLOCATE(int_mask)
       ALLOCATE(shine_through(nbeams))
       ALLOCATE(shine_port(nbeams))
       ALLOCATE(nlost(nbeams))
       ALLOCATE(norbit(nbeams))
       ALLOCATE(tlow(nbeams))
       ALLOCATE(thigh(nbeams))
+      ALLOCATE(dbl_help(mystart:myend))
+      ALLOCATE(int_mask(mystart:myend))
 #if defined(MPI_OPT)
       CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)
       IF (ierr_mpi /= 0) CALL handle_err(MPI_BARRIER_ERR, 'beams3d_follow', ierr_mpi)
@@ -85,14 +87,20 @@ SUBROUTINE beams3d_diagnostics
 
       ! Calculate shinethrough and loss
       shine_through = 0
+      dbl_help(mystart:myend) = t_end(mystart:myend)
+      int_mask(mystart:myend) = beam(mystart:myend)
+
+
       DO i = 1, nbeams
          norbit(i)         =      SUM(weight(mystart:myend), MASK = (end_state(mystart:myend) == 0 .and. (beam(mystart:myend)==i)))
          nlost(i)          =      SUM(weight(mystart:myend), MASK = (end_state(mystart:myend) == 2 .and. (beam(mystart:myend)==i)))
          shine_through(i)  = 100.*SUM(weight(mystart:myend), MASK = (end_state(mystart:myend) == 3 .and. (beam(mystart:myend)==i)))/SUM(weight,MASK=(beam==i))
          shine_port(i)     = 100.*SUM(weight(mystart:myend), MASK = (end_state(mystart:myend) == 4 .and. (beam(mystart:myend)==i)))/SUM(weight,MASK=(beam==i))
-         tlow(i)           = MINVAL(t_end(mystart:myend), MASK = (beam(mystart:myend)==i))
-         thigh(i)          = MAXVAL(t_end(mystart:myend), MASK = (beam(mystart:myend)==i))
+         tlow(i)           = MINVAL(dbl_help, MASK = (int_mask==i))
+         thigh(i)          = MAXVAL(dbl_help, MASK = (int_mask==i))
       END DO
+      DEALLOCATE(dbl_help)
+      DEALLOCATE(int_mask)
 
 #if defined(MPI_OPT)
       IF (myworkid == master) THEN

BEAMS3D/Sources/beams3d_follow.f90

@@ -307,7 +307,7 @@ SUBROUTINE beams3d_follow
 #endif
 
     ! Adjust T_END back to values of T_last
-    t_end(mystart:myend) = t_last(mystart:myend)
+    t_end(mystart_save:myend_save) = t_last(mystart_save:myend_save)
     IF (ALLOCATED(t_last)) DEALLOCATE(t_last)
 
     RETURN

BEAMS3D/Sources/beams3d_init_vmec.f90

@@ -307,6 +307,8 @@ SUBROUTINE beams3d_init_vmec
       CALL MPI_CALC_MYRANGE(MPI_COMM_LOCAL, 1, nr*nphi*nz, mystart, myend)
       ALLOCATE(lsmooth(mystart:myend))
       lsmooth = .false.
+
+      uflx = 0.0
 
       IF (mylocalid == mylocalmaster) THEN
          TE = 0; NE = 0; TI=0; S_ARR=scaleup*scaleup; U_ARR=0; POT_ARR=0; ZEFF_ARR = 1;

BEAMS3D/Sources/beams3d_physics_mod.f90

@@ -281,7 +281,7 @@ SUBROUTINE beams3d_physics_gc(t, q)
             CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
                             hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),&
                             S4D(1,1,1,1),nr,nphi,nz)
-            s_temp = fval(1)
+            s_temp = max(fval(1),zero)
 
             !-----------------------------------------------------------
             !  Helpers
@@ -488,7 +488,7 @@ SUBROUTINE beams3d_physics_fo(t, q)
             CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
                             hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),&
                             S4D(1,1,1,1),nr,nphi,nz)
-            s_temp = fval(1)
+            s_temp = max(fval(1),zero)
             CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
                             hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),&
                             BR4D(1,1,1,1),nr,nphi,nz)
@@ -867,7 +867,7 @@ SUBROUTINE beams3d_follow_neut(t, q)
             DO l = 1, num_depo
                nelocal(l)  = MAX(MIN(nelocal(l),1E21),1E18)
                telocal(l)  = MAX(MIN(telocal(l),energy(l)*0.5),energy(l)*0.01)
-               ni_in = nilocal(:,l)
+               ni_in = MAX(MIN(nilocal(:,l),1E21),1E18)
                CALL suzuki_sigma(NION,energy(l),nelocal(l),telocal(l),ni_in,A_in,Z_in,tau_inv(l))
             END DO
             tau_inv = tau_inv*nelocal*ABS(q(4))*1E-4 !cm^2 to m^2 for sigma
@@ -1943,8 +1943,8 @@ SUBROUTINE beams3d_suv2rzp(s,u,v,r_out,z_out,phi_out)
          x0 = s * COS(u)
          y0 = s * SIN(U)
 
-         fnorm = x0*x0+y0*y0
-         fnorm = MIN(1./fnorm,1E5)
+         fnorm = MAX(x0*x0+y0*y0,1E-5)
+         fnorm = 1./fnorm
          n = 1
 
          ! Loop Basically a NEWTON's METHOD

LIBSTELL/Sources/Lsode/dlsode_aux2.f

@@ -359,8 +359,9 @@ DOUBLE PRECISION FUNCTION DNRM23 (N, DX, INCX)
 C
       DATA CUTLO, CUTHI /8.232D-11,  1.304D19/
 C***FIRST EXECUTABLE STATEMENT  DNRM2
+      DNRM23 = ZERO
       IF (N .GT. 0) GO TO 10
-         DNRM2  = ZERO
+         DNRM23  = ZERO
          GO TO 300
 C
    10 ASSIGN 30 TO NEXT
@@ -423,7 +424,7 @@ DOUBLE PRECISION FUNCTION DNRM23 (N, DX, INCX)
       DO 95 J = I,NN,INCX
       IF (ABS(DX(J)) .GE. HITEST) GO TO 100
    95    SUM = SUM + DX(J)**2
-      DNRM2 = SQRT(SUM)
+      DNRM23 = SQRT(SUM)
       GO TO 300
 C
   200 CONTINUE
@@ -434,7 +435,7 @@ DOUBLE PRECISION FUNCTION DNRM23 (N, DX, INCX)
 C
 C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
 C
-      DNRM2 = XMAX * SQRT(SUM)
+      DNRM23 = XMAX * SQRT(SUM)
   300 CONTINUE
       RETURN
       END

LIBSTELL/Sources/Modules/ez_hdf5.f90

@@ -232,7 +232,7 @@ SUBROUTINE write_scalar_hdf5(file_id,var,ierr,BOOVAR,INTVAR,FLTVAR,DBLVAR,STRVAR
       INTEGER        :: arank = 1
       INTEGER(HID_T) :: dset_id
       INTEGER(HID_T) :: type
-      INTEGER(HID_T) :: strlen
+      INTEGER(HSIZE_T) :: strlen
       INTEGER(HID_T) :: attr_id  
       INTEGER(HID_T) :: dspace_id
       INTEGER(HID_T) :: aspace_id
@@ -329,7 +329,7 @@ SUBROUTINE write_vector_hdf5(file_id,var,n,ier,BOOVAR,INTVAR,FLTVAR,DBLVAR,STRVA
       INTEGER        :: arank = 1
       INTEGER        :: boo_temp(n)
       INTEGER(HID_T) :: type
-      INTEGER(HID_T) :: strlen
+      INTEGER(HSIZE_T) :: strlen
       INTEGER(HID_T) :: dset_id
       INTEGER(HID_T) :: attr_id  
       INTEGER(HID_T) :: dspace_id

LIBSTELL/Sources/Modules/vmec_utils.f

@@ -90,13 +90,14 @@ SUBROUTINE GetBcyl_WOUT(R1, Phi, Z1, Br, Bphi, Bz,
       CALL cyl2flx(rzl_local, r_cyl, c_flx, ns_w, ntor_w, mpol_w, 
      1     ntmax_w, lthreed_w, lasym_w, info_loc, nfe, fmin, 
      2     RU=Ru1, ZU=Zu1, RV=Rv1, ZV=Zv1, RS=Rs1, ZS=Zs1)
-      Rv1 = nfp*Rv1;  Zv1 = nfp*Zv1
 
       IF (info_loc.eq.-1 .and. (fmin .le. fmin_acceptable)) info_loc = 0
 
       IF (PRESENT(info)) info = info_loc
       IF (info_loc .ne. 0) RETURN
 
+      Rv1 = nfp*Rv1;  Zv1 = nfp*Zv1
+
       IF (PRESENT(sflx)) sflx = c_flx(1)  
       IF (PRESENT(uflx)) uflx = c_flx(2)
 
@@ -733,24 +734,43 @@ SUBROUTINE flx2cyl(rzl_array, c_flux, r_cyl, ns, ntor,
       whi  = (si - slo)/hs1 ! x
       wlo  = one - whi ! (1-x)
 
-      ! Odd modes we include a factor of sqrt(s)=rho
-      wlo_odd = wlo*rho/rholo
-      whi_odd = whi*rho/rhohi
-
       ! Derivative values
       dwlo = -DBLE(ns-1)
       dwhi =  DBLE(ns-1)
-      dwlo_odd = -(3*rho-shi/rho)/(2*hs1*rholo)
-      dwhi_odd =  (3*rho-slo/rho)/(2*hs1*rhohi)
-
+
       ! Adjust near axis
       IF (jslo .eq. 1) THEN
          wlo_odd = 0
          whi_odd = rho/rhohi
          dwlo_odd = 0
          dwhi_odd = 0.5/(rho*rhohi)
+      ELSE ! Otherwise
+         ! Odd modes we include a factor of sqrt(s)=rho
+         wlo_odd = wlo*rho/rholo
+         whi_odd = whi*rho/rhohi
+         dwlo_odd = -(3*rho-shi/rho)/(2*hs1*rholo)
+         dwhi_odd =  (3*rho-slo/rho)/(2*hs1*rhohi)
       END IF
 
+!     Old way with floating point exception when rholo=0
+!      ! Odd modes we include a factor of sqrt(s)=rho
+!      wlo_odd = wlo*rho/rholo
+!      whi_odd = whi*rho/rhohi
+!
+!      ! Derivative values
+!      dwlo = -DBLE(ns-1)
+!      dwhi =  DBLE(ns-1)
+!      dwlo_odd = -(3*rho-shi/rho)/(2*hs1*rholo)
+!      dwhi_odd =  (3*rho-slo/rho)/(2*hs1*rhohi)
+!
+!      ! Adjust near axis
+!      IF (jslo .eq. 1) THEN
+!         wlo_odd = 0
+!         whi_odd = rho/rhohi
+!         dwlo_odd = 0
+!         dwhi_odd = 0.5/(rho*rhohi)
+!      END IF
+
       mpol1 = mpol-1
 
       wmins(:,0:mpol1:2) = wlo
@@ -1138,7 +1158,7 @@ SUBROUTINE newt2d(xc_opt, fmin, ftol, nfe, nvar, iflag)
       REAL(rprec) :: c_flx(3), r_cyl_out(3), fvec(nvar), sflux, 
      1               uflux, eps0, eps, epu, xc_min(2), factor
       REAL(rprec) :: x0(3), xs(3), xu(3), dels, delu, tau, fmin0,
-     1               ru1, zu1, edge_value, snew, rs1, zs1
+     1               ru1, zu1, edge_value, snew, rs1, zs1, z_small
 C-----------------------------------------------
 !
 !     INPUT/OUTPUT:
@@ -1165,6 +1185,7 @@ SUBROUTINE newt2d(xc_opt, fmin, ftol, nfe, nvar, iflag)
 !     without yielding a lower value of F.
 !
       iflag = -1      
+      z_small=TINY(1.0_rprec)
       eps0 = SQRT(EPSILON(eps))
       xc_min = xc_opt
 
@@ -1235,7 +1256,8 @@ SUBROUTINE newt2d(xc_opt, fmin, ftol, nfe, nvar, iflag)
 
 !           NEWTON STEP
             tau = xu(1)*xs(3) - xu(3)*xs(1)
-            IF (ABS(tau) .le. ABS(eps)*r_target**2) THEN
+            !IF (ABS(tau) .le. ABS(eps)*r_target**2) THEN
+            IF (ABS(tau) .le. ABS(z_small)*r_target**2) THEN
                iflag = -2
                EXIT
             END IF

LIBSTELL/Sources/Modules/wall_mod.f90

@@ -1083,17 +1083,12 @@ SUBROUTINE collide_double(x0,y0,z0,x1,y1,z1,xw,yw,zw,lhit)
       r0(1) = x0
       r0(2) = y0
       r0(3) = z0
-      DO i=1,3
-         outlow(i) = .false.
-         outhigh(i) = .false.
-         tcomp(i) = zero
-         ! check direction line to determine in which direction to step through the blocks
-         IF (dr(i) < zero) THEN
-            step(i) = -1
-         ELSE
-            step(i) = 1
-         END IF
-      END DO
+      outlow  = .FALSE.
+      outhigh = .FALSE.
+      tcomp   = zero
+      ! check direction line to determine in which direction to step through the blocks
+      step = 1
+      WHERE(dr < zero) step = -1
 
       k1 = 1; k2 = wall%nblocks
       b_found = -1
@@ -1171,6 +1166,7 @@ SUBROUTINE collide_double(x0,y0,z0,x1,y1,z1,xw,yw,zw,lhit)
             alphal = FN(ik,1)*dr(1) + FN(ik,2)*dr(2) + FN(ik,3)*dr(3)
             betal = FN(ik,1)*r0(1) + FN(ik,2)*r0(2) + FN(ik,3)*r0(3)
             ! tloc indicated when hit. If hit between r0 and r1, tloc between 0 and 1
+            IF (alphal == zero) CYCLE
             tloc = (d(ik)-betal)/alphal
             IF (tloc > one) CYCLE
             IF (tloc <= zero) CYCLE
@@ -1193,14 +1189,18 @@ SUBROUTINE collide_double(x0,y0,z0,x1,y1,z1,xw,yw,zw,lhit)
          END DO
 
          ! Check when the line will leave the block
-         DO i=1,3
-            IF (step(i) .eq. 1) THEN
-               tDelta(i) = (b%rmax(i) - r0(i)) / dr(i)
-            ELSE
-               tDelta(i) = (b%rmin(i) - r0(i)) / dr(i)
-            END IF
-         END DO
-
+         WHERE(step == 1)
+            tDelta = (b%rmax - r0)
+         ELSEWHERE
+            tDelta = (b%rmin - r0)
+         END WHERE
+         ! Handle division by zero
+         WHERE (dr == zero)
+            tDelta = 1.0E20
+         ELSEWHERE
+            tDelta = tDelta / dr
+         END WHERE
+
          ! Check if leaves block before hits. 
          ! Compare with tcomp to make sure that ray always moves in positive time direction 
          IF (ANY(tDelta < tmin .and. tDelta > tcomp)) THEN
@@ -1735,9 +1735,26 @@ SUBROUTINE LINE_BOX_INTERSECTION(x1i, y1i, z1i, x2i, y2i, z2i, xmin, ymin, zmin,
       z2 = MAX(z1i,z2i)
 
       ! Helpers
-      dx = one/(x2-x1)
-      dy = one/(y2-y1)
-      dz = one/(z2-z1)
+      dx = x2-x1
+      dy = y2-y1
+      dz = z2-z1
+
+      ! IF statemnt for parallel lines
+      IF (dx == zero) THEN
+         dx = 1.0D20
+      ELSE
+         dx = one/dx
+      END IF
+      IF (dy == zero) THEN
+         dy = 1.0D20
+      ELSE
+         dy = one/dy
+      END IF
+      IF (dz == zero) THEN
+         dz = 1.0D20
+      ELSE
+         dz = one/dz
+      ENDIF
 
       ! Calculate the minimum and maximum values of t for each axis
       tmin = zero

SHARE/make_cobra.inc

@@ -14,10 +14,12 @@
 #######################################################################
 #            Define Compiler Flags
 #######################################################################
-  FLAGS_R = -I$(MKL_HOME)/include -O2 -fp-model strict -ip \
-            -assume noold_unit_star -xCORE-AVX512 -qopt-zmm-usage=high -g -traceback
-  FLAGS_D = -I$(MKL_HOME)/include -O2 -fp-model strict -ip \
-            -assume noold_unit_star -g -traceback -xCORE-AVX512 -qopt-zmm-usage=high
+  FLAGS_R = -I$(MKL_HOME)/include -O2 -assume noold_unit_star \
+            -xCORE-AVX512 -qopt-zmm-usage=high -fp-model strict -ip 
+  FLAGS_D = -I$(MKL_HOME)/include -O0 -assume noold_unit_star \
+            -xCORE-AVX512 -qopt-zmm-usage=high -fp-model strict -ip \
+            -check all,noarg_temp_created -ftrapuv -g -traceback
+
   LIBS    = -Wl,-rpath,$(MKL_HOME)/lib/intel64 \
             -L$(MKL_HOME)/lib/intel64 -lmkl_scalapack_lp64 \
             -lmkl_intel_lp64 -lmkl_core -lmkl_sequential \

SHARE/make_raven.inc

@@ -14,10 +14,13 @@
 #######################################################################
 #            Define Compiler Flags
 #######################################################################
-  FLAGS_R = -I${MKLROOT}/include/intel64/lp64 -I$(MKL_HOME)/include -O2 -fp-model strict -ip \
-            -assume noold_unit_star -xCORE-AVX512 -qopt-zmm-usage=high
-  FLAGS_D = -I${MKLROOT}/include/intel64/lp64 -I$(MKL_HOME)/include -O2 -fp-model strict -ip \
-            -assume noold_unit_star -g -traceback -check all -ftrapuv -xCORE-AVX512 -qopt-zmm-usage=high
+  FLAGS_R = -I${MKLROOT}/include/intel64/lp64 -I$(MKL_HOME)/include \
+            -O2 -assume noold_unit_star \
+            -xCORE-AVX512 -qopt-zmm-usage=high -fp-model strict -ip
+  FLAGS_D = -I${MKLROOT}/include/intel64/lp64 -I$(MKL_HOME)/include \
+            -O0 -assume noold_unit_star \
+            -xCORE-AVX512 -qopt-zmm-usage=high -fp-model strict -ip \
+            -check all,noarg_temp_created -ftrapuv -g -traceback
   LIBS    =  ${MKLROOT}/lib/intel64/libmkl_blas95_lp64.a \
              ${MKLROOT}/lib/intel64/libmkl_lapack95_lp64.a \
              ${MKLROOT}/lib/intel64/libmkl_scalapack_lp64.a \

WALL_ACCELERATE/Sources/wall_accelerate_main.f90

@@ -125,8 +125,7 @@ PROGRAM WALL_ACCELERATE
 
 #if defined(MPI_OPT)
          CALL wall_load_txt(vessel_string, ier, lverb, MPI_COMM_WALLACC)  
-         WRITE(6,'(A,I2.2)') 'WALL ERRORCODE: ',ier
-         CALL MPI_BARRIER(MPI_COMM_WALLACC,ierr_mpi)
+         WRITE(6,'(A,I5)') 'WALL ERRORCODE: ',ier
          IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_BARRIER_ERR,'wall_acc_main',ierr_mpi)
 #else
          CALL wall_load_txt(vessel_string, ier, lverb)