diff --git a/BEAMS3D/Sources/beams3d_distnorm.f90 b/BEAMS3D/Sources/beams3d_distnorm.f90 index 84740244..c5362879 100644 --- a/BEAMS3D/Sources/beams3d_distnorm.f90 +++ b/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, & diff --git a/BEAMS3D/Sources/beams3d_fix_poloidal.f90 b/BEAMS3D/Sources/beams3d_fix_poloidal.f90 index fbdb20f9..cf29426f 100644 --- a/BEAMS3D/Sources/beams3d_fix_poloidal.f90 +++ b/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 \ No newline at end of file + END SUBROUTINE beams3d_fix_poloidal diff --git a/BEAMS3D/Sources/beams3d_follow.f90 b/BEAMS3D/Sources/beams3d_follow.f90 index 79620f8a..d34f4a13 100644 --- a/BEAMS3D/Sources/beams3d_follow.f90 +++ b/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) diff --git a/BEAMS3D/Sources/beams3d_follow_fo.f90 b/BEAMS3D/Sources/beams3d_follow_fo.f90 index 544148c6..d8ec4866 100644 --- a/BEAMS3D/Sources/beams3d_follow_fo.f90 +++ b/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 diff --git a/BEAMS3D/Sources/beams3d_follow_gc.f90 b/BEAMS3D/Sources/beams3d_follow_gc.f90 index c6a29ab5..07af522a 100644 --- a/BEAMS3D/Sources/beams3d_follow_gc.f90 +++ b/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 diff --git a/BEAMS3D/Sources/beams3d_free.f90 b/BEAMS3D/Sources/beams3d_free.f90 index a15339a1..52a6a3aa 100644 --- a/BEAMS3D/Sources/beams3d_free.f90 +++ b/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) diff --git a/BEAMS3D/Sources/beams3d_grid.f90 b/BEAMS3D/Sources/beams3d_grid.f90 index 683c62e7..0df67c50 100644 --- a/BEAMS3D/Sources/beams3d_grid.f90 +++ b/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 diff --git a/BEAMS3D/Sources/beams3d_init.f90 b/BEAMS3D/Sources/beams3d_init.f90 index 69908e3a..57d84800 100644 --- a/BEAMS3D/Sources/beams3d_init.f90 +++ b/BEAMS3D/Sources/beams3d_init.f90 @@ -342,77 +342,75 @@ SUBROUTINE beams3d_init CALL mpialloc(energy_fida, nenergy_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_energy_fida) CALL mpialloc(pitch_fida, npitch_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_pitch_fida) END IF - !IF (lcontinue_grid) THEN - !CALL beams3d_init_restart - !ELSE - ! Create the background grid - CALL mpialloc(raxis, nr, myid_sharmem, 0, MPI_COMM_SHARMEM, win_raxis) - CALL mpialloc(phiaxis, nphi, myid_sharmem, 0, MPI_COMM_SHARMEM, win_phiaxis) - CALL mpialloc(zaxis, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_zaxis) - CALL mpialloc(hr, nr-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hr) - CALL mpialloc(hp, nphi-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hp) - CALL mpialloc(hz, nz-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hz) - CALL mpialloc(hri, nr-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hri) - CALL mpialloc(hpi, nphi-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hpi) - CALL mpialloc(hzi, nz-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hzi) - CALL mpialloc(B_R, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_B_R) - CALL mpialloc(B_PHI, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_B_PHI) - CALL mpialloc(B_Z, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_B_Z) - CALL mpialloc(MODB, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_MODB) - CALL mpialloc(TE, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_TE) - CALL mpialloc(NE, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_NE) - CALL mpialloc(TI, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_TI) - CALL mpialloc(ZEFF_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_ZEFF_ARR) - CALL mpialloc(POT_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_POT_ARR) - CALL mpialloc(S_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_S_ARR) - CALL mpialloc(U_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_U_ARR) - CALL mpialloc(X_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_X_ARR) - CALL mpialloc(Y_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_Y_ARR) - CALL mpialloc(NI, NION, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_NI) - IF (lfidasim) THEN - CALL mpialloc(raxis_fida, nr_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_raxis_fida) - CALL mpialloc(phiaxis_fida, nphi_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_phiaxis_fida) - CALL mpialloc(zaxis_fida, nz_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_zaxis_fida) - CALL mpialloc(energy_fida, nenergy_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_energy_fida) - CALL mpialloc(pitch_fida, npitch_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_pitch_fida) - END IF - IF (myid_sharmem == 0) THEN - FORALL(i = 1:nr) raxis(i) = (i-1)*(rmax-rmin)/(nr-1) + rmin - FORALL(i = 1:nz) zaxis(i) = (i-1)*(zmax-zmin)/(nz-1) + zmin - FORALL(i = 1:nphi) phiaxis(i) = (i-1)*(phimax-phimin)/(nphi-1) + phimin - S_ARR = 1.5 - X_ARR = 1.5 - Y_ARR = 1.5 - POT_ARR = 0 - NI = 0 - ! Setup grid helpers - ! Note: All helpers are defined in terms of differences on half grid - ! so values are indexed from 1 to n-1. Which we store at n - ! i = MIN(MAX(COUNT(raxis < r_temp),1),nr-1) - ! hr(i) = raxis(i+1) - raxis(i) - ! hri = one / hr - FORALL(i = 1:nr-1) hr(i) = raxis(i+1) - raxis(i) - FORALL(i = 1:nz-1) hz(i) = zaxis(i+1) - zaxis(i) - FORALL(i = 1:nphi-1) hp(i) = phiaxis(i+1) - phiaxis(i) - hri = one / hr - hpi = one / hp - hzi = one / hz - ! Do this here so EQDSK vac RMP works. - B_R = 0 - B_PHI = 0 - B_Z = 0 - MODB = 0 + ! Create the background grid + CALL mpialloc(raxis, nr, myid_sharmem, 0, MPI_COMM_SHARMEM, win_raxis) + CALL mpialloc(phiaxis, nphi, myid_sharmem, 0, MPI_COMM_SHARMEM, win_phiaxis) + CALL mpialloc(zaxis, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_zaxis) + CALL mpialloc(hr, nr-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hr) + CALL mpialloc(hp, nphi-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hp) + CALL mpialloc(hz, nz-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hz) + CALL mpialloc(hri, nr-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hri) + CALL mpialloc(hpi, nphi-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hpi) + CALL mpialloc(hzi, nz-1, myid_sharmem, 0, MPI_COMM_SHARMEM, win_hzi) + CALL mpialloc(B_R, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_B_R) + CALL mpialloc(B_PHI, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_B_PHI) + CALL mpialloc(B_Z, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_B_Z) + CALL mpialloc(MODB, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_MODB) + CALL mpialloc(TE, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_TE) + CALL mpialloc(NE, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_NE) + CALL mpialloc(TI, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_TI) + CALL mpialloc(ZEFF_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_ZEFF_ARR) + CALL mpialloc(POT_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_POT_ARR) + CALL mpialloc(S_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_S_ARR) + CALL mpialloc(RHO_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_RHO_ARR) + CALL mpialloc(U_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_U_ARR) + CALL mpialloc(XRHO_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_XRHO_ARR) + CALL mpialloc(YRHO_ARR, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_YRHO_ARR) + CALL mpialloc(NI, NION, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_NI) + IF (lfidasim) THEN + CALL mpialloc(raxis_fida, nr_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_raxis_fida) + CALL mpialloc(phiaxis_fida, nphi_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_phiaxis_fida) + CALL mpialloc(zaxis_fida, nz_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_zaxis_fida) + CALL mpialloc(energy_fida, nenergy_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_energy_fida) + CALL mpialloc(pitch_fida, npitch_fida, myid_sharmem, 0, MPI_COMM_SHARMEM, win_pitch_fida) + END IF + IF (myid_sharmem == 0) THEN + FORALL(i = 1:nr) raxis(i) = (i-1)*(rmax-rmin)/(nr-1) + rmin + FORALL(i = 1:nz) zaxis(i) = (i-1)*(zmax-zmin)/(nz-1) + zmin + FORALL(i = 1:nphi) phiaxis(i) = (i-1)*(phimax-phimin)/(nphi-1) + phimin + S_ARR = 1.5 + RHO_ARR = 1.5 + XRHO_ARR = 1.5 + YRHO_ARR = 1.5 + POT_ARR = 0 + NI = 0 + ! Setup grid helpers + ! Note: All helpers are defined in terms of differences on half grid + ! so values are indexed from 1 to n-1. Which we store at n + ! i = MIN(MAX(COUNT(raxis < r_temp),1),nr-1) + ! hr(i) = raxis(i+1) - raxis(i) + ! hri = one / hr + FORALL(i = 1:nr-1) hr(i) = raxis(i+1) - raxis(i) + FORALL(i = 1:nz-1) hz(i) = zaxis(i+1) - zaxis(i) + FORALL(i = 1:nphi-1) hp(i) = phiaxis(i+1) - phiaxis(i) + hri = one / hr + hpi = one / hp + hzi = one / hz + ! Do this here so EQDSK vac RMP works. + B_R = 0 + B_PHI = 0 + B_Z = 0 + MODB = 0 - END IF - CALL MPI_BARRIER(MPI_COMM_SHARMEM, ier) + END IF + CALL MPI_BARRIER(MPI_COMM_SHARMEM, ier) - ! Put the vacuum field on the background grid - IF (lmgrid) THEN - CALL beams3d_init_mgrid - ELSE IF (lcoil) THEN - CALL beams3d_init_coil - END IF - !END IF + ! Put the vacuum field on the background grid + IF (lmgrid) THEN + CALL beams3d_init_mgrid + ELSE IF (lcoil) THEN + CALL beams3d_init_coil + END IF ! Put the plasma field on the background grid IF (lvmec .and. .not.lvac) THEN @@ -552,8 +550,6 @@ SUBROUTINE beams3d_init ! Construct MODB IF (myid_sharmem == master) MODB = SQRT(B_R*B_R+B_PHI*B_PHI+B_Z*B_Z) - - ! Construct Splines on shared memory master nodes IF (myid_sharmem == master) THEN bcs1=(/ 0, 0/) @@ -567,16 +563,16 @@ SUBROUTINE beams3d_init IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:BZ_spl',ier) CALL EZspline_init(MODB_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:MODB_spl',ier) - CALL EZspline_init(S_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) - IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:S_spl',ier) CALL EZspline_init(U_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:U_spl',ier) - CALL EZspline_init(X_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) - IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:X_spl',ier) - CALL EZspline_init(Y_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) - IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:Y_spl',ier) CALL EZspline_init(POT_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:POT_spl',ier) + CALL EZspline_init(RHO_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) + IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:RHO_spl',ier) + CALL EZspline_init(XRHO_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) + IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:XRHO_spl',ier) + CALL EZspline_init(YRHO_spl,nr,nphi,nz,bcs1,bcs2,bcs3,ier) + IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:YRHO_spl',ier) BR_spl%isHermite = 1 BR_spl%x1 = raxis BR_spl%x2 = phiaxis @@ -593,26 +589,30 @@ SUBROUTINE beams3d_init MODB_spl%x1 = raxis MODB_spl%x2 = phiaxis MODB_spl%x3 = zaxis - S_spl%isHermite = 1 - S_spl%x1 = raxis - S_spl%x2 = phiaxis - S_spl%x3 = zaxis U_spl%isHermite = 1 U_spl%x1 = raxis U_spl%x2 = phiaxis U_spl%x3 = zaxis - X_spl%isHermite = 1 - X_spl%x1 = raxis - X_spl%x2 = phiaxis - X_spl%x3 = zaxis - Y_spl%isHermite = 1 - Y_spl%x1 = raxis - Y_spl%x2 = phiaxis - Y_spl%x3 = zaxis POT_spl%isHermite = 1 POT_spl%x1 = raxis POT_spl%x2 = phiaxis POT_spl%x3 = zaxis + RHO_spl%isHermite = 1 + RHO_spl%x1 = raxis + RHO_spl%x2 = phiaxis + RHO_spl%x3 = zaxis + XRHO_spl%isHermite = 1 + XRHO_spl%x1 = raxis + XRHO_spl%x2 = phiaxis + XRHO_spl%x3 = zaxis + YRHO_spl%isHermite = 1 + YRHO_spl%x1 = raxis + YRHO_spl%x2 = phiaxis + YRHO_spl%x3 = zaxis + ! Make sure we setup RHO from S + RHO_ARR = SQRT(S_ARR) + XRHO_ARR = RHO_ARR * COS(U_ARR) + YRHO_ARR = RHO_ARR * SIN(U_ARR) CALL EZspline_setup(BR_spl,B_R,ier,EXACT_DIM=.true.) IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:BR_spl',ier) CALL EZspline_setup(BPHI_spl,B_PHI,ier,EXACT_DIM=.true.) @@ -621,12 +621,16 @@ SUBROUTINE beams3d_init IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:BZ_spl',ier) CALL EZspline_setup(MODB_spl,MODB,ier,EXACT_DIM=.true.) IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:MODB_spl',ier) - CALL EZspline_setup(S_spl,S_ARR,ier,EXACT_DIM=.true.) - IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:S_spl',ier) CALL EZspline_setup(U_spl,U_ARR,ier,EXACT_DIM=.true.) IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:U_spl',ier) CALL EZspline_setup(POT_spl,POT_ARR,ier,EXACT_DIM=.true.) IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:POT_spl',ier) + CALL EZspline_setup(RHO_spl,RHO_ARR,ier,EXACT_DIM=.true.) + IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:RHO_spl',ier) + CALL EZspline_setup(XRHO_spl,XRHO_ARR,ier,EXACT_DIM=.true.) + IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:XRHO_spl',ier) + CALL EZspline_setup(YRHO_spl,YRHO_ARR,ier,EXACT_DIM=.true.) + IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:YRHO_spl',ier) END IF ! Allocate Shared memory space @@ -635,10 +639,10 @@ SUBROUTINE beams3d_init CALL mpialloc(BPHI4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_BPHI4D) CALL mpialloc(BZ4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_BZ4D) CALL mpialloc(MODB4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_MODB4D) - CALL mpialloc(S4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_S4D) + CALL mpialloc(RHO4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_RHO4D) CALL mpialloc(U4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_U4D) - CALL mpialloc(X4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_X4D) - CALL mpialloc(Y4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_Y4D) + CALL mpialloc(XRHO4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_XRHO4D) + CALL mpialloc(YRHO4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_YRHO4D) CALL mpialloc(POT4D, 8, nr, nphi, nz, myid_sharmem, 0, MPI_COMM_SHARMEM, win_POT4D) ! Copy Spline info to shared memory and Free IF (myid_sharmem == master) THEN @@ -646,29 +650,20 @@ SUBROUTINE beams3d_init BPHI4D = BPHI_SPL%fspl BZ4D = BZ_SPL%fspl MODB4D = MODB_SPL%fspl - S4D = S_SPL%fspl U4D = U_SPL%fspl POT4D = POT_SPL%fspl - - X_ARR = S4D(1,:,:,:) * COS(U4D(1,:,:,:)) - Y_ARR = S4D(1,:,:,:) * SIN(U4D(1,:,:,:)) - CALL EZspline_setup(X_spl,X_ARR,ier,EXACT_DIM=.true.) - IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:X_spl',ier) - CALL EZspline_setup(Y_spl,Y_ARR,ier,EXACT_DIM=.true.) - IF (ier /=0) CALL handle_err(EZSPLINE_ERR,'beams3d_init:Y_spl',ier) - X4D = X_SPL%fspl - Y4D = Y_SPL%fspl - + RHO4D = RHO_SPL%fspl + XRHO4D = XRHO_SPL%fspl + YRHO4D = YRHO_SPL%fspl CALL EZspline_free(BR_spl,ier) CALL EZspline_free(BPHI_spl,ier) CALL EZspline_free(BZ_spl,ier) CALL EZspline_free(MODB_spl,ier) - CALL EZspline_free(S_spl,ier) CALL EZspline_free(U_spl,ier) - CALL EZspline_free(X_spl,ier) - CALL EZspline_free(Y_spl,ier) - CALL EZspline_free(POT_spl,ier) + CALL EZspline_free(RHO_spl,ier) + CALL EZspline_free(XRHO_spl,ier) + CALL EZspline_free(YRHO_spl,ier) END IF ! These are helpers for range eps1 = (rmax-rmin)*small @@ -696,9 +691,10 @@ SUBROUTINE beams3d_init CALL mpidealloc(MODB,win_MODB) CALL mpidealloc(S_ARR,win_S_ARR) CALL mpidealloc(U_ARR,win_U_ARR) - CALL mpidealloc(X_ARR,win_X_ARR) - CALL mpidealloc(Y_ARR,win_Y_ARR) CALL mpidealloc(POT_ARR,win_POT_ARR) + CALL mpidealloc(RHO_ARR,win_RHO_ARR) + CALL mpidealloc(XRHO_ARR,win_XRHO_ARR) + CALL mpidealloc(YRHO_ARR,win_YRHO_ARR) IF (.not. lvac) THEN CALL mpidealloc(TE,win_TE) CALL mpidealloc(NE,win_NE) diff --git a/BEAMS3D/Sources/beams3d_init_fusion.f90 b/BEAMS3D/Sources/beams3d_init_fusion.f90 index 4f5d47ce..2f0a8f6a 100644 --- a/BEAMS3D/Sources/beams3d_init_fusion.f90 +++ b/BEAMS3D/Sources/beams3d_init_fusion.f90 @@ -12,7 +12,7 @@ SUBROUTINE beams3d_init_fusion USE stel_kinds, ONLY: rprec USE beams3d_runtime USE beams3d_grid, ONLY: nr, nphi, nz, hr, hp, hz, raxis, zaxis, & - phiaxis, S4D, U4D, dexionHe3, & + phiaxis, RHO4D, U4D, dexionHe3, & NEUTRONS_ARR, win_NEUTRONS, & E_NEUTRONS, win_E_NEUTRONS USE beams3d_lines, ONLY: nparticles, partvmax @@ -335,9 +335,9 @@ SUBROUTINE beams3d_init_fusion ! Calc u = [0,2*pi] utemp = Z1_rand*pi2 ! Calc sval =[0,1] - rtemp = ABS(X1_rand*X1_rand-roffset) + rtemp = ABS(X1_rand-roffset) ! Now we need to find the proper point - f2d = ((U4D(1,:,j,:)-utemp)**2)*((S4D(1,:,j,:)-rtemp)**2) + f2d = ((U4D(1,:,j,:)-utemp)**2)*((RHO4D(1,:,j,:)-rtemp)**2) minik = MINLOC(f2d) i = MIN(MAX(minik(1),2),nr1); k = MIN(MAX(minik(2),2),nz1) IF (l3d(i,j,k)) EXIT diff --git a/BEAMS3D/Sources/beams3d_interface_mod.f90 b/BEAMS3D/Sources/beams3d_interface_mod.f90 index b64b5900..9aa04e52 100644 --- a/BEAMS3D/Sources/beams3d_interface_mod.f90 +++ b/BEAMS3D/Sources/beams3d_interface_mod.f90 @@ -76,10 +76,12 @@ SUBROUTINE beams3d_init_pointers IMPLICIT NONE ! Nullify pointers NULLIFY(raxis,phiaxis,zaxis,hr,hp,hz,hri,hpi,hzi,B_R,B_PHI,B_Z, & - MODB,TE,NE,TI,ZEFF_ARR,POT_ARR,S_ARR,U_ARR,X_ARR,Y_ARR,NI, & + MODB,TE,NE,TI,ZEFF_ARR,POT_ARR,S_ARR,U_ARR,NI, & + RHO_ARR,XRHO_ARR,YRHO_ARR, & raxis_fida,zaxis_fida,phiaxis_fida,energy_fida,pitch_fida, & req_axis,zeq_axis,TE4D,NE4D,TI4D,ZEFF4D,NI5D,BR4D,BPHI4D, & - BZ4D,MODB4D,S4D,U4D,X4D,Y4D,POT4D,dist5d_prof,dist5d_fida, & + BZ4D,MODB4D,U4D,POT4D,dist5d_prof,dist5d_fida, & + RHO4D,XRHO4D,YRHO4D, & BEAM_DENSITY,wall_load,wall_shine, ndot_prof, epower_prof, & ipower_prof,j_prof,dense_prof) END SUBROUTINE diff --git a/BEAMS3D/Sources/beams3d_physics_mod.f90 b/BEAMS3D/Sources/beams3d_physics_mod.f90 index 933e08a9..dca93a61 100644 --- a/BEAMS3D/Sources/beams3d_physics_mod.f90 +++ b/BEAMS3D/Sources/beams3d_physics_mod.f90 @@ -29,7 +29,8 @@ MODULE beams3d_physics_mod ns_prof5, my_end, h1_prof, fact_crit_legacy USE beams3d_grid, ONLY: BR_spl, BZ_spl, delta_t, BPHI_spl, & MODB_spl, MODB4D, & - phimax, S4D, X4D, Y4D, TE4D, NE4D, TI4D, ZEFF4D, & + phimax, TE4D, NE4D, TI4D, ZEFF4D, & + RHO4D, XRHO4D, YRHO4D, & nr, nphi, nz, rmax, rmin, zmax, zmin, & phimin, eps1, eps2, eps3, raxis, phiaxis,& zaxis, U4D,nzeff, dexionT, dexionD, dexionHe3, & @@ -237,6 +238,7 @@ SUBROUTINE beams3d_physics_gc(t, q) zeta, sigma, zeta_mean, zeta_o, v_s, tau_inv, tau_spit_inv, & reduction, dve,dvi, tau_spit, v_crit, coulomb_log, te_cube, & inv_mymass, speed_cube, vcrit_cube, vfrac, modb, s_temp, & + rho_temp, & vc3_tauinv, vbeta, zeff_temp,& !omega_p2, Omega_p, bmax, mu_ip, u_ip2, bmin_c, bmin_q, bmin sm,omega2,vrel2,bmax,bmincl,bminqu,bmin @@ -305,8 +307,8 @@ SUBROUTINE beams3d_physics_gc(t, q) zeff_temp = max(fval(1),one) 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 = max(fval(1),zero) + RHO4D(1,1,1,1),nr,nphi,nz) + rho_temp = max(fval(1),zero) DO l = 1, NION CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,& hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),& @@ -373,7 +375,7 @@ SUBROUTINE beams3d_physics_gc(t, q) q(4) = vll RETURN END IF - l = MAX(MIN(CEILING(SQRT(s_temp)*h1_prof),ns_prof1),1) + l = MAX(MIN(CEILING(rho_temp*h1_prof),ns_prof1),1) epower_prof(mybeam,l) = epower_prof(mybeam,l) + mymass*dve*dt*speed*weight(myline) ipower_prof(mybeam,l) = ipower_prof(mybeam,l) + mymass*dvi*dt*speed*weight(myline) vll = vfrac*vll @@ -450,6 +452,7 @@ SUBROUTINE beams3d_physics_fo(t, q) zeta, sigma, zeta_mean, zeta_o, v_s, tau_inv, tau_spit_inv, & reduction, dve,dvi, tau_spit, v_crit, coulomb_log, te_cube, & inv_mymass, speed_cube, vcrit_cube, vfrac, modb, s_temp, & + rho_temp, & vc3_tauinv, vbeta, zeff_temp, br_temp, bphi_temp, bz_temp, vperp, & sm,omega2,vrel2,bmax,bmincl,bminqu,bmin, binv DOUBLE PRECISION :: Ebench ! for ASCOT Benchmark @@ -512,8 +515,8 @@ SUBROUTINE beams3d_physics_fo(t, q) zeff_temp = max(fval(1),one) 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 = max(fval(1),zero) + RHO4D(1,1,1,1),nr,nphi,nz) + rho_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) @@ -596,7 +599,7 @@ SUBROUTINE beams3d_physics_fo(t, q) q(6) = q(6) + vll*bz_temp RETURN END IF - l = MAX(MIN(CEILING(SQRT(s_temp)*h1_prof),ns_prof1),1) + l = MAX(MIN(CEILING(rho_temp*h1_prof),ns_prof1),1) epower_prof(mybeam,l) = epower_prof(mybeam,l) + mymass*dve*dt*speed*weight(myline) ipower_prof(mybeam,l) = ipower_prof(mybeam,l) + mymass*dvi*dt*speed*weight(myline) vll = vfrac*vll @@ -684,7 +687,8 @@ SUBROUTINE beams3d_follow_neut(t, q) LOGICAL :: ltest INTEGER :: ier, l, m DOUBLE PRECISION :: rinv, phi_temp, dt_local, ti_temp, ne_temp,& - s_temp, x0, y0, z0, xw, yw, zw, te_temp, Zeff_temp + s_temp, x0, y0, z0, xw, yw, zw, te_temp, Zeff_temp, & + rho_temp DOUBLE PRECISION :: qf(3),qs(3),qe(3) DOUBLE PRECISION :: rlocal(num_depo), plocal(num_depo), zlocal(num_depo) DOUBLE PRECISION :: tilocal(num_depo), telocal(num_depo), nelocal(num_depo) @@ -737,8 +741,6 @@ SUBROUTINE beams3d_follow_neut(t, q) t = t + dt_local phi_temp = MODULO(q(2), phimax) IF (phi_temp < 0) phi_temp = phi_temp + phimax - !CALL EZspline_isInDomain(S_spl,q(1),phi_temp,q(3),ier) - !IF (ier==0) THEN IF ((q(1) >= rmin-eps1) .and. (q(1) <= rmax+eps1) .and. & (phi_temp >= phimin-eps2) .and. (phi_temp <= phimax+eps2) .and. & (q(3) >= zmin-eps3) .and. (q(3) <= zmax+eps3)) THEN @@ -748,13 +750,12 @@ SUBROUTINE beams3d_follow_neut(t, q) xparam = (q(1) - raxis(i)) * hri(i) yparam = (phi_temp - phiaxis(j)) * hpi(j) zparam = (q(3) - zaxis(k)) * hzi(k) - s_temp =1.5 - !CALL EZspline_interp(S_spl,q(1),phi_temp,q(3),s_temp,ier) + rho_temp =1.5 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) - IF (s_temp < one) EXIT + RHO4D(1,1,1,1),nr,nphi,nz) + rho_temp = fval(1) + IF (rho_temp < one) EXIT END IF IF ((q(1) > 5*rmax) .or. (q(1) < rmin)) THEN t = my_end+dt_local @@ -810,9 +811,9 @@ SUBROUTINE beams3d_follow_neut(t, q) zparam = (q(3) - zaxis(k)) * hzi(k) 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) - IF (s_temp > one) EXIT INNER + RHO4D(1,1,1,1),nr,nphi,nz) + rho_temp = fval(1) + IF (rho_temp > one) EXIT INNER ELSE EXIT INNER END IF @@ -944,11 +945,10 @@ SUBROUTINE beams3d_follow_neut(t, q) xparam = (rlocal(l) - raxis(i)) * hri(i) yparam = (plocal(l) - phiaxis(j)) * hpi(j) zparam = (zlocal(l) - zaxis(k)) * hzi(k) - !CALL EZspline_interp(S_spl,rlocal(l),plocal(l),zlocal(l),s_temp,ier) 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) + RHO4D(1,1,1,1),nr,nphi,nz) + rho_temp = fval(1) lneut=.false. xlast = qf(1) ylast = qf(2) @@ -1887,6 +1887,7 @@ SUBROUTINE beams3d_SFLX(q,S) ! fval Spline output array !-------------------------------------------------------------- DOUBLE PRECISION :: r_temp, z_temp, phi_temp + DOUBLE PRECISION :: RHO ! For splines INTEGER :: i,j,k REAL*8 :: xparam, yparam, zparam @@ -1919,8 +1920,9 @@ SUBROUTINE beams3d_SFLX(q,S) ! Evaluate the Splines 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 = max(fval(1),zero) + RHO4D(1,1,1,1),nr,nphi,nz) + RHO = max(fval(1),zero) + S = RHO*RHO END IF RETURN @@ -1955,6 +1957,7 @@ SUBROUTINE beams3d_SUFLX(q,S,U) ! fval Spline output array !-------------------------------------------------------------- DOUBLE PRECISION :: r_temp, z_temp, phi_temp + DOUBLE PRECISION :: RHO ! For splines INTEGER :: i,j,k REAL*8 :: xparam, yparam, zparam @@ -1987,8 +1990,9 @@ SUBROUTINE beams3d_SUFLX(q,S,U) ! Evaluate the Splines 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 = max(fval(1),zero) + RHO4D(1,1,1,1),nr,nphi,nz) + RHO = max(fval(1),zero) + S = RHO*RHO CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,& hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),& U4D(1,1,1,1),nr,nphi,nz) @@ -2022,7 +2026,6 @@ SUBROUTINE beams3d_suv2rzp(s,u,v,r_out,z_out,phi_out) DOUBLE PRECISION, INTENT(inout) :: r_out DOUBLE PRECISION, INTENT(inout) :: z_out DOUBLE PRECISION, INTENT(out) :: phi_out - !REAL(rprec), POINTER, DIMENSION(:,:,:,:), INTENT(inout) :: X4D, Y4D !-------------------------------------------------------------- ! Local Variables @@ -2031,12 +2034,15 @@ SUBROUTINE beams3d_suv2rzp(s,u,v,r_out,z_out,phi_out) INTEGER :: n DOUBLE PRECISION :: s0, u0, residual, detJ, delR, delZ, fnorm, & factor, x, y, x0, y0, x_term, y_term, dxdR, dxdZ, dydR, dydZ + DOUBLE PRECISION :: rho, rho0, xrho0, yrho0 + DOUBLE PRECISION :: xrho_term, yrho_term ! For splines INTEGER :: i,j,k, ier REAL*8 :: xparam, yparam, zparam INTEGER, parameter :: ict(8)=(/1,1,1,1,0,0,0,0/) REAL*8 :: fvalx(1,4),fvaly(1,4) !(f,df/fR,df/dphi,dfdZ) + REAL*8 :: fvalxrho(1,4),fvalyrho(1,4) !(f,df/fR,df/dphi,dfdZ) !-------------------------------------------------------------- @@ -2055,9 +2061,12 @@ SUBROUTINE beams3d_suv2rzp(s,u,v,r_out,z_out,phi_out) ! Adjust u u = MOD(u,pi2) + rho = sqrt(s) x0 = s * COS(u) y0 = s * SIN(U) + xrho0 = rho * COS(u) + yrho0 = rho * SIN(u) fnorm = MAX(x0*x0+y0*y0,1E-5) fnorm = 1./fnorm @@ -2073,26 +2082,26 @@ SUBROUTINE beams3d_suv2rzp(s,u,v,r_out,z_out,phi_out) xparam = (r_out - raxis(i)) * hri(i) zparam = (z_out - zaxis(k)) * hzi(k) ! Evaluate the Splines - CALL R8HERM3FCN(ict,1,1,fvalx,i,j,k,xparam,yparam,zparam,& + CALL R8HERM3FCN(ict,1,1,fvalxrho,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) - CALL R8HERM3FCN(ict,1,1,fvaly,i,j,k,xparam,yparam,zparam,& + XRHO4D(1,1,1,1),nr,nphi,nz) + CALL R8HERM3FCN(ict,1,1,fvalyrho,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) + YRHO4D(1,1,1,1),nr,nphi,nz) - x_term = x0 - fvalx(1,1) - y_term = y0 - fvaly(1,1) + xrho_term = xrho0 - fvalxrho(1,1) + yrho_term = yrho0 - fvalyrho(1,1) - detJ = fvalx(1,2) * fvaly(1,4) - fvaly(1,2) * fvalx(1,4) + detJ = fvalxrho(1,2) * fvalyrho(1,4) - fvalyrho(1,2) * fvalxrho(1,4) detJ = MAX(detJ,0.0001) !Upper bound for step size as detJ enters in denominator - delR = -(-x_term*fvaly(1,4) + y_term*fvalx(1,4))/detJ - delZ = -( x_term*fvaly(1,2) - y_term*fvalx(1,2))/detJ + delR = -(-xrho_term*fvalyrho(1,4) + yrho_term*fvalxrho(1,4))/detJ + delZ = -( xrho_term*fvalyrho(1,2) - yrho_term*fvalxrho(1,2))/detJ delR = MIN(MAX(delR,-hr(1)),hr(1)) delZ = MIN(MAX(delZ,-hz(1)),hz(1)) - residual = (x_term*x_term+y_term*y_term)*fnorm - !WRITE(6,*) '----- ',s,u,s0,u0,r_out,z_out,residual,tau,delR,delZ + residual = (xrho_term*xrho_term+yrho_term*yrho_term)*fnorm + !WRITE(6,*) '----- ',s,rho,u,s0,rho0,u0,r_out,z_out,residual,tau,delR,delZ IF (residual < 0.01) THEN !"Damping" of oscillation delR = delR*0.5 diff --git a/BEAMS3D/Sources/beams3d_volume.f90 b/BEAMS3D/Sources/beams3d_volume.f90 index 98bdb4dc..f24b2c14 100644 --- a/BEAMS3D/Sources/beams3d_volume.f90 +++ b/BEAMS3D/Sources/beams3d_volume.f90 @@ -149,4 +149,4 @@ SUBROUTINE beams3d_volume !----------------------------------------------------------------------- ! END SUBROUTINE !----------------------------------------------------------------------- - END SUBROUTINE beams3d_volume \ No newline at end of file + END SUBROUTINE beams3d_volume diff --git a/BEAMS3D/Sources/beams3d_write.f90 b/BEAMS3D/Sources/beams3d_write.f90 index bd25fcaa..104a1bc8 100644 --- a/BEAMS3D/Sources/beams3d_write.f90 +++ b/BEAMS3D/Sources/beams3d_write.f90 @@ -16,6 +16,7 @@ SUBROUTINE beams3d_write(write_type) USE beams3d_lines USE beams3d_grid, ONLY: nr, nphi, nz, B_R, B_PHI, B_Z, raxis, & zaxis, phiaxis, S_ARR, U_ARR, POT_ARR, & + RHO_ARR, & ZEFF_ARR, TE, TI, NE, wall_load, wall_shine, & plasma_mass, plasma_Zmean, & B_kick_min, B_kick_max, freq_kick, & @@ -114,6 +115,8 @@ SUBROUTINE beams3d_write(write_type) IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'B_PHI',ier) CALL write_var_hdf5(fid,'S_ARR',nr,nphi,nz,ier,DBLVAR=S_ARR,ATT='Normalized Toroidal Flux',ATT_NAME='description') IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'S_ARR',ier) + CALL write_var_hdf5(fid,'RHO_ARR',nr,nphi,nz,ier,DBLVAR=RHO_ARR,ATT='Normalized Radius',ATT_NAME='description') + IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'RHO_ARR',ier) CALL write_var_hdf5(fid,'U_ARR',nr,nphi,nz,ier,DBLVAR=U_ARR,ATT='Equilibrium Poloidal Angle [rad]',ATT_NAME='description') IF (ier /= 0) CALL handle_err(HDF5_WRITE_ERR,'U_ARR',ier) CALL write_var_hdf5(fid,'POT_ARR',nr,nphi,nz,ier,DBLVAR=POT_ARR,ATT='Electrostatic Potential [V]',ATT_NAME='description') diff --git a/BEAMS3D/Sources/fidasim_input_mod.f90 b/BEAMS3D/Sources/fidasim_input_mod.f90 index e6f281bb..f80b8d90 100644 --- a/BEAMS3D/Sources/fidasim_input_mod.f90 +++ b/BEAMS3D/Sources/fidasim_input_mod.f90 @@ -20,7 +20,7 @@ MODULE fidasim_input_mod TE, TI, NE, npot, nte, nti, & POT4D, NE4D, TE4D, TI4D, ZEFF4D, & BR4D, BPHI4D, BZ4D, & - hr, hp, hz, hri, hpi, hzi, S4D, U4D, & + hr, hp, hz, hri, hpi, hzi, U4D, & rmin, rmax, phimin, phimax, & rmin_fida, rmax_fida, zmin_fida, zmax_fida, phimin_fida, phimax_fida, & raxis_fida, zaxis_fida, phiaxis_fida, nr_fida, nphi_fida, nz_fida, & diff --git a/BEAMS3D/Sources/out_beams3d_gc.f90 b/BEAMS3D/Sources/out_beams3d_gc.f90 index 662a0a5c..85a8e50a 100644 --- a/BEAMS3D/Sources/out_beams3d_gc.f90 +++ b/BEAMS3D/Sources/out_beams3d_gc.f90 @@ -62,14 +62,10 @@ SUBROUTINE out_beams3d_gc(t, q) neut_lines(mytdex,myline) = lneut x0 = MOD(q(2), phimax) IF (x0 < 0) x0 = x0 + phimax - !CALL EZspline_isInDomain(S_spl,q(1),x0,q(3),ier) rho_help = 0 ! If we're out of domain then don't worry about collisions - !IF (ier==0) THEN IF ((q(1) >= rmin-eps1) .and. (q(1) <= rmax+eps1) .and. & (x0 >= phimin-eps2) .and. (x0 <= phimax+eps2) .and. & (q(3) >= zmin-eps3) .and. (q(3) <= zmax+eps3)) THEN - - i = MIN(MAX(COUNT(raxis < q(1)),1),nr-1) j = MIN(MAX(COUNT(phiaxis < x0),1),nphi-1) k = MIN(MAX(COUNT(zaxis < q(3)),1),nz-1) @@ -78,12 +74,12 @@ SUBROUTINE out_beams3d_gc(t, q) zparam = (q(3) - zaxis(k)) * hzi(k) 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) + XRHO4D(1,1,1,1),nr,nphi,nz) CALL R8HERM3FCN(ict,1,1,fval2,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) - y0 = SQRT(fval(1)*fval(1) + fval2(1) * fval2(1)) - rho_help = sqrt(y0) + YRHO4D(1,1,1,1),nr,nphi,nz) + y0 = fval(1)*fval(1)+fval2(1)*fval2(1) ! s = rho*rho + rho_help = SQRT(y0) z0 = ATAN2(fval2(1),fval(1)) S_lines(mytdex, myline) = y0 U_lines(mytdex, myline) = z0 diff --git a/BEAMS3D/Sources/out_beams3d_part.f90 b/BEAMS3D/Sources/out_beams3d_part.f90 index 36cfaf4d..63fd9d10 100644 --- a/BEAMS3D/Sources/out_beams3d_part.f90 +++ b/BEAMS3D/Sources/out_beams3d_part.f90 @@ -45,7 +45,8 @@ SUBROUTINE out_beams3d_part(t, q) INTEGER :: ier, d1, d2, d3, d4, d5 DOUBLE PRECISION :: x0,y0,z0,x1,y1,z1,xw,yw,zw,vperp, & br_temp, bphi_temp, bz_temp, & - v_total, binv, vll_temp + v_total, binv, vll_temp, & + rho0 DOUBLE PRECISION :: q2(6),qdot(6), q4(4) ! For splines INTEGER :: i,j,k,l @@ -65,7 +66,7 @@ SUBROUTINE out_beams3d_part(t, q) neut_lines(mytdex,myline) = lneut x0 = MOD(q(2), phimax) IF (x0 < 0) x0 = x0 + phimax - y0 = 0 + rho_help = 0 ! If we're out of domain then don't worry about collisions IF ((q(1) >= rmin-eps1) .and. (q(1) <= rmax+eps1) .and. & (x0 >= phimin-eps2) .and. (x0 <= phimax+eps2) .and. & @@ -78,18 +79,16 @@ SUBROUTINE out_beams3d_part(t, q) zparam = (q(3) - zaxis(k)) * hzi(k) 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) + XRHO4D(1,1,1,1),nr,nphi,nz) CALL R8HERM3FCN(ict,1,1,fval2,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) - - y0 = SQRT(fval(1)*fval(1) + fval2(1) * fval2(1)) - !z0 = fval(1) + YRHO4D(1,1,1,1),nr,nphi,nz) + y0 = fval(1)*fval(1)+fval2(1)*fval2(1) ! s = rho*rho + rho_help = SQRT(y0) z0 = ATAN2(fval2(1),fval(1)) S_lines(mytdex, myline) = y0 U_lines(mytdex, myline) = z0 - ! Now We get Br, Bphi, Bz to calc B (and vll) CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,& hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),& @@ -122,7 +121,7 @@ SUBROUTINE out_beams3d_part(t, q) ! Calc dist func bins x0 = MOD(q(2),pi2) IF (x0 < 0) x0 = x0 + pi2 - d1 = MAX(MIN(CEILING(SQRT(y0)*h1_prof ), ns_prof1), 1) ! Rho Bin + d1 = MAX(MIN(CEILING(rho_help*h1_prof ), ns_prof1), 1) ! Rho Bin d2 = MAX(MIN(CEILING( z0*h2_prof ), ns_prof2), 1) ! U Bin d3 = MAX(MIN(CEILING( x0*h3_prof ), ns_prof3), 1) ! V Bin d4 = MAX(MIN(1+nsh_prof4+FLOOR(h4_prof*vll_temp), ns_prof4), 1) ! vll @@ -131,17 +130,15 @@ SUBROUTINE out_beams3d_part(t, q) !CALL MPI_WIN_LOCK(MPI_LOCK_EXCLUSIVE,myworkid,0,win_dist5d,ier) dist5d_prof(mybeam,d1,d2,d3,d4,d5) = dist5d_prof(mybeam,d1,d2,d3,d4,d5) + xw !CALL MPI_WIN_UNLOCK(myworkid,win_dist5d,ier) - IF (lfidasim_cyl) THEN - !x0 = MOD(q(2), phimax_fida) - !IF (x0 < 0) x0 = x0 + phimax_fida - d1 = MIN(MAX(CEILING((q(1)-rmin_fida)*r_h),1),nr_fida) - d2 = MIN(MAX(CEILING((x0-phimin_fida)*p_h),1),nphi_fida) - d3 = MIN(MAX(CEILING((q(3)-zmin_fida)*z_h),1),nz_fida) - y0 = (q(4)**2+vperp**2) - d4 = MIN(MAX(CEILING((y0*E_by_v-emin_fida)*e_h),1),nenergy_fida) - d5 = MIN(MAX(CEILING((q(4)/SQRT(y0)-pimin_fida)*pi_h),1),npitch_fida) - dist5d_fida(d1,d3,d2,d4,d5) = dist5d_fida(d1,d3,d2,d4,d5) + xw - END IF + IF (lfidasim_cyl) THEN + d1 = MIN(MAX(CEILING((q(1)-rmin_fida)*r_h),1),nr_fida) + d2 = MIN(MAX(CEILING((x0-phimin_fida)*p_h),1),nphi_fida) + d3 = MIN(MAX(CEILING((q(3)-zmin_fida)*z_h),1),nz_fida) + y0 = (q(4)**2+vperp**2) + d4 = MIN(MAX(CEILING((y0*E_by_v-emin_fida)*e_h),1),nenergy_fida) + d5 = MIN(MAX(CEILING((q(4)/SQRT(y0)-pimin_fida)*pi_h),1),npitch_fida) + dist5d_fida(d1,d3,d2,d4,d5) = dist5d_fida(d1,d3,d2,d4,d5) + xw + END IF IF (lcollision) CALL beams3d_physics_fo(t,q) IF (ltherm) THEN ndot_prof(mybeam,d1) = ndot_prof(mybeam,d1) + weight(myline) @@ -153,7 +150,7 @@ SUBROUTINE out_beams3d_part(t, q) v_total = SUM(q(4:6)*q(4:6)) vll_temp = SQRT(v_total) ! Vll=Vtotal END IF - IF (lvessel .and. mytdex > 0 .and. y0 > 0.5) THEN + IF (lvessel .and. mytdex > 0 .and. rho_help > 0.5) THEN lhit = .false. x0 = xlast y0 = ylast