Merge branch 'develop' into bugfix/BEAMS3D_minor

BEAMS3D/BEAMS3D.dep

@@ -1,6 +1,14 @@
 # Microsoft Developer Studio Generated Dependency File, included by BEAMS3D.mak
 
 
+beams3d_beam_density.o: \
+      beams3d_runtime.o \
+      beams3d_grid.o \
+      beams3d_lines.o \
+      ../../LIBSTELL/$(LOCTYPE)/mpi_inc.o \
+      ../../LIBSTELL/$(LOCTYPE)/mpi_sharmem.o \
+      ../../LIBSTELL/$(LOCTYPE)/mpi_params.o 
+
 beams3d_interface_mod.o: \
       beams3d_input_mod.o \
       ../../LIBSTELL/$(LOCTYPE)/wall_mod.o \

BEAMS3D/ObjectList

@@ -1,4 +1,5 @@
 ObjectFiles = \
+beams3d_beam_density.o \
 beams3d_duplicate_part.o \
 outpart_beams3d_nag.o \
 beams3d_follow.o \

BEAMS3D/Sources/beams3d_beam_density.f90

@@ -0,0 +1,116 @@
+!-----------------------------------------------------------------------
+!     Module:        beams3d_beam_density
+!     Authors:       S. Lazerson (lazerson@pppl.gov)
+!     Date:          04/03/2023
+!     Description:   This subroutine computes the neutral beam density
+!                    on a grid.
+!-----------------------------------------------------------------------
+      SUBROUTINE beams3d_beam_density
+!-----------------------------------------------------------------------
+!     Libraries
+!-----------------------------------------------------------------------
+      USE beams3d_runtime
+      USE beams3d_grid
+      USE beams3d_lines
+      USE mpi_params ! MPI
+      USE mpi_inc
+      USE mpi_sharmem
+!-----------------------------------------------------------------------
+!     Local Variables
+!          ier          Error Flag
+!          iunit        File ID
+!          ndist        Number of Vll divisions for dist function
+!          ns           Number of flux divisions for current calculation
+!-----------------------------------------------------------------------
+      IMPLICIT NONE
+      LOGICAL :: lline_in_box
+      INTEGER :: ier, l, nl, i, j, k, m, s
+      DOUBLE PRECISION :: x0, y0, z0, x1, y1, z1, hr2, hz2, hp2, d, &
+                          denbeam, dl, dV
+      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xl,yl,zl,rl,sl,pl
+#if defined(MPI_OPT)
+      INTEGER :: MPI_COMM_LOCAL
+      INTEGER :: mystart, mypace
+#endif
+!-----------------------------------------------------------------------
+!     Begin Subroutine
+!-----------------------------------------------------------------------
+
+      ! Allocate the Grid
+      CALL mpialloc(BEAM_DENSITY, nbeams, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_BEAM_DENSITY)
+      IF (myid_sharmem == 0) BEAM_DENSITY = 0
+      nl = nr + nr/2 + 1
+      ALLOCATE(xl(nl),yl(nl),zl(nl),rl(nl),pl(nl),sl(nl))
+      FORALL(s=1:nl) sl(s) = DBLE(s-0.5)/DBLE(nl) ! half grid
+
+
+      IF (lverb) WRITE(6,'(A)') '----- BEAM Density Calc. -----'
+      hr2 = 0.5*hr(1)
+      hz2 = 0.5*hz(1)
+      hp2 = 0.5*hp(1)
+      ! Loop over particles
+      DO l = mystart_save, myend_save
+         x0 = R_lines(0,l)*cos(PHI_lines(0,l))
+         y0 = R_lines(0,l)*sin(PHI_lines(0,l))
+         z0 = Z_lines(0,l)
+         x1 = R_lines(1,l)*cos(PHI_lines(1,l))
+         y1 = R_lines(1,l)*sin(PHI_lines(1,l))
+         z1 = Z_lines(1,l)
+         m  = Beam(l)
+         ! Because we're in cylindrical space we just need to brute force this
+         ! We have a chord with particles/s (weight) going along it.
+         ! Since the velocity is constant then the TOF is just L/v = t
+         ! The the total number of particles is just n=W*t = W*L/v
+         ! then we just need to calc the volume of the voxel to get particles/m^3
+         d  = SQRT((x1-x0)**2 + (y1-y0)**2 + (z1-z0)**2)/DBLE(nl) ! so it's dl
+         denbeam = weight(l)*d/vll_lines(0,l) ! This is the total number of particles per step
+         xl = x0 + sl*(x1-x0)
+         yl = y0 + sl*(y1-y0)
+         zl = z0 + sl*(z1-z0)
+         rl = sqrt(xl*xl+yl*yl)
+         pl = atan2(yl,xl)
+         pl = MODULO(pl,phimax)
+         DO s = 1, nl
+            ! Find gridpoint but use half grid
+            i = MIN(MAX(COUNT(raxis-hr2 < rl(s)),1),nr)
+            j = MIN(MAX(COUNT(phiaxis-hp2 < pl(s)),1),nphi)
+            k = MIN(MAX(COUNT(zaxis-hz2 < zl(s)),1),nz)
+            BEAM_DENSITY(m,i,j,k) = BEAM_DENSITY(m,i,j,k) + denbeam
+         END DO
+      END DO
+      DEALLOCATE(xl,yl,zl,rl,pl,sl)
+
+      ! Fix edges which are double counts
+      IF (myid_sharmem == 0) THEN
+         BEAM_DENSITY(:,1,:,:) = 0
+         BEAM_DENSITY(:,nr,:,:) = 0
+         BEAM_DENSITY(:,:,:,1) = 0
+         BEAM_DENSITY(:,1,:,nz) = 0
+      ENDIF
+
+
+      ! Barrier so  calculation is done
+#if defined(MPI_OPT)
+      CALL MPI_BARRIER(MPI_COMM_SHARMEM,ierr_mpi)
+#endif
+      ! Now Divide by the grid volumes
+      CALL MPI_CALC_MYRANGE(MPI_COMM_SHARMEM, 1, (nr-1)*(nphi-1)*(nz-1), mystart, myend)
+      DO s = mystart,myend
+         i = MOD(s-1,nr-1)+1
+         j = MOD(s-1,(nr-1)*(nphi-1))
+         j = FLOOR(REAL(j) / REAL(nr-1))+1
+         k = CEILING(REAL(s) / REAL((nr-1)*(nphi-1)))
+         dV = raxis(i)*hr(i)*hp(j)*hz(k)
+         BEAM_DENSITY(:,i,j,k) = BEAM_DENSITY(:,i,j,k) / dV
+      END DO
+
+#if defined(MPI_OPT)
+      CALL MPI_BARRIER(MPI_COMM_SHARMEM,ierr_mpi)
+#endif
+
+      RETURN
+
+!-----------------------------------------------------------------------
+!     End Subroutine
+!----------------------------------------------------------------------- 
+      END SUBROUTINE beams3d_beam_density

