Merge pull request #274 from PrincetonUniversity/feature/XYRHO4D

Interpolate using xrho and yrho, not S

BEAMS3D/Sources/beams3d_distnorm.f90

@@ -12,7 +12,7 @@ SUBROUTINE beams3d_distnorm
       USE stel_kinds, ONLY: rprec
       USE beams3d_runtime, ONLY: lfidasim, lfidasim_cyl, nbeams, pi2, &
                               EZSPLINE_ERR, MPI_BARRIER_ERR
-      USE beams3d_grid, ONLY: raxis, phiaxis, zaxis, X4D, Y4D, S4D, &
+      USE beams3d_grid, ONLY: raxis, phiaxis, zaxis, &
                               nr, nphi, nz
       USE beams3d_lines, ONLY: ns_prof1, ns_prof2, ns_prof3, ns_prof4, &
                                ns_prof4, ns_prof5, dist5d_prof, &

BEAMS3D/Sources/beams3d_fix_poloidal.f90

@@ -26,11 +26,8 @@ SUBROUTINE beams3d_fix_poloidal
 !-----------------------------------------------------------------------
 !     Local Variables
 !-----------------------------------------------------------------------
-      INTEGER :: mystart, mystep,  ier !win_X4D, win_Y4D,
+      INTEGER :: mystart, mystep,  ier 
       INTEGER :: bcs1(2), bcs2(2), bcs3(2)
-      !TYPE(EZspline3_r8) :: X_spl, Y_spl
-      !REAL(rprec), POINTER, DIMENSION(:,:,:) :: X_ARR, Y_ARR
-      !REAL(rprec), POINTER, DIMENSION(:,:,:,:) :: X4D, Y4D
 
       ! For splines
       INTEGER :: i,j,k,l
@@ -42,10 +39,9 @@ SUBROUTINE beams3d_fix_poloidal
 !-----------------------------------------------------------------------
 
       ! Recalculate U for all particles
-      !mystart = LBOUND(R_lines,2)
       DO myline = mystart_save, myend_save
          DO mystep = 0, npoinc
-            s1 = S_lines(mystep,myline)
+            !s1 = S_lines(mystep,myline)
             r1 = R_lines(mystep,myline)
             p1 = MOD(PHI_lines(mystep,myline),phimax)
             IF (p1 < 0) p1 = p1 + phimax
@@ -56,14 +52,16 @@ SUBROUTINE beams3d_fix_poloidal
             xparam = (r1 - raxis(i)) * hri(i)
             yparam = (p1 - phiaxis(j)) * hpi(j)
             zparam = (z1 - zaxis(k)) * hzi(k)
+            ! Even taking X4D, Y4D to XRHO4D, YRHO4D it's still 
+            ! true that U = atan2(yrho / xrho)
             CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
                             hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),&
-                            X4D(1,1,1,1),nr,nphi,nz)
-            xt = fval(1)/s1
+                            XRHO4D(1,1,1,1),nr,nphi,nz)
+            xt = fval(1)
             CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
                             hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),&
-                            Y4D(1,1,1,1),nr,nphi,nz)
-            yt = fval(1)/s1
+                            YRHO4D(1,1,1,1),nr,nphi,nz)
+            yt = fval(1)
             U_lines(mystep,myline) = ATAN2(yt,xt)
          END DO
       END DO
@@ -73,4 +71,4 @@ SUBROUTINE beams3d_fix_poloidal
 !-----------------------------------------------------------------------
 !     END SUBROUTINE
 !-----------------------------------------------------------------------
-      END SUBROUTINE beams3d_fix_poloidal
+      END SUBROUTINE beams3d_fix_poloidal

BEAMS3D/Sources/beams3d_follow.f90

@@ -13,9 +13,7 @@ SUBROUTINE beams3d_follow
     USE stel_kinds, ONLY: rprec
     USE beams3d_runtime
     USE beams3d_lines
