Skip to content

Commit

Permalink
creating ACC routines for aorsa fields
Browse files Browse the repository at this point in the history
  • Loading branch information
@mbeidler3
mbeidler3 committed Aug 3, 2023
1 parent 21cf596 commit e4a5356
Show file tree
Hide file tree
Showing 8 changed files with 370 additions and 78 deletions.
1 change: 1 addition & 0 deletions src/korc_fields.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 2 additions & 0 deletions src/korc_initialize.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
11 changes: 7 additions & 4 deletions src/korc_input.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand All @@ -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, &
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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'
Expand Down
101 changes: 93 additions & 8 deletions src/korc_interp.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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\).
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit e4a5356

Please sign in to comment.