BEAMS3D/Sources/beams3d_fix_poloidal.f90

@@ -42,8 +42,8 @@ SUBROUTINE beams3d_fix_poloidal
 !-----------------------------------------------------------------------
 
       ! Recalculate U for all particles
-      mystart = LBOUND(R_lines,2)
-      DO myline = mystart,myend
+      !mystart = LBOUND(R_lines,2)
+      DO myline = mystart_save, myend_save
          DO mystep = 0, npoinc
             s1 = S_lines(mystep,myline)
             r1 = R_lines(mystep,myline)
@@ -69,9 +69,6 @@ SUBROUTINE beams3d_fix_poloidal
       END DO
       WHERE ( R_lines < 0 ) U_lines = 0
 
-      ! Wait to deallocate till everyone is done with the memory
-      CALL MPI_BARRIER(MPI_COMM_SHARMEM, ier)
-
       RETURN
 !-----------------------------------------------------------------------
 !     END SUBROUTINE

BEAMS3D/Sources/beams3d_follow.f90

@@ -17,7 +17,8 @@ SUBROUTINE beams3d_follow
                             MODB_spl, S_spl, U_spl, TE_spl, NE_spl, TI_spl, &
                             TE_spl, TI_spl, wall_load, wall_shine, rho_fullorbit, &
                             plasma_mass, plasma_Zmean, therm_factor, &