-    USE beams3d_grid, ONLY: tmin, tmax, delta_t, BR_spl, BZ_spl, BPHI_spl, &
-                            MODB_spl, S_spl, U_spl, TE_spl, NE_spl, TI_spl, &
-                            TE_spl, TI_spl, wall_load, wall_shine, rho_fullorbit, &
+    USE beams3d_grid, ONLY: tmin, tmax, delta_t, wall_load, wall_shine, rho_fullorbit, &
                             plasma_mass, plasma_Zmean, therm_factor, &
                             nr_fida, nphi_fida, nz_fida, nenergy_fida, &
                             npitch_fida, BEAM_DENSITY
@@ -156,7 +154,7 @@ SUBROUTINE beams3d_follow
        mycharge = charge(i)
        myZ = Zatom(i)
        mymass = mass(i)
-	    E_by_v=mymass*0.5d-3/e_charge
+       E_by_v=mymass*0.5d-3/e_charge
        mybeam = Beam(i)
        moment = mu_start(i)
        fact_pa   = plasma_mass/(mymass*plasma_Zmean)
@@ -188,7 +186,7 @@ SUBROUTINE beams3d_follow
           mycharge = charge(i)
           myZ = Zatom(i)
           mymass = mass(i)
-		    E_by_v=mymass*0.5d-3/e_charge
+          E_by_v=mymass*0.5d-3/e_charge
           mybeam = Beam(i)
           moment = mu_start(i)
           my_end = t_end(i)

BEAMS3D/Sources/beams3d_follow_fo.f90

@@ -19,11 +19,8 @@ SUBROUTINE beams3d_follow_fo
     USE stel_kinds, ONLY: rprec
     USE beams3d_runtime
     USE beams3d_lines
-    USE beams3d_grid, ONLY: tmin, tmax, delta_t, BR_spl, BZ_spl, BPHI_spl, &
-                            MODB_spl, S_spl, U_spl, TE_spl, NE_spl, TI_spl, &
-                            TE_spl, TI_spl, wall_load, wall_shine, &
-                            plasma_mass, plasma_Zmean, therm_factor, &
-                            rho_fullorbit
+    USE beams3d_grid, ONLY: plasma_mass, plasma_Zmean, &
+                            rho_fullorbit, rho_help, E_kick, freq_kick
     USE mpi_params ! MPI
     USE beams3d_write_par
     USE beams3d_physics_mod, ONLY: beams3d_gc2fo, beams3d_calc_dt

BEAMS3D/Sources/beams3d_follow_gc.f90

@@ -19,9 +19,7 @@ SUBROUTINE beams3d_follow_gc
     USE stel_kinds, ONLY: rprec
     USE beams3d_runtime
     USE beams3d_lines
-    USE beams3d_grid, ONLY: tmin, tmax, delta_t, BR_spl, BZ_spl, BPHI_spl, &
-                            MODB_spl, S_spl, U_spl, TE_spl, NE_spl, TI_spl, &
-                            TE_spl, TI_spl, wall_load, wall_shine, &
+    USE beams3d_grid, ONLY: tmin, tmax, delta_t, wall_load, wall_shine, &
                             plasma_mass, plasma_Zmean, therm_factor, &
                             rho_fullorbit, rho_help, E_kick, freq_kick, &
                             nr_fida, nphi_fida, nz_fida, nenergy_fida, npitch_fida,raxis

BEAMS3D/Sources/beams3d_free.f90

@@ -41,7 +41,7 @@ SUBROUTINE beams3d_free(IN_COMM)
       IF (EZspline_allocated(BZ_spl))   CALL EZspline_free(BZ_spl,ier)
       IF (EZspline_allocated(BPHI_spl)) CALL EZspline_free(BPHI_spl,ier)
       IF (EZspline_allocated(MODB_spl)) CALL EZspline_free(MODB_spl,ier)
