From 6b02250bdf0b9115a98d92d96da0ef719298155a Mon Sep 17 00:00:00 2001 From: lazersos Date: Sun, 27 Aug 2023 11:09:28 +0200 Subject: [PATCH 1/9] BEAMS3D: Implemented beam density calculation. --- BEAMS3D/BEAMS3D.dep | 8 ++ BEAMS3D/ObjectList | 1 + BEAMS3D/Sources/beams3d_beam_density.f90 | 117 +++++++++++++++++++++++ BEAMS3D/Sources/beams3d_fix_poloidal.f90 | 7 +- BEAMS3D/Sources/beams3d_follow.f90 | 43 +++++---- BEAMS3D/Sources/beams3d_grid.f90 | 2 +- BEAMS3D/Sources/beams3d_write.f90 | 8 +- 7 files changed, 163 insertions(+), 23 deletions(-) create mode 100644 BEAMS3D/Sources/beams3d_beam_density.f90 diff --git a/BEAMS3D/BEAMS3D.dep b/BEAMS3D/BEAMS3D.dep index 8a893719a..03bea0f4f 100644 --- a/BEAMS3D/BEAMS3D.dep +++ b/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 \ diff --git a/BEAMS3D/ObjectList b/BEAMS3D/ObjectList index a5c0e44ad..1dc3af942 100644 --- a/BEAMS3D/ObjectList +++ b/BEAMS3D/ObjectList @@ -1,4 +1,5 @@ ObjectFiles = \ +beams3d_beam_density.o \ beams3d_duplicate_part.o \ outpart_beams3d_nag.o \ beams3d_follow.o \ diff --git a/BEAMS3D/Sources/beams3d_beam_density.f90 b/BEAMS3D/Sources/beams3d_beam_density.f90 new file mode 100644 index 000000000..bc1e5b352 --- /dev/null +++ b/BEAMS3D/Sources/beams3d_beam_density.f90 @@ -0,0 +1,117 @@ +!----------------------------------------------------------------------- +! 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, xt, yt, zt, rt, pt, dV, & + dx, dy, dz , dgrid + LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE :: lgrid +#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) + ALLOCATE(lgrid(nbeams,nr,nphi,nz)) + + IF (lverb) WRITE(6,'(A)') '----- BEAM Density Calc. -----' + + dgrid = sqrt(hr(1)**2+hz(1)**2) + 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 jsut need to calc the volume of the voxel to get particles/m^3 + d = SQRT((x1-x0)**2 + (y1-y0)**2 + (z1-z0)**2) + denbeam = weight(l)*L/vll_lines(0,l) ! This is the total number of particles + dl = d/dgrid + nl = NINT(L/dl)+1 + dl = d/nl + dx = (x1-x0)/(nl-1) + dy = (y1-y0)/(nl-1) + dz = (z1-z0)/(nl-1) + LGRID = .FALSE. + DO s = 1, nl + xt = x0 + (i-1)*dx + yt = y0 + (i-1)*dy + zt = z0 + (i-1)*dz + rt = SQRT(xt*xt+yt*yt) + IF ( (rt > raxis(nr-1)) .or. (rt < raxis(2)) .or. & + (zt > zaxis(nz-1)) .or. (zt < zaxis(2)) ) CYCLE + pt = ATAN2(yt,xt) + pt = MODULO(pt, phimax) + ! Find gridpoint but use half grid + i = COUNT(raxis-hr2 < rt) + j = COUNT(phiaxis-hp2 < pt) + k = COUNT(zaxis-hz2 < zt) + LGRID(m,i,j,k) = .true. + END DO + WHERE(lgrid) BEAM_DENSITY = BEAM_DENSITY + denbeam + END DO + DEALLOCATE(lgrid) + + ! Barrier so calculation is done +#if defined(MPI_OPT) + CALL MPI_BARRIER(MPI_COMM_LOCAL,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_LOCAL,ierr_mpi) +#endif + + RETURN + +!----------------------------------------------------------------------- +! End Subroutine +!----------------------------------------------------------------------- + END SUBROUTINE beams3d_beam_density diff --git a/BEAMS3D/Sources/beams3d_fix_poloidal.f90 b/BEAMS3D/Sources/beams3d_fix_poloidal.f90 index d8e7b13a1..fbdb20f9b 100644 --- a/BEAMS3D/Sources/beams3d_fix_poloidal.f90 +++ b/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 diff --git a/BEAMS3D/Sources/beams3d_follow.f90 b/BEAMS3D/Sources/beams3d_follow.f90 index 4933244d0..46aadfb7b 100644 --- a/BEAMS3D/Sources/beams3d_follow.f90 +++ b/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 @@ -169,6 +170,7 @@ SUBROUTINE beams3d_follow ! Beam Deposition IF (lbeam) THEN + IF (lverb) WRITE(6,'(A)') '----- BEAM Deposition -----' DO i = mystart, myend lneut = lbeam ltherm = .FALSE. @@ -220,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) CALL beams3d_beam_density + ! Follow Trajectories IF (.not.ldepo) THEN CALL beams3d_follow_gc @@ -269,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) diff --git a/BEAMS3D/Sources/beams3d_grid.f90 b/BEAMS3D/Sources/beams3d_grid.f90 index 8710d2774..8e99d88c1 100644 --- a/BEAMS3D/Sources/beams3d_grid.f90 +++ b/BEAMS3D/Sources/beams3d_grid.f90 @@ -63,7 +63,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, & diff --git a/BEAMS3D/Sources/beams3d_write.f90 b/BEAMS3D/Sources/beams3d_write.f90 index 7777537f2..35e2c8fba 100644 --- a/BEAMS3D/Sources/beams3d_write.f90 +++ b/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, & @@ -210,6 +211,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) From c6da6c7e99608c5730d025a3556b0f4d8cb5d6d7 Mon Sep 17 00:00:00 2001 From: lazersos Date: Sun, 27 Aug 2023 13:37:44 +0200 Subject: [PATCH 2/9] BEAMS3D: Attempts to speed up beam density calculation. --- BEAMS3D/Sources/beams3d_beam_density.f90 | 49 +++++++++++++----------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/BEAMS3D/Sources/beams3d_beam_density.f90 b/BEAMS3D/Sources/beams3d_beam_density.f90 index bc1e5b352..67efcd2ee 100644 --- a/BEAMS3D/Sources/beams3d_beam_density.f90 +++ b/BEAMS3D/Sources/beams3d_beam_density.f90 @@ -26,9 +26,9 @@ SUBROUTINE beams3d_beam_density 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, xt, yt, zt, rt, pt, dV, & - dx, dy, dz , dgrid + denbeam, dl, dV LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE :: lgrid + DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xl,yl,zl,rl,sl,pl #if defined(MPI_OPT) INTEGER :: MPI_COMM_LOCAL INTEGER :: mystart, mypace @@ -39,11 +39,14 @@ SUBROUTINE beams3d_beam_density ! 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 ALLOCATE(lgrid(nbeams,nr,nphi,nz)) + 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-1)/DBLE(nl-1) - IF (lverb) WRITE(6,'(A)') '----- BEAM Density Calc. -----' - dgrid = sqrt(hr(1)**2+hz(1)**2) + IF (lverb) WRITE(6,'(A)') '----- BEAM Density Calc. -----' hr2 = 0.5*hr(1) hz2 = 0.5*hz(1) hp2 = 0.5*hp(1) @@ -62,32 +65,34 @@ SUBROUTINE beams3d_beam_density ! The the total number of particles is just n=W*t = W*L/v ! then we jsut need to calc the volume of the voxel to get particles/m^3 d = SQRT((x1-x0)**2 + (y1-y0)**2 + (z1-z0)**2) - denbeam = weight(l)*L/vll_lines(0,l) ! This is the total number of particles - dl = d/dgrid - nl = NINT(L/dl)+1 - dl = d/nl - dx = (x1-x0)/(nl-1) - dy = (y1-y0)/(nl-1) - dz = (z1-z0)/(nl-1) + denbeam = weight(l)*d/vll_lines(0,l) ! This is the total number of particles + 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) LGRID = .FALSE. DO s = 1, nl - xt = x0 + (i-1)*dx - yt = y0 + (i-1)*dy - zt = z0 + (i-1)*dz - rt = SQRT(xt*xt+yt*yt) - IF ( (rt > raxis(nr-1)) .or. (rt < raxis(2)) .or. & - (zt > zaxis(nz-1)) .or. (zt < zaxis(2)) ) CYCLE - pt = ATAN2(yt,xt) - pt = MODULO(pt, phimax) ! Find gridpoint but use half grid - i = COUNT(raxis-hr2 < rt) - j = COUNT(phiaxis-hp2 < pt) - k = COUNT(zaxis-hz2 < zt) + 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) LGRID(m,i,j,k) = .true. END DO WHERE(lgrid) BEAM_DENSITY = BEAM_DENSITY + denbeam END DO DEALLOCATE(lgrid) + 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) From f594889954e4d82ca4315a44979e2e42bc53f7cf Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 6 Sep 2023 21:03:51 +0200 Subject: [PATCH 3/9] BEAMS3D: Added option for beam density calculation to be optional. --- BEAMS3D/Sources/beams3d_follow.f90 | 2 +- BEAMS3D/Sources/beams3d_interface_mod.f90 | 6 ++++++ BEAMS3D/Sources/beams3d_runtime.f90 | 2 +- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/BEAMS3D/Sources/beams3d_follow.f90 b/BEAMS3D/Sources/beams3d_follow.f90 index 46aadfb7b..d0f7c03b1 100644 --- a/BEAMS3D/Sources/beams3d_follow.f90 +++ b/BEAMS3D/Sources/beams3d_follow.f90 @@ -226,7 +226,7 @@ SUBROUTINE beams3d_follow #if defined(MPI_OPT) CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi) #endif - IF (lbeam) CALL beams3d_beam_density + IF (lbeam .and. lbeamdensity) CALL beams3d_beam_density ! Follow Trajectories IF (.not.ldepo) THEN diff --git a/BEAMS3D/Sources/beams3d_interface_mod.f90 b/BEAMS3D/Sources/beams3d_interface_mod.f90 index e8173b37c..6ec8529df 100644 --- a/BEAMS3D/Sources/beams3d_interface_mod.f90 +++ b/BEAMS3D/Sources/beams3d_interface_mod.f90 @@ -157,6 +157,7 @@ SUBROUTINE beams3d_init_commandline lfusion_alpha = .false. lboxsim = .false. lfieldlines = .false. + lbeamdensity = .true. id_string = '' coil_string = '' mgrid_string = '' @@ -264,6 +265,8 @@ SUBROUTINE beams3d_init_commandline lfusion_alpha = .true. case ("-boxsim") lboxsim = .true. + case ("-nobeamdensity") + lbeamdensity = .false. case ("-help", "-h") ! Output Help message write(6, *) ' Beam MC Code' write(6, *) ' Usage: xbeams3d ' @@ -293,6 +296,7 @@ SUBROUTINE beams3d_init_commandline write(6, *) ' -fusion: Fusion Reaction Rates for birth' write(6, *) ' -fusion_alpha: Fusion Reaction Rates for birth (alphas 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 @@ -378,6 +382,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 diff --git a/BEAMS3D/Sources/beams3d_runtime.f90 b/BEAMS3D/Sources/beams3d_runtime.f90 index 663ff9ee6..48f865b24 100644 --- a/BEAMS3D/Sources/beams3d_runtime.f90 +++ b/BEAMS3D/Sources/beams3d_runtime.f90 @@ -134,7 +134,7 @@ MODULE beams3d_runtime ldepo, lbeam_simple, ldebug, lcollision, lw7x, lsuzuki, & lascot, lascot4, lbbnbi, lfidasim, lfidasim2, lsplit, lvessel_beam, lascotfl, lrandomize, & lfusion, lfusion_alpha, leqdsk, lhint, lkick, lgcsim, & - lboxsim, limas, lfieldlines + lboxsim, limas, lfieldlines, lbeamdensity INTEGER :: nextcur, npoinc, nbeams, nparticles_start, nprocs_beams, & ndt, ndt_max, duplicate_factor INTEGER, DIMENSION(MAXBEAMS) :: Dex_beams From e7c674dda4e20bedc838aa6ae416f19ba6f6978c Mon Sep 17 00:00:00 2001 From: lazersos Date: Wed, 6 Sep 2023 21:04:27 +0200 Subject: [PATCH 4/9] BENCHMARKS: Turned density calcualtion off for depo runs. --- BENCHMARKS/makefile | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/BENCHMARKS/makefile b/BENCHMARKS/makefile index dd611956b..12c244cd9 100644 --- a/BENCHMARKS/makefile +++ b/BENCHMARKS/makefile @@ -136,20 +136,22 @@ test_beams3d_slow: test_beams3d_depo: @cd $(BEAMS3D_TEST_DIR); rm -f *.h5 *.txt; - @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_depo -plasma -depo -suzuki - @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_depo -plasma -depo -suzuki + @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_depo -plasma -depo -suzuki -nobeamdensity + @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_depo -plasma -depo -suzuki -nobeamdensity @cd $(BEAMS3D_TEST_DIR); ./test_ORBITS_depo.py @echo "" test_beams3d_depo_adas: @cd $(BEAMS3D_TEST_DIR); rm -f *.h5 *.txt; - @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_depo_adas -plasma -depo - @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_depo_adas -plasma -depo + @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_depo_adas -plasma -depo -nobeamdensity + @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_depo_adas -plasma -depo -nobeamdensity @cd $(BEAMS3D_TEST_DIR); ./test_ORBITS_depo_adas.py @echo "" -test_beams3d_restart: +beams3d_ORBITS.h5: $(MAKE) test_beams3d_depo + +test_beams3d_restart: beams3d_ORBITS.h5 @cd $(BEAMS3D_TEST_DIR); @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_restart -plasma -restart beams3d_ORBITS_depo.h5 @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_restart -plasma -restart beams3d_ORBITS_depo.h5 @@ -172,15 +174,15 @@ test_beams3d_fusion: test_beams3d_eqdsk: @cd $(BEAMS3D_TEST_DIR); rm -f *.h5 *.txt; - @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -eqdsk ORBITS_eqdsk g164723.03059 -depo -suzuki - @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -eqdsk ORBITS_eqdsk g164723.03059 -depo -suzuki + @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -eqdsk ORBITS_eqdsk g164723.03059 -depo -suzuki -nobeamdensity + @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -eqdsk ORBITS_eqdsk g164723.03059 -depo -suzuki -nobeamdensity @cd $(BEAMS3D_TEST_DIR); ./test_ORBITS_eqdsk.py @echo "" test_beams3d_multiion: @cd $(BEAMS3D_TEST_DIR); rm -f *.h5 *.txt; - @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_multiion -plasma -depo -suzuki - @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_multiion -plasma -depo -suzuki + @echo $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_multiion -plasma -depo -suzuki -nobeamdensity + @cd $(BEAMS3D_TEST_DIR); $(MPI_RUN) $(MPI_RUN_OPTS_SM) $(MYHOME)/xbeams3d -vmec ORBITS_multiion -plasma -depo -suzuki -nobeamdensity @cd $(BEAMS3D_TEST_DIR); ./test_ORBITS_multiion.py @echo "" From 6f61d0177f529af5a654bd642d335414eebd1089 Mon Sep 17 00:00:00 2001 From: kudav Date: Thu, 7 Sep 2023 15:27:03 +0200 Subject: [PATCH 5/9] BEAMS3D: MPI comm fix --- BEAMS3D/Sources/beams3d_beam_density.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/BEAMS3D/Sources/beams3d_beam_density.f90 b/BEAMS3D/Sources/beams3d_beam_density.f90 index 67efcd2ee..931f4bc69 100644 --- a/BEAMS3D/Sources/beams3d_beam_density.f90 +++ b/BEAMS3D/Sources/beams3d_beam_density.f90 @@ -96,9 +96,8 @@ SUBROUTINE beams3d_beam_density ! Barrier so calculation is done #if defined(MPI_OPT) - CALL MPI_BARRIER(MPI_COMM_LOCAL,ierr_mpi) + 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 @@ -111,7 +110,7 @@ SUBROUTINE beams3d_beam_density END DO #if defined(MPI_OPT) - CALL MPI_BARRIER(MPI_COMM_LOCAL,ierr_mpi) + CALL MPI_BARRIER(MPI_COMM_SHARMEM,ierr_mpi) #endif RETURN From bd914a7dfd4f7eee6a039e05dd3a5f3dbbbed80b Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 6 Oct 2023 11:12:16 +0200 Subject: [PATCH 6/9] BEAMS3D: Added error message about grid size --- BEAMS3D/Sources/beams3d_init.f90 | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/BEAMS3D/Sources/beams3d_init.f90 b/BEAMS3D/Sources/beams3d_init.f90 index 2b8151cb2..ce6ab90a1 100644 --- a/BEAMS3D/Sources/beams3d_init.f90 +++ b/BEAMS3D/Sources/beams3d_init.f90 @@ -737,6 +737,18 @@ 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 + print *,phiaxis(nphi)-phiaxis(1),phitemp + 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 From afaedb121b27304392004af094afed608aed32cf Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 6 Oct 2023 11:29:14 +0200 Subject: [PATCH 7/9] BEAMS3D: Fixed issue where negative potential may not be printed. --- BEAMS3D/Sources/beams3d_init.f90 | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/BEAMS3D/Sources/beams3d_init.f90 b/BEAMS3D/Sources/beams3d_init.f90 index ce6ab90a1..464a66c5c 100644 --- a/BEAMS3D/Sources/beams3d_init.f90 +++ b/BEAMS3D/Sources/beams3d_init.f90 @@ -254,8 +254,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 From a2b81cc836ee373f38a1fcaaf388248110159a5d Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 6 Oct 2023 11:30:03 +0200 Subject: [PATCH 8/9] BEAMS3D: Small cleanup of screen output. --- BEAMS3D/Sources/beams3d_init.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/BEAMS3D/Sources/beams3d_init.f90 b/BEAMS3D/Sources/beams3d_init.f90 index 464a66c5c..d283a070f 100644 --- a/BEAMS3D/Sources/beams3d_init.f90 +++ b/BEAMS3D/Sources/beams3d_init.f90 @@ -751,7 +751,6 @@ SUBROUTINE beams3d_init IF (lverb) THEN phitemp = MAXVAL(phi_start)-MINVAL(phi_start) IF ((phitemp < hp(1)) .and. (phitemp > 0) .and. lbeamdensity) THEN - print *,phiaxis(nphi)-phiaxis(1),phitemp i = (phiaxis(nphi)-phiaxis(1))/phitemp WRITE(6,'(A)') ' WARNING: The range of PHI_START > dphi.' WRITE(6,'(A,I4)') ' Consider increasing nphi >=',i From 8b3392ef492d974c622885c646bf5bb59d2a5691 Mon Sep 17 00:00:00 2001 From: lazersos Date: Fri, 6 Oct 2023 11:51:19 +0200 Subject: [PATCH 9/9] BEAMS3D: Fixed density calculation so better desnity is calculated. --- BEAMS3D/Sources/beams3d_beam_density.f90 | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/BEAMS3D/Sources/beams3d_beam_density.f90 b/BEAMS3D/Sources/beams3d_beam_density.f90 index 931f4bc69..44202129e 100644 --- a/BEAMS3D/Sources/beams3d_beam_density.f90 +++ b/BEAMS3D/Sources/beams3d_beam_density.f90 @@ -27,7 +27,6 @@ SUBROUTINE beams3d_beam_density INTEGER :: ier, l, nl, i, j, k, m, s DOUBLE PRECISION :: x0, y0, z0, x1, y1, z1, hr2, hz2, hp2, d, & denbeam, dl, dV - LOGICAL, DIMENSION(:,:,:,:), ALLOCATABLE :: lgrid DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xl,yl,zl,rl,sl,pl #if defined(MPI_OPT) INTEGER :: MPI_COMM_LOCAL @@ -40,10 +39,9 @@ SUBROUTINE beams3d_beam_density ! 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 - ALLOCATE(lgrid(nbeams,nr,nphi,nz)) 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-1)/DBLE(nl-1) + FORALL(s=1:nl) sl(s) = DBLE(s-0.5)/DBLE(nl) ! half grid IF (lverb) WRITE(6,'(A)') '----- BEAM Density Calc. -----' @@ -63,26 +61,23 @@ SUBROUTINE beams3d_beam_density ! 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 jsut need to calc the volume of the voxel to get particles/m^3 - d = SQRT((x1-x0)**2 + (y1-y0)**2 + (z1-z0)**2) - denbeam = weight(l)*d/vll_lines(0,l) ! This is the total number of particles + ! 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) - LGRID = .FALSE. 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) - LGRID(m,i,j,k) = .true. + BEAM_DENSITY(m,i,j,k) = BEAM_DENSITY(m,i,j,k) + denbeam END DO - WHERE(lgrid) BEAM_DENSITY = BEAM_DENSITY + denbeam END DO - DEALLOCATE(lgrid) DEALLOCATE(xl,yl,zl,rl,pl,sl) ! Fix edges which are double counts