-                            nr_fida, nphi_fida, nz_fida, nenergy_fida, npitch_fida
+                            nr_fida, nphi_fida, nz_fida, nenergy_fida, &
+                            npitch_fida, BEAM_DENSITY
     USE mpi_params ! MPI
     USE beams3d_physics_mod
     USE beams3d_write_par
@@ -221,6 +222,12 @@ SUBROUTINE beams3d_follow
 
     DEALLOCATE(q)
 
+    ! Calcualte Beam Density
+#if defined(MPI_OPT)
+    CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)
+#endif
+    IF (lbeam .and. lbeamdensity) CALL beams3d_beam_density
+
     ! Follow Trajectories
     IF (.not.ldepo) THEN
         CALL beams3d_follow_gc
@@ -270,25 +277,28 @@ SUBROUTINE beams3d_follow
        IF (ASSOCIATED(wall_shine)) THEN
           CALL MPI_ALLREDUCE(MPI_IN_PLACE,wall_shine,nface*nbeams,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_LOCAL,ierr_mpi)
        END IF
+       IF (ASSOCIATED(BEAM_DENSITY)) THEN
+          CALL MPI_ALLREDUCE(MPI_IN_PLACE, BEAM_DENSITY, nbeams*nr*nphi*nz, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_LOCAL, ierr_mpi)
+       END IF
        CALL MPI_COMM_FREE(MPI_COMM_LOCAL,ierr_mpi)
     END IF
     CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)
 
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'R_lines', DBLVAR=R_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,    'PHI_lines', DBLVAR=PHI_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'Z_lines', DBLVAR=Z_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,    'vll_lines', DBLVAR=vll_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend, 'moment_lines', DBLVAR=moment_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,     'vr_lines', DBLVAR=vr_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,   'vphi_lines', DBLVAR=vphi_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,     'vz_lines', DBLVAR=vz_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'S_lines', DBLVAR=S_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'U_lines', DBLVAR=U_lines)
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,      'B_lines', DBLVAR=B_lines)
-    CALL beams3d_write1d_parhdf5(         1, nparticles, mystart, myend,      't_end',   DBLVAR=t_last,FILENAME='beams3d_'//TRIM(id_string))
-    ALLOCATE(itemp(0:npoinc,mystart:myend))
-    itemp = 0; WHERE(neut_lines) itemp=1;
-    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart, myend,   'neut_lines', INTVAR=itemp)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,      'R_lines', DBLVAR=R_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,    'PHI_lines', DBLVAR=PHI_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,      'Z_lines', DBLVAR=Z_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,    'vll_lines', DBLVAR=vll_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save, 'moment_lines', DBLVAR=moment_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,     'vr_lines', DBLVAR=vr_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,   'vphi_lines', DBLVAR=vphi_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,     'vz_lines', DBLVAR=vz_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,      'S_lines', DBLVAR=S_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,      'U_lines', DBLVAR=U_lines)
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,      'B_lines', DBLVAR=B_lines)
+    CALL beams3d_write1d_parhdf5(         1, nparticles, mystart_save, myend_save,      't_end',   DBLVAR=t_last,FILENAME='beams3d_'//TRIM(id_string))
+    ALLOCATE(itemp(0:npoinc,mystart_save:myend_save))
+    itemp = 0; WHERE(neut_lines) itemp=1
+    CALL beams3d_write_parhdf5(0, npoinc, 1, nparticles, mystart_save, myend_save,   'neut_lines', INTVAR=itemp)
     DEALLOCATE(itemp)
     IF (ALLOCATED(mnum)) DEALLOCATE(mnum)
     IF (ALLOCATED(moffsets)) DEALLOCATE(moffsets)

BEAMS3D/Sources/beams3d_grid.f90

@@ -64,7 +64,7 @@ MODULE beams3d_grid
       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(:,:,:)
-      REAL(rprec), POINTER :: NI(:,:,:,:), BEAM_DENSE(:,:,:,:)
+      REAL(rprec), POINTER :: NI(:,:,:,:), BEAM_DENSITY(:,:,:,:)
       REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: X_BEAMLET, Y_BEAMLET, Z_BEAMLET, &
                                                   NX_BEAMLET, NY_BEAMLET, NZ_BEAMLET
       REAL(rprec), DIMENSION(:,:,:,:), POINTER :: BR4D, BPHI4D, BZ4D, MODB4D, &

BEAMS3D/Sources/beams3d_init.f90

@@ -255,8 +255,18 @@ SUBROUTINE beams3d_init
             POT_spl_s%isHermite   = 0
             CALL EZspline_setup(POT_spl_s,POT_AUX_F(1:npot),ier,EXACT_DIM=.true.)
             IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init10',ier)
-         IF (lverb) WRITE(6,'(A,F9.5,A,F9.5,A,I4,A,F8.5)') '   V    = [', &
-            MINVAL(POT_AUX_F(1:npot))*1E-3,',',MAXVAL(POT_AUX_F(1:npot))*1E-3,'] kV;  NPOT: ',npot, ';  S_MAX_POT: ',s_max_pot
+            IF (lverb) THEN
+               IF (MAXVAL(ABS(POT_AUX_F(1:npot))) > 1E6) THEN
+                  WRITE(6,'(A,F9.4,A,F9.4,A,I4,A,F8.5)') '   V    = [', &
+                     MINVAL(POT_AUX_F(1:npot))*1E-6,',',MAXVAL(POT_AUX_F(1:npot))*1E-6,'] MV;  NPOT: ',npot, ';  S_MAX_POT: ',s_max_pot
+               ELSEIF (MAXVAL(ABS(POT_AUX_F(1:npot))) > 1E3) THEN
+                  WRITE(6,'(A,F9.4,A,F9.4,A,I4,A,F8.5)') '   V    = [', &
+                     MINVAL(POT_AUX_F(1:npot))*1E-3,',',MAXVAL(POT_AUX_F(1:npot))*1E-3,'] kV;  NPOT: ',npot, ';  S_MAX_POT: ',s_max_pot
+               ELSE
+                  WRITE(6,'(A,F9.4,A,F9.4,A,I4,A,F8.5)') '   V    = [', &
+                     MINVAL(POT_AUX_F(1:npot)),',',MAXVAL(POT_AUX_F(1:npot)),'] V;  NPOT: ',npot, ';  S_MAX_POT: ',s_max_pot
+               END IF
+            END IF
          END IF
 
          IF (lverb) THEN
@@ -738,6 +748,17 @@ SUBROUTINE beams3d_init
       ! Duplicate particles if requested
       IF (duplicate_factor > 1) CALL beams3d_duplicate_part
 
+      ! Print a warning if the range of phi_start > hphi
+      IF (lverb) THEN
+         phitemp = MAXVAL(phi_start)-MINVAL(phi_start)
+         IF ((phitemp < hp(1)) .and. (phitemp > 0) .and. lbeamdensity) THEN
+            i = (phiaxis(nphi)-phiaxis(1))/phitemp
+            WRITE(6,'(A)')    ' WARNING: The range of PHI_START > dphi.'
+            WRITE(6,'(A,I4)') '          Consider increasing nphi >=',i
+            CALL FLUSH(6)
+         END IF
+      END IF
+
       ! In all cases create an end_state array
       ALLOCATE(end_state(nparticles))
       end_state=0

BEAMS3D/Sources/beams3d_interface_mod.f90

@@ -163,6 +163,7 @@ SUBROUTINE beams3d_init_commandline
          lfusion_He3     = .false.
          lboxsim = .false.
          lfieldlines = .false.
+         lbeamdensity = .true.
          id_string = ''
          coil_string = ''
          mgrid_string = ''
@@ -284,6 +285,8 @@ SUBROUTINE beams3d_init_commandline
                 lfusion_He3 = .true.
             case ("-boxsim")
                 lboxsim = .true.
+            case ("-nobeamdensity")
+                lbeamdensity = .false.
             case ("-help", "-h") ! Output Help message
                 write(6, *) ' Beam MC Code'
                 write(6, *) ' Usage: xbeams3d <options>'
@@ -316,6 +319,7 @@ SUBROUTINE beams3d_init_commandline
                 write(6, *) '     -fusion_proton:  Thermal Fusion Births Rates  (proton only)'
                 write(6, *) '     -fusion_he3:     Fusion Reaction Rates for birth (He3 only)'
                 write(6, *) '     -boxsim:         Inject charged particles for box modeling'
+                write(6, *) '     -nobeamdensity:  Supress beam density calc.'
                 write(6, *) '     -noverb:         Supress all screen output'
                 write(6, *) '     -help:           Output help message'
             end select
@@ -408,6 +412,8 @@ SUBROUTINE beams3d_init_commandline
       IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_BCAST_ERR,'beams3d_main',ierr_mpi)
       CALL MPI_BCAST(lboxsim,1,MPI_LOGICAL, master, MPI_COMM_BEAMS,ierr_mpi)
       IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_BCAST_ERR,'beams3d_main',ierr_mpi)