-      IF (EZspline_allocated(S_spl)) CALL EZspline_free(S_spl,ier)
+      IF (EZspline_allocated(RHO_spl)) CALL EZspline_free(RHO_spl,ier)
       IF (EZspline_allocated(U_spl)) CALL EZspline_free(U_spl,ier)
       IF (EZspline_allocated(TE_spl))   CALL EZspline_free(TE_spl,ier)
       IF (EZspline_allocated(NE_spl))   CALL EZspline_free(NE_spl,ier)
@@ -88,9 +88,10 @@ SUBROUTINE beams3d_free(IN_COMM)
          IF (ASSOCIATED(B_Z))      CALL mpidealloc(B_Z,win_B_Z)
          IF (ASSOCIATED(MODB))     CALL mpidealloc(MODB,win_MODB)
          IF (ASSOCIATED(S_ARR))    CALL mpidealloc(S_ARR,win_S_ARR)
+         IF (ASSOCIATED(RHO_ARR))  CALL mpidealloc(RHO_ARR,win_RHO_ARR)
          IF (ASSOCIATED(U_ARR))    CALL mpidealloc(U_ARR,win_U_ARR)
-         IF (ASSOCIATED(X_ARR))    CALL mpidealloc(X_ARR,win_X_ARR)
-         IF (ASSOCIATED(Y_ARR))    CALL mpidealloc(Y_ARR,win_Y_ARR)
+         IF (ASSOCIATED(XRHO_ARR)) CALL mpidealloc(XRHO_ARR,win_XRHO_ARR)
+         IF (ASSOCIATED(YRHO_ARR)) CALL mpidealloc(YRHO_ARR,win_YRHO_ARR)
          IF (ASSOCIATED(TE))       CALL mpidealloc(TE,win_TE)
          IF (ASSOCIATED(TI))       CALL mpidealloc(TI,win_TI)
          IF (ASSOCIATED(NE))       CALL mpidealloc(NE,win_NE)
@@ -106,10 +107,10 @@ SUBROUTINE beams3d_free(IN_COMM)
          IF (ASSOCIATED(NI5D))     CALL mpidealloc(NI5D,win_NI5D)
          IF (ASSOCIATED(TI4D))     CALL mpidealloc(TI4D,win_TI4D)
          IF (ASSOCIATED(ZEFF4D))   CALL mpidealloc(ZEFF4D,win_ZEFF4D)
-         IF (ASSOCIATED(S4D))      CALL mpidealloc(S4D,win_S4D)
+         IF (ASSOCIATED(RHO4D))    CALL mpidealloc(RHO4D,win_RHO4D)
          IF (ASSOCIATED(U4D))      CALL mpidealloc(U4D,win_U4D)
-         IF (ASSOCIATED(X4D))      CALL mpidealloc(X4D,win_X4D)
-         IF (ASSOCIATED(Y4D))      CALL mpidealloc(Y4D,win_Y4D)
+         IF (ASSOCIATED(XRHO4D))   CALL mpidealloc(XRHO4D,win_XRHO4D)
+         IF (ASSOCIATED(YRHO4D))   CALL mpidealloc(YRHO4D,win_YRHO4D)
          IF (ASSOCIATED(POT4D))    CALL mpidealloc(POT4D,win_POT4D)
          IF (ASSOCIATED(wall_load))   CALL mpidealloc(wall_load,win_wall_load)
          IF (ASSOCIATED(wall_shine))  CALL mpidealloc(wall_shine,win_wall_shine)
@@ -166,7 +167,7 @@ SUBROUTINE beams3d_free(IN_COMM)
          IF (ASSOCIATED(NI5D))     DEALLOCATE(NI5D)
          IF (ASSOCIATED(TI4D))     DEALLOCATE(TI4D)
          IF (ASSOCIATED(ZEFF4D))   DEALLOCATE(ZEFF4D)
-         IF (ASSOCIATED(S4D))      DEALLOCATE(S4D)
+         IF (ASSOCIATED(RHO4D))      DEALLOCATE(RHO4D)
          IF (ASSOCIATED(U4D))      DEALLOCATE(U4D)
          IF (ASSOCIATED(POT4D))    DEALLOCATE(POT4D)
          IF (ASSOCIATED(wall_load))    DEALLOCATE(wall_load)

