From 12dddc53cb3f361625ccc63349231d9da7acd52a Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Fri, 12 Jan 2024 18:30:58 +0100 Subject: [PATCH 01/16] LIBSTELL: Type fix for ez_hdf5 --- LIBSTELL/Sources/Modules/ez_hdf5.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/LIBSTELL/Sources/Modules/ez_hdf5.f90 b/LIBSTELL/Sources/Modules/ez_hdf5.f90 index 3077daee8..650650b37 100644 --- a/LIBSTELL/Sources/Modules/ez_hdf5.f90 +++ b/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 From de53d02f8ef4c5ac17778e3d3f0b1216797c7882 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Fri, 12 Jan 2024 18:45:21 +0100 Subject: [PATCH 02/16] LIBSTELL: fixed DNRM23 not assigned bug. --- LIBSTELL/Sources/Lsode/dlsode_aux2.f | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/LIBSTELL/Sources/Lsode/dlsode_aux2.f b/LIBSTELL/Sources/Lsode/dlsode_aux2.f index 858b15684..e6e99b777 100644 --- a/LIBSTELL/Sources/Lsode/dlsode_aux2.f +++ b/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 From e95ebfcc9fdc2f2468d6941818861d8998dd43f8 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 08:44:40 +0100 Subject: [PATCH 03/16] LIBSTELL: Fixed possible divide by zero error --- LIBSTELL/Sources/Modules/vmec_utils.f | 40 +++++++++++++++++++++------ 1 file changed, 31 insertions(+), 9 deletions(-) diff --git a/LIBSTELL/Sources/Modules/vmec_utils.f b/LIBSTELL/Sources/Modules/vmec_utils.f index 9f8284e74..9d3b7906f 100644 --- a/LIBSTELL/Sources/Modules/vmec_utils.f +++ b/LIBSTELL/Sources/Modules/vmec_utils.f @@ -90,6 +90,7 @@ 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) + PRINT *,nfp,c_flx,r_cyl,Rv1,Zv1 Rv1 = nfp*Rv1; Zv1 = nfp*Zv1 IF (info_loc.eq.-1 .and. (fmin .le. fmin_acceptable)) info_loc = 0 @@ -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 From b379742ca4f06b2312b8ff7278dabd5e30ccf285 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 08:57:19 +0100 Subject: [PATCH 04/16] LIBSTELL: Fixed an issue of a floating point exception when error is caught. --- LIBSTELL/Sources/Modules/vmec_utils.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/LIBSTELL/Sources/Modules/vmec_utils.f b/LIBSTELL/Sources/Modules/vmec_utils.f index 9d3b7906f..d361f0f10 100644 --- a/LIBSTELL/Sources/Modules/vmec_utils.f +++ b/LIBSTELL/Sources/Modules/vmec_utils.f @@ -90,14 +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) - PRINT *,nfp,c_flx,r_cyl,Rv1,Zv1 - 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) From 6d842ca534a35bf00d215504915d58c63d6af29c Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 09:04:39 +0100 Subject: [PATCH 05/16] BEAMS3D: Fix for u not being initialized before call --- BEAMS3D/Sources/beams3d_physics_mod.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/BEAMS3D/Sources/beams3d_physics_mod.f90 b/BEAMS3D/Sources/beams3d_physics_mod.f90 index 225cbceb1..3fff779e5 100644 --- a/BEAMS3D/Sources/beams3d_physics_mod.f90 +++ b/BEAMS3D/Sources/beams3d_physics_mod.f90 @@ -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 From a0e1d8b146b3304f68b10a6fada694fa1f88a4d9 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 09:04:58 +0100 Subject: [PATCH 06/16] BEAMS3D: Fix for u not being initialized before call --- BEAMS3D/Sources/beams3d_init_vmec.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/BEAMS3D/Sources/beams3d_init_vmec.f90 b/BEAMS3D/Sources/beams3d_init_vmec.f90 index a0c605199..3cea9edda 100644 --- a/BEAMS3D/Sources/beams3d_init_vmec.f90 +++ b/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; From 1d78afa6361000ccecdf9dc370303b200037cd50 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 10:52:24 +0100 Subject: [PATCH 07/16] SHARE: Improvement to debug line for cobra --- SHARE/make_cobra.inc | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/SHARE/make_cobra.inc b/SHARE/make_cobra.inc index 79c5e3022..ab62bf0b0 100644 --- a/SHARE/make_cobra.inc +++ b/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 \ From d59652e1fe5026db73314eb39d098fad57f88ee8 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 11:00:42 +0100 Subject: [PATCH 08/16] WALL_ACCELERATE: Minor screen bug fix --- WALL_ACCELERATE/Sources/wall_accelerate_main.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/WALL_ACCELERATE/Sources/wall_accelerate_main.f90 b/WALL_ACCELERATE/Sources/wall_accelerate_main.f90 index 0c4c93509..e2ac9ceb0 100644 --- a/WALL_ACCELERATE/Sources/wall_accelerate_main.f90 +++ b/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) From c413288f00a84634d4fcca17c9a6bc0e9a0d4f49 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 11:07:45 +0100 Subject: [PATCH 09/16] BEAMS3D: Minor bugfix for Suzuki model --- BEAMS3D/Sources/beams3d_physics_mod.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BEAMS3D/Sources/beams3d_physics_mod.f90 b/BEAMS3D/Sources/beams3d_physics_mod.f90 index 3fff779e5..08684e80b 100644 --- a/BEAMS3D/Sources/beams3d_physics_mod.f90 +++ b/BEAMS3D/Sources/beams3d_physics_mod.f90 @@ -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 From b0491d6c717cdd1a1b2169bd69d523e0e2c5fbaf Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 11:34:12 +0100 Subject: [PATCH 10/16] BEAMS3D: Fixed debug issue with diagnostics --- BEAMS3D/Sources/beams3d_diagnostics.f90 | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/BEAMS3D/Sources/beams3d_diagnostics.f90 b/BEAMS3D/Sources/beams3d_diagnostics.f90 index 8fb273076..7129e01aa 100644 --- a/BEAMS3D/Sources/beams3d_diagnostics.f90 +++ b/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(rprec), 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,18 @@ 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 From ba928a82ebbbbfe800a255897a6076f721fd4550 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 13:41:39 +0100 Subject: [PATCH 11/16] BEASM3D: Fixed type issue in diagnostics --- BEAMS3D/Sources/beams3d_diagnostics.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/BEAMS3D/Sources/beams3d_diagnostics.f90 b/BEAMS3D/Sources/beams3d_diagnostics.f90 index 7129e01aa..60b5dec0a 100644 --- a/BEAMS3D/Sources/beams3d_diagnostics.f90 +++ b/BEAMS3D/Sources/beams3d_diagnostics.f90 @@ -39,7 +39,7 @@ SUBROUTINE beams3d_diagnostics INTEGER, ALLOCATABLE :: dist_func(:,:,:) REAL, ALLOCATABLE :: vllaxis(:),vperpaxis(:), nlost(:), norbit(:), tlow(:), thigh(:) REAL, ALLOCATABLE :: help3d(:,:,:) - REAL(rprec), ALLOCATABLE :: dbl_help(:) + REAL, ALLOCATABLE :: dbl_help(:) INTEGER :: mystart !----------------------------------------------------------------------- ! Begin Subroutine @@ -89,6 +89,8 @@ SUBROUTINE beams3d_diagnostics 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))) From f6c67295ae790afaa38db9ac8ac035f5dcd9bf72 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 13:55:26 +0100 Subject: [PATCH 12/16] LIBSTELL: Debug fixes for wall_mod.f90 --- LIBSTELL/Sources/Modules/wall_mod.f90 | 61 +++++++++++++++++---------- 1 file changed, 39 insertions(+), 22 deletions(-) diff --git a/LIBSTELL/Sources/Modules/wall_mod.f90 b/LIBSTELL/Sources/Modules/wall_mod.f90 index f06ef53fa..8270de154 100644 --- a/LIBSTELL/Sources/Modules/wall_mod.f90 +++ b/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 From 15898237847465f847be0ae24b2b8cd544fdfb3b Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 13:56:44 +0100 Subject: [PATCH 13/16] SHARE: Cleanup of make_cobra.inc --- SHARE/make_cobra.inc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SHARE/make_cobra.inc b/SHARE/make_cobra.inc index ab62bf0b0..6a13b29ff 100644 --- a/SHARE/make_cobra.inc +++ b/SHARE/make_cobra.inc @@ -19,7 +19,7 @@ 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 \ From 90c4647b4edb965ae7195f01c0e95ef0d7d0f836 Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 14:27:33 +0100 Subject: [PATCH 14/16] SHARE: Updated compiler flags for Raven --- SHARE/make_raven.inc | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/SHARE/make_raven.inc b/SHARE/make_raven.inc index 2adcfd3fa..1ec780c7c 100644 --- a/SHARE/make_raven.inc +++ b/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 \ From 9809943c9650b069ba32dfa25c878b0ebedd648a Mon Sep 17 00:00:00 2001 From: Samuel Lazerson Date: Sat, 13 Jan 2024 15:02:44 +0100 Subject: [PATCH 15/16] BEAMS3D: Fixed segfault bug when beam density is requested. --- BEAMS3D/Sources/beams3d_follow.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/BEAMS3D/Sources/beams3d_follow.f90 b/BEAMS3D/Sources/beams3d_follow.f90 index 79ac31c42..070a060dd 100644 --- a/BEAMS3D/Sources/beams3d_follow.f90 +++ b/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 From b0d4eb379044b3c2e8202e39078f8a1d93c1dfa7 Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 17 Jan 2024 14:04:31 +0100 Subject: [PATCH 16/16] BEAMS3D: Bugfix for possible s<0 in physics mod. --- BEAMS3D/Sources/beams3d_physics_mod.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/BEAMS3D/Sources/beams3d_physics_mod.f90 b/BEAMS3D/Sources/beams3d_physics_mod.f90 index 08684e80b..d5bf68b0c 100644 --- a/BEAMS3D/Sources/beams3d_physics_mod.f90 +++ b/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)