+      CALL MPI_BCAST(lbeamdensity,1,MPI_LOGICAL, master, MPI_COMM_BEAMS,ierr_mpi)
+      IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_BCAST_ERR,'beams3d_main',ierr_mpi)
       CALL MPI_BCAST(rminor_norm,1,MPI_DOUBLE_PRECISION, master, MPI_COMM_BEAMS,ierr_mpi)
       IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_BCAST_ERR,'beams3d_main',ierr_mpi)
 #endif

BEAMS3D/Sources/beams3d_runtime.f90

@@ -135,7 +135,7 @@ MODULE beams3d_runtime
                ldepo, lbeam_simple, ldebug, lcollision, lw7x, lsuzuki, &
                lascot, lascot4, lbbnbi, lfidasim, lfidasim_cyl, lsplit, lvessel_beam, lascotfl, lrandomize, &
                lfusion, lfusion_alpha, leqdsk, lhint, lkick, lgcsim, &
-               lboxsim, limas, lfieldlines, lfusion_tritium, lfusion_proton, lfusion_He3
+               lboxsim, limas, lfieldlines, lfusion_tritium, lfusion_proton, lfusion_He3, lbeamdensity
     INTEGER :: nextcur, npoinc, nbeams, nparticles_start, nprocs_beams, &
                ndt, ndt_max, duplicate_factor
     INTEGER, DIMENSION(MAXBEAMS) :: Dex_beams