BEAMS3D/Sources/beams3d_grid.f90

@@ -56,8 +56,10 @@ MODULE beams3d_grid
       INTEGER  ::    win_raxis, win_phiaxis, win_zaxis, win_B_R, win_B_PHI, win_B_Z,&
                      win_MODB, win_TE, win_NE, win_TI, win_ZEFF_ARR,&
                      win_S_ARR, win_U_ARR,win_X_ARR,win_Y_ARR, win_POT_ARR, win_BR4D, win_BPHI4D, &
+                     win_RHO_ARR, win_XRHO_ARR, win_YRHO_ARR, &
                      win_BZ4D, win_MODB4D, win_TE4D, win_NE4D, win_TI4D, win_ZEFF4D, &
-                     win_S4D, win_U4D,  win_X4D, win_Y4D, win_POT4D, win_req_axis, win_zeq_axis, &
+                     win_U4D,  win_POT4D, win_req_axis, win_zeq_axis, &
+                     win_RHO4D, win_XRHO4D, win_YRHO4D, &
                      win_wall_load, win_wall_shine, win_hr, win_hp, win_hz, &
                      win_hri, win_hpi, win_hzi, win_NI5D, win_NI, &
                      win_raxis_fida, win_phiaxis_fida, win_zaxis_fida, win_energy_fida, win_pitch_fida, &
@@ -71,20 +73,23 @@ MODULE beams3d_grid
       REAL(rprec), POINTER :: wall_load(:,:), wall_shine(:,:)
       REAL(rprec), POINTER :: B_R(:,:,:),B_PHI(:,:,:), B_Z(:,:,:), MODB(:,:,:),&
                                   TE(:,:,:), NE(:,:,:), TI(:,:,:), ZEFF_ARR(:,:,:), &
-                                  S_ARR(:,:,:), U_ARR(:,:,:), X_ARR(:,:,:), Y_ARR(:,:,:), POT_ARR(:,:,:)
+                                  RHO_ARR(:,:,:), XRHO_ARR(:,:,:), YRHO_ARR(:,:,:), &
+                                  S_ARR(:,:,:), U_ARR(:,:,:), POT_ARR(:,:,:)
       REAL(rprec), POINTER :: NI(:,:,:,:), BEAM_DENSITY(:,:,:,:), NEUTRONS_ARR(:,:,:,:)
       REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: X_BEAMLET, Y_BEAMLET, Z_BEAMLET, &
                                                   NX_BEAMLET, NY_BEAMLET, NZ_BEAMLET
       REAL(rprec), DIMENSION(:,:,:,:), POINTER :: BR4D, BPHI4D, BZ4D, MODB4D, &
                                   TE4D, NE4D, TI4D, ZEFF4D, &
-                                  S4D, U4D, X4D, Y4D, POT4D
+                                  U4D, POT4D, &
+                                  RHO4D, XRHO4D, YRHO4D
       REAL(rprec), DIMENSION(:,:,:,:,:), POINTER :: NI5D
       REAL*8 ::      eps1, eps2, eps3
       REAL*8, parameter :: small = 1.e-10_ezspline_r8
       REAL*8, POINTER :: hr(:), hp(:), hz(:)
       REAL*8, POINTER :: hri(:), hpi(:), hzi(:)
       TYPE(EZspline3_r8) :: BR_spl, BPHI_spl, BZ_spl, MODB_spl, TE_spl, NE_spl, &
-                            TI_spl, ZEFF_spl, S_spl, U_spl, X_spl, Y_spl,POT_spl
+                            RHO_spl, XRHO_spl, YRHO_spl, &
+                            TI_spl, ZEFF_spl, U_spl, POT_spl
       TYPE(EZspline1_r8) :: TE_spl_s, NE_spl_s, TI_spl_S, ZEFF_spl_s, Vp_spl_s, POT_spl_s
       TYPE(EZspline1_r8), DIMENSION(4) :: NI_spl_s