From e4a53561a67591ae3821051be0c2e95afa724cf5 Mon Sep 17 00:00:00 2001 From: mbeidler3 Date: Thu, 3 Aug 2023 15:47:14 -0400 Subject: [PATCH] creating ACC routines for aorsa fields --- src/korc_fields.f90 | 1 + src/korc_initialize.f90 | 2 + src/korc_input.f90 | 11 +- src/korc_interp.f90 | 101 +++++++++- src/korc_ppusher.f90 | 303 ++++++++++++++++++++++++------ src/korc_spatial_distribution.f90 | 15 +- src/korc_types.f90 | 4 +- src/main.f90 | 11 +- 8 files changed, 370 insertions(+), 78 deletions(-) diff --git a/src/korc_fields.f90 b/src/korc_fields.f90 index 62b00eff..12e02919 100755 --- a/src/korc_fields.f90 +++ b/src/korc_fields.f90 @@ -1291,6 +1291,7 @@ subroutine initialize_fields(params,F) F%AORSA_AMP_Scale=AORSA_AMP_Scale F%AORSA_freq=AORSA_freq F%AORSA_nmode=AORSA_nmode + F%AORSA_mmode=AORSA_mmode F%useLCFS = useLCFS F%useDiMES = useDiMES F%DiMESloc = DiMESloc diff --git a/src/korc_initialize.f90 b/src/korc_initialize.f90 index 9b3f7a76..e831415c 100755 --- a/src/korc_initialize.f90 +++ b/src/korc_initialize.f90 @@ -509,6 +509,7 @@ subroutine initialize_particles(params,F,P,spp) spp(ii)%sigmaZ = sigmaZ(ii) spp(ii)%theta_gauss = theta_gauss(ii) spp(ii)%psi_max = psi_max(ii) + spp(ii)%psi_min = psi_min(ii) spp(ii)%Spong_w = Spong_w(ii) spp(ii)%Spong_b = Spong_b(ii) spp(ii)%Spong_dlam = Spong_dlam(ii) @@ -659,6 +660,7 @@ subroutine initialize_particles(params,F,P,spp) DEALLOCATE(sigmaZ) DEALLOCATE(theta_gauss) DEALLOCATE(psi_max) + DEALLOCATE(psi_min) DEALLOCATE(Spong_b) DEALLOCATE(Spong_w) DEALLOCATE(Spong_dlam) diff --git a/src/korc_input.f90 b/src/korc_input.f90 index 1ececdda..aae6093a 100755 --- a/src/korc_input.f90 +++ b/src/korc_input.f90 @@ -114,7 +114,7 @@ module korc_input REAL(rp), DIMENSION(:), ALLOCATABLE :: sigmaR REAL(rp), DIMENSION(:), ALLOCATABLE :: sigmaZ REAL(rp), DIMENSION(:), ALLOCATABLE :: theta_gauss - REAL(rp), DIMENSION(:), ALLOCATABLE :: psi_max + REAL(rp), DIMENSION(:), ALLOCATABLE :: psi_max,psi_min !! Maximum value of the argument of the 2D gaussian exponential, used for an !! indicator function that limits the region of MH sampling ! goes as R^2 for HOLLMANN-3D, is psiN_max for HOLLMANN-3D-PSI @@ -211,6 +211,7 @@ module korc_input REAL(rp) :: AORSA_AMP_Scale=1.0 REAL(rp) :: AORSA_freq=0.0 REAL(rp) :: AORSA_nmode=0.0 + REAL(rp) :: AORSA_mmode=0.0 LOGICAL :: Analytic_D3D_IWL=.FALSE. INTEGER :: ntiles=42 REAL(rp) :: circumradius=1.016 @@ -450,7 +451,7 @@ subroutine read_namelist(params,infile,echo_in,outdir) NAMELIST /plasma_species/ ppp,q,m,Eno,etao,Eo_lims,etao_lims,runaway, & spatial_distribution,energy_distribution,pitch_distribution,Ro, & PHIo,Zo,r_inner,r_outter,falloff_rate,shear_factor,sigmaR,sigmaZ, & - theta_gauss,psi_max,Xtrace,Spong_b,Spong_w,Spong_dlam,dth,dR,dZ,dgam,& + theta_gauss,psi_max,psi_min,Xtrace,Spong_b,Spong_w,Spong_dlam,dth,dR,dZ,dgam,& pinit NAMELIST /analytical_fields_params/ Bo,minor_radius,major_radius,& qa,qo,Eo,current_direction,nR,nZ,nPHI,dim_1D,dt_E_SC,Ip_exp, & @@ -459,7 +460,7 @@ subroutine read_namelist(params,infile,echo_in,outdir) axisymmetric_fields, Eo,E_dyn,E_pulse,E_width,res_double, & dim_1D,dt_E_SC,Ip_exp,PSIp_lim,Dim2x1t,t0_2x1t,E_2x1t,ReInterp_2x1t, & ind0_2x1t,PSIp_0,B1field,psip_conv,MARS_AMP_Scale,Analytic_D3D_IWL, & - ntiles,circumradius,AORSA_AMP_Scale,AORSA_freq,AORSA_nmode,E1field, & + ntiles,circumradius,AORSA_AMP_Scale,AORSA_freq,AORSA_nmode,AORSA_mmode,E1field, & useLCFS,useDiMES,DiMESloc,DiMESdims,MARS_phase NAMELIST /plasmaProfiles/ radius_profile,ne_profile,neo,n_ne,a_ne, & Te_profile,Teo,n_Te,a_Te,n_REr0,n_tauion,n_lamfront,n_lamback, & @@ -552,6 +553,7 @@ subroutine read_namelist(params,infile,echo_in,outdir) ALLOCATE(sigmaZ(num_species)) ALLOCATE(theta_gauss(num_species)) ALLOCATE(psi_max(num_species)) + ALLOCATE(psi_min(num_species)) ALLOCATE(falloff_rate(num_species)) ALLOCATE(energy_distribution(num_species)) ALLOCATE(pitch_distribution(num_species)) @@ -584,7 +586,8 @@ subroutine read_namelist(params,infile,echo_in,outdir) sigmaR = 1.e6 sigmaZ = 0.2 theta_gauss = 0.0 - psi_max=.8446 + psi_max=1. + psi_min=0. falloff_rate = 0.0 energy_distribution = 'MONOENERGETIC' pitch_distribution = 'MONOPITCH' diff --git a/src/korc_interp.f90 b/src/korc_interp.f90 index 317c30bb..61d4a3c5 100755 --- a/src/korc_interp.f90 +++ b/src/korc_interp.f90 @@ -87,7 +87,7 @@ module korc_interp !! ends of the \(Z\) direction. END TYPE KORC_2X1T_FIELDS_INTERPOLANT -TYPE, PRIVATE :: KORC_2DX_FIELDS_INTERPOLANT +TYPE :: KORC_2DX_FIELDS_INTERPOLANT !! @note Derived type containing 2-D PSPLINE interpolants for !! cylindrical components of vector fields \(\mathbf{F}(R,Z) = !! F_R\hat{e}_R + F_\phi\hat{e}_phi+ F_Z\hat{e}_Z\). @@ -292,7 +292,7 @@ module korc_interp TYPE(KORC_2D_FIELDS_INTERPOLANT) :: b1Imfield_2d TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: b1Refield_2dx TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: b1Imfield_2dx -TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: e1Refield_2dx +TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: e1Refield_2dx TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: e1Imfield_2dx TYPE(KORC_2X1T_FIELDS_INTERPOLANT) :: bfield_2X1T TYPE(KORC_2D_FIELDS_INTERPOLANT) :: efield_2d @@ -1611,7 +1611,7 @@ subroutine check_if_in_fields_domain_2D_p_ACC(fields_domain_local,bfield_2d_loca !! to the vertical position of the particles. REAL(rp) :: Rwall,lscale,xtmp,ytmp,ztmp - REAL(rp) :: DiMESsurf,DiMESrad + REAL(rp) :: DiMESsurf,DiMESrad,DiMESloc_deg REAL(rp),DIMENSION(3) :: DiMESloc_cart REAL(rp),DIMENSION(3),INTENT(IN) :: DiMESloc_cyl REAL(rp),DIMENSION(2),INTENT(IN) :: DiMESdims @@ -1654,18 +1654,18 @@ subroutine check_if_in_fields_domain_2D_p_ACC(fields_domain_local,bfield_2d_loca if (useDiMES) THEN - DiMESloc_cyl(2)=C_PI*DiMESloc_cyl(2)/180._rp + DiMESloc_deg=C_PI*DiMESloc_cyl(2)/180._rp if ((abs(Y_R-DiMESloc_cyl(1)).le.DiMESdims(1)).and. & (abs(Y_Z-DiMESloc_cyl(3)).le.DiMESdims(2)).and.& - (abs(Y_PHI-DiMESloc_cyl(2)).le.asin(DiMESdims(1)/DiMESloc_cyl(1)))) THEN + (abs(Y_PHI-DiMESloc_deg).le.asin(DiMESdims(1)/DiMESloc_cyl(1)))) THEN xtmp=Y_R*cos(Y_PHI) ytmp=Y_R*sin(Y_PHI) ztmp=Y_Z - DiMESloc_cart(1)=DiMESloc_cyl(1)*cos(DiMESloc_cyl(2)) - DiMESloc_cart(2)=DiMESloc_cyl(1)*sin(DiMESloc_cyl(2)) + DiMESloc_cart(1)=DiMESloc_cyl(1)*cos(DiMESloc_deg) + DiMESloc_cart(2)=DiMESloc_cyl(1)*sin(DiMESloc_deg) DiMESloc_cart(3)=DiMESloc_cyl(3) DiMESrad=DiMESdims(1)**2-(xtmp-DiMESloc_cart(1))**2-(ytmp-DiMESloc_cart(2))**2 @@ -2426,7 +2426,7 @@ subroutine interp_FOfields_mars(prtcls, F, params) end subroutine interp_FOfields_mars -subroutine provide_ezspline_mars_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local,fields_domain_local +subroutine provide_ezspline_mars_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local,fields_domain_local) TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT) :: bfield_2d_local TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT) :: b1Refield_2d_local TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT) :: b1Imfield_2d_local @@ -2551,6 +2551,24 @@ subroutine interp_FOfields_mars_p_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfie end subroutine interp_FOfields_mars_p_ACC +subroutine provide_ezspline_aorsa_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & + e1Refield_2d_local,e1Imfield_2d_local,fields_domain_local) + TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(OUT) :: bfield_2d_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT),INTENT(OUT) :: b1Refield_2d_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT),INTENT(OUT) :: b1Imfield_2d_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT),INTENT(OUT) :: e1Refield_2d_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT),INTENT(OUT) :: e1Imfield_2d_local + TYPE(KORC_INTERPOLANT_DOMAIN),INTENT(OUT) :: fields_domain_local + + bfield_2d_local=bfield_2d + b1Refield_2d_local=b1Refield_2dx + b1Imfield_2d_local=b1Imfield_2dx + e1Refield_2d_local=e1Refield_2dx + e1Imfield_2d_local=e1Imfield_2dx + fields_domain_local=fields_domain + +end subroutine provide_ezspline_aorsa_ACC + subroutine interp_FOfields_aorsa(prtcls, F, params) TYPE(KORC_PARAMS), INTENT(IN) :: params TYPE(PARTICLES), INTENT(INOUT) :: prtcls @@ -2651,6 +2669,73 @@ subroutine interp_FOfields_aorsa(prtcls, F, params) end subroutine interp_FOfields_aorsa +subroutine interp_FOfields_aorsa_p_ACC(time,bfield_2d_local,b1Refield_2dx_local,b1Imfield_2dx_local, & + e1Refield_2dx_local,e1Imfield_2dx_local,psip_conv,amp,nmode,omega,Bo,Ro, & + Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp) + !$acc routine seq + REAL(rp), INTENT(IN) :: time + REAL(rp),INTENT(IN) :: Y_R,Y_PHI,Y_Z + REAL(rp),INTENT(OUT) :: B_X,B_Y,B_Z + REAL(rp),INTENT(OUT) :: E_X,E_Y,E_Z + REAL(rp) :: B0_R,B0_PHI,B0_Z,B0_X,B0_Y + REAL(rp) :: B1_X,B1_Y,B1_Z + REAL(rp) :: B1Re_X,B1Re_Y,B1Re_Z + REAL(rp) :: B1Im_X,B1Im_Y,B1Im_Z + REAL(rp) :: E1Re_X,E1Re_Y,E1Re_Z + REAL(rp) :: E1Im_X,E1Im_Y,E1Im_Z + REAL(rp),INTENT(OUT) :: PSIp + REAL(rp) :: cP,sP,cnP,snP + REAL(rp), DIMENSION(3) :: A + INTEGER :: ezerr_local + REAL(rp),INTENT(IN) :: psip_conv,amp,nmode,omega,Bo,Ro + TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN) :: bfield_2d_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT),INTENT(IN) :: b1Refield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT),INTENT(IN) :: b1Imfield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT),INTENT(IN) :: e1Refield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT),INTENT(IN) :: e1Imfield_2dx_local + + !$acc routine (EZspline_interp2_FOaorsa) seq + !$acc routine (EZspline_error) seq + + call EZspline_interp2_FOaorsa(bfield_2d_local%A,b1Refield_2dx_local%X,b1Refield_2dx_local%Y, & + b1Refield_2dx_local%Z,b1Imfield_2dx_local%X,b1Imfield_2dx_local%Y,b1Imfield_2dx_local%Z, & + e1Refield_2dx_local%X,e1Refield_2dx_local%Y,e1Refield_2dx_local%Z, & + e1Imfield_2dx_local%X,e1Imfield_2dx_local%Y,e1Imfield_2dx_local%Z, & + Y_R,Y_Z,A,B1Re_X,B1Re_Y,B1Re_Z,B1Im_X,B1Im_Y,B1Im_Z, & + E1Re_X,E1Re_Y,E1Re_Z,E1Im_X,E1Im_Y,E1Im_Z,ezerr) + call EZspline_error(ezerr_local) + + + PSIp=A(1) + + B0_R = psip_conv*A(3)/Y_R + B0_PHI = -Bo*Ro/Y_R + B0_Z = -psip_conv*A(2)/Y_R + + cP=cos(Y_PHI) + sP=sin(Y_PHI) + + B0_X = B0_R*cP - B0_PHI*sP + B0_Y = B0_R*sP + B0_PHI*cP + + cnP=cos(omega*time+nmode*Y_PHI) + snP=sin(omega*time+nmode*Y_PHI) + + B1_X = amp*(B1Re_X*cnP-B1Im_X*snP) + B1_Y = amp*(B1Re_Y*cnP-B1Im_Y*snP) + B1_Z = amp*(B1Re_Z*cnP-B1Im_Z*snP) + + B_X = B0_X+B1_X + B_Y = B0_Y+B1_Y + B_Z = B0_Z+B1_Z + + E_X = amp*(E1Re_X*cnP-E1Im_X*snP) + E_Y = amp*(E1Re_Y*cnP-E1Im_Y*snP) + E_Z = amp*(E1Re_Z*cnP-E1Im_Z*snP) + + +end subroutine interp_FOfields_aorsa_p_ACC + subroutine interp_FOfields_aorsa_p(time,params,pchunk,F,Y_R,Y_PHI,Y_Z, & B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp,flag_cache) TYPE(KORC_PARAMS), INTENT(IN) :: params diff --git a/src/korc_ppusher.f90 b/src/korc_ppusher.f90 index 75841034..e893e8f2 100755 --- a/src/korc_ppusher.f90 +++ b/src/korc_ppusher.f90 @@ -473,7 +473,7 @@ subroutine FO_init_ACC(params,F,spp,output,step) REAL(rp) :: B_X,B_Y,B_Z REAL(rp) :: E_X,E_Y,E_Z REAL(rp) :: PSIp - REAL(rp) :: m_cache,q_cache,psip_conv,amp,phase,Ro,Bo,circumradius,ntiles + REAL(rp) :: m_cache,q_cache,psip_conv,amp,phase,Ro,Bo,circumradius,ntiles,nmode,mmode,omega INTEGER(is) :: flagCon,flagCol LOGICAL :: Analytic_D3D_IWL,useDiMES,Dim2x1t REAL(rp),DIMENSION(2) :: DiMESdims @@ -481,6 +481,10 @@ subroutine FO_init_ACC(params,F,spp,output,step) TYPE(KORC_2D_FIELDS_INTERPOLANT) :: bfield_2d_local TYPE(KORC_2D_FIELDS_INTERPOLANT) :: b1Refield_2d_local TYPE(KORC_2D_FIELDS_INTERPOLANT) :: b1Imfield_2d_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: b1Refield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: b1Imfield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: e1Refield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: e1Imfield_2dx_local TYPE(KORC_INTERPOLANT_DOMAIN) :: fields_domain_local !$acc routine (cart_to_cyl_p_ACC) seq @@ -495,6 +499,9 @@ subroutine FO_init_ACC(params,F,spp,output,step) psip_conv=F%psip_conv amp=F%AMP phase=F%MARS_phase + nmode=F%AORSA_nmode + omega=2*C_PI*F%AORSA_freq + mmode=F%AORSA_mmode Ro=F%Ro Bo=F%Bo @@ -508,26 +515,19 @@ subroutine FO_init_ACC(params,F,spp,output,step) if(output) then - call provide_ezspline_mars_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & - fields_domain_local) - - !$acc enter data copyin(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local,fields_domain_local) - - !$acc parallel loop !& - !!$acc& firstprivate(E0,m_cache,q_cache,psip_conv,amp,phase,Ro,Bo, & - !!$acc& Analytic_D3D_IWL,circumradius,ntiles,useDiMES,DiMESloc_cyl(3), & - !!$acc& DiMESdims(2),Dim2x1t) & - !!$acc& copyin(ii,spp(ii)%ppp,bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & - !!$acc& fields_domain_local) & - !!$acc& copy(spp(ii)%vars%X(1:spp(ii)%ppp,1:3),spp(ii)%vars%V(1:spp(ii)%ppp,1:3), & - !!$acc& spp(ii)%vars%flagCon(1:spp(ii)%ppp),spp(ii)%vars%flagCol(1:spp(ii)%ppp)) & - !!$acc& copyout(spp(ii)%vars%B(1:spp(ii)%ppp,1:3),spp(ii)%vars%E(1:spp(ii)%ppp,1:3), & - !!$acc& spp(ii)%vars%PSI_P(1:spp(ii)%ppp),spp(ii)%vars%eta(1:spp(ii)%ppp), & - !!$acc& spp(ii)%vars%mu(1:spp(ii)%pppb),spp(ii)%vars%Pin(1:spp(ii)%ppp),spp(ii)%vars%Prad(1:spp(ii)%ppp)) & - !!$acc& PRIVATE(pp,X_X,X_Y,X_Z,B_X,B_Y,B_Z,V_X,V_Y,V_Z, & - !!$acc& E_X,E_Y,E_Z,Y_R,Y_PHI,Y_Z,flagCon,flagCol,PSIp,Bmag, & - !!$acc& b_unit_X,b_unit_Y,b_unit_Z,v,vpar,vperp,tmp, & - !!$acc& cross_X,cross_Y,cross_Z,vec_X,vec_Y,vec_Z,g) + if (params%field_model(10:13).eq.'MARS') then + call provide_ezspline_mars_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & + fields_domain_local) + elseif (params%field_model(10:13).eq.'AORSA') then + call provide_ezspline_aorsa_ACC(bfield_2d_local,b1Refield_2dx_local,b1Imfield_2dx_local, & + e1Refield_2dx_local,e1Imfield_2dx_local,fields_domain_local) + endif + + !$acc enter data copyin(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & + !$acc& b1Refield_2dx_local,b1Imfield_2dx_local,e1Refield_2dx_local,e1Imfield_2dx_local, & + !$acc& fields_domain_local) + + !$acc parallel loop do pp=1_idef,spp(ii)%ppp X_X=spp(ii)%vars%X(pp,1) @@ -563,8 +563,14 @@ subroutine FO_init_ACC(params,F,spp,output,step) Dim2x1t,Analytic_D3D_IWL,circumradius, & ntiles,useDiMES,DiMESloc_cyl,DiMESdims,Y_R,Y_PHI,Y_Z,flagCon) - call interp_FOfields_mars_p_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & - psip_conv,amp,phase,Bo,Ro,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,PSIp) + if (params%field_model(10:13).eq.'MARS') then + call interp_FOfields_mars_p_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & + psip_conv,amp,phase,Bo,Ro,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,PSIp) + elseif (params%field_model(10:13).eq.'AORSA') then + call interp_FOfields_aorsa_p_ACC(0._rp,bfield_2d_local,b1Refield_2dx_local,b1Imfield_2dx_local, & + e1Refield_2dx_local,e1Imfield_2dx_local,psip_conv,amp,nmode,omega,Bo,Ro, & + Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp) + endif #endif PSPLINE @@ -646,18 +652,16 @@ subroutine FO_init_ACC(params,F,spp,output,step) enddo !$acc end parallel loop - !$acc exit data delete(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local,fields_domain_local) + !$acc exit data delete(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & + !$acc& b1Refield_2dx_local,b1Imfield_2dx_local,e1Refield_2dx_local,e1Imfield_2dx_local, & + !$acc& fields_domain_local) endif !(if output) if(step.and.(.not.params%FokPlan)) then dt=0.5_rp*params%dt - !$acc parallel loop !& - !!$acc& firstprivate(dt) & - !!$acc& private(pp) & - !!$acc& copyin(ii,spp(ii)%ppp) & - !!$acc& copyout(spp(ii)%vars%X(1:pp,1:3)) + !$acc parallel loop do pp=1_idef,spp(ii)%ppp spp(ii)%vars%X(pp,1) = spp(ii)%vars%X(pp,1) + & @@ -1923,8 +1927,8 @@ subroutine adv_FOinterp_mars_top_ACC(params,F,P,spp) REAL(rp) :: B_X,B_Y,B_Z REAL(rp) :: E_X,E_Y,E_Z REAL(rp) :: PSIp - INTEGER(is) :: flagCon,flagCol - REAL(rp) :: a,m_cache,q_cache,tskip,psip_conv,amp,phase + INTEGER(is) :: flagCon,flagCol,tskip + REAL(rp) :: a,m_cache,q_cache,psip_conv,amp,phase REAL(rp) :: Ro,Bo,circumradius,ntiles,dt INTEGER :: ii,pp,ss,tt,ppp LOGICAL :: Analytic_D3D_IWL,useDiMES,Dim2x1t @@ -1968,28 +1972,7 @@ subroutine adv_FOinterp_mars_top_ACC(params,F,P,spp) call provide_ezspline_mars_ACC(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local, & fields_domain_local) - !$acc parallel loop !& - !!$acc firstprivate(E0,m_cache,q_cache,psip_conv,amp,phase,Ro,Bo, & - !!$acc& Analytic_D3D_IWL,circumradius,ntiles,useDiMES,DiMESloc_cyl, & - !!$acc& DiMESdims,tskip,ii,ppp,dt,Dim2x1t) & - !!$acc& PRIVATE(pp,tt,Bmag,X_X,X_Y,X_Z,V_X,V_Y,V_Z,B_X,B_Y,B_Z, & - !!$acc& E_X,E_Y,E_Z,b_unit_X,b_unit_Y,b_unit_Z,v,vpar,vperp,tmp, & - !!$acc& cross_X,cross_Y,cross_Z,vec_X,vec_Y,vec_Z,g, & - !!$acc& Y_R,Y_PHI,Y_Z,flagCon,flagCol,PSIp) & - !!$acc& copyin(bfield_2d_local,b1Refield_2d_local,b1Imfield_2d_local,fields_domain_local) & - !!$acc& copy(spp(ii)%vars%X(1:spp(ii)%ppp,1:3), & - !!$acc& spp(ii)%vars%V(1:spp(ii)%ppp,1:3), & - !!$acc& spp(ii)%vars%B(1:spp(ii)%ppp,1:3), & - !!$acc& spp(ii)%vars%E(1:spp(ii)%ppp,1:3), & - !!$acc& spp(ii)%vars%PSI_P(1:spp(ii)%ppp), & - !!$acc& spp(ii)%vars%g(1:spp(ii)%ppp), & - !!$acc& spp(ii)%vars%flagCon(1:spp(ii)%ppp), & - !!$acc& spp(ii)%vars%flagCol(1:spp(ii)%ppp)) & - !!$acc& copyout(spp(ii)%vars%eta(1:spp(ii)%ppp), & - !!$acc& spp(ii)%vars%mu(1:spp(ii)%ppp), & - !!$acc& spp(ii)%vars%Prad(1:spp(ii)%ppp), & - !!$acc& spp(ii)%vars%Pin(1:spp(ii)%ppp)) - + !$acc parallel loop do pp=1_idef,ppp X_X=spp(ii)%vars%X(pp,1) @@ -2333,6 +2316,220 @@ subroutine adv_FOinterp_aorsa_top(params,F,P,spp) end subroutine adv_FOinterp_aorsa_top +subroutine adv_FOinterp_aorsa_top_ACC(params,F,P,spp) + TYPE(KORC_PARAMS), INTENT(INOUT) :: params + !! Core KORC simulation parameters. + TYPE(FIELDS), INTENT(IN) :: F + !! An instance of the KORC derived type FIELDS. + TYPE(PROFILES), INTENT(IN) :: P + !! An instance of the KORC derived type PROFILES. + TYPE(SPECIES), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: spp + !! An instance of the derived type SPECIES containing all the parameters + !! and simulation variables of the different species in the simulation. + REAL(rp) :: Bmag + REAL(rp) :: b_unit_X,b_unit_Y,b_unit_Z + REAL(rp) :: v,vpar,vperp + REAL(rp) :: tmp + REAL(rp) :: g + REAL(rp) :: cross_X,cross_Y,cross_Z + REAL(rp) :: vec_X,vec_Y,vec_Z + REAL(rp) :: X_X,X_Y,X_Z + REAL(rp) :: Y_R,Y_PHI,Y_Z + REAL(rp) :: V_X,V_Y,V_Z + REAL(rp) :: B_X,B_Y,B_Z + REAL(rp) :: E_X,E_Y,E_Z + REAL(rp) :: PSIp + INTEGER(is) :: flagCon,flagCol,tskip + REAL(rp) :: a,m_cache,q_cache,psip_conv,amp,phase,nmode,mmode,omega + REAL(rp) :: Ro,Bo,circumradius,ntiles,dt,time,t0,tnorm + INTEGER :: ii,pp,ss,tt,ppp + LOGICAL :: Analytic_D3D_IWL,useDiMES,Dim2x1t + REAL(rp),DIMENSION(2) :: DiMESdims + REAL(rp),DIMENSION(3) :: DiMESloc_cyl + TYPE(KORC_2D_FIELDS_INTERPOLANT) :: bfield_2d_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: b1Refield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: b1Imfield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: e1Refield_2dx_local + TYPE(KORC_2DX_FIELDS_INTERPOLANT) :: e1Imfield_2dx_local + TYPE(KORC_INTERPOLANT_DOMAIN) :: fields_domain_local + + !$acc routine (cart_to_cyl_p_ACC) seq + !$acc routine (check_if_in_fields_domain_2D_p_ACC) seq + !$acc routine (interp_FOfields_aorsa_p_ACC) seq + !$acc routine (advance_FOinterp_vars_ACC) seq + + ss=params%num_species + + do ii = 1_idef,ss + + m_cache=spp(ii)%m + q_cache=spp(ii)%q + a =params%dt + tskip=params%t_skip + ppp=spp(ii)%ppp + dt=params%dt + t0=params%time + tnorm=params%cpp%time + + + psip_conv=F%psip_conv + amp=F%AMP + phase=F%MARS_phase + nmode=F%AORSA_nmode + omega=2*C_PI*F%AORSA_freq + mmode=F%AORSA_mmode + Ro=F%Ro + Bo=F%Bo + + Dim2x1t=F%Dim2x1t + Analytic_D3D_IWL=F%Analytic_D3D_IWL + circumradius=F%circumradius + ntiles=F%ntiles + useDiMES=F%useDiMES + DiMESloc_cyl=F%DiMESloc + DiMESdims=F%DiMESdims + + call provide_ezspline_aorsa_ACC(bfield_2d_local,b1Refield_2dx_local,b1Imfield_2dx_local, & + e1Refield_2dx_local,e1Imfield_2dx_local,fields_domain_local) + + !$acc parallel loop + do pp=1_idef,ppp + + X_X=spp(ii)%vars%X(pp,1) + X_Y=spp(ii)%vars%X(pp,2) + X_Z=spp(ii)%vars%X(pp,3) + + Y_R=0._rp + Y_PHI=0._rp + Y_Z=0._rp + + V_X=spp(ii)%vars%V(pp,1) + V_Y=spp(ii)%vars%V(pp,2) + V_Z=spp(ii)%vars%V(pp,3) + + B_X=spp(ii)%vars%B(pp,1) + B_Y=spp(ii)%vars%B(pp,2) + B_Z=spp(ii)%vars%B(pp,3) + + E_X=spp(ii)%vars%E(pp,1) + E_Y=spp(ii)%vars%E(pp,2) + E_Z=spp(ii)%vars%E(pp,3) + + PSIp=spp(ii)%vars%PSI_P(pp) + + g=spp(ii)%vars%g(pp) + + flagCon=spp(ii)%vars%flagCon(pp) + flagCol=spp(ii)%vars%flagCol(pp) + + do tt=1_ip,tskip + + time=(t0+dt*tt)*tnorm + + call cart_to_cyl_p_ACC(X_X,X_Y,X_Z,Y_R,Y_PHI,Y_Z) + + call check_if_in_fields_domain_2D_p_ACC(fields_domain_local,bfield_2d_local, & + Dim2x1t,Analytic_D3D_IWL,circumradius, & + ntiles,useDiMES,DiMESloc_cyl,DiMESdims,Y_R,Y_PHI,Y_Z,flagCon) + + call interp_FOfields_aorsa_p_ACC(time,bfield_2d_local,b1Refield_2dx_local,b1Imfield_2dx_local, & + e1Refield_2dx_local,e1Imfield_2dx_local,psip_conv,amp,nmode,omega,Bo,Ro, & + Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z,PSIp) + + call advance_FOinterp_vars_ACC(dt,tt,a,q_cache,m_cache, & + X_X,X_Y,X_Z,V_X,V_Y,V_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z, & + g,flagCon,flagCol,PSIp) + end do !timestep iterator + + + spp(ii)%vars%X(pp,1)=X_X + spp(ii)%vars%X(pp,2)=X_Y + spp(ii)%vars%X(pp,3)=X_Z + + spp(ii)%vars%V(pp,1)=V_X + spp(ii)%vars%V(pp,2)=V_Y + spp(ii)%vars%V(pp,3)=V_Z + + spp(ii)%vars%g(pp) = g + spp(ii)%vars%flagCon(pp) = flagCon + spp(ii)%vars%flagCol(pp) = flagCol + + spp(ii)%vars%B(pp,1) = B_X + spp(ii)%vars%B(pp,2) = B_Y + spp(ii)%vars%B(pp,3) = B_Z + + spp(ii)%vars%E(pp,1) = E_X + spp(ii)%vars%E(pp,2) = E_Y + spp(ii)%vars%E(pp,3) = E_Z + + spp(ii)%vars%PSI_P(pp) = PSIp + + !Derived output data + Bmag = SQRT(B_X*B_X+B_Y*B_Y+B_Z*B_Z) + + ! Parallel unit vector + b_unit_X = B_X/Bmag + b_unit_Y = B_Y/Bmag + b_unit_Z = B_Z/Bmag + + v = SQRT(V_X*V_X+V_Y*V_Y+V_Z*V_Z) + if (v.GT.korc_zero) then + ! Parallel and perpendicular components of velocity + vpar = (V_X*b_unit_X+V_Y*b_unit_Y+ & + V_Z*b_unit_Z) + + vperp = v**2 - vpar**2 + if ( vperp .GE. korc_zero ) then + vperp = SQRT( vperp ) + else + vperp = 0.0_rp + end if + + ! Pitch angle + spp(ii)%vars%eta(pp) = 180.0_rp* & + MODULO(ATAN2(vperp,vpar),2.0_rp*C_PI)/C_PI + + ! Magnetic moment + spp(ii)%vars%mu(pp) = 0.5_rp*m_cache* & + g**2*vperp**2/Bmag + ! See Northrop's book (The adiabatic motion of charged + ! particles) + + ! Radiated power + tmp = q_cache**4/(6.0_rp*C_PI*E0*m_cache**2) + + cross_X = V_Y*B_Z-V_Z*B_Y + cross_Y = V_Z*B_X-V_X*B_Z + cross_Z = V_X*B_Y-V_Y*B_X + + vec_X = E_X + cross_X + vec_Y = E_Y + cross_Y + vec_Z = E_Z + cross_Z + + spp(ii)%vars%Prad(pp) = tmp* & + ( E_X*E_X+E_Y*E_Y+E_Z*E_Z + & + cross_X*E_X+cross_Y*E_Y+ & + cross_Z*E_Z + g**2* & + ((E_X*V_X+E_Y*V_Y+E_Z*V_Z)**2 & + - vec_X*vec_X+vec_Y*vec_Y+ & + vec_Z*vec_Z) ) + + ! Input power due to electric field + spp(ii)%vars%Pin(pp) = q_cache*(E_X*V_X+ & + E_Y*V_Y+E_Z*V_Z) + else + spp(ii)%vars%eta(pp) = 0.0_rp + spp(ii)%vars%mu(pp) = 0.0_rp + spp(ii)%vars%Prad(pp) = 0.0_rp + spp(ii)%vars%Pin(pp) = 0.0_rp + end if + + end do !particle iterator + + end do !species iterator + +end subroutine adv_FOinterp_aorsa_top_ACC + subroutine advance_FOinterp_vars(tt,a,q_cache,m_cache,params,X_X,X_Y,X_Z, & V_X,V_Y,V_Z,B_X,B_Y,B_Z,E_X,E_Y,E_Z,g,flagCon,flagCol,P,F,PSIp) TYPE(KORC_PARAMS), INTENT(IN) :: params diff --git a/src/korc_spatial_distribution.f90 b/src/korc_spatial_distribution.f90 index a544ba4e..54b8e4a6 100755 --- a/src/korc_spatial_distribution.f90 +++ b/src/korc_spatial_distribution.f90 @@ -746,7 +746,6 @@ FUNCTION indicator(psi,psi_max) END FUNCTION indicator - FUNCTION random_norm(mean,sigma) REAL(rp), INTENT(IN) :: mean REAL(rp), INTENT(IN) :: sigma @@ -1322,7 +1321,7 @@ subroutine MH_psi(params,spp,F) !! Present sample of Z location REAL(rp) :: PHI_test - REAL(rp) :: psi_max,psi_max_buff + REAL(rp) :: psi_max,psi_min,psi_max_buff REAL(rp) :: PSIp_lim,PSIP0,PSIN,PSIN0,PSIN1,sigma,psi0,psi1 REAL(rp) :: rand_unif @@ -1364,7 +1363,9 @@ subroutine MH_psi(params,spp,F) max_Z=maxval(F%X%Z) PSIp0=F%PSIP_min + psi_max = spp%psi_max + psi_min = spp%psi_min psi_max_buff = spp%psi_max*2._rp end if @@ -1417,10 +1418,6 @@ subroutine MH_psi(params,spp,F) !Z_test = Z_buffer + get_random_mkl_N(0.0_rp,spp%dZ) Z_test = Z_buffer + get_random_N()*spp%dZ - - - - do while ((R_test.GT.max_R).OR.(R_test .LT. min_R)) !eta_test = eta_buffer + random_norm(0.0_rp,spp%dth) !eta_test = eta_buffer + get_random_mkl_N(0.0_rp,spp%dth) @@ -1517,7 +1514,7 @@ subroutine MH_psi(params,spp,F) if (.not.F%useLCFS) spp%vars%initLCFS(1)=1_is ratio = real(spp%vars%flagCon(1))*real(spp%vars%initLCFS(1))* & - indicator(PSIN1,psi_max)* & + indicator(PSIN1,psi_max)*indicator(psi_min,PSIN1)* & R_test*EXP(-PSIN1/sigma)/ & (R_buffer*EXP(-PSIN0/sigma)) @@ -1625,7 +1622,7 @@ subroutine MH_psi(params,spp,F) if (.not.F%useLCFS) spp%vars%initLCFS(1)=1_is ratio = real(spp%vars%flagCon(1))*real(spp%vars%initLCFS(1))* & - indicator(PSIN1,psi_max)* & + indicator(PSIN1,psi_max)*indicator(psi_min,PSIN1)* & R_test*EXP(-PSIN1/sigma)/ & (R_buffer*EXP(-PSIN0/sigma)) @@ -1656,7 +1653,7 @@ subroutine MH_psi(params,spp,F) ! add to MC above if within buffer. This helps make the boundary ! more defined. IF ((INT(indicator(PSIN1,psi_max)).EQ.1).AND. & - ACCEPTED) THEN + (INT(indicator(psi_min,PSIN1)).EQ.1).AND.ACCEPTED) THEN R_samples(ii) = R_buffer Z_samples(ii) = Z_buffer PHI_samples(ii) = PHI_buffer diff --git a/src/korc_types.f90 b/src/korc_types.f90 index d5f90a9e..a9381eb9 100755 --- a/src/korc_types.f90 +++ b/src/korc_types.f90 @@ -481,7 +481,7 @@ module korc_types REAL(rp) :: theta_gauss !! Angle of counter-clockwise rotation (in degrees) of 2D Gaussian !! distribution relative to R,Z - REAL(rp) :: psi_max + REAL(rp) :: psi_max,psi_min !! Maximum value of the argument of the 2D gaussian exponential, used !! for an indicator function that limits the region of MH sampling REAL(rp) :: Spong_b @@ -621,7 +621,7 @@ module korc_types REAL(rp) :: MARS_phase REAL(rp) :: AORSA_AMP_Scale REAL(rp) :: AORSA_freq - REAL(rp) :: AORSA_nmode + REAL(rp) :: AORSA_nmode,AORSA_mmode !! interpolated E field LOGICAL :: Analytic_D3D_IWL INTEGER :: ntiles diff --git a/src/main.f90 b/src/main.f90 index 13d54346..6b9487a4 100755 --- a/src/main.f90 +++ b/src/main.f90 @@ -492,12 +492,19 @@ program main if (params%orbit_model(1:2).eq.'FO'.and. & params%field_model(10:14).eq.'AORSA') then - call FO_init(params,F,spp,.false.,.true.) +#ifdef ACC + call FO_init_ACC(params,F,spp,.false.,.true.) +#else + call FO_init(params,F,spp,.false.,.true.) +#endif ACC ! Initial half-time particle push do it=params%ito,params%t_steps,params%t_skip +#ifdef ACC + call adv_FOinterp_aorsa_top_ACC(params,F,P,spp) +#else call adv_FOinterp_aorsa_top(params,F,P,spp) - +#endif ACC params%time = params%init_time & +REAL(it-1_ip+params%t_skip,rp)*params%dt params%it = it-1_ip+params%t_skip