BEAMS3D/Sources/beams3d_write.f90

@@ -18,7 +18,8 @@ SUBROUTINE beams3d_write(write_type)
                                  zaxis, phiaxis, S_ARR, U_ARR, POT_ARR, &
                                  ZEFF_ARR, TE, TI, NE, wall_load, wall_shine, &
                                  plasma_mass, plasma_Zmean, &
-                                 B_kick_min, B_kick_max, freq_kick, E_kick, NI
+                                 B_kick_min, B_kick_max, freq_kick, &
+                                 E_kick, NI, beam_density
       USE beams3d_runtime, ONLY: id_string, npoinc, nbeams, beam, t_end, lverb, &
                                     lvmec, lpies, lspec, lcoil, lmgrid, lbeam, lascot, &
                                     lvessel, lvac, lbeam_simple, handle_err, nparticles_start, &
@@ -211,6 +212,11 @@ SUBROUTINE beams3d_write(write_type)
                                    ATT='Neutral Beam Shine-through [W/m^2]',ATT_NAME='description')
                   IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'wall_shine',ier)
                END IF
+               IF (ASSOCIATED(BEAM_DENSITY)) THEN
+                  CALL write_var_hdf5(fid,'beam_density',nbeams,nr,nphi,nz,ier,DBLVAR=beam_density,&
+                                   ATT='Neutral Beam Density [1/m^3]',ATT_NAME='description')
+                  IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'beam_density',ier)
+               END IF
             CASE('TRAJECTORY_FULL')
                CALL open_hdf5('beams3d_'//TRIM(id_string)//'.h5',fid,ier,LCREATE=.false.)
                IF (ier /= 0) CALL handle_err(HDF5_OPEN_ERR,'beams3d_'//TRIM(id_string)//'.h5',ier)