diff --git a/build.sh b/build.sh index b434d4e0..4c38edb3 100755 --- a/build.sh +++ b/build.sh @@ -15,8 +15,9 @@ cmake -DCMAKE_BUILD_TYPE:String=$BUILD_TYPE \ -DCMAKE_Fortran_FLAGS="-malign-double -fconvert='big-endian'" \ -DCMAKE_C_FLAGS="-malign-double" \ -DCMAKE_CXX_FLAGS="-malign-double" \ - -DCMAKE_Fortran_FLAGS_DEBUG="-g3 -ffpe-trap=zero,overflow -fbacktrace" \ - -DCMAKE_C_FLAGS_DEBUG="-g3" \ - -DCMAKE_CXX_FLAGS_DEBUG="-g3" ../ + -DCMAKE_Fortran_FLAGS_DEBUG="-g -ffpe-trap=zero,overflow -fbacktrace" \ + -DCMAKE_C_FLAGS_DEBUG="-g" \ + -DCMAKE_CXX_FLAGS_DEBUG="-g" \ + ../ make -j VERBOSE=1 diff --git a/src/korc_HDF5.f90 b/src/korc_HDF5.f90 index ba73e21b..cbb4c725 100755 --- a/src/korc_HDF5.f90 +++ b/src/korc_HDF5.f90 @@ -1547,25 +1547,84 @@ subroutine save_simulation_parameters(params,spp,F,P) call save_1d_array_to_hdf5(h5file_id,dset, & P%X%Z*params%cpp%length,attr_array) - dset = TRIM(gname) // "/ne" - units = params%cpp%density - call rsave_2d_array_to_hdf5(h5file_id, dset, & - units*P%ne_2D) - - dset = TRIM(gname) // "/Te" - units = params%cpp%temperature - call rsave_2d_array_to_hdf5(h5file_id, dset, & - units*P%Te_2D) - - dset = TRIM(gname) // "/Zeff" - call rsave_2d_array_to_hdf5(h5file_id, dset, & - P%Zeff_2D) + if (.not.(params%profile_model(10:12).eq.'Hdt')) then + dset = TRIM(gname) // "/ne" + units = params%cpp%density + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%ne_2D) + + dset = TRIM(gname) // "/Te" + units = params%cpp%temperature + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%Te_2D) + + dset = TRIM(gname) // "/Zeff" + call rsave_2d_array_to_hdf5(h5file_id, dset, & + P%Zeff_2D) + else + dset = TRIM(gname) // "/ne" + units = params%cpp%density + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%ne_3D) + + dset = TRIM(gname) // "/Te" + units = params%cpp%temperature + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%Te_3D) + + dset = TRIM(gname) // "/Zeff" + call rsave_3d_array_to_hdf5(h5file_id, dset, & + P%Zeff_3D) + endif - if (params%profile_model(10:10).eq.'H') then + if (params%profile_model(10:11).eq.'H0') then + + dset = TRIM(gname) // "/RHON" + call rsave_2d_array_to_hdf5(h5file_id, dset, & + P%RHON_2D) + + dset = TRIM(gname) // "/nRE" + units = params%cpp%density + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%nRE_2D) + + dset = TRIM(gname) // "/nAr0" + units = params%cpp%density + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%nAr0_2D) + + dset = TRIM(gname) // "/nAr1" + units = params%cpp%density + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%nAr1_2D) + + dset = TRIM(gname) // "/nAr2" + units = params%cpp%density + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%nAr2_2D) + + dset = TRIM(gname) // "/nAr3" + units = params%cpp%density + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%nAr3_2D) + + dset = TRIM(gname) // "/nD" + units = params%cpp%density + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%nD_2D) + + dset = TRIM(gname) // "/nD1" + units = params%cpp%density + call rsave_2d_array_to_hdf5(h5file_id, dset, & + units*P%nD1_2D) + + endif + + if (params%profile_model(10:12).eq.'Hdt') then dset = TRIM(gname) // "/RHON" - call rsave_2d_array_to_hdf5(h5file_id, dset, & - P%RHON) + call rsave_3d_array_to_hdf5(h5file_id, dset, & + P%RHON_3D) dset = TRIM(gname) // "/nRE" units = params%cpp%density @@ -1574,33 +1633,38 @@ subroutine save_simulation_parameters(params,spp,F,P) dset = TRIM(gname) // "/nAr0" units = params%cpp%density - call rsave_2d_array_to_hdf5(h5file_id, dset, & - units*P%nAr0_2D) + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%nAr0_3D) dset = TRIM(gname) // "/nAr1" units = params%cpp%density - call rsave_2d_array_to_hdf5(h5file_id, dset, & - units*P%nAr1_2D) + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%nAr1_3D) dset = TRIM(gname) // "/nAr2" units = params%cpp%density - call rsave_2d_array_to_hdf5(h5file_id, dset, & - units*P%nAr2_2D) + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%nAr2_3D) dset = TRIM(gname) // "/nAr3" units = params%cpp%density - call rsave_2d_array_to_hdf5(h5file_id, dset, & - units*P%nAr3_2D) + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%nAr3_3D) + + dset = TRIM(gname) // "/nAr4" + units = params%cpp%density + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%nAr4_3D) dset = TRIM(gname) // "/nD" units = params%cpp%density - call rsave_2d_array_to_hdf5(h5file_id, dset, & - units*P%nD_2D) + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%nD_3D) dset = TRIM(gname) // "/nD1" units = params%cpp%density - call rsave_2d_array_to_hdf5(h5file_id, dset, & - units*P%nD1_2D) + call rsave_3d_array_to_hdf5(h5file_id, dset, & + units*P%nD1_3D) endif @@ -2707,20 +2771,33 @@ subroutine save_restart_variables(params,spp,F) !! Iterator for reading all the entried of params::outputs_list. INTEGER :: mpierr !! MPI error status. - INTEGER :: numel_send, numel_receive + INTEGER :: numel_send, numel_receive,exei !! Variable used by MPI to count the amount of data sent by each MPI !! procces. !! Variable used by MPI to count the amount of data received by the main !! MPI procces. + LOGICAL :: exist ! if ( MODULO(params%it,params%restart_output_cadence) .EQ. 0_ip ) then if (params%mpi_params%rank.EQ.0_idef) then - write(output_unit_write,'("Saving restart: ",I15)') & - params%it/(params%t_skip*params%t_it_SC) + write(output_unit_write,'("Saving restart: ",I15)') & + params%it/(params%t_skip*params%t_it_SC) filename = TRIM(params%path_to_outputs) // "restart_file.h5" + + inquire(FILE=filename,exist=exist) + + if (exist) then + call execute_command_line("rm " // TRIM(params%path_to_outputs) // "restart_file.h5",exitstat=exei) + IF (exei/=0) then + write(6,*) 'Error removing restart_file.h5' + call korc_abort(28) + end if + end if + + call h5fcreate_f(TRIM(filename), H5F_ACC_TRUNC_F, h5file_id, h5error) dset = "it" @@ -3123,7 +3200,7 @@ subroutine load_prev_iter(params) call h5fclose_f(h5file_id, h5error) end if - CALL MPI_BCAST(params%prev_iter_2x1t,1,MPI_REAL8,0,MPI_COMM_WORLD,mpierr) + CALL MPI_BCAST(params%prev_iter_2x1t,1,MPI_INTEGER8,0,MPI_COMM_WORLD,mpierr) end subroutine load_prev_iter diff --git a/src/korc_collisions.f90 b/src/korc_collisions.f90 index cde3e84f..1015497b 100755 --- a/src/korc_collisions.f90 +++ b/src/korc_collisions.f90 @@ -226,8 +226,8 @@ subroutine load_params_ms(params) !close(output_unit_write) if (params%profile_model(10:10).eq.'H') then - cparams_ms%num_impurity_species = 6 - params%num_impurity_species = 6 + cparams_ms%num_impurity_species = 7 + params%num_impurity_species = 7 else cparams_ms%num_impurity_species = num_impurity_species params%num_impurity_species = num_impurity_species @@ -263,10 +263,12 @@ subroutine load_params_ms(params) cparams_ms%Zo(3)=18 cparams_ms%Zj(4)=3 cparams_ms%Zo(4)=18 - cparams_ms%Zj(5)=0 - cparams_ms%Zo(5)=1 - cparams_ms%Zj(6)=1 + cparams_ms%Zj(5)=4 + cparams_ms%Zo(5)=18 + cparams_ms%Zj(6)=0 cparams_ms%Zo(6)=1 + cparams_ms%Zj(7)=1 + cparams_ms%Zo(7)=1 cparams_ms%nz(1) = nz_mult(1) @@ -307,10 +309,12 @@ subroutine load_params_ms(params) cparams_ms%aZj(3) = cparams_ms%aAr(3) cparams_ms%IZj(4) = C_E*cparams_ms%IAr(4) cparams_ms%aZj(4) = cparams_ms%aAr(4) - cparams_ms%IZj(5) = C_E*cparams_ms%IH(1) - cparams_ms%aZj(5) = cparams_ms%aH(1) - cparams_ms%IZj(6) = C_E*cparams_ms%IH(2) - cparams_ms%aZj(6) = cparams_ms%aH(2) + cparams_ms%IZj(5) = C_E*cparams_ms%IAr(5) + cparams_ms%aZj(5) = cparams_ms%aAr(5) + cparams_ms%IZj(6) = C_E*cparams_ms%IH(1) + cparams_ms%aZj(6) = cparams_ms%aH(1) + cparams_ms%IZj(7) = C_E*cparams_ms%IH(2) + cparams_ms%aZj(7) = cparams_ms%aH(2) endif cparams_ms%nef = ne_mult + sum(cparams_ms%Zj*cparams_ms%nz) @@ -1086,6 +1090,7 @@ end function VTe_wu function VTe(Te) + !$acc routine seq !! Dimensionless temperature REAL(rp), INTENT(IN) :: Te REAL(rp) :: VTe @@ -1106,12 +1111,15 @@ end function Gammac_wu function Gammacee(v,ne,Te) + !$acc routine seq !! Dimensionless ne and Te REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne REAL(rp), INTENT(IN) :: Te REAL(rp) :: Gammacee + !$acc routine (CLogee) seq + Gammacee = ne*CLogee(v,ne,Te)*cparams_ss%Gammaco end function Gammacee @@ -1194,6 +1202,7 @@ function CLog(ne,Te) ! Dimensionless ne and Te end function CLog function CLog0(ne,Te) ! Dimensionless ne and Te + !$acc routine seq REAL(rp), INTENT(IN) :: ne REAL(rp), INTENT(IN) :: Te REAL(rp) :: CLog0 @@ -1203,7 +1212,7 @@ function CLog0(ne,Te) ! Dimensionless ne and Te end function CLog0 function CLogee(v,ne,Te) - + !$acc routine seq REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne !! ne is in m^-3 and below is converted to cm^-3 @@ -1213,6 +1222,8 @@ function CLogee(v,ne,Te) REAL(rp) :: gam REAL(rp) :: gam_min + !$acc routine (CLog0) seq + gam=1/sqrt(1-v**2) gam_min=cparams_ss%gam_min @@ -1233,7 +1244,7 @@ function CLogee(v,ne,Te) end function CLogee function CLogei(v,ne,Te) - + !$acc routine seq REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne !! ne is in m^-3 and below is converted to cm^-3 @@ -1242,6 +1253,8 @@ function CLogei(v,ne,Te) REAL(rp) :: k=5._rp REAL(rp) :: gam,p + !$acc routine (CLog0) seq + gam=1/sqrt(1-v**2) p=gam*v @@ -1254,14 +1267,18 @@ function CLogei(v,ne,Te) end function CLogei function delta(Te) + !$acc routine seq REAL(rp), INTENT(IN) :: Te REAL(rp) :: delta + !$acc routine (VTe) seq + delta = VTe(Te)*cparams_ss%deltao end function delta function psi(x) + !$acc routine seq REAL(rp), INTENT(IN) :: x REAL(rp) :: psi @@ -1285,12 +1302,14 @@ function CA(v) end function CA function CA_SD(v,ne,Te) + !$acc routine seq REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne REAL(rp), INTENT(IN) :: Te REAL(rp) :: CA_SD REAL(rp) :: x + !$acc routine (Gammacee) seq ! write(6,*) ne,Te x = v/VTe(Te) @@ -1306,6 +1325,7 @@ function CA_SD(v,ne,Te) end function CA_SD function dCA_SD(v,me,ne,Te) + !$acc routine seq REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: me REAL(rp), INTENT(IN) :: ne @@ -1314,6 +1334,8 @@ function dCA_SD(v,me,ne,Te) REAL(rp) :: x real(rp) :: gam + !$acc routine (Gammacee) seq + gam=1/sqrt(1-v**2) x = v/VTe(Te) dCA_SD = Gammacee(v,ne,Te)*((2*(gam*v)**2-1)*psi(x)+ & @@ -1395,6 +1417,7 @@ function CF_FIO(params,v) end function CF_FIO function CF_SD(params,v,ne,Te,P,Y_R,Y_Z) + !$acc routine seq TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne @@ -1407,6 +1430,12 @@ function CF_SD(params,v,ne,Te,P,Y_R,Y_Z) TYPE(PROFILES), INTENT(IN) :: P REAL(rp), INTENT(IN) :: Y_R,Y_Z + !$acc routine (VTe) seq + !$acc routine (Gammacee) seq + !$acc routine (psi) seq + !$acc routine (h_j) seq + !$acc routine (CLogee) seq + x = v/VTe(Te) CF_SD = Gammacee(v,ne,Te)*psi(x)/Te @@ -1458,6 +1487,7 @@ function CF_SD(params,v,ne,Te,P,Y_R,Y_Z) end function CF_SD function CF_SD_FIO(params,v,ne,Te,nimp) + !$acc routine seq TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne @@ -1469,6 +1499,10 @@ function CF_SD_FIO(params,v,ne,Te,nimp) INTEGER :: i REAL(rp) :: k=5._rp + !$acc routine (Gammacee) seq + !$acc routine (CLogee) seq + !$acc routine (h_j) seq + x = v/VTe(Te) CF_SD_FIO = Gammacee(v,ne,Te)*psi(x)/Te @@ -1579,6 +1613,7 @@ function CB_ei_FIO(params,v) end function CB_ei_FIO function CB_ee_SD(v,ne,Te,Zeff) + !$acc routine seq REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne REAL(rp), INTENT(IN) :: Te @@ -1586,6 +1621,10 @@ function CB_ee_SD(v,ne,Te,Zeff) REAL(rp) :: CB_ee_SD REAL(rp) :: x + !$acc routine (Gammacee) seq + !$acc routine (psi) seq + !$acc routine (delta) seq + x = v/VTe(Te) CB_ee_SD = (0.5_rp*Gammacee(v,ne,Te)/v)* & (ERF(x) - psi(x) + & @@ -1593,6 +1632,7 @@ function CB_ee_SD(v,ne,Te,Zeff) end function CB_ee_SD function CB_ei_SD(params,v,ne,Te,Zeff,P,Y_R,Y_Z) + !$acc routine seq TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne @@ -1605,6 +1645,12 @@ function CB_ei_SD(params,v,ne,Te,Zeff,P,Y_R,Y_Z) TYPE(PROFILES), INTENT(IN) :: P REAL(rp), INTENT(IN) :: Y_R,Y_Z + !$acc routine (VTe) seq + !$acc routine (Gammacee) seq + !$acc routine (CLogei) seq + !$acc routine (CLogee) seq + !$acc routine (g_j) seq + x = v/VTe(Te) CB_ei_SD = (0.5_rp*Gammacee(v,ne,Te)/v)* & (Zeff*CLogei(v,ne,Te)/CLogee(v,ne,Te)) @@ -1648,6 +1694,7 @@ function CB_ei_SD(params,v,ne,Te,Zeff,P,Y_R,Y_Z) end function CB_ei_SD function CB_ei_SD_FIO(params,v,ne,Te,nimp,Zeff) + !$acc routine seq TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp), INTENT(IN) :: v REAL(rp), INTENT(IN) :: ne @@ -1659,6 +1706,11 @@ function CB_ei_SD_FIO(params,v,ne,Te,nimp,Zeff) REAL(rp) :: x INTEGER :: i + !$acc routine (Gammacee) seq + !$acc routine (CLogee) seq + !$acc routine (CLogei) seq + !$acc routine (g_j) seq + x = v/VTe(Te) CB_ei_SD_FIO = (0.5_rp*Gammacee(v,ne,Te)/v)* & (Zeff*CLogei(v,ne,Te)/CLogee(v,ne,Te)) @@ -1704,6 +1756,7 @@ function nu_S_FIO(params,v) end function nu_S_FIO function h_j(i,v) + !$acc routine seq INTEGER, INTENT(IN) :: i REAL(rp), INTENT(IN) :: v REAL(rp) :: gam @@ -1718,6 +1771,7 @@ function h_j(i,v) end function h_j function g_j(i,v) + !$acc routine seq INTEGER, INTENT(IN) :: i REAL(rp), INTENT(IN) :: v REAL(rp) :: gam @@ -2293,7 +2347,7 @@ end subroutine include_CoulombCollisions_FOfio_p #endif FIO subroutine include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & - Ppll,Pmu,me,flagCon,flagCol,F,P,E_PHI,ne,PSIp) + Ppll,Pmu,me,flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) TYPE(PROFILES), INTENT(IN) :: P TYPE(FIELDS), INTENT(IN) :: F @@ -2306,13 +2360,12 @@ subroutine include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & REAL(rp), DIMENSION(params%pchunk) :: curlb_R,curlb_PHI,curlb_Z REAL(rp), DIMENSION(params%pchunk) :: gradB_R,gradB_PHI,gradB_Z REAL(rp), DIMENSION(params%pchunk) :: E_R,E_Z - REAL(rp), DIMENSION(params%pchunk), INTENT(OUT) :: E_PHI,ne,PSIp + REAL(rp), DIMENSION(params%pchunk), INTENT(OUT) :: E_PHI,ne,PSIp,Te,Zeff REAL(rp), DIMENSION(params%pchunk), INTENT(IN) :: Y_R,Y_PHI,Y_Z INTEGER(is), DIMENSION(params%pchunk), INTENT(INOUT) :: flagCol INTEGER(is), DIMENSION(params%pchunk), INTENT(INOUT) :: flagCon REAL(rp), INTENT(IN) :: me - REAL(rp), DIMENSION(params%pchunk) :: Te,Zeff - REAL(rp), DIMENSION(params%pchunk) :: nAr0,nAr1,nAr2,nAr3 + REAL(rp), DIMENSION(params%pchunk) :: nAr0,nAr1,nAr2,nAr3,nAr4 REAL(rp), DIMENSION(params%pchunk) :: nD,nD1 REAL(rp), DIMENSION(params%pchunk,2) :: dW REAL(rp), DIMENSION(params%pchunk,2) :: rnd1 @@ -2331,7 +2384,7 @@ subroutine include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & REAL(rp) :: kappa integer :: cc,pchunk integer(ip),INTENT(IN) :: tt - REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species) :: nimp + REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species), INTENT(OUT) :: nimp REAL(rp), DIMENSION(params%pchunk) :: E_PHI_tmp pchunk=params%pchunk @@ -2372,16 +2425,20 @@ subroutine include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & else if (params%profile_model(1:8).eq.'EXTERNAL') then #ifdef PSPLINE if (params%profile_model(10:10).eq.'H') then - call interp_Hcollision_p(pchunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, & - nAr0,nAr1,nAr2,nAr3,nD,nD1,flagCon) + call interp_Hcollision_p(params,pchunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,flagCon) do cc=1_idef,pchunk nimp(cc,1)=nAr0(cc) nimp(cc,2)=nAr1(cc) nimp(cc,3)=nAr2(cc) nimp(cc,4)=nAr3(cc) - nimp(cc,5)=nD(cc) - nimp(cc,6)=nD1(cc) + nimp(cc,5)=nAr4(cc) + nimp(cc,6)=nD(cc) + nimp(cc,7)=nD1(cc) end do + + !write(6,*) 'collisions ne',ne(1)*params%cpp%density + else call interp_FOcollision_p(pchunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff,flagCon) endif @@ -2577,8 +2634,220 @@ subroutine include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & end subroutine include_CoulombCollisions_GC_p +subroutine include_CoulombCollisions_GC_p_ACC(spp,pp,tt,params,random, & + Y_R,Y_PHI,Y_Z,Ppll,Pmu,me,flagCon,flagCol,F,P,E_PHI, & + ne,Te,Zeff,nimp,PSIp) + !$acc routine seq + TYPE(SPECIES), INTENT(INOUT) :: spp + TYPE(PROFILES), INTENT(IN) :: P + TYPE(FIELDS), INTENT(IN) :: F + TYPE(KORC_PARAMS), INTENT(INOUT) :: params + CLASS(random_context), INTENT(INOUT) :: random + REAL(rp), INTENT(INOUT) :: Ppll + REAL(rp), INTENT(INOUT) :: Pmu + REAL(rp) :: Bmag + REAL(rp) :: B_R,B_PHI,B_Z + REAL(rp) :: curlb_R,curlb_PHI,curlb_Z + REAL(rp) :: gradB_R,gradB_PHI,gradB_Z + REAL(rp) :: E_R,E_Z,E_PHI_LAC + REAL(rp), INTENT(OUT) :: E_PHI,ne,PSIp,Te,Zeff + REAL(rp), INTENT(IN) :: Y_R,Y_PHI,Y_Z + INTEGER(is), INTENT(INOUT) :: flagCol + INTEGER(is), INTENT(INOUT) :: flagCon + REAL(rp),INTENT(IN) :: me + REAL(rp) :: nAr0,nAr1,nAr2,nAr3,nAr4 + REAL(rp) :: nD,nD1 + REAL(rp),DIMENSION(2) :: dW + REAL(rp),DIMENSION(2) :: rnd1 + REAL(rp) :: dt,time + REAL(rp) :: pm + REAL(rp) :: dp + REAL(rp) :: xi + REAL(rp) :: dxi + REAL(rp) :: v,gam + !! speed of particle + REAL(rp) :: CAL + REAL(rp) :: dCAL + REAL(rp) :: CFL + REAL(rp) :: CBL + REAL(rp) :: SC_p,SC_mu,BREM_p + REAL(rp) :: kappa + integer(ip),INTENT(IN) :: tt + INTEGER,INTENT(IN) :: pp + REAL(rp), DIMENSION(params%num_impurity_species), INTENT(OUT) :: nimp + REAL(rp) :: E_PHI_tmp,ntot,ra + INTEGER :: ii + + !$acc routine (calculate_GCfieldswE_p_ACC) seq + !$acc routine (analytical_profiles_p_ACC) seq + !$acc routine (interp_Hcollision_p_ACC) seq + !$acc routine (CA_SD) seq + !$acc routine (dCA_SD) seq + !$acc routine (CF_SD_FIO) seq + !$acc routine (CF_SD) seq + !$acc routine (CB_ee_SD) seq + !$acc routine (CB_ei_SD_FIO) seq + !$acc routine (CB_ei_SD) seq + !$acc routine (get_random_U) seq + !$acc routine (large_angle_source_ACC) seq + + if (MODULO(params%it+tt,cparams_ss%subcycling_iterations) .EQ. 0_ip) then + dt = REAL(cparams_ss%subcycling_iterations,rp)*params%dt + time=params%init_time+(params%it-1+tt)*params%dt + +#ifdef PSPLINE + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z, & + B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp) +#endif + + if (params%profile_model(1:10).eq.'ANALYTICAL') then + call analytical_profiles_p_ACC(time,params,Y_R,Y_Z, & + P,F,ne,Te,Zeff,PSIp) +#ifdef PSPLINE + else if (params%profile_model(10:10).eq.'H') then + call interp_Hcollision_p_ACC(params,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,flagCon) + + nimp(1)=nAr0 + nimp(2)=nAr1 + nimp(3)=nAr2 + nimp(4)=nAr3 + nimp(5)=nAr4 + nimp(6)=nD + nimp(7)=nD1 + +#endif PSPLINE + endif + + + E_PHI_LAC=E_PHI + E_PHI_tmp=E_PHI + if (.not.params%FokPlan) E_PHI=0._rp + + Bmag=sqrt(B_R*B_R+B_PHI*B_PHI+B_Z*B_Z) + ! Transform p_pll,mu to P,eta + pm = SQRT(Ppll*Ppll+2*me*Bmag*Pmu) + xi = Ppll/pm + + gam = sqrt(1+pm*pm) + + v = pm/gam + ! normalized speed (v_K=v_P/c) + + CALL random%uniform%set(0.0_rp, 1.0_rp) + + rnd1(1) = random%uniform%get() + rnd1(2) = random%uniform%get() + ! rnd1(:,1) = get_random_mkl() + ! rnd1(:,2) = get_random_mkl() + + dW(1) = SQRT(3*dt)*(-1+2*rnd1(1)) + dW(2) = SQRT(3*dt)*(-1+2*rnd1(2)) + + ! write(output_unit_write,'("dW1: ",E17.10)') dW(cc,1) + ! write(output_unit_write,'("dW2: ",E17.10)') dW(cc,2) + + if (params%profile_model(10:10).eq.'H') then + CAL = CA_SD(v,ne,Te) + dCAL= dCA_SD(v,me,ne,Te) + CFL = CF_SD_FIO(params,v,ne,Te,nimp) + CBL = (CB_ee_SD(v,ne,Te,Zeff)+ & + CB_ei_SD_FIO(params,v,ne,Te,nimp,Zeff)) + else + CAL = CA_SD(v,ne,Te) + dCAL= dCA_SD(v,me,ne,Te) + CFL = CF_SD(params,v,ne,Te,P,Y_R,Y_Z) + CBL = (CB_ee_SD(v,ne,Te,Zeff)+ & + CB_ei_SD(params,v,ne,Te,Zeff,P,Y_R,Y_Z)) + endif + + if (.not.cparams_ss%slowing_down) CFL=0._rp + if (.not.cparams_ss%pitch_diffusion) CBL=0._rp + if (.not.cparams_ss%energy_diffusion) THEN + CAL=0._rp + dCAL=0._rp + ENDIF + + dp=REAL(flagCol)*REAL(flagCon)* & + ((-CFL+dCAL+E_PHI*xi)*dt+ & + sqrt(2.0_rp*CAL)*dW(1)) + + dxi=REAL(flagCol)*REAL(flagCon)* & + ((-2*xi*CBL/(pm*pm)+ & + E_PHI*(1-xi*xi)/pm)*dt- & + sqrt(2.0_rp*CBL*(1-xi*xi))/pm*dW(2)) + + + pm=pm+dp + xi=xi+dxi + + if (xi>1) then + + xi=1-mod(xi,1._rp) + else if (xi<-1) then + xi=-1-mod(xi,-1._rp) + + endif + + ! Transform P,xi to p_pll,mu + Ppll=pm*xi + Pmu=(pm*pm-Ppll*Ppll)/(2*me*Bmag) + + if ((pm.lt.min(cparams_ss%p_min*cparams_ss%pmin_scale, & + p_therm)).and.flagCol.eq.1_ip) then + flagCol=0_ip + end if + + if (cparams_ss%avalanche) then + + ntot=ne + if (params%bound_electron_model.eq.'HESSLOW') then + + if (params%profile_model(1:10).eq.'ANALYTICAL') then + if ((cparams_ms%Zj(1).eq.0.0).and. & + (neut_prof.eq.'UNIFORM')) then + ntot=ntot+cparams_ms%nz(1)* & + (cparams_ms%Zo(1)-cparams_ms%Zj(1)) + else if ((cparams_ms%Zj(1).eq.0.0).and. & + (neut_prof.eq.'HOLLOW')) then + ntot=ntot+max(cparams_ms%nz(1)-ne,0._rp)* & + (cparams_ms%Zo(1)-cparams_ms%Zj(1)) + else if ((cparams_ms%Zj(1).eq.0.0).and. & + (neut_prof.eq.'EDGE')) then + ra=sqrt((Y_R-P%R0)**2+(Y_Z-P%Z0)**2)/P%a + ntot=ntot+cparams_ms%nz(1)*ra**cparams_ms%neut_edge_fac* & + (cparams_ms%Zo(1)-cparams_ms%Zj(1)) + else + ntot=ntot+ne*cparams_ms%nz(1)/cparams_ms%ne* & + (cparams_ms%Zo(1)-cparams_ms%Zj(1)) + endif + + do ii=2,cparams_ms%num_impurity_species + ntot=ntot+ne*cparams_ms%nz(ii)/cparams_ms%ne* & + (cparams_ms%Zo(ii)-cparams_ms%Zj(ii)) + end do + else if (params%profile_model(10:10).eq.'H') then + ntot=ntot+nimp(1)*18+nimp(2)*17+nimp(3)*16+nimp(4)*15+nimp(5)*14+ & + nimp(6) + end if + + end if + + call large_angle_source_ACC(spp,params,random,F,Y_R,Y_PHI,Y_Z, & + pm,xi,ne,ntot,Te,Bmag,E_PHI_LAC,me,flagCol,flagCon,B_R,B_PHI,B_Z) + + end if + + if (.not.params%FokPlan) E_PHI=E_PHI_tmp + + end if + +end subroutine include_CoulombCollisions_GC_p_ACC + subroutine include_CoulombCollisionsLA_GC_p(spp,achunk,tt,params,random, & Y_R,Y_PHI,Y_Z,Ppll,Pmu,me,flagCon,flagCol,F,P,E_PHI,ne,Te,PSIp) + TYPE(SPECIES), INTENT(INOUT) :: spp TYPE(PROFILES), INTENT(IN) :: P TYPE(FIELDS), INTENT(IN) :: F @@ -2598,7 +2867,7 @@ subroutine include_CoulombCollisionsLA_GC_p(spp,achunk,tt,params,random, & INTEGER(is), DIMENSION(achunk), INTENT(INOUT) :: flagCon REAL(rp), INTENT(IN) :: me REAL(rp), DIMENSION(achunk) :: Zeff - REAL(rp), DIMENSION(achunk) :: nAr0,nAr1,nAr2,nAr3 + REAL(rp), DIMENSION(achunk) :: nAr0,nAr1,nAr2,nAr3,nAr4 REAL(rp), DIMENSION(achunk) :: nD,nD1 REAL(rp), DIMENSION(achunk,2) :: dW REAL(rp), DIMENSION(achunk,2) :: rnd1 @@ -2654,16 +2923,17 @@ subroutine include_CoulombCollisionsLA_GC_p(spp,achunk,tt,params,random, & ne,Te,Zeff,PSIp) else if (params%profile_model(1:8).eq.'EXTERNAL') then #ifdef PSPLINE - if (params%profile_model(10:10).eq.'H') then - call interp_Hcollision_p(achunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, & - nAr0,nAr1,nAr2,nAr3,nD,nD1,flagCon) + if (params%profile_model(10:11).eq.'H') then + call interp_Hcollision_p(params,achunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,flagCon) do cc=1_idef,achunk nimp(cc,1)=nAr0(cc) nimp(cc,2)=nAr1(cc) nimp(cc,3)=nAr2(cc) nimp(cc,4)=nAr3(cc) - nimp(cc,5)=nD(cc) - nimp(cc,6)=nD1(cc) + nimp(cc,5)=nAr4(cc) + nimp(cc,6)=nD(cc) + nimp(cc,7)=nD1(cc) end do else call interp_FOcollision_p(achunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff,flagCon) @@ -3334,19 +3604,6 @@ subroutine large_angle_source(spp,params,random,achunk,F,Y_R,Y_PHI,Y_Z, & CALL random%uniform%set(0.0_rp, 1.0_rp) - !$OMP SIMD - do cc=1_idef,achunk - pm0(cc)=pm(cc) - xi0(cc)=xi(cc) - - gam(cc) = sqrt(1+pm(cc)*pm(cc)) - gam0(cc)=gam(cc) - - prob0(cc) = random%uniform%get() - - end do - !$OMP END SIMD - vmin=1/sqrt(1+1/(cparams_ss%p_min*cparams_ss%pmin_scale)**2) !! For each primary RE, calculating probability to generate a secondary RE @@ -3361,7 +3618,13 @@ subroutine large_angle_source(spp,params,random,achunk,F,Y_R,Y_PHI,Y_Z, & E_C=Gammacee(vmin,ne(cc),Te(cc)) end if - + pm0(cc)=pm(cc) + xi0(cc)=xi(cc) + + gam(cc) = sqrt(1+pm(cc)*pm(cc)) + gam0(cc)=gam(cc) + + prob0(cc) = random%uniform%get() !write(6,*) 'E',E_PHI*params%cpp%Eo !write(6,*) 'E_C',E_C*params%cpp%Eo @@ -3718,6 +3981,384 @@ subroutine large_angle_source(spp,params,random,achunk,F,Y_R,Y_PHI,Y_Z, & end subroutine large_angle_source + subroutine large_angle_source_ACC(spp,params,random,F,Y_R,Y_PHI,Y_Z, & + pm,xi,ne,netot,Te,Bmag,E_PHI,me,flagCol,flagCon,B_R,B_PHI,B_Z) + !$acc routine seq + TYPE(SPECIES), INTENT(INOUT) :: spp + TYPE(KORC_PARAMS), INTENT(IN) :: params + CLASS(random_context), INTENT(INOUT) :: random + TYPE(FIELDS), INTENT(IN) :: F + REAL(rp), INTENT(INOUT) :: pm,xi + REAL(rp), INTENT(IN) :: Y_R,Y_PHI,Y_Z + REAL(rp), INTENT(IN) :: B_R,B_PHI,B_Z + REAL(rp), INTENT(IN) :: ne,netot,Te + REAL(rp), INTENT(IN) :: Bmag,E_PHI + INTEGER(is), INTENT(IN) :: flagCol,flagCon + REAL(rp), INTENT(IN) :: me + REAL(rp) :: gam,prob0,prob1,pm0,xi0,gam0 + REAL(rp) :: gam_min,p_min,gammax,dt,gamsecmax,psecmax,ptrial + REAL(rp) :: gamtrial,cosgam1,xirad,xip,xim,xitrial,sinsq1,cossq1,pitchprob1 + REAL(rp) :: dsigdgam1,S_LAmax,S_LA1,tmppm,gamvth,vmin,E_C,p_c,gam_c,pRE + INTEGER :: ngam1,neta1 + INTEGER :: ii,jj,cc,seciter + REAL(rp), DIMENSION(cparams_ss%ngrid1) :: gam1,pm1,tmpgam1,tmpcosgam,tmpdsigdgam,tmpsecthreshgam,probtmp,intpitchprob + REAL(rp), DIMENSION(cparams_ss%ngrid1-1) :: dpm1 + REAL(rp), DIMENSION(cparams_ss%ngrid1) :: eta1,tmpsinsq,tmpcossq + REAL(rp), DIMENSION(cparams_ss%ngrid1-1) :: deta1 + REAL(rp), DIMENSION(cparams_ss%ngrid1,cparams_ss%ngrid1) :: cosgam,sinsq,cossq,tmpcossq1,pitchprob,dsigdgam + REAL(rp), DIMENSION(cparams_ss%ngrid1,cparams_ss%ngrid1) :: secthreshgam,pm11,gam11,eta11,S_LA,pitchrad + LOGICAL :: accepted + + !$acc routine (Gammacee) seq + !$acc routine (get_random) seq + + CALL random%uniform%set(0.0_rp, 1.0_rp) + + + dt=cparams_ss%coll_per_dump_dt*params%cpp%time + ngam1=cparams_ss%ngrid1 + neta1=cparams_ss%ngrid1 + + + pm0=pm + xi0=xi + + gam = sqrt(1+pm*pm) + gam0=gam + + + prob0 = random%uniform%get() + + + + vmin=1/sqrt(1+1/(cparams_ss%p_min*cparams_ss%pmin_scale)**2) + + !! For each primary RE, calculating probability to generate a secondary RE + + if ((flagCol.lt.1).or.(flagCon.lt.1)) return + + if (.not.(cparams_ms%lowKE_REs)) then + E_C=netot/ne*Gammacee(vmin,ne,Te) + else + E_C=Gammacee(vmin,ne,Te) + end if + + !write(6,*) 'E',E_PHI*params%cpp%Eo + !write(6,*) 'E_C',E_C*params%cpp%Eo + !write(6,*) 'E_c,min',cparams_ms%Ec_min*params%cpp%Eo + !write(6,*) 'ne',ne*params%cpp%density + !write(6,*) 'netot',netot*params%cpp%density + !write(6,*) 'Te',Te*params%cpp%temperature + !write(6,*) 'Clog',CLogee_wu(params,ne*params%cpp%density,Te*params%cpp%temperature) + + if (.not.(cparams_ss%always_aval)) then + if (E_C.gt.abs(E_PHI)) return + + p_c=cparams_ss%pmin_scale/sqrt(abs(E_PHI)/E_C-1) + gam_c=sqrt(1+p_c**2) + + if(cparams_ss%min_secRE.eq.'THERM') then + gam_min=(gam_c+1)/2 + p_min=sqrt(gam_min**2-1) + else + gam_min=gam_c + p_min=p_c + end if + else + p_min=cparams_ss%p_min*cparams_ss%pmin_scale + gam_min=sqrt(1+p_min**2) + endif + + gammax=(gam+1._rp)/2._rp + + if (gam_min.eq.1._rp) then + write(6,*) 'R',Y_R*params%cpp%length,'Z',Y_Z*params%cpp%length + write(6,*) 'vmin',vmin,'netot',netot*params%cpp%density,'Te',Te*params%cpp%temperature/C_E + write(6,*) 'E',E_PHI*params%cpp%Eo,'E_c',E_C*params%cpp%Eo + write(6,*) 'p_c',p_c,'gam_c',gam_c + write(6,*) 'p_min',p_min,'gam_min',gam_min + !write(6,*) 'LAC_gam_resolution: ',TRIM(cparams_ss%LAC_gam_resolution) + !write(6,*) 'gam_min,gammax',gam_min,gammax + end if + + !! Generating 1D and 2D ranges for secondary RE distribution + + if (TRIM(cparams_ss%LAC_gam_resolution).eq.'LIN') then + do ii=1,ngam1 + gam1(ii)=gam_min+(gammax-gam_min)* & + REAL(ii-1)/REAL(ngam1-1) + end do + elseif (TRIM(cparams_ss%LAC_gam_resolution).eq.'EXP') then + do ii=1,ngam1 + tmpgam1(ii)=log10(gam_min)+ & + (log10(gammax)-log10(gam_min))* & + REAL(ii-1)/REAL(ngam1-1) + end do + gam1=10**(tmpgam1) + elseif (TRIM(cparams_ss%LAC_gam_resolution).eq.'2EXP') then + do ii=1,ngam1 + tmpgam1(ii)=log10(log10(gam_min))+ & + (log10(log10(gammax))-log10(log10(gam_min)))* & + REAL(ii-1)/REAL(ngam1-1) + end do + gam1=10**(10**(tmpgam1)) + end if + + !write(6,*) 'tmpgam1',tmpgam1 + !write(6,*) 'gam1',gam1 + + + pm1=sqrt(gam1**2-1) + + do ii=1,ngam1-1 + dpm1(ii)=pm1(ii+1)-pm1(ii) + end do + + do ii=1,neta1 + eta1(ii)=C_PI*(ii-1)/(neta1-1) + end do + + !write(6,*) 'eta1',eta1 + + do ii=1,neta1-1 + deta1(ii)=eta1(ii+1)-eta1(ii) + end do + + do ii=1,neta1 + pm11(:,ii)=pm1 + gam11(:,ii)=gam1 + end do + + do ii=1,ngam1 + eta11(ii,:)=eta1 + end do + + + tmpcosgam=sqrt(((gam+1)*(gam1-1))/((gam-1)*(gam1+1))) + tmpdsigdgam=2*C_PI*C_RE**2/(gam**2-1)* & + (((gam-1)**2*gam**2)/((gam1-1)**2*(gam-gam1)**2)- & + (2*gam**2+2*gam-1)/((gam1-1)*(gam-gam1))+1) + tmpsecthreshgam=1._rp + where(gam1.gt.(gam+1)/2._rp) tmpsecthreshgam=0._rp + + !write(6,*) 'tmpcosgam',tmpcosgam + !write(6,*) 'tmpdsigdgam',tmpdsigdgam + !write(6,*) 'tmpsecthreshgam',tmpsecthreshgam + + + do ii=1,neta1 + cosgam(:,ii)=tmpcosgam + dsigdgam(:,ii)=tmpdsigdgam + secthreshgam(:,ii)=tmpsecthreshgam + end do + + !if (cc.eq.1) then + !write(6,*) cosgam + !end if + + tmpsinsq=(1-xi**2)*sin(eta1)**2 + tmpcossq=xi*cos(eta1) + + do ii=1,ngam1 + sinsq(ii,:)=tmpsinsq + tmpcossq1(ii,:)=tmpcossq + end do + + !if (cc.eq.1) then + !write(6,*) sinsq + !end if + + cossq=(cosgam-tmpcossq1)**2 + + pitchrad=sinsq-cossq + where(pitchrad.lt.0) pitchrad=tiny(0._rp) + + !if (cc.eq.1) then + !write(6,*) cossq + !write(6,*) 'sinsq-cossq',sqrt(sinsq-cossq) + !write(6,*) 'pitchrad',pitchrad + !end if + + pitchprob=1/(C_PI*sqrt(pitchrad)) + + where(pitchprob.eq.1/(C_PI*sqrt(tiny(0._rp)))) pitchprob=0._rp + + !if (cc.eq.1) then + !write(6,*) 'pitchprob',pitchprob + !end if + + S_LA=netot*params%cpp%density*C_C/(2*C_PI)* & + (pm11/gam11)*(pm/gam)* & + pitchprob*dsigdgam*secthreshgam + + !! Saving maximum secondary RE source for use in rejection-acceptance + !! sampling algorithm + S_LAmax=maxval(S_LA) + + S_LA=S_LA*sin(eta11) + + !! Trapezoidal integration of secondary RE source to find probabilty + + do ii=1,ngam1 + probtmp(ii)=S_LA(ii,1)*deta1(1)/2+S_LA(ii,neta1)*deta1(neta1-1)/2 + !intpitchprob(ii)=pitchprob(ii,1)*sin(eta1(1))*deta1(1)/2+ & + ! pitchprob(ii,neta1)*sin(eta1(neta1))*deta1(neta1-1)/2 + do jj=2,neta1-1 + probtmp(ii)=probtmp(ii)+S_LA(ii,jj)*(deta1(jj)+deta1(jj-1))/2 + ! intpitchprob(ii)=intpitchprob(ii)+ & + ! pitchprob(ii,jj)*sin(eta1(jj))*(deta1(jj)+deta1(jj-1))/2 + end do + end do + + prob1=probtmp(1)*dpm1(1)/2+probtmp(ngam1)*dpm1(ngam1-1)/2 + do jj=2,ngam1-1 + prob1=prob1+probtmp(jj)*(dpm1(jj)+dpm1(jj-1))/2 + end do + + !write(6,*) 'prob1pre',prob1,'flagCol',flagCol,'flagCon',flagCon,'dt',dt + + !write(6,*) 'intpitchprob',intpitchprob + + prob1=prob1*dt*2*C_PI + +#ifdef __NVCOMPILER + if (IEEE_IS_NAN(prob1)) then + +#else + if (ISNAN(prob1)) then +#endif + write(6,*) 'NaN probability from secondary RE source' + write(6,*) 'p,xi',pm,xi + write(6,*) 'gam_min,gammax',gam_min,gammax + write(6,*) 'E',E_PHI*params%cpp%Eo + write(6,*) 'E_C',E_C*params%cpp%Eo + !write(6,*) 'pitchprob',pitchprob + !write(6,*) 'S_LA',S_LA + call korc_abort(24) + end if + + if (prob1.gt.1._rp) then + write(6,*) 'Multiple secondary REs generated in a collisional time step' + write(6,*) 'p,xi',pm,xi + write(6,*) 'gam_min,gammax',gam_min,gammax + write(6,*) 'E',E_PHI*params%cpp%Eo + write(6,*) 'E_C',E_C*params%cpp%Eo + call korc_abort(24) + end if + + !write(6,*) 'gam',gam,'xi',xi + !write(6,*) 'prob1',prob1,'prob0',prob0 + + if (prob1.gt.prob0) then + + !! If secondary RE generated, begin acceptance-rejection sampling + !! algorithm + + !write(6,*) 'secondary RE from ',prob1,prob0 + + accepted=.false. + + gamsecmax=(gam+1)/2 + psecmax=sqrt(gamsecmax**2-1) + + seciter=0 + do while (.not.accepted) + + ptrial=p_min+(psecmax-p_min)*random%uniform%get() + gamtrial=sqrt(1+ptrial*ptrial) + + cosgam1=sqrt(((gam+1)*(gamtrial-1))/((gam-1)*(gamtrial+1))) + + xitrial=-1+2*random%uniform%get() + + sinsq1=(1-xi*xi)*(1-xitrial*xitrial) + cossq1=(cosgam1-xi*xitrial)**2 + + if ((sinsq1-cossq1).le.0._rp) cycle + + pitchprob1=1/(C_PI*sqrt(sinsq1-cossq1)) + + !if (isnan(pitchprob1)) cycle + + dsigdgam1=2*C_PI*C_RE**2/(gam**2-1)* & + (((gam-1)**2*gam**2)/ & + ((gamtrial-1)**2*(gam-gamtrial)**2)- & + (2*gam**2+2*gam-1)/ & + ((gamtrial-1)*(gam-gamtrial))+1) + + S_LA1=netot*params%cpp%density*C_C/(2*C_PI)* & + (ptrial/gamtrial)*(pm/gam)* & + pitchprob1*dsigdgam1 + + if (S_LA1/S_LAmax.gt.random%uniform%get()) accepted=.true. + + seciter=seciter+1 + + !if (mod(seciter,100).eq.0) then + ! write(6,*) 'iteration',seciter + !end if + + end do + + !! Write secondary RE degrees of freedom to particle derived type + + + !$acc ATOMIC UPDATE + spp%pRE=spp%pRE+1 + !$acc ATOMIC WRITE + spp%vars%flagRE(spp%pRE)=1 + !$acc ATOMIC WRITE + spp%vars%flagCon(spp%pRE)=1 + !$acc ATOMIC WRITE + spp%vars%flagCol(spp%pRE)=1 + !$acc ATOMIC WRITE + spp%vars%V(spp%pRE,1)=ptrial*xitrial + !$acc ATOMIC WRITE + spp%vars%V(spp%pRE,2)=ptrial*ptrial*(1-xitrial*xitrial)/ & + (2*me*Bmag) + !$acc ATOMIC WRITE + spp%vars%Y(spp%pRE,1)=Y_R + !$acc ATOMIC WRITE + spp%vars%Y(spp%pRE,2)=Y_PHI + !$acc ATOMIC WRITE + spp%vars%Y(spp%pRE,3)=Y_Z + !$acc ATOMIC WRITE + spp%vars%Yborn(spp%pRE,1)=Y_R + !$acc ATOMIC WRITE + spp%vars%Yborn(spp%pRE,2)=Y_PHI + !$acc ATOMIC WRITE + spp%vars%Yborn(spp%pRE,3)=Y_Z + !$acc ATOMIC WRITE + spp%vars%B(spp%pRE,1)=B_R + !$acc ATOMIC WRITE + spp%vars%B(spp%pRE,2)=B_PHI + !$acc ATOMIC WRITE + spp%vars%B(spp%pRE,3)=B_Z + + + !! Write changes to primary RE degrees of freedom to temporary + !! arrays for passing back out to particle derived type + + if (cparams_ss%ConserveLA) then + tmppm=pm + gamvth=1/sqrt(1-2*Te) + gam=gam-gamtrial+gamvth + pm=sqrt(gam*gam-1) + xi=(tmppm*xi-ptrial*xitrial)/pm + end if + + !$acc ATOMIC READ + pRE=spp%pRE + + if (pRE.eq.spp%ppp) then + write(6,*) 'All REs allocated on proc',params%mpi_params%rank + call korc_abort(24) + end if + + end if + +end subroutine large_angle_source_ACC + subroutine save_params_ms(params) TYPE(KORC_PARAMS), INTENT(IN) :: params diff --git a/src/korc_experimental_pdf.f90 b/src/korc_experimental_pdf.f90 index b3e4453f..83120d15 100755 --- a/src/korc_experimental_pdf.f90 +++ b/src/korc_experimental_pdf.f90 @@ -590,11 +590,11 @@ FUNCTION fRE_H_3D(params,F,eta,g,R,Z,R0,Z0,EPHI,rho_ind) END FUNCTION fRE_H_3D - FUNCTION fRE_H_pitch(params,eta,g,EPHI,ne,Te,nAr0,nAr1,nAr2,nAr3,nD,nD1) + FUNCTION fRE_H_pitch(params,eta,g,EPHI,ne,Te,nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1) TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp), INTENT(IN) :: eta ! pitch angle in degrees REAL(rp), INTENT(IN) :: g ! Relativistic gamma factor - REAL(rp), INTENT(IN) :: EPHI,ne,Te,nAr0,nAr1,nAr2,nAr3,nD,nD1 + REAL(rp), INTENT(IN) :: EPHI,ne,Te,nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1 REAL(rp) :: fRE_H_pitch REAL(rp) :: A REAL(rp) :: E_G,nf,pdm,Z_brac @@ -611,9 +611,9 @@ FUNCTION fRE_H_pitch(params,eta,g,EPHI,ne,Te,nAr0,nAr1,nAr2,nAr3,nD,nD1) !write(6,*) 'VTe',VTe,'pdm',pdm !write(6,*) 'CLog0',CLog0,'CLogee',CLogee,'CLogei',CLogei - nf=nAr0*18+nAr1*17+nAr2*16+nAr3*15+nD + nf=nAr0*18+nAr1*17+nAr2*16+nAr3*15+nAr4*14+nD - Z_brac=((nAr0+nAr1+nAr2+nAr3)*18**2+(nD+nD1))/nf & + Z_brac=((nAr0+nAr1+nAr2+nAr3+nAr4)*18**2+(nD+nD1))/nf & *(CLogei/CLogee) E_CH=nf*C_E**3*CLogee/(4*C_PI*C_E0**2*C_ME*C_C**2) @@ -1968,8 +1968,7 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) !! mpi error indicator REAL(rp) :: dgmin,dgmax,deta LOGICAL :: accepted - REAL(rp) :: EPHI,fRE_out,nAr0,nAr1,nAr2,nAr3,nD,nD1,ne,Te,Zeff,nRE - + REAL(rp) :: EPHI,fRE_out,nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,ne,Te,Zeff,nRE nsamples = spp%ppp*params%mpi_params%nmpi @@ -2072,7 +2071,7 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) !write(6,*) 'G_buffer',G_buffer call interp_nRE(params,R_buffer,0._rp,Z_buffer,psi0,EPHI,ne,Te,nRE, & - nAr0,nAr1,nAr2,nAr3,nD,nD1,G_buffer,fRE_out, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,G_buffer,fRE_out, & rho1D=h_params%rho_axis(h_params%rho_ind)) !write(6,*) 'after first interp_nRE' @@ -2083,7 +2082,7 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) f0=nRE*fRE_out* & fRE_H_pitch(params,eta_buffer,G_buffer,EPHI,ne,Te, & - nAr0,nAr1,nAr2,nAr3,nD,nD1) + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1) !write(6,*) 'nRE',nRE*params%cpp%density !write(6,*) 'fRE_out',fRE_out @@ -2096,8 +2095,8 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) do while (ii .LE. 1000_idef) !if (modulo(ii,100).eq.0) then - ! write(output_unit_write,'("Burn: ",I10)') ii - ! write(6,'("Burn: ",I10)') ii + !write(output_unit_write,'("Burn: ",I10)') ii + !write(6,'("Burn: ",I10)') ii !end if ! Test that pitch angle and momentum are within chosen boundary @@ -2141,7 +2140,7 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) !write(6,*) 'G_test',G_test call interp_nRE(params,R_test,0._rp,Z_test,psi1,EPHI,ne,Te,nRE, & - nAr0,nAr1,nAr2,nAr3,nD,nD1,G_test,fRE_out, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,G_test,fRE_out, & rho1D=h_params%rho_axis(h_params%rho_ind)) PSIN1=(psi1-PSIP0)/(PSIp_lim-PSIP0) @@ -2156,7 +2155,7 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) f1=nRE*fRE_out* & fRE_H_pitch(params,eta_test,G_test,EPHI,ne,Te, & - nAr0,nAr1,nAr2,nAr3,nD,nD1) + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1) ! f1=fRE_H(eta_test,G_test) !write(6,*) 'nRE',nRE*params%cpp%density @@ -2222,9 +2221,9 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) ! write(output_unit_write,'("sample:",I15)') ii - if (modulo(ii,nsamples/100).eq.0) then + if (modulo(ii,nsamples/10).eq.0) then write(output_unit_write,'("Sample: ",I10)') ii - write(6,'("Sample: ",I10)') ii + !write(6,'("Sample: ",I10)') ii end if ! Test that pitch angle and momentum are within chosen boundary @@ -2262,7 +2261,7 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) ! spp%sigmaZ,theta_rad) call interp_nRE(params,R_test,0._rp,Z_test,psi1,EPHI,ne,Te,nRE, & - nAr0,nAr1,nAr2,nAr3,nD,nD1,g_test,fRE_out, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,g_test,fRE_out, & rho1D=h_params%rho_axis(h_params%rho_ind)) PSIN1=(psi1-PSIP0)/(PSIp_lim-PSIP0) @@ -2280,7 +2279,7 @@ subroutine sample_Hollmann_distribution_1Dtransport(params,random,spp,F) f1=nRE*fRE_out* & fRE_H_pitch(params,eta_test,G_test,EPHI,ne,Te, & - nAr0,nAr1,nAr2,nAr3,nD,nD1) + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1) ! f1=fRE_H(eta_test,G_test) ! ratio = indicator_exp(PSIN,psi_max_buff)* & diff --git a/src/korc_fields.f90 b/src/korc_fields.f90 index bc4cec50..6a4439db 100755 --- a/src/korc_fields.f90 +++ b/src/korc_fields.f90 @@ -1194,6 +1194,7 @@ subroutine unitVectors(params,Xo,F,b1,b2,b3,flag,cart,hint,Bo) DEALLOCATE( vars%curlb ) DEALLOCATE( vars%E ) DEALLOCATE( vars%flagCon ) + DEALLOCATE( vars%initLCFS ) #ifdef FIO DEALLOCATE( vars%hint) #endif @@ -1907,7 +1908,7 @@ subroutine which_fields_in_file(params,Bfield,Efield,Bflux,B1field, & gname = "BR" call h5lexists_f(h5file_id,TRIM(gname),Bfield,h5error) - gname = "ER" + gname = "EPHI" call h5lexists_f(h5file_id,TRIM(gname),Efield,h5error) gname = "PSIp" @@ -1950,7 +1951,7 @@ subroutine load_field_data_from_hdf5(params,F) INTEGER(HID_T) :: group_id INTEGER(HID_T) :: subgroup_id INTEGER :: h5error - LOGICAL :: Efield + LOGICAL :: Efield,Bflux,Bflux3D filename = TRIM(params%magnetic_field_filename) call h5fopen_f(filename, H5F_ACC_RDONLY_F, h5file_id, h5error) @@ -2054,9 +2055,9 @@ subroutine load_field_data_from_hdf5(params,F) dset = "/PSIp" gname = 'PSIp' - call h5lexists_f(h5file_id,TRIM(gname),Efield,h5error) + call h5lexists_f(h5file_id,TRIM(gname),Bflux,h5error) - if (Efield) then + if (Bflux) then call load_array_from_hdf5(h5file_id,dset,F%PSIp) else F%PSIp = 0.0_rp @@ -2070,9 +2071,9 @@ subroutine load_field_data_from_hdf5(params,F) dset = "/PSIp" gname = 'PSIp' - call h5lexists_f(h5file_id,TRIM(gname),Efield,h5error) + call h5lexists_f(h5file_id,TRIM(gname),Bflux3D,h5error) - if (Efield) then + if (Bflux3D) then call load_array_from_hdf5(h5file_id,dset,F%PSIp3D) else F%PSIp3D = 0.0_rp diff --git a/src/korc_hpc.f90 b/src/korc_hpc.f90 index 52b3db4a..ee471b5c 100755 --- a/src/korc_hpc.f90 +++ b/src/korc_hpc.f90 @@ -45,6 +45,7 @@ subroutine korc_abort(errorcode) !! 25: korc_ppusher:adv_GCinterp_psiwE_top !! 26: korc_collisions:define_collisions_time_step !! 27: main:initialize_fields_interpolant + !! 28: korc_HDF5:save_restart_variables flush(output_unit_write) @@ -332,7 +333,7 @@ subroutine initialize_communications(params) !$OMP PARALLEL !$OMP MASTER write(output_unit_write,'(/,"OMP threads per MPI process: ",I3)') params%num_omp_threads - write(output_unit_write,'(/,"Cores available per MPI process: ",I3)') get_thread_number() + write(output_unit_write,'(/,"Cores available per MPI process: ",I3)') get_max_threads() !$OMP END MASTER !$OMP END PARALLEL #ifdef GNU @@ -352,9 +353,9 @@ function get_num_threads() integer :: get_num_threads #if defined(_OPENMP) -!$OMP PARALLEL +!$omp PARALLEL get_num_threads = OMP_GET_NUM_THREADS() -!$OMP END PARALLEL +!$omp END PARALLEL #else get_num_threads = 1 #endif diff --git a/src/korc_input.f90 b/src/korc_input.f90 index ba871d1b..5c9c9bc8 100755 --- a/src/korc_input.f90 +++ b/src/korc_input.f90 @@ -201,7 +201,7 @@ module korc_input LOGICAL :: Dim2x1t =.FALSE. LOGICAL :: E_2x1t = .FALSE. REAL(rp) :: t0_2x1t = 1.405 - INTEGER :: ind0_2x1t = 11 + INTEGER :: ind0_2x1t = 1 LOGICAL :: ReInterp_2x1t = .FALSE. INTEGER :: res_double=0 INTEGER :: dim_1D=50 diff --git a/src/korc_interp.f90 b/src/korc_interp.f90 index f43b1de2..a39a980b 100755 --- a/src/korc_interp.f90 +++ b/src/korc_interp.f90 @@ -208,6 +208,7 @@ module korc_interp TYPE(EZspline2) :: nAr1 TYPE(EZspline2) :: nAr2 TYPE(EZspline2) :: nAr3 + TYPE(EZspline2) :: nAr4 TYPE(EZspline2) :: nD TYPE(EZspline2) :: nD1 @@ -1582,9 +1583,108 @@ subroutine check_if_in_fields_domain_p(pchunk,F,Y_R,Y_PHI,Y_Z,flag) end if end subroutine check_if_in_fields_domain_p -subroutine 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,flag) + subroutine 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,flag) + !$acc routine seq + !! @note Subrotuine that checks if particles in the simulation are within + !! the spatial domain where interpolants and fields are known. @endnote + !! External fields and interpolants can have different spatial domains where + !! they are defined. Therefore, it is necessary to + !! check if a given particle has left these spatial domains to + !! stop following it, otherwise this will cause an error in the simulation. + REAL(rp), INTENT(IN) :: Y_R,Y_PHI,Y_Z + INTEGER(is), INTENT(INOUT) :: flag + !! Flag that determines whether particles are followed in the + !! simulation (flag=1), or not (flag=0). + INTEGER :: IR + !! Variable used to localize the grid cell in the \((R,\phi,Z)\) + !! or \((R,Z)\) grid containing the fields data that corresponds + !! to the radial position of the particles. + INTEGER :: IPHI + !! Variable used to localize the grid cell in the \((R,\phi,Z)\) + !! or \((R,Z)\) grid containing the fields data that corresponds + !! to the azimuthal position of the particles. + INTEGER :: IZ + !! Variable used to localize the grid cell in the \((R,\phi,Z)\) + !! or \((R,Z)\) grid containing the fields data that corresponds + !! to the vertical position of the particles. + + REAL(rp) :: Rwall,lscale,xtmp,ytmp,ztmp + 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 + REAL(rp),INTENT(IN) :: circumradius,ntiles + LOGICAL,INTENT(IN) :: Dim2x1t,Analytic_D3D_IWL,useDiMES + TYPE(KORC_INTERPOLANT_DOMAIN),INTENT(IN) :: fields_domain_local + TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN) :: bfield_2d_local + + IR = INT(FLOOR((Y_R - fields_domain_local%Ro + & + 0.5_rp*fields_domain_local%DR)/fields_domain_local%DR) + 1.0_rp,idef) + IZ = INT(FLOOR((Y_Z + ABS(fields_domain_local%Zo) + & + 0.5_rp*fields_domain_local%DZ)/fields_domain_local%DZ) + 1.0_rp,idef) + + if ((IR.lt.1).or.(IZ.lt.1).or. & + (IR.GT.bfield_2d_local%NR).OR.(IZ.GT.bfield_2d_local%NZ).or. & + (fields_domain_local%FLAG2D(IR,IZ).NE.1_is)) then + + if (.not.Analytic_D3D_IWL) then + flag = 0_is + else if (Analytic_D3D_IWL) then + if ((IR.lt.floor(bfield_2d_local%NR/6._rp)).and. & + (IZ.gt.floor(bfield_2d_local%NZ/5._rp)).and. & + (IZ.lt.floor(4._rp*bfield_2d_local%NZ/5._rp))) then + + Rwall=circumradius*cos(C_PI/ntiles)/ & + (cos(modulo(Y_PHI,2*C_PI/ntiles)-C_PI/ntiles)) + + !write(6,*) 'Rc,nt',F%circumradius*lscale,F%ntiles + !write(6,*) 'Rwall',Rwall*lscale + !write(6,*) 'mod',modulo(Y_PHI,2*C_PI/F%ntiles) + + if (Y_R.lt.Rwall) flag = 0_is + + else + flag = 0_is + endif + endif + + end if + + if (useDiMES) THEN + + 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_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_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 + + if (DiMESrad.le.0._rp) THEN + return + end if + + !DiMESsurf=DiMESloc_cart(3)+(DiMESdims(2)/DiMESdims(1))*sqrt(DiMESrad) + DiMESsurf=(DiMESloc_cart(3)-(DiMESdims(1)-DiMESdims(2)))+sqrt(DiMESrad) + + if (ztmp.le.DiMESsurf) flag=0_is + + end if + endif !DiMES + + end subroutine check_if_in_fields_domain_2D_p_ACC + +subroutine check_if_in_fields_domain_2D_ACC(F,Y_R,Y_PHI,Y_Z,flag) !$acc routine seq !! @note Subrotuine that checks if particles in the simulation are within !! the spatial domain where interpolants and fields are known. @endnote @@ -1594,6 +1694,7 @@ subroutine check_if_in_fields_domain_2D_p_ACC(fields_domain_local,bfield_2d_loca !! stop following it, otherwise this will cause an error in the simulation. REAL(rp), INTENT(IN) :: Y_R,Y_PHI,Y_Z INTEGER(is), INTENT(INOUT) :: flag + TYPE(FIELDS), INTENT(IN) :: F !! Flag that determines whether particles are followed in the !! simulation (flag=1), or not (flag=0). INTEGER :: IR @@ -1611,32 +1712,26 @@ subroutine check_if_in_fields_domain_2D_p_ACC(fields_domain_local,bfield_2d_loca REAL(rp) :: Rwall,lscale,xtmp,ytmp,ztmp 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 - REAL(rp),INTENT(IN) :: circumradius,ntiles - LOGICAL,INTENT(IN) :: Dim2x1t,Analytic_D3D_IWL,useDiMES - TYPE(KORC_INTERPOLANT_DOMAIN),INTENT(IN) :: fields_domain_local - TYPE(KORC_2D_FIELDS_INTERPOLANT),INTENT(IN) :: bfield_2d_local + REAL(rp),DIMENSION(3) :: DiMESloc_cart,DiMESloc_cyl - IR = INT(FLOOR((Y_R - fields_domain_local%Ro + & - 0.5_rp*fields_domain_local%DR)/fields_domain_local%DR) + 1.0_rp,idef) - IZ = INT(FLOOR((Y_Z + ABS(fields_domain_local%Zo) + & - 0.5_rp*fields_domain_local%DZ)/fields_domain_local%DZ) + 1.0_rp,idef) + IR = INT(FLOOR((Y_R - fields_domain%Ro + & + 0.5_rp*fields_domain%DR)/fields_domain%DR) + 1.0_rp,idef) + IZ = INT(FLOOR((Y_Z + ABS(fields_domain%Zo) + & + 0.5_rp*fields_domain%DZ)/fields_domain%DZ) + 1.0_rp,idef) if ((IR.lt.1).or.(IZ.lt.1).or. & - (IR.GT.bfield_2d_local%NR).OR.(IZ.GT.bfield_2d_local%NZ).or. & - (fields_domain_local%FLAG2D(IR,IZ).NE.1_is)) then + (IR.GT.bfield_2d%NR).OR.(IZ.GT.bfield_2d%NZ).or. & + (fields_domain%FLAG2D(IR,IZ).NE.1_is)) then - if (.not.Analytic_D3D_IWL) then + if (.not.F%Analytic_D3D_IWL) then flag = 0_is - else if (Analytic_D3D_IWL) then - if ((IR.lt.floor(bfield_2d_local%NR/6._rp)).and. & - (IZ.gt.floor(bfield_2d_local%NZ/5._rp)).and. & - (IZ.lt.floor(4._rp*bfield_2d_local%NZ/5._rp))) then + else if (F%Analytic_D3D_IWL) then + if ((IR.lt.floor(bfield_2d%NR/6._rp)).and. & + (IZ.gt.floor(bfield_2d%NZ/5._rp)).and. & + (IZ.lt.floor(4._rp*bfield_2d%NZ/5._rp))) then - Rwall=circumradius*cos(C_PI/ntiles)/ & - (cos(modulo(Y_PHI,2*C_PI/ntiles)-C_PI/ntiles)) + Rwall=F%circumradius*cos(C_PI/F%ntiles)/ & + (cos(modulo(Y_PHI,2*C_PI/F%ntiles)-C_PI/F%ntiles)) !write(6,*) 'Rc,nt',F%circumradius*lscale,F%ntiles !write(6,*) 'Rwall',Rwall*lscale @@ -1651,37 +1746,7 @@ subroutine check_if_in_fields_domain_2D_p_ACC(fields_domain_local,bfield_2d_loca end if - if (useDiMES) THEN - - 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_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_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 - - if (DiMESrad.le.0._rp) THEN - return - end if - - !DiMESsurf=DiMESloc_cart(3)+(DiMESdims(2)/DiMESdims(1))*sqrt(DiMESrad) - DiMESsurf=(DiMESloc_cart(3)-(DiMESdims(1)-DiMESdims(2)))+sqrt(DiMESrad) - - if (ztmp.le.DiMESsurf) flag=0_is - - end if - endif !DiMES - -end subroutine check_if_in_fields_domain_2D_p_ACC +end subroutine check_if_in_fields_domain_2D_ACC subroutine check_if_in_LCFS(F,Y,inLCFS) !! @note Subrotuine that checks if particles in the simulation are within @@ -1796,56 +1861,85 @@ subroutine initialize_profiles_interpolant(params,P) write(output_unit_write,'("* * * * INITIALIZING PROFILES INTERPOLANT * * * *")') end if - if (P%axisymmetric) then + if (P%axisymmetric.and.P%ReInterp_2x1t) then if (params%mpi_params%rank .EQ. 0) then write(output_unit_write,*) '2D ne, Te, Zeff' flush(output_unit_write) end if - profiles_2d%NR = P%dims(1) - profiles_2d%NZ = P%dims(3) + if (.not.(EZspline_allocated2(profiles_2d%ne))) then + profiles_2d%NR = P%dims(1) + profiles_2d%NZ = P%dims(3) - !write(output_unit_write,'("NR",I15)') profiles_2d%NR - !write(output_unit_write,'("NZ",I15)') profiles_2d%NR + !write(output_unit_write,'("NR",I15)') profiles_2d%NR + !write(output_unit_write,'("NZ",I15)') profiles_2d%NR - ! Initializing ne - call EZspline_init2(profiles_2d%ne,profiles_2d%NR,profiles_2d%NZ, & - profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) + ! Initializing ne + call EZspline_init2(profiles_2d%ne,profiles_2d%NR,profiles_2d%NZ, & + profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%ne%x1 = P%X%R + profiles_2d%ne%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%ne, P%ne_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) - profiles_2d%ne%x1 = P%X%R - profiles_2d%ne%x2 = P%X%Z + + !write(6,*) 'ne',P%ne_3D(12,P%ind_2x1t,15)*params%cpp%density + !write(6,*) 'ne_spline',profiles_2d%ne%fspl(1,12,15)*params%cpp%density - call EZspline_setup2(profiles_2d%ne, P%ne_2D, ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%ne, P%ne_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) + + ! Initializing Te - call EZspline_init2(profiles_2d%Te,profiles_2d%NR,profiles_2d%NZ, & - profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) + if (.not.(EZspline_allocated2(profiles_2d%Te))) then + call EZspline_init2(profiles_2d%Te,profiles_2d%NR,profiles_2d%NZ, & + profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) - profiles_2d%Te%x1 = P%X%R - profiles_2d%Te%x2 = P%X%Z + profiles_2d%Te%x1 = P%X%R + profiles_2d%Te%x2 = P%X%Z + endif !write(output_unit_write,'("Te_interp_R",E17.10)') profiles_2d%Te%x1 !write(output_unit_write,'("Te_interp_Z",E17.10)') profiles_2d%Te%x2 !write(output_unit_write,'("Te",E17.10)') P%Te_2D(10,:) - call EZspline_setup2(profiles_2d%Te, P%Te_2D, ezerr, .TRUE.) + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%Te, P%Te_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + + !write(6,*) 'Te',P%Te_3D(:,P%ind_2x1t,:)*params%cpp%temperature/C_E + !write(6,*) 'Te_spline',profiles_2d%Te%fspl(1,:,:)*params%cpp%temperature/C_E + + else + call EZspline_setup2(profiles_2d%Te, P%Te_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) ! Initializing Zeff - call EZspline_init2(profiles_2d%Zeff,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) + if (.not.(EZspline_allocated2(profiles_2d%Zeff))) then + call EZspline_init2(profiles_2d%Zeff,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) - profiles_2d%Zeff%x1 = P%X%R - profiles_2d%Zeff%x2 = P%X%Z + profiles_2d%Zeff%x1 = P%X%R + profiles_2d%Zeff%x2 = P%X%Z + endif - call EZspline_setup2(profiles_2d%Zeff, P%Zeff_2D, ezerr, .TRUE.) + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%Zeff, P%Zeff_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%Zeff, P%Zeff_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) if (params%profile_model(10:10) .EQ. 'H') then @@ -1856,102 +1950,170 @@ subroutine initialize_profiles_interpolant(params,P) end if ! Initializing RHON - call EZspline_init2(profiles_2d%RHON,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) - - profiles_2d%RHON%x1 = P%X%R - profiles_2d%RHON%x2 = P%X%Z - - call EZspline_setup2(profiles_2d%RHON, P%RHON, ezerr, .TRUE.) + if (.not.(EZspline_allocated2(profiles_2d%RHON))) then + call EZspline_init2(profiles_2d%RHON,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%RHON%x1 = P%X%R + profiles_2d%RHON%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%RHON, P%RHON_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%RHON, P%RHON_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) !write(output_unit_write,'("profiles_2d%RHON: ",E17.10)') profiles_2d%RHON%fspl(1,:,:) ! Initializing nRE - call EZspline_init2(profiles_2d%nRE,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) + if (.not.(EZspline_allocated2(profiles_2d%nRE))) then + call EZspline_init2(profiles_2d%nRE,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) - profiles_2d%nRE%x1 = P%X%R - profiles_2d%nRE%x2 = P%X%Z + profiles_2d%nRE%x1 = P%X%R + profiles_2d%nRE%x2 = P%X%Z - call EZspline_setup2(profiles_2d%nRE, P%nRE_2D, ezerr, .TRUE.) - call EZspline_error(ezerr) + call EZspline_setup2(profiles_2d%nRE, P%nRE_2D, ezerr, .TRUE.) + call EZspline_error(ezerr) + endif ! Initializing nAr0 - call EZspline_init2(profiles_2d%nAr0,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) - - profiles_2d%nAr0%x1 = P%X%R - profiles_2d%nAr0%x2 = P%X%Z - - call EZspline_setup2(profiles_2d%nAr0, P%nAr0_2D, ezerr, .TRUE.) + if (.not.(EZspline_allocated2(profiles_2d%nAr0))) then + call EZspline_init2(profiles_2d%nAr0,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%nAr0%x1 = P%X%R + profiles_2d%nAr0%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%nAr0, P%nAr0_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%nAr0, P%nAr0_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) ! Initializing nAr1 - call EZspline_init2(profiles_2d%nAr1,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) - - profiles_2d%nAr1%x1 = P%X%R - profiles_2d%nAr1%x2 = P%X%Z - - call EZspline_setup2(profiles_2d%nAr1, P%nAr1_2D, ezerr, .TRUE.) + if (.not.(EZspline_allocated2(profiles_2d%nAr1))) then + call EZspline_init2(profiles_2d%nAr1,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%nAr1%x1 = P%X%R + profiles_2d%nAr1%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%nAr1, P%nAr1_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%nAr1, P%nAr1_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) ! Initializing nAr2 - call EZspline_init2(profiles_2d%nAr2,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) - - profiles_2d%nAr2%x1 = P%X%R - profiles_2d%nAr2%x2 = P%X%Z - - call EZspline_setup2(profiles_2d%nAr2, P%nAr2_2D, ezerr, .TRUE.) + if (.not.(EZspline_allocated2(profiles_2d%nAr2))) then + call EZspline_init2(profiles_2d%nAr2,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%nAr2%x1 = P%X%R + profiles_2d%nAr2%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%nAr2, P%nAr2_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%nAr2, P%nAr2_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) ! Initializing nAr3 - call EZspline_init2(profiles_2d%nAr3,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + if (.not.(EZspline_allocated2(profiles_2d%nAr3))) then + call EZspline_init2(profiles_2d%nAr3,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%nAr3%x1 = P%X%R + profiles_2d%nAr3%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%nAr3, P%nAr3_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%nAr3, P%nAr3_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) - profiles_2d%nAr3%x1 = P%X%R - profiles_2d%nAr3%x2 = P%X%Z - - call EZspline_setup2(profiles_2d%nAr3, P%nAr3_2D, ezerr, .TRUE.) + ! Initializing nAr4 + if (.not.(EZspline_allocated2(profiles_2d%nAr4))) then + call EZspline_init2(profiles_2d%nAr4,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%nAr4%x1 = P%X%R + profiles_2d%nAr4%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%nAr4, P%nAr4_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%nAr4, P%nAr4_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) ! Initializing nD - call EZspline_init2(profiles_2d%nD,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) - - profiles_2d%nD%x1 = P%X%R - profiles_2d%nD%x2 = P%X%Z - - call EZspline_setup2(profiles_2d%nD, P%nD_2D, ezerr, .TRUE.) + if (.not.(EZspline_allocated2(profiles_2d%nD))) then + call EZspline_init2(profiles_2d%nD,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%nD%x1 = P%X%R + profiles_2d%nD%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%nD, P%nD_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%nD, P%nD_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) ! Initializing nD1 - call EZspline_init2(profiles_2d%nD1,profiles_2d%NR, & - profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) - call EZspline_error(ezerr) - - profiles_2d%nD1%x1 = P%X%R - profiles_2d%nD1%x2 = P%X%Z - - call EZspline_setup2(profiles_2d%nD1, P%nD1_2D, ezerr, .TRUE.) + if (.not.(EZspline_allocated2(profiles_2d%nD1))) then + call EZspline_init2(profiles_2d%nD1,profiles_2d%NR, & + profiles_2d%NZ,profiles_2d%BCSR,profiles_2d%BCSZ,ezerr) + call EZspline_error(ezerr) + + profiles_2d%nD1%x1 = P%X%R + profiles_2d%nD1%x2 = P%X%Z + endif + + if (P%ReInterp_2x1t) then + call EZspline_setup2(profiles_2d%nD1, P%nD1_3D(:,P%ind_2x1t,:), ezerr, .TRUE.) + else + call EZspline_setup2(profiles_2d%nD1, P%nD1_2D, ezerr, .TRUE.) + endif call EZspline_error(ezerr) end if - ALLOCATE(profiles_domain%FLAG2D(profiles_2d%NR,profiles_2d%NZ)) - profiles_domain%FLAG2D = P%FLAG2D - - profiles_domain%DR = ABS(P%X%R(2) - P%X%R(1)) - profiles_domain%DZ = ABS(P%X%Z(2) - P%X%Z(1)) + if (.not.ALLOCATED(profiles_domain%FLAG2D)) then + ALLOCATE(profiles_domain%FLAG2D(profiles_2d%NR,profiles_2d%NZ)) + + if (P%ReInterp_2x1t) then + profiles_domain%FLAG2D = P%FLAG3D(:,1,:) + else + profiles_domain%FLAG2D = P%FLAG2D + endif + + profiles_domain%DR = ABS(P%X%R(2) - P%X%R(1)) + profiles_domain%DZ = ABS(P%X%Z(2) - P%X%Z(1)) + endif else if (params%mpi_params%rank .EQ. 0) then @@ -3050,7 +3212,7 @@ subroutine calculate_magnetic_field(params,Y,F,B,E,PSI_P,flag) A=0._rp if(F%Dim2x1t.and.(.not.F%ReInterp_2x1t)) then - !$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,ezerr) & + !$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,ezerr,Atmp) & !$OMP& SHARED(F,Y,A,B,flag,bfield_2X1T,PSI_P) do pp=1_idef,ss @@ -3119,7 +3281,7 @@ subroutine calculate_magnetic_field(params,Y,F,B,E,PSI_P,flag) else - !$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,ezerr) & + !$OMP PARALLEL DO FIRSTPRIVATE(ss) PRIVATE(pp,ezerr,Atmp) & !$OMP& SHARED(F,Y,A,B,flag,bfield_2d,PSI_P) do pp=1_idef,ss @@ -3170,10 +3332,14 @@ subroutine calculate_magnetic_field(params,Y,F,B,E,PSI_P,flag) B(pp,1) = A(pp,1)*COS(Y(pp,2)) - A(pp,2)*SIN(Y(pp,2)) B(pp,2) = A(pp,1)*SIN(Y(pp,2)) + A(pp,2)*COS(Y(pp,2)) B(pp,3) = A(pp,3) + + !write(6,*) 'R',Y(pp,1)*params%cpp%length,'Z',Y(pp,3)*params%cpp%length,'PSIP', & + ! PSI_P(pp)*(params%cpp%Bo*params%cpp%length**2) else B(pp,1) = A(pp,1) B(pp,2) = A(pp,2) B(pp,3) = A(pp,3) + end if @@ -3182,6 +3348,12 @@ subroutine calculate_magnetic_field(params,Y,F,B,E,PSI_P,flag) !$OMP END PARALLEL DO end if + !if (params%GC_coords) then + ! write(6,*) 'R',Y(:,1)*params%cpp%length,'Z',Y(:,3)*params%cpp%length,'PSIP', & + ! PSI_P(:)*(params%cpp%Bo*params%cpp%length**2),'BR',B(:,1)*params%cpp%Bo, & + ! 'BZ',B(:,3)*params%cpp%Bo + !endif + ! write(output_unit_write,'("calculate_fields")') ! write(output_unit_write,'("B_R: ",E17.10)') A(:,1) @@ -3256,9 +3428,9 @@ subroutine calculate_GCfieldswE_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI call check_if_in_fields_domain_p(pchunk,F,Y_R,Y_PHI,Y_Z,flag_cache) - !$OMP SIMD - ! !$OMP& aligned(PSIp,A,B_R,Y_R,B_PHI,B_Z,Bmag,gradB_R,gradB_PHI,gradB_Z, & - ! !$OMP& curlb_R,curlb_PHI,curlb_Z,E_R,E_PHI,E_Z) + !!$OMP SIMD + !! !$OMP& aligned(PSIp,A,B_R,Y_R,B_PHI,B_Z,Bmag,gradB_R,gradB_PHI,gradB_Z, & + !! !$OMP& curlb_R,curlb_PHI,curlb_Z,E_R,E_PHI,E_Z) do cc=1_idef,pchunk call EZspline_interp2_GCvarswE(bfield_2d%A, efield_2d%PHI, Y_R(cc), Y_Z(cc), A, & EPHI, ezerr) @@ -3311,6 +3483,75 @@ subroutine calculate_GCfieldswE_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI end subroutine calculate_GCfieldswE_p +subroutine calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,flag_cache,PSIp) + !$acc routine seq + REAL(rp), INTENT(IN) :: Y_R,Y_PHI,Y_Z + TYPE(FIELDS), INTENT(IN) :: F + REAL(rp), INTENT(OUT) :: B_R,B_PHI,B_Z + REAL(rp), INTENT(OUT) :: gradB_R,gradB_PHI,gradB_Z + REAL(rp), INTENT(OUT) :: curlb_R,curlb_PHI,curlb_Z + REAL(rp), INTENT(OUT) :: E_R,E_PHI,E_Z + REAL(rp) :: Bmag,EPHI + INTEGER :: cc + REAL(rp),INTENT(OUT) :: PSIp + REAL(rp), DIMENSION(6) :: A + INTEGER(is),INTENT(INOUT) :: flag_cache + REAL(rp) :: psip_conv + + !$acc routine (EZspline_interp2_GCvarswE) seq + !$acc routine (EZspline_error) seq + + psip_conv=F%psip_conv + + call EZspline_interp2_GCvarswE(bfield_2d%A, efield_2d%PHI, Y_R, Y_Z, A, & + EPHI, ezerr) + call EZspline_error(ezerr) + + !A(:,1) = PSIp + !A(:,2) = dPSIp/dR + !A(:,3) = dPSIp/dZ + !A(:,4) = d^2PSIp/dR^2 + !A(:,5) = d^2PSIp/dZ^2 + !A(:,6) = d^2PSIp/dRdZ + + PSIp=A(1) + + B_R = psip_conv*A(3)/Y_R + + ! BR = (dA/dZ)/R + B_PHI = -F%Bo*F%Ro/Y_R + ! BPHI = Fo*Ro/R + B_Z = -psip_conv*A(2)/Y_R + ! BR = -(dA/dR)/R + + Bmag=sqrt(B_R*B_R+B_PHI*B_PHI+B_Z*B_Z) + + gradB_R=(B_R*psip_conv*A(6)-B_Z*psip_conv*A(4)- & + Bmag*Bmag)/(Y_R*Bmag) + gradB_PHI=0._rp + gradB_Z=(B_R*psip_conv*A(5)-B_Z*psip_conv*A(6))/ & + (Y_R*Bmag) + + curlb_R=B_PHI*gradB_Z/(Bmag*Bmag) + curlb_PHI=(Bmag/Y_R*(B_Z+psip_conv*A(4)+ & + psip_conv*A(5))-B_R*gradB_Z+B_Z*gradB_R)/ & + (Bmag*Bmag) + curlb_Z=-B_PHI*gradB_R/(Bmag*Bmag) + + if (F%E_2x1t) then + E_R = 0._rp + E_PHI = EPHI + E_Z = 0._rp + else + E_R = 0._rp + E_PHI = 0._rp + E_Z = 0._rp + end if + + +end subroutine calculate_GCfieldswE_p_ACC + subroutine calculate_GCfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, & E_R,E_PHI,E_Z, & curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,flag_cache,PSIp) @@ -3722,39 +3963,71 @@ subroutine interp_fields(params,prtcls,F) end subroutine interp_fields #ifdef PSPLINE -subroutine interp_Hcollision_p(pchunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, & - nAr0,nAr1,nAr2,nAr3,nD,nD1,flag_cache) - INTEGER, INTENT(IN) :: pchunk - REAL(rp),DIMENSION(pchunk),INTENT(IN) :: Y_R,Y_PHI,Y_Z - REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: ne,Te,Zeff - REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: nAr0,nAr1,nAr2,nAr3,nD,nD1 - INTEGER(is),DIMENSION(pchunk),INTENT(INOUT) :: flag_cache - INTEGER :: cc +subroutine interp_Hcollision_p(params,pchunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,flag_cache) + TYPE(KORC_PARAMS), INTENT(IN) :: params + INTEGER, INTENT(IN) :: pchunk + REAL(rp),DIMENSION(pchunk),INTENT(IN) :: Y_R,Y_PHI,Y_Z + REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: ne,Te,Zeff + REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1 + INTEGER(is),DIMENSION(pchunk),INTENT(INOUT) :: flag_cache + INTEGER :: cc + + call check_if_in_profiles_domain_p(pchunk,Y_R,Y_PHI,Y_Z,flag_cache) + ! write(output_unit_write,'("YR: ",E17.10)') Y_R(1) + ! write(output_unit_write,'("YPHI: ",E17.10)') Y_PHI(1) + ! write(output_unit_write,'("YZ: ",E17.10)') Y_Z(1) + + ! write(output_unit_write,'("Te_interp_R",E17.10)') profiles_2d%Te%x1 + ! write(output_unit_write,'("Te_interp_Z",E17.10)') profiles_2d%Te%x2 + + do cc=1_idef,pchunk + call EZspline_interp2_Hcollision(profiles_2d%ne,profiles_2d%Te,profiles_2d%Zeff, & + profiles_2d%nAr0,profiles_2d%nAr1,profiles_2d%nAr2,profiles_2d%nAr3,profiles_2d%nAr4, & + profiles_2d%nD,profiles_2d%nD1,Y_R(cc),Y_Z(cc),ne(cc),Te(cc),Zeff(cc), & + nAr0(cc),nAr1(cc),nAr2(cc),nAr3(cc),nAr4(cc),nD(cc),nD1(cc),ezerr) + call EZspline_error(ezerr) + end do + + !write(6,*) 'interp ne%fspl',profiles_2d%ne%fspl(1,12,15)*params%cpp%density + + !write(6,*) 'interp ne',ne(1) + + end subroutine interp_Hcollision_p + +subroutine interp_Hcollision_p_ACC(params,Y_R,Y_PHI,Y_Z, & + ne,Te,Zeff,nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,flag_cache) + !$acc routine seq + TYPE(KORC_PARAMS), INTENT(IN) :: params + REAL(rp),INTENT(IN) :: Y_R,Y_PHI,Y_Z + REAL(rp),INTENT(OUT) :: ne,Te,Zeff + REAL(rp),INTENT(OUT) :: nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1 + INTEGER(is),INTENT(INOUT) :: flag_cache + - call check_if_in_profiles_domain_p(pchunk,Y_R,Y_PHI,Y_Z,flag_cache) - ! write(output_unit_write,'("YR: ",E17.10)') Y_R(1) - ! write(output_unit_write,'("YPHI: ",E17.10)') Y_PHI(1) - ! write(output_unit_write,'("YZ: ",E17.10)') Y_Z(1) + !$acc routine (EZspline_interp2_Hcollision) seq + !$acc routine (EZspline_error) seq - ! write(output_unit_write,'("Te_interp_R",E17.10)') profiles_2d%Te%x1 - ! write(output_unit_write,'("Te_interp_Z",E17.10)') profiles_2d%Te%x2 - do cc=1_idef,pchunk - call EZspline_interp2_collision(profiles_2d%ne,profiles_2d%Te,profiles_2d%Zeff, & - profiles_2d%nAr0,profiles_2d%nAr1,profiles_2d%nAr2,profiles_2d%nAr3, & - profiles_2d%nD,profiles_2d%nD1,Y_R(cc),Y_Z(cc),ne(cc),Te(cc),Zeff(cc), & - nAr0(cc),nAr1(cc),nAr2(cc),nAr3(cc),nD(cc),nD1(cc),ezerr) + call EZspline_interp2_Hcollision(profiles_2d%ne,profiles_2d%Te,profiles_2d%Zeff, & + profiles_2d%nAr0,profiles_2d%nAr1,profiles_2d%nAr2,profiles_2d%nAr3,profiles_2d%nAr4, & + profiles_2d%nD,profiles_2d%nD1,Y_R,Y_Z,ne,Te,Zeff, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,ezerr) call EZspline_error(ezerr) - end do -end subroutine interp_Hcollision_p + + !write(6,*) 'interp ne%fspl',profiles_2d%ne%fspl(1,12,15)*params%cpp%density + + !write(6,*) 'interp ne',ne(1) + +end subroutine interp_Hcollision_p_ACC subroutine interp_nRE(params,Y_R,Y_PHI,Y_Z,PSIp,EPHI,ne,Te,nRE, & - nAr0,nAr1,nAr2,nAr3,nD,nD1,g_test,fRE_out,rho1D) + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,g_test,fRE_out,rho1D) TYPE(KORC_PARAMS), INTENT(IN) :: params REAL(rp),INTENT(IN) :: Y_R,Y_PHI,Y_Z,g_test REAL(rp),INTENT(OUT) :: PSIp,EPHI,ne,Te,nRE,fRE_out - REAL(rp),INTENT(OUT) :: nAr0,nAr1,nAr2,nAr3,nD,nD1 + REAL(rp),INTENT(OUT) :: nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1 REAL(rp) :: RHON,dummy INTEGER(is),DIMENSION(1) :: flag_cache REAL(rp), INTENT(IN),optional :: rho1D @@ -3768,11 +4041,13 @@ subroutine interp_nRE(params,Y_R,Y_PHI,Y_Z,PSIp,EPHI,ne,Te,nRE, & call EZspline_interp2_FOaorsa(bfield_2d%A,efield_2d%PHI, & profiles_2d%ne,profiles_2d%Te, & profiles_2d%nRE,profiles_2d%nAr0,profiles_2d%nAr1, & - profiles_2d%nAr2,profiles_2d%nAr3,profiles_2d%nD,profiles_2d%nD1, & - profiles_2d%RHON,bfield_2d%A,Y_R,Y_Z, & - A,EPHI,ne,Te,nRE,nAr0,nAr1,nAr2,nAr3,nD,nD1,RHON,dummy,ezerr) + profiles_2d%nAr2,profiles_2d%nAr3,profiles_2d%nAr4,profiles_2d%nD,profiles_2d%nD1, & + profiles_2d%RHON,Y_R,Y_Z, & + A,EPHI,ne,Te,nRE,nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,RHON,ezerr) call EZspline_error(ezerr) + !write(6,*) 'EPHI',EPHI + PSIp=A(1) !write(6,*) 'RHON',RHON diff --git a/src/korc_ppusher.f90 b/src/korc_ppusher.f90 index 1903a09d..042f536d 100755 --- a/src/korc_ppusher.f90 +++ b/src/korc_ppusher.f90 @@ -9,10 +9,7 @@ module korc_ppusher use korc_interp use korc_collisions use korc_hpc - -#ifdef PARALLEL_RANDOM use korc_random -#endif IMPLICIT NONE @@ -4139,8 +4136,371 @@ subroutine advance_FP3Dinterp_vars(params,random,X_X,X_Y,X_Z,V_X,V_Y,V_Z,g, & end subroutine advance_FP3Dinterp_vars +subroutine GC_init(params,random,F,P,spp) + !! @note Subroutine to advance GC variables \(({\bf X},p_\parallel)\) + !! @endnote + !! Comment this section further with evolution equations, numerical + !! methods, and descriptions of both. + TYPE(KORC_PARAMS), INTENT(INOUT) :: params + CLASS(random_context), POINTER, INTENT(INOUT) :: random + !! Core KORC simulation parameters. + TYPE(FIELDS), INTENT(INOUT) :: F + !! An instance of the KORC derived type FIELDS. + TYPE(PROFILES), INTENT(IN) :: P + 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. + + INTEGER :: ii + !! Species iterator. + INTEGER :: pp + !! Particles iterator. + INTEGER :: cc,pchunk + !! Chunk iterator. + REAL(rp) :: Bmag1,pmag + REAL(rp) :: Bmagc + REAL(rp) :: rm + REAL(rp),DIMENSION(:,:),ALLOCATABLE :: RAphi + REAL(rp), DIMENSION(3) :: bhat + REAL(rp), DIMENSION(3) :: bhatc + REAL(rp), DIMENSION(params%pchunk) :: E_PHI + REAL(rp),DIMENSION(:),ALLOCATABLE :: RVphi + + REAL(rp),DIMENSION(params%pchunk) :: rm8,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,vpll,gam,ne,Te,Zeff,PSIp + REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species) :: nimp + INTEGER(is),DIMENSION(params%pchunk) :: flagCon,flagCol + real(rp),dimension(F%dim_1D) :: Vpart,Vpartave,VpartOMP + real(rp) :: dr,m_cache + integer :: rind + + ! write(output_unit_write,'("eta",E17.10)') spp(ii)%vars%eta(pp) + ! write(output_unit_write,'("gam",E17.10)') spp(ii)%vars%g(pp) + + do ii = 1_idef,params%num_species + + pchunk=params%pchunk + + if (spp(ii)%spatial_distribution.eq.'TRACER'.and. & + params%FO_GC_compare) then + call get_fields(params,spp(ii)%vars,F) + !! Calls [[get_fields]] in [[korc_fields]]. + ! Interpolates fields at local particles' position and keeps in + ! spp%vars. Fields in (R,\(\phi\),Z) coordinates. + + ALLOCATE(RAphi(spp(ii)%ppp,2)) + ALLOCATE(RVphi(spp(ii)%ppp)) + RAphi=0.0_rp + + call cart_to_cyl(spp(ii)%vars%X,spp(ii)%vars%Y) + + !$OMP PARALLEL DO SHARED(params,ii,spp,F,RAphi,RVphi) & + !$OMP& PRIVATE(pp,Bmag1,bhat,rm) + ! Call OpenMP to calculate p_par and mu for each particle and + ! put into spp%vars%V + do pp=1_idef,spp(ii)%ppp + if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then + + RVphi(pp)=(-sin(spp(ii)%vars%Y(pp,2))*spp(ii)%vars%V(pp,1)+ & + cos(spp(ii)%vars%Y(pp,2))*spp(ii)%vars%V(pp,2))* & + spp(ii)%vars%Y(pp,1) + + Bmag1 = SQRT(spp(ii)%vars%B(pp,1)*spp(ii)%vars%B(pp,1)+ & + spp(ii)%vars%B(pp,2)*spp(ii)%vars%B(pp,2)+ & + spp(ii)%vars%B(pp,3)*spp(ii)%vars%B(pp,3)) + + ! write(output_unit_write,'("pp: ",I16)') pp + ! write(output_unit_write,'("Bmag: ",E17.10)') Bmag + + + bhat = spp(ii)%vars%B(pp,:)/Bmag1 + + if (params%field_model(1:10).eq.'ANALYTICAL') then + rm=sqrt((spp(ii)%vars%Y(pp,1)-F%AB%Ro)**2+ & + (spp(ii)%vars%Y(pp,3))**2) + + RAphi(pp,1)=-F%AB%lambda**2*F%AB%Bo/(2*F%AB%qo)* & + log(1+(rm/F%AB%lambda)**2) + + else if (params%field_model(1:8).eq.'EXTERNAL') then + + RAphi(pp,1)=spp(ii)%vars%PSI_P(pp)/(2*C_PI) + + end if + + ! write(output_unit_write,'("bhat: ",E17.10)') bhat + ! write(output_unit_write,'("V: ",E17.10)') spp(ii)%vars%V(pp,:) + + + spp(ii)%vars%X(pp,:)=spp(ii)%vars%X(pp,:)- & + spp(ii)%m*spp(ii)%vars%g(pp)* & + cross(bhat,spp(ii)%vars%V(pp,:))/(spp(ii)%q*Bmag1) + + ! transforming from particle location to associated + ! GC location + + end if ! if particle in domain, i.e. spp%vars%flagCon==1 + end do ! loop over particles on an mpi process + !$OMP END PARALLEL DO + + call cart_to_cyl(spp(ii)%vars%X,spp(ii)%vars%Y) + call get_fields(params,spp(ii)%vars,F) + + !$OMP PARALLEL DO SHARED(params,ii,spp,F,RAphi,RVphi) & + !$OMP& PRIVATE(pp,rm) + ! Call OpenMP to calculate p_par and mu for each particle and + ! put into spp%vars%V + do pp=1_idef,spp(ii)%ppp + if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then + + if (params%field_model(1:10).eq.'ANALYTICAL') then + rm=sqrt((spp(ii)%vars%Y(pp,1)-F%AB%Ro)**2+ & + (spp(ii)%vars%Y(pp,3))**2) + RAphi(pp,2)=-F%AB%lambda**2*F%AB%Bo/(2*F%AB%qo)* & + log(1+(rm/F%AB%lambda)**2) + + else if (params%field_model(1:8).eq.'EXTERNAL') then + + RAphi(pp,2)=spp(ii)%vars%PSI_P(pp)/(2*C_PI) + + end if + + write(output_unit_write,'("RAphi1: ",E17.10)') RAphi(pp,1) + write(output_unit_write,'("RAphi2: ",E17.10)') RAphi(pp,2) + + spp(ii)%vars%V(pp,1)=(spp(ii)%m*spp(ii)%vars%g(pp)* & + RVphi(pp)+spp(ii)%q*(RAphi(pp,1)-RAphi(pp,2)))/ & + spp(ii)%vars%Y(pp,1) + !GC ppar + + end if ! if particle in domain, i.e. spp%vars%flagCon==1 + end do ! loop over particles on an mpi process + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO SHARED(ii,spp) PRIVATE(pp,Bmagc,bhatc) + ! Call OpenMP to calculate p_par and mu for each particle and + ! put into spp%vars%V + do pp=1_idef,spp(ii)%ppp + if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then + + Bmagc = SQRT( DOT_PRODUCT(spp(ii)%vars%B(pp,:), & + spp(ii)%vars%B(pp,:))) + + bhatc = spp(ii)%vars%B(pp,:)/Bmagc + + spp(ii)%vars%V(pp,1)=spp(ii)%vars%V(pp,1)/ & + bhatc(2) + !GC ppar + + spp(ii)%vars%V(pp,2)=spp(ii)%m/(2*Bmagc)* & + (spp(ii)%vars%g(pp)**2- & + (1+(spp(ii)%vars%V(pp,1)/spp(ii)%m)**2)) + !GC mu + + + end if ! if particle in domain, i.e. spp%vars%flagCon==1 + end do ! loop over particles on an mpi process + !$OMP END PARALLEL DO + + params%GC_coords=.TRUE. + DEALLOCATE(RAphi) + DEALLOCATE(RVphi) + + !Preparing Output Data + call get_fields(params,spp(ii)%vars,F) + + !$OMP PARALLEL DO shared(F,params,spp) & + !$OMP& PRIVATE(cc,pp,E_PHI,Y_R) firstprivate(pchunk) + do pp=1_idef,spp(ii)%ppp,pchunk + + !$OMP SIMD + do cc=1_idef,pchunk + E_PHI(cc)=spp(ii)%vars%E(pp-1+cc,2) + Y_R(cc)=spp(ii)%vars%Y(pp-1+cc,1) + end do + !$OMP END SIMD + + call add_analytical_E_p(params,int8(0),F,E_PHI,Y_R,Y_Z) + + + !$OMP SIMD + do cc=1_idef,pchunk + spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) + end do + !$OMP END SIMD + + end do + !$OMP END PARALLEL DO + + + !$OMP PARALLEL DO SHARED(ii,spp) PRIVATE(pp,Bmag1) + ! Call OpenMP to calculate p_par and mu for each particle and + ! put into spp%vars%V + do pp=1_idef,spp(ii)%ppp + if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then + + Bmag1 = SQRT( DOT_PRODUCT(spp(ii)%vars%B(pp,:), & + spp(ii)%vars%B(pp,:))) + + spp(ii)%vars%g(pp)=sqrt(1+(spp(ii)%vars%V(pp,1))**2+ & + 2*spp(ii)%vars%V(pp,2)*Bmag1) + + ! write(output_unit_write,'("Bmag:",E17.10)') Bmag1 + ! write(output_unit_write,'("PPLL:",E17.10)') spp(ii)%vars%V(pp,1) + ! write(output_unit_write,'("MU:",E17.10)') spp(ii)%vars%V(pp,2) + + spp(ii)%vars%eta(pp) = atan2(sqrt(2*spp(ii)%m*Bmag1* & + spp(ii)%vars%V(pp,2)),spp(ii)%vars%V(pp,1))*180.0_rp/C_PI + + ! write(output_unit_write,'("BR",E17.10)') spp(ii)%vars%B(pp,1) + ! write(output_unit_write,'("BPHI",E17.10)') spp(ii)%vars%B(pp,2) + ! write(output_unit_write,'("BZ",E17.10)') spp(ii)%vars%B(pp,3) + + ! write(output_unit_write,'("ppll",E17.10)') spp(ii)%vars%V(pp,1) + ! write(output_unit_write,'("pperp",E17.10)') sqrt(2*spp(ii)%m*Bmag1* & + ! spp(ii)%vars%V(pp,2)) + + ! write(output_unit_write,'("eta GCinit",E17.10)') spp(ii)%vars%eta(pp) + ! write(output_unit_write,'("gam",E17.10)') spp(ii)%vars%g(pp) + + + end if ! if particle in domain, i.e. spp%vars%flagCon==1 + end do ! loop over particles on an mpi process + !$OMP END PARALLEL DO + else + + if ((spp(ii)%spatial_distribution.eq.'TRACER').or. & + (spp(ii)%spatial_distribution.eq.'TORUS').or. & + (spp(ii)%spatial_distribution.eq.'DISK').or. & + (spp(ii)%spatial_distribution.eq. & + '2D-GAUSSIAN-ELLIPTIC-TORUS-MH')) & + call cart_to_cyl(spp(ii)%vars%X,spp(ii)%vars%Y) + + params%GC_coords=.TRUE. + + do pp=1_idef,spp(ii)%ppp + spp(ii)%vars%E(pp,1)=0._rp + spp(ii)%vars%E(pp,2)=0._rp + spp(ii)%vars%E(pp,3)=0._rp + end do + + + !write(6,*) 'before second get fields' + call get_fields(params,spp(ii)%vars,F) + !write(6,*) 'after second get fields' + + !write(6,*) 'R',spp(ii)%vars%Y(:,1)*params%cpp%length,'Z',spp(ii)%vars%Y(:,3)*params%cpp%length + !write(6,*) 'BR',spp(ii)%vars%B(:,1)*params%cpp%Bo,'BZ',spp(ii)%vars%B(:,3)*params%cpp%Bo + + !write(output_unit_write,*) spp(1)%vars%PSI_P + + !$OMP PARALLEL DO SHARED(ii,spp) PRIVATE(pp,Bmag1) + + do pp=1_idef,spp(ii)%ppp + + ! if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then + + ! write(output_unit_write,'("BR: ",E17.10)') spp(ii)%vars%B(pp,1) + ! write(output_unit_write,'("BPHI: ",E17.10)') spp(ii)%vars%B(pp,2) + ! write(output_unit_write,'("BZ: ",E17.10)') spp(ii)%vars%B(pp,3) + + Bmag1 = SQRT( DOT_PRODUCT(spp(ii)%vars%B(pp,:), & + spp(ii)%vars%B(pp,:))) + + pmag=sqrt(spp(ii)%vars%g(pp)**2-1) + + spp(ii)%vars%V(pp,1)=pmag*cos(deg2rad(spp(ii)%vars%eta(pp))) + + !write(6,*) 'GC_init' + !write(6,*) spp(ii)%m,Bmag1 + + spp(ii)%vars%V(pp,2)=(pmag* & + sin(deg2rad(spp(ii)%vars%eta(pp))))**2/ & + (2*spp(ii)%m*Bmag1) + + ! write(output_unit_write,'("BR",E17.10)') spp(ii)%vars%B(pp,1) + ! write(output_unit_write,'("BPHI",E17.10)') spp(ii)%vars%B(pp,2) + ! write(output_unit_write,'("BZ",E17.10)') spp(ii)%vars%B(pp,3) + + !write(output_unit_write,'("ppll",E17.10)') spp(ii)%vars%V(pp,1) + !write(output_unit_write,'("mu",E17.10)') spp(ii)%vars%V(pp,2) + + ! write(output_unit_write,'("eta",E17.10)') spp(ii)%vars%eta(pp) + ! write(output_unit_write,'("gam",E17.10)') spp(ii)%vars%g(pp) + + + ! end if ! if particle in domain, i.e. spp%vars%flagCon==1 + end do ! loop over particles on an mpi process + !$OMP END PARALLEL DO + + m_cache=spp(ii)%m + + !$OMP PARALLEL DO shared(F,params,spp) & + !$OMP& PRIVATE(pp,cc,E_PHI,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,flagCon,flagCol,ne,Te,Zeff,nimp) & + !$OMP& firstprivate(pchunk,m_cache) + do pp=1_idef,spp(ii)%ppp,pchunk + + !$OMP SIMD + do cc=1_idef,pchunk + E_PHI(cc)=spp(ii)%vars%E(pp-1+cc,2) + Y_R(cc)=spp(ii)%vars%Y(pp-1+cc,1) + Y_PHI(cc)=spp(ii)%vars%Y(pp-1+cc,2) + Y_Z(cc)=spp(ii)%vars%Y(pp-1+cc,3) + + V_PLL(cc)=spp(ii)%vars%V(pp-1+cc,1) + V_MU(cc)=spp(ii)%vars%V(pp-1+cc,2) + + flagCon(cc)=spp(ii)%vars%flagCon(pp-1+cc) + flagCol(cc)=spp(ii)%vars%flagCol(pp-1+cc) + end do + !$OMP END SIMD + + if (params%field_model(1:8).eq.'EXTERNAL') then + call add_analytical_E_p(params,int8(0),F,E_PHI,Y_R,Y_Z) + end if + + !$OMP SIMD + do cc=1_idef,pchunk + + spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) + end do + !$OMP END SIMD + + if (params%collisions) then + + call include_CoulombCollisions_GC_p(0_ip,params,random, & + Y_R,Y_PHI,Y_Z, V_PLL,V_MU,m_cache, & + flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) + + !$OMP SIMD + do cc=1_idef,pchunk + spp(ii)%vars%ne(pp-1+cc) = ne(cc) + spp(ii)%vars%Te(pp-1+cc) = Te(cc) + spp(ii)%vars%Zeff(pp-1+cc) = Zeff(cc) + spp(ii)%vars%nimp(pp-1+cc,:) = nimp(cc,:) + end do + !$OMP END SIMD + + end if + + end do + !$OMP END PARALLEL DO + + + + end if + + spp(ii)%vars%Yborn=spp(ii)%vars%Y + + end do ! loop over particle species + + !write(6,*) 'R',spp(1)%vars%Y(1,1)*params%cpp%length + !write(6,*) 'Z',spp(1)%vars%Y(1,3)*params%cpp%length + !write(6,*) 'ppusher ne',spp(1)%vars%ne(1)*params%cpp%density + + + end subroutine GC_init -subroutine GC_init(params,F,spp) +subroutine GC_init_ACC(params,F,P,spp) !! @note Subroutine to advance GC variables \(({\bf X},p_\parallel)\) !! @endnote !! Comment this section further with evolution equations, numerical @@ -4149,7 +4509,7 @@ subroutine GC_init(params,F,spp) !! Core KORC simulation parameters. TYPE(FIELDS), INTENT(INOUT) :: F !! An instance of the KORC derived type FIELDS. - + TYPE(PROFILES), INTENT(IN) :: P 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. @@ -4158,309 +4518,100 @@ subroutine GC_init(params,F,spp) !! Species iterator. INTEGER :: pp !! Particles iterator. - INTEGER :: cc,pchunk - !! Chunk iterator. - REAL(rp) :: Bmag1,pmag - REAL(rp) :: Bmagc - REAL(rp) :: rm - REAL(rp),DIMENSION(:,:),ALLOCATABLE :: RAphi - REAL(rp), DIMENSION(3) :: bhat - REAL(rp), DIMENSION(3) :: bhatc - REAL(rp), DIMENSION(params%pchunk) :: E_PHI - REAL(rp),DIMENSION(:),ALLOCATABLE :: RVphi - - REAL(rp),DIMENSION(params%pchunk) :: rm8,Y_R,Y_Z,V_PLL,vpll,gam - real(rp),dimension(F%dim_1D) :: Vpart,Vpartave,VpartOMP - real(rp) :: dr - integer :: rind - - ! write(output_unit_write,'("eta",E17.10)') spp(ii)%vars%eta(pp) - ! write(output_unit_write,'("gam",E17.10)') spp(ii)%vars%g(pp) + REAL(rp) :: Bmag,pmag + REAL(rp) :: Y_R,Y_PHI,Y_Z,V_PLL,V_MU,vpll,gam,ne,Te,Zeff,PSIp + REAL(rp) :: B_R,B_PHI,B_Z,E_R,E_PHI,E_Z + REAL(rp) :: curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z + REAL(rp) :: nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1 + INTEGER(is) :: flagCon,flagCol + real(rp) :: m_cache + + !$acc routine (calculate_GCfieldswE_p_ACC) seq + !$acc routine (cart_to_cyl_p_ACC) seq + !$acc routine (interp_Hcollision_p_ACC) seq do ii = 1_idef,params%num_species - pchunk=params%pchunk - - if (spp(ii)%spatial_distribution.eq.'TRACER'.and. & - params%FO_GC_compare) then - call get_fields(params,spp(ii)%vars,F) - !! Calls [[get_fields]] in [[korc_fields]]. - ! Interpolates fields at local particles' position and keeps in - ! spp%vars. Fields in (R,\(\phi\),Z) coordinates. - - ALLOCATE(RAphi(spp(ii)%ppp,2)) - ALLOCATE(RVphi(spp(ii)%ppp)) - RAphi=0.0_rp - - call cart_to_cyl(spp(ii)%vars%X,spp(ii)%vars%Y) - - !$OMP PARALLEL DO SHARED(params,ii,spp,F,RAphi,RVphi) & - !$OMP& PRIVATE(pp,Bmag1,bhat,rm) - ! Call OpenMP to calculate p_par and mu for each particle and - ! put into spp%vars%V - do pp=1_idef,spp(ii)%ppp - if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then - - RVphi(pp)=(-sin(spp(ii)%vars%Y(pp,2))*spp(ii)%vars%V(pp,1)+ & - cos(spp(ii)%vars%Y(pp,2))*spp(ii)%vars%V(pp,2))* & - spp(ii)%vars%Y(pp,1) - - Bmag1 = SQRT(spp(ii)%vars%B(pp,1)*spp(ii)%vars%B(pp,1)+ & - spp(ii)%vars%B(pp,2)*spp(ii)%vars%B(pp,2)+ & - spp(ii)%vars%B(pp,3)*spp(ii)%vars%B(pp,3)) - - ! write(output_unit_write,'("pp: ",I16)') pp - ! write(output_unit_write,'("Bmag: ",E17.10)') Bmag - - - bhat = spp(ii)%vars%B(pp,:)/Bmag1 - - if (params%field_model(1:10).eq.'ANALYTICAL') then - rm=sqrt((spp(ii)%vars%Y(pp,1)-F%AB%Ro)**2+ & - (spp(ii)%vars%Y(pp,3))**2) - - RAphi(pp,1)=-F%AB%lambda**2*F%AB%Bo/(2*F%AB%qo)* & - log(1+(rm/F%AB%lambda)**2) - - else if (params%field_model(1:8).eq.'EXTERNAL') then - - RAphi(pp,1)=spp(ii)%vars%PSI_P(pp)/(2*C_PI) - - end if - - ! write(output_unit_write,'("bhat: ",E17.10)') bhat - ! write(output_unit_write,'("V: ",E17.10)') spp(ii)%vars%V(pp,:) - - - spp(ii)%vars%X(pp,:)=spp(ii)%vars%X(pp,:)- & - spp(ii)%m*spp(ii)%vars%g(pp)* & - cross(bhat,spp(ii)%vars%V(pp,:))/(spp(ii)%q*Bmag1) - - ! transforming from particle location to associated - ! GC location - - end if ! if particle in domain, i.e. spp%vars%flagCon==1 - end do ! loop over particles on an mpi process - !$OMP END PARALLEL DO - - call cart_to_cyl(spp(ii)%vars%X,spp(ii)%vars%Y) - call get_fields(params,spp(ii)%vars,F) - - !$OMP PARALLEL DO SHARED(params,ii,spp,F,RAphi,RVphi) & - !$OMP& PRIVATE(pp,rm) - ! Call OpenMP to calculate p_par and mu for each particle and - ! put into spp%vars%V - do pp=1_idef,spp(ii)%ppp - if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then - - if (params%field_model(1:10).eq.'ANALYTICAL') then - rm=sqrt((spp(ii)%vars%Y(pp,1)-F%AB%Ro)**2+ & - (spp(ii)%vars%Y(pp,3))**2) - RAphi(pp,2)=-F%AB%lambda**2*F%AB%Bo/(2*F%AB%qo)* & - log(1+(rm/F%AB%lambda)**2) - - else if (params%field_model(1:8).eq.'EXTERNAL') then - - RAphi(pp,2)=spp(ii)%vars%PSI_P(pp)/(2*C_PI) - - end if - - write(output_unit_write,'("RAphi1: ",E17.10)') RAphi(pp,1) - write(output_unit_write,'("RAphi2: ",E17.10)') RAphi(pp,2) - - spp(ii)%vars%V(pp,1)=(spp(ii)%m*spp(ii)%vars%g(pp)* & - RVphi(pp)+spp(ii)%q*(RAphi(pp,1)-RAphi(pp,2)))/ & - spp(ii)%vars%Y(pp,1) - !GC ppar - - end if ! if particle in domain, i.e. spp%vars%flagCon==1 - end do ! loop over particles on an mpi process - !$OMP END PARALLEL DO - - !$OMP PARALLEL DO SHARED(ii,spp) PRIVATE(pp,Bmagc,bhatc) - ! Call OpenMP to calculate p_par and mu for each particle and - ! put into spp%vars%V - do pp=1_idef,spp(ii)%ppp - if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then - - Bmagc = SQRT( DOT_PRODUCT(spp(ii)%vars%B(pp,:), & - spp(ii)%vars%B(pp,:))) - - bhatc = spp(ii)%vars%B(pp,:)/Bmagc - - spp(ii)%vars%V(pp,1)=spp(ii)%vars%V(pp,1)/ & - bhatc(2) - !GC ppar - - spp(ii)%vars%V(pp,2)=spp(ii)%m/(2*Bmagc)* & - (spp(ii)%vars%g(pp)**2- & - (1+(spp(ii)%vars%V(pp,1)/spp(ii)%m)**2)) - !GC mu - - - end if ! if particle in domain, i.e. spp%vars%flagCon==1 - end do ! loop over particles on an mpi process - !$OMP END PARALLEL DO + m_cache=spp(ii)%m params%GC_coords=.TRUE. - DEALLOCATE(RAphi) - DEALLOCATE(RVphi) - - !Preparing Output Data - call get_fields(params,spp(ii)%vars,F) - - !$OMP PARALLEL DO shared(F,params,spp) & - !$OMP& PRIVATE(cc,pp,E_PHI,Y_R) firstprivate(pchunk) - do pp=1_idef,spp(ii)%ppp,pchunk - - !$OMP SIMD - do cc=1_idef,pchunk - E_PHI(cc)=spp(ii)%vars%E(pp-1+cc,2) - Y_R(cc)=spp(ii)%vars%Y(pp-1+cc,1) - end do - !$OMP END SIMD - - call add_analytical_E_p(params,int8(0),F,E_PHI,Y_R,Y_Z) - - - !$OMP SIMD - do cc=1_idef,pchunk - spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) - end do - !$OMP END SIMD - - end do - !$OMP END PARALLEL DO - - !$OMP PARALLEL DO SHARED(ii,spp) PRIVATE(pp,Bmag1) - ! Call OpenMP to calculate p_par and mu for each particle and - ! put into spp%vars%V + !$acc parallel loop private(B_R,B_PHI,B_Z) do pp=1_idef,spp(ii)%ppp - if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then - Bmag1 = SQRT( DOT_PRODUCT(spp(ii)%vars%B(pp,:), & - spp(ii)%vars%B(pp,:))) - - spp(ii)%vars%g(pp)=sqrt(1+(spp(ii)%vars%V(pp,1))**2+ & - 2*spp(ii)%vars%V(pp,2)*Bmag1) - - ! write(output_unit_write,'("Bmag:",E17.10)') Bmag1 - ! write(output_unit_write,'("PPLL:",E17.10)') spp(ii)%vars%V(pp,1) - ! write(output_unit_write,'("MU:",E17.10)') spp(ii)%vars%V(pp,2) - - spp(ii)%vars%eta(pp) = atan2(sqrt(2*spp(ii)%m*Bmag1* & - spp(ii)%vars%V(pp,2)),spp(ii)%vars%V(pp,1))*180.0_rp/C_PI - - ! write(output_unit_write,'("BR",E17.10)') spp(ii)%vars%B(pp,1) - ! write(output_unit_write,'("BPHI",E17.10)') spp(ii)%vars%B(pp,2) - ! write(output_unit_write,'("BZ",E17.10)') spp(ii)%vars%B(pp,3) - - ! write(output_unit_write,'("ppll",E17.10)') spp(ii)%vars%V(pp,1) - ! write(output_unit_write,'("pperp",E17.10)') sqrt(2*spp(ii)%m*Bmag1* & - ! spp(ii)%vars%V(pp,2)) - - ! write(output_unit_write,'("eta GCinit",E17.10)') spp(ii)%vars%eta(pp) - ! write(output_unit_write,'("gam",E17.10)') spp(ii)%vars%g(pp) - - - end if ! if particle in domain, i.e. spp%vars%flagCon==1 - end do ! loop over particles on an mpi process - !$OMP END PARALLEL DO - else - - if ((spp(ii)%spatial_distribution.eq.'TRACER').or. & + if ((spp(ii)%spatial_distribution.eq.'TRACER').or. & (spp(ii)%spatial_distribution.eq.'TORUS').or. & (spp(ii)%spatial_distribution.eq.'DISK').or. & (spp(ii)%spatial_distribution.eq. & '2D-GAUSSIAN-ELLIPTIC-TORUS-MH')) & call cart_to_cyl(spp(ii)%vars%X,spp(ii)%vars%Y) - params%GC_coords=.TRUE. - - do pp=1_idef,spp(ii)%ppp - spp(ii)%vars%E(pp,1)=0._rp - spp(ii)%vars%E(pp,2)=0._rp - spp(ii)%vars%E(pp,3)=0._rp - end do - - - !write(6,*) 'before second get fields' - call get_fields(params,spp(ii)%vars,F) - !write(6,*) 'after second get fields' + Y_R=spp(ii)%vars%Y(pp,1) + Y_PHI=spp(ii)%vars%Y(pp,2) + Y_Z=spp(ii)%vars%Y(pp,3) - !write(output_unit_write,*) spp(1)%vars%PSI_P - - !$OMP PARALLEL DO SHARED(ii,spp) PRIVATE(pp,Bmag1) - - do pp=1_idef,spp(ii)%ppp - ! if ( spp(ii)%vars%flagCon(pp) .EQ. 1_is ) then + flagCon=spp(ii)%vars%flagCon(pp) + flagCol=spp(ii)%vars%flagCol(pp) - ! write(output_unit_write,'("BR: ",E17.10)') spp(ii)%vars%B(pp,1) - ! write(output_unit_write,'("BPHI: ",E17.10)') spp(ii)%vars%B(pp,2) - ! write(output_unit_write,'("BZ: ",E17.10)') spp(ii)%vars%B(pp,3) + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R, & + E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + flagCon,PSIp) - Bmag1 = SQRT( DOT_PRODUCT(spp(ii)%vars%B(pp,:), & - spp(ii)%vars%B(pp,:))) + Bmag=sqrt(B_R*B_R+B_PHI*B_PHI+B_Z*B_Z) - pmag=spp(ii)%m*sqrt(spp(ii)%vars%g(pp)**2-1) + pmag=sqrt(spp(ii)%vars%g(pp)**2-1) - spp(ii)%vars%V(pp,1)=pmag*cos(deg2rad(spp(ii)%vars%eta(pp))) + spp(ii)%vars%V(pp,1)=pmag*cos(deg2rad(spp(ii)%vars%eta(pp))) - !write(6,*) 'GC_init' - !write(6,*) spp(ii)%m,Bmag1 + spp(ii)%vars%V(pp,2)=(pmag* & + sin(deg2rad(spp(ii)%vars%eta(pp))))**2/ & + (2*spp(ii)%m*Bmag) - spp(ii)%vars%V(pp,2)=(pmag* & - sin(deg2rad(spp(ii)%vars%eta(pp))))**2/ & - (2*spp(ii)%m*Bmag1) + spp(ii)%vars%B(pp,1) = B_R + spp(ii)%vars%B(pp,2) = B_PHI + spp(ii)%vars%B(pp,3) = B_Z - ! write(output_unit_write,'("BR",E17.10)') spp(ii)%vars%B(pp,1) - ! write(output_unit_write,'("BPHI",E17.10)') spp(ii)%vars%B(pp,2) - ! write(output_unit_write,'("BZ",E17.10)') spp(ii)%vars%B(pp,3) + spp(ii)%vars%E(pp,2) = E_PHI + spp(ii)%vars%PSI_P(pp) = PSIp - !write(output_unit_write,'("ppll",E17.10)') spp(ii)%vars%V(pp,1) - !write(output_unit_write,'("mu",E17.10)') spp(ii)%vars%V(pp,2) + V_PLL=spp(ii)%vars%V(pp,1) + V_MU=spp(ii)%vars%V(pp,2) - ! write(output_unit_write,'("eta",E17.10)') spp(ii)%vars%eta(pp) - ! write(output_unit_write,'("gam",E17.10)') spp(ii)%vars%g(pp) + if (params%collisions) then + if (params%profile_model(1:10).eq.'ANALYTICAL') then + call analytical_profiles_p_ACC(0._rp,params,Y_R,Y_Z, & + P,F,ne,Te,Zeff,PSIp) + else if (params%profile_model(10:10).eq.'H') then + call interp_Hcollision_p_ACC(params,Y_R,Y_PHI,Y_Z,ne,Te,Zeff, & + nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1,flagCon) - ! end if ! if particle in domain, i.e. spp%vars%flagCon==1 - end do ! loop over particles on an mpi process - !$OMP END PARALLEL DO - - !$OMP PARALLEL DO shared(F,params,spp) & - !$OMP& PRIVATE(pp,cc,E_PHI,Y_R) firstprivate(pchunk) - do pp=1_idef,spp(ii)%ppp,pchunk - - !$OMP SIMD - do cc=1_idef,pchunk - E_PHI(cc)=spp(ii)%vars%E(pp-1+cc,2) - Y_R(cc)=spp(ii)%vars%Y(pp-1+cc,1) - end do - !$OMP END SIMD - - if (params%field_model(1:8).eq.'EXTERNAL') then - call add_analytical_E_p(params,int8(0),F,E_PHI,Y_R,Y_Z) - end if + spp(ii)%vars%ne(pp) = ne + spp(ii)%vars%Te(pp) = Te + spp(ii)%vars%Zeff(pp) = Zeff + spp(ii)%vars%nimp(pp,1) = nAr0 + spp(ii)%vars%nimp(pp,1) = nAr1 + spp(ii)%vars%nimp(pp,1) = nAr2 + spp(ii)%vars%nimp(pp,1) = nAr3 + spp(ii)%vars%nimp(pp,1) = nAr4 + spp(ii)%vars%nimp(pp,1) = nD + spp(ii)%vars%nimp(pp,1) = nD1 + endif - !$OMP SIMD - do cc=1_idef,pchunk + end if - spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) - end do - !$OMP END SIMD + spp(ii)%vars%flagCon(pp)=flagCon + spp(ii)%vars%flagCol(pp)=flagCol end do - !$OMP END PARALLEL DO + !$acc end parallel loop - end if + spp(ii)%vars%Yborn=spp(ii)%vars%Y - spp(ii)%vars%Yborn=spp(ii)%vars%Y + end do ! loop over particle species - end do ! loop over particle species -end subroutine GC_init +end subroutine GC_init_ACC FUNCTION deg2rad(x) REAL(rp), INTENT(IN) :: x @@ -4514,6 +4665,8 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !! time iterator. real(rp),dimension(F%dim_1D) :: Vden,Vdenave,VdenOMP INTEGER :: newREs + REAL(rp),DIMENSION(params%pchunk) :: Zeff + REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species) :: nimp do ii = 1_idef,params%num_species @@ -4535,7 +4688,7 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !$OMP& PRIVATE(pp,tt,ttt,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & !$OMP& flagCon,flagCol,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,PSIp,ne, & !$OMP& Vden,Vdenave,gradB_R,gradB_PHI,gradB_Z, & - !$OMP& curlb_R,curlb_PHI,curlb_Z) & + !$OMP& curlb_R,curlb_PHI,curlb_Z,Te,Zeff,nimp) & !$OMP& REDUCTION(+:VdenOMP) do pp=1_idef,spp(ii)%ppp,pchunk @@ -4574,7 +4727,7 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) call include_CoulombCollisions_GC_p(tt+params%t_skip*(ttt-1),params, & random,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,m_cache,flagCon,flagCol, & - F,P,E_PHI,ne,PSIp) + F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) end if @@ -4619,7 +4772,7 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) call include_CoulombCollisions_GC_p(tt,params,random, & Y_R,Y_PHI,Y_Z,V_PLL,V_MU,m_cache,flagCon,flagCol, & - F,P,E_PHI,ne,PSIp) + F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) end do @@ -5263,12 +5416,13 @@ subroutine advance_FPeqn_vars(params,random,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,flagCon,fla REAL(rp),DIMENSION(params%pchunk) :: E_PHI INTEGER(is),DIMENSION(params%pchunk), INTENT(INOUT) :: flagCon,flagCol REAL(rp),intent(in) :: m_cache - REAL(rp),DIMENSION(params%pchunk) :: ne + REAL(rp),DIMENSION(params%pchunk) :: ne,Te,Zeff + REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species) :: nimp do tt=1_ip,params%t_skip call include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & - V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,PSIp) + V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) ! write(output_unit_write,'("Collision Loop in FP")') @@ -5300,6 +5454,7 @@ subroutine adv_GCinterp_psi_top(params,random,spp,P,F) REAL(rp),DIMENSION(params%pchunk) :: gradB_R,gradB_PHI,gradB_Z INTEGER(is),DIMENSION(params%pchunk) :: flagCon,flagCol REAL(rp) :: m_cache,q_cache,B0,EF0,R0,q0,lam,ar + REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species) :: nimp INTEGER :: ii @@ -5327,7 +5482,7 @@ subroutine adv_GCinterp_psi_top(params,random,spp,P,F) !$OMP& SHARED(params,ii,spp,P,F,random) & !$OMP& PRIVATE(pp,tt,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,B_R,B_PHI,B_Z, & !$OMP& flagCon,flagCol,E_PHI,PSIp,curlb_R,curlb_PHI,curlb_Z, & - !$OMP& gradB_R,gradB_PHI,gradB_Z,ne,E_R,E_Z, & + !$OMP& gradB_R,gradB_PHI,gradB_Z,ne,Te,Zeff,nimp,E_R,E_Z, & !$OMP& Y_R0,Y_PHI0,Y_Z0) do pp=1_idef,spp(ii)%ppp,pchunk @@ -5365,7 +5520,7 @@ subroutine adv_GCinterp_psi_top(params,random,spp,P,F) if (params%collisions) then call include_CoulombCollisions_GC_p(tt,params,random, & Y_R,Y_PHI,Y_Z,V_PLL,V_MU,m_cache,flagCon,flagCol, & - F,P,E_PHI,ne,PSIp) + F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) end if end do !timestep iterator @@ -5786,6 +5941,238 @@ end subroutine adv_GCinterp_fio_top #ifdef PSPLINE +subroutine adv_GCinterp_psiwE_top_ACC(params,random,spp,P,F) + + IMPLICIT NONE + + TYPE(KORC_PARAMS), INTENT(INOUT) :: params + !! Core KORC simulation parameters. + CLASS(random_context), POINTER, INTENT(INOUT) :: random + TYPE(PROFILES), INTENT(IN) :: P + TYPE(FIELDS), INTENT(INOUT) :: F + 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) :: Y_R,Y_PHI,Y_Z + REAL(rp) :: Y_R0,Y_PHI0,Y_Z0 + REAL(rp) :: Y_R1,Y_PHI1,Y_Z1 + REAL(rp) :: B_R,B_PHI,B_Z + REAL(rp) :: E_R,E_PHI,E_Z + REAL(rp) :: ne,Te,Zeff + REAL(rp) :: V_PLL,V_MU + REAL(rp) :: PSIp + REAL(rp) :: curlb_R,curlb_PHI,curlb_Z + REAL(rp) :: gradB_R,gradB_PHI,gradB_Z + INTEGER(is) :: flagCon,flagCol + REAL(rp) :: m_cache,q_cache,B0,EF0,R0,q0,lam,ar + REAL(rp), DIMENSION(params%num_impurity_species) :: nimp + INTEGER :: ii + !! Species iterator. + INTEGER :: pp + !! Particles iterator. + INTEGER(ip) :: tt,ttt + !! time iterator. + INTEGER :: thread_num + + !$acc routine (advance_GCinterp_psiwE_vars_ACC) seq + !$acc routine (include_CoulombCollisions_GC_p_ACC) seq + + do ii = 1_idef,params%num_species + + q_cache=spp(ii)%q + m_cache=spp(ii)%m + + !write(6,*) q_cache,m_cache + + if (.not.params%LargeCollisions) then + + !$acc parallel loop firstprivate(m_cache) & + !$acc& private(Y_R,Y_PHI,Y_Z,Y_R0,Y_PHI0,Y_Z0, & + !$acc& Y_R1,Y_PHI1,Y_Z1,V_PLL,V_MU,flagCon,flagCol, & + !$acc& B_R,B_PHI,B_Z,PSIp,E_PHI,ne,Te,Zeff,nimp) + do pp=1_idef,spp(ii)%ppp + + Y_R=spp(ii)%vars%Y(pp,1) + Y_PHI=spp(ii)%vars%Y(pp,2) + Y_Z=spp(ii)%vars%Y(pp,3) + + Y_R0=spp(ii)%vars%Y0(pp,1) + Y_PHI0=spp(ii)%vars%Y0(pp,2) + Y_Z0=spp(ii)%vars%Y0(pp,3) + Y_R1=spp(ii)%vars%Y1(pp,1) + Y_PHI1=spp(ii)%vars%Y1(pp,2) + Y_Z1=spp(ii)%vars%Y1(pp,3) + + V_PLL=spp(ii)%vars%V(pp,1) + V_MU=spp(ii)%vars%V(pp,2) + + PSIp=spp(ii)%vars%PSI_P(pp) + + flagCon=spp(ii)%vars%flagCon(pp) + flagCol=spp(ii)%vars%flagCol(pp) + + !$acc loop seq + do tt=1_ip,params%t_skip + + call advance_GCinterp_psiwE_vars_ACC(spp(ii),pp,tt, & + params, & + Y_R,Y_PHI,Y_Z,V_PLL,V_MU,flagCon,flagCol, & + F,P,B_R,B_PHI,B_Z,E_PHI,PSIp,curlb_R,curlb_PHI, & + curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + Y_R0,Y_PHI0,Y_Z0,Y_R1,Y_PHI1,Y_Z1) + + if (params%collisions) then + + call include_CoulombCollisions_GC_p_ACC(spp(ii),pp,tt,params, & + random,Y_R,Y_PHI,Y_Z, V_PLL,V_MU,m_cache, & + flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) + + end if + + + end do !timestep iterator + + spp(ii)%vars%Y(pp,1)=Y_R + spp(ii)%vars%Y(pp,2)=Y_PHI + spp(ii)%vars%Y(pp,3)=Y_Z + spp(ii)%vars%V(pp,1)=V_PLL + spp(ii)%vars%V(pp,2)=V_MU + + spp(ii)%vars%Y0(pp,1)=Y_R0 + spp(ii)%vars%Y0(pp,2)=Y_PHI0 + spp(ii)%vars%Y0(pp,3)=Y_Z0 + spp(ii)%vars%Y1(pp,1)=Y_R1 + spp(ii)%vars%Y1(pp,2)=Y_PHI1 + spp(ii)%vars%Y1(pp,3)=Y_Z1 + + spp(ii)%vars%flagCon(pp)=flagCon + spp(ii)%vars%flagCol(pp)=flagCol + + spp(ii)%vars%B(pp,1) = B_R + spp(ii)%vars%B(pp,2) = B_PHI + spp(ii)%vars%B(pp,3) = B_Z + + spp(ii)%vars%E(pp,2) = E_PHI + spp(ii)%vars%PSI_P(pp) = PSIp + + spp(ii)%vars%ne(pp) = ne + spp(ii)%vars%Te(pp) = Te + spp(ii)%vars%Zeff(pp) = Zeff + spp(ii)%vars%nimp(pp,:) = nimp + + Bmag=sqrt(B_R*B_R+B_PHI*B_PHI+ & + B_Z*B_Z) + + spp(ii)%vars%g(pp)=sqrt(1+V_PLL**2+ & + 2*V_MU*Bmag) + + spp(ii)%vars%eta(pp) = atan2(sqrt(2*m_cache*Bmag* & + spp(ii)%vars%V(pp,2)),spp(ii)%vars%V(pp,1))* & + 180.0_rp/C_PI + + end do + !$acc end parallel loopCLASS(random_context), POINTER, INTENT(INOUT) :: random + + else if (params%LargeCollisions) then + + do tt=1_ip,params%coll_per_dump + + !$acc parallel loop firstprivate(m_cache,tt) & + !$acc& private(ttt,Y_R,Y_PHI,Y_Z,Y_R0,Y_PHI0,Y_Z0, & + !$acc& Y_R1,Y_PHI1,Y_Z1,V_PLL,V_MU,flagCon,flagCol, & + !$acc& B_R,B_PHI,B_Z,PSIp,E_PHI,ne,Te,Zeff,nimp) + do pp=1_idef,spp(ii)%pRE + + Y_R=spp(ii)%vars%Y(pp,1) + Y_PHI=spp(ii)%vars%Y(pp,2) + Y_Z=spp(ii)%vars%Y(pp,3) + + Y_R0=spp(ii)%vars%Y0(pp,1) + Y_PHI0=spp(ii)%vars%Y0(pp,2) + Y_Z0=spp(ii)%vars%Y0(pp,3) + Y_R1=spp(ii)%vars%Y1(pp,1) + Y_PHI1=spp(ii)%vars%Y1(pp,2) + Y_Z1=spp(ii)%vars%Y1(pp,3) + + V_PLL=spp(ii)%vars%V(pp,1) + V_MU=spp(ii)%vars%V(pp,2) + + PSIp=spp(ii)%vars%PSI_P(pp) + + flagCon=spp(ii)%vars%flagCon(pp) + flagCol=spp(ii)%vars%flagCol(pp) + + !$acc loop seq + do ttt=1_ip,params%orbits_per_coll + + call advance_GCinterp_psiwE_vars_ACC(spp(ii),pp,tt, & + params, & + Y_R,Y_PHI,Y_Z,V_PLL,V_MU,flagCon,flagCol, & + F,P,B_R,B_PHI,B_Z,E_PHI,PSIp,curlb_R,curlb_PHI, & + curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + Y_R0,Y_PHI0,Y_Z0,Y_R1,Y_PHI1,Y_Z1) + + if (params%collisions) then + + call include_CoulombCollisions_GC_p_ACC(spp(ii),pp,tt,params, & + random,Y_R,Y_PHI,Y_Z, V_PLL,V_MU,m_cache, & + flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) + + end if + + + end do !timestep iterator + + spp(ii)%vars%Y(pp,1)=Y_R + spp(ii)%vars%Y(pp,2)=Y_PHI + spp(ii)%vars%Y(pp,3)=Y_Z + spp(ii)%vars%V(pp,1)=V_PLL + spp(ii)%vars%V(pp,2)=V_MU + + spp(ii)%vars%Y0(pp,1)=Y_R0 + spp(ii)%vars%Y0(pp,2)=Y_PHI0 + spp(ii)%vars%Y0(pp,3)=Y_Z0 + spp(ii)%vars%Y1(pp,1)=Y_R1 + spp(ii)%vars%Y1(pp,2)=Y_PHI1 + spp(ii)%vars%Y1(pp,3)=Y_Z1 + + spp(ii)%vars%flagCon(pp)=flagCon + spp(ii)%vars%flagCol(pp)=flagCol + + spp(ii)%vars%B(pp,1) = B_R + spp(ii)%vars%B(pp,2) = B_PHI + spp(ii)%vars%B(pp,3) = B_Z + + spp(ii)%vars%E(pp,2) = E_PHI + spp(ii)%vars%PSI_P(pp) = PSIp + + spp(ii)%vars%ne(pp) = ne + spp(ii)%vars%Te(pp) = Te + spp(ii)%vars%Zeff(pp) = Zeff + spp(ii)%vars%nimp(pp,:) = nimp + + Bmag=sqrt(B_R*B_R+B_PHI*B_PHI+ & + B_Z*B_Z) + + spp(ii)%vars%g(pp)=sqrt(1+V_PLL**2+ & + 2*V_MU*Bmag) + + spp(ii)%vars%eta(pp) = atan2(sqrt(2*m_cache*Bmag* & + spp(ii)%vars%V(pp,2)),spp(ii)%vars%V(pp,1))* & + 180.0_rp/C_PI + + end do + !$acc end parallel loop + + end do + + endif + + end do !species iterator + + end subroutine adv_GCinterp_psiwE_top_ACC + subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) IMPLICIT NONE @@ -5811,6 +6198,7 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) REAL(rp),DIMENSION(params%pchunk) :: gradB_R,gradB_PHI,gradB_Z INTEGER(is),DIMENSION(params%pchunk) :: flagCon,flagCol REAL(rp) :: m_cache,q_cache,B0,EF0,R0,q0,lam,ar + REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species) :: nimp INTEGER :: ii @@ -5836,7 +6224,7 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) !$OMP& SHARED(params,ii,spp,P,F,random) & !$OMP& PRIVATE(pp,tt,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,B_R,B_PHI,B_Z, & !$OMP& flagCon,flagCol,E_PHI,PSIp,curlb_R,curlb_PHI,curlb_Z, & - !$OMP& gradB_R,gradB_PHI,gradB_Z,ne,E_R,E_Z,thread_num, & + !$OMP& gradB_R,gradB_PHI,gradB_Z,ne,Te,Zeff,nimp,E_R,E_Z,thread_num, & !$OMP& Y_R0,Y_PHI0,Y_Z0,Y_R1,Y_PHI1,Y_Z1) do pp=1_idef,spp(ii)%ppp,pchunk @@ -5891,7 +6279,7 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) call include_CoulombCollisions_GC_p(tt,params,random, & Y_R,Y_PHI,Y_Z, V_PLL,V_MU,m_cache, & - flagCon,flagCol,F,P,E_PHI,ne,PSIp) + flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) end if @@ -5921,18 +6309,13 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) spp(ii)%vars%B(pp-1+cc,2) = B_PHI(cc) spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc) - spp(ii)%vars%gradB(pp-1+cc,1) = gradB_R(cc) - spp(ii)%vars%gradB(pp-1+cc,2) = gradB_PHI(cc) - spp(ii)%vars%gradB(pp-1+cc,3) = gradB_Z(cc) - - spp(ii)%vars%curlb(pp-1+cc,1) = curlb_R(cc) - spp(ii)%vars%curlb(pp-1+cc,2) = curlb_PHI(cc) - spp(ii)%vars%curlb(pp-1+cc,3) = curlb_Z(cc) - spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc) spp(ii)%vars%ne(pp-1+cc) = ne(cc) + spp(ii)%vars%Te(pp-1+cc) = Te(cc) + spp(ii)%vars%Zeff(pp-1+cc) = Zeff(cc) + spp(ii)%vars%nimp(pp-1+cc,:) = nimp(cc,:) end do !$OMP END SIMD @@ -5940,7 +6323,7 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) do tt=1_ip,params%t_skip call include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & - V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,PSIp) + V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) end do @@ -6023,7 +6406,7 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) !$OMP& PRIVATE(pp,ttt,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & !$OMP& B_R,B_PHI,B_Z,achunk, & !$OMP& flagCon,flagCol,E_PHI,PSIp,curlb_R,curlb_PHI,curlb_Z, & - !$OMP& gradB_R,gradB_PHI,gradB_Z,ne,Te,E_R,E_Z,thread_num, & + !$OMP& gradB_R,gradB_PHI,gradB_Z,ne,Te,Zeff,nimp,E_R,E_Z,thread_num, & !$OMP& Y_R0,Y_PHI0,Y_Z0,Y_R1,Y_PHI1,Y_Z1) do pp=1_idef,spp(ii)%pRE,pchunk @@ -6179,6 +6562,9 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc) spp(ii)%vars%ne(pp-1+cc) = ne(cc) + spp(ii)%vars%Te(pp-1+cc) = Te(cc) + spp(ii)%vars%Zeff(pp-1+cc) = Zeff(cc) + spp(ii)%vars%nimp(pp-1+cc,:) = nimp(cc,:) end do !$OMP END SIMD @@ -7698,6 +8084,305 @@ subroutine advance_GCinterp_psiwE_vars(spp,pchunk,pp,tt,params,random,Y_R,Y_PHI, #endif DBG_CHECK end subroutine advance_GCinterp_psiwE_vars +subroutine advance_GCinterp_psiwE_vars_ACC(spp,pp,tt,params,Y_R,Y_PHI,Y_Z, & + V_PLL,V_MU,flagCon,flagCol,F,P,B_R,B_PHI,B_Z,E_PHI,PSIp, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + Y_R0,Y_PHI0,Y_Z0,Y_R1,Y_PHI1,Y_Z1) + !$acc routine seq + !! @note Subroutine to advance GC variables \(({\bf X},p_\parallel)\) + !! @endnote + !! Comment this section further with evolution equations, numerical + !! methods, and descriptions of both. + TYPE(KORC_PARAMS), INTENT(INOUT) :: params + !! Core KORC simulation parameters. + TYPE(SPECIES), INTENT(INOUT) :: spp + TYPE(PROFILES), INTENT(IN) :: P + TYPE(FIELDS), INTENT(IN) :: F + REAL(rp) :: dt + !! Time step used in the leapfrog step (\(\Delta t\)). + + INTEGER :: cc,ii + INTEGER(ip),intent(in) :: tt + !! time iterator. + INTEGER,intent(in) :: pp + REAL(rp) :: Bmag + REAL(rp) :: a1 = 1./5._rp + REAL(rp) :: a21 = 3./40._rp,a22=9./40._rp + REAL(rp) :: a31 = 3./10._rp,a32=-9./10._rp,a33=6./5._rp + REAL(rp) :: a41 = -11./54._rp,a42=5./2._rp,a43=-70./27._rp,a44=35./27._rp + REAL(rp) :: a51 = 1631./55296._rp,a52=175./512._rp,a53=575./13824._rp,a54=44275./110592._rp,a55=253./4096._rp + REAL(rp) :: b1=37./378._rp,b2=0._rp,b3=250./621._rp,b4=125./594._rp,b5=0._rp,b6=512./1771._rp + + REAL(rp) :: k1_R,k1_PHI,k1_Z,k1_PLL,k1_MU + REAL(rp):: k2_R,k2_PHI,k2_Z,k2_PLL,k2_MU + REAL(rp) :: k3_R,k3_PHI,k3_Z,k3_PLL,k3_MU + REAL(rp) :: k4_R,k4_PHI,k4_Z,k4_PLL,k4_MU + REAL(rp) :: k5_R,k5_PHI,k5_Z,k5_PLL,k5_MU + REAL(rp) :: k6_R,k6_PHI,k6_Z,k6_PLL,k6_MU + REAL(rp) :: Y0_R,Y0_PHI,Y0_Z + REAL(rp),INTENT(INOUT) :: Y_R,Y_PHI,Y_Z + REAL(rp),INTENT(INOUT) :: Y_R0,Y_PHI0,Y_Z0 + REAL(rp),INTENT(INOUT) :: Y_R1,Y_PHI1,Y_Z1 + REAL(rp),INTENT(OUT) :: B_R,B_PHI,B_Z + REAL(rp) :: E_R,E_Z + REAL(rp),INTENT(OUT) :: E_PHI + REAL(rp),INTENT(OUT) :: PSIp + REAL(rp),INTENT(OUT) :: curlb_R,curlb_PHI,curlb_Z + REAL(rp),INTENT(OUT) :: gradB_R,gradB_PHI,gradB_Z + REAL(rp),INTENT(INOUT) :: V_PLL,V_MU + REAL(rp) :: RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU + REAL(rp) :: V0_PLL,V0_MU + REAL(rp) :: Te,Zeff + + INTEGER(is),intent(INOUT) :: flagCon,flagCol + INTEGER(is) :: flagCon0 + REAL(rp) :: q_cache,m_cache + LOGICAL :: accepted + REAL(rp) :: Rmin,Rmax,Zmin,Zmax,Rtrial,Ztrial,rm_trial,pmag0,Bmag0,maxRnRE + REAL(rp),dimension(spp%BMC_Nra) :: RnRE + + !$acc routine (check_if_in_fields_domain_2D_ACC) seq + !$acc routine (calculate_GCfieldswE_p_ACC) seq + !$acc routine (GCEoM1_p_ACC) seq + + dt=params%dt + q_cache=spp%q + m_cache=spp%m + + + + Y0_R=Y_R + Y0_PHI=Y_PHI + Y0_Z=Y_Z + V0_PLL=V_PLL + V0_MU=V_MU + + flagCon0=flagCon + + + ! write(output_unit_write,*) 'R0',Y_R(1) + ! write(output_unit_write,*) 'PHI0',Y_PHI(1) + ! write(output_unit_write,*) 'Z0',Y_Z(1) + ! write(output_unit_write,*) 'PPLL0',V_PLL(1) + ! write(output_unit_write,*) 'MU0',V_MU(1) + + call check_if_in_fields_domain_2D_ACC(F,Y_R,Y_PHI,Y_Z,flagCon) + + ! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + flagCon,PSIp) + + call GCEoM1_p_ACC(pp,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU,B_R,B_PHI, & + B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & + gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp,flagCon) + + ! write(output_unit_write,*) 'R0',Y_R(1) + ! write(output_unit_write,*) 'PHI0',Y_PHI(1) + ! write(output_unit_write,*) 'Z0',Y_Z(1) + ! write(output_unit_write,*) 'PPLL0',V_PLL(1) + ! write(output_unit_write,*) 'MU0',V_MU(1) + + ! write(output_unit_write,*) 'BR',B_R(1) + ! write(output_unit_write,*) 'BPHI',B_PHI(1) + ! write(output_unit_write,*) 'BZ',B_Z(1) + + ! write(output_unit_write,*) 'gradBR',gradB_R(1) + ! write(output_unit_write,*) 'gradBPHI',gradB_PHI(1) + ! write(output_unit_write,*) 'gradBZ',gradB_Z(1) + + ! write(output_unit_write,*) 'curlBR',curlB_R(1) + ! write(output_unit_write,*) 'curlBPHI',curlB_PHI(1) + ! write(output_unit_write,*) 'curlBZ',curlB_Z(1) + + ! write(output_unit_write,*) 'dt',params%dt + ! write(output_unit_write,*) 'RHS_R',RHS_R(1) + ! write(output_unit_write,*) 'RHS_PHI',RHS_PHI(1) + ! write(output_unit_write,*) 'RHS_Z',RHS_Z(1) + ! write(output_unit_write,*) 'RHS_PLL',RHS_PLL(1) + ! write(output_unit_write,*) 'RHS_MU',RHS_MU(1) + + + k1_R=dt*RHS_R + k1_PHI=dt*RHS_PHI + k1_Z=dt*RHS_Z + k1_PLL=dt*RHS_PLL + k1_MU=dt*RHS_MU + + !write(6,*) 'advance1_gceom',pp,q_cache,m_cache + !flush(6) + + Y_R=Y0_R+a1*k1_R + Y_PHI=Y0_PHI+a1*k1_PHI + Y_Z=Y0_Z+a1*k1_Z + V_PLL=V0_PLL +a1*k1_PLL + V_MU=V0_MU +a1*k1_MU + + + ! write(output_unit_write,*) 'R1',Y_R(1) + ! write(output_unit_write,*) 'PHI1',Y_PHI(1) + !! write(output_unit_write,*) 'Z1',Y_Z(1) + ! write(output_unit_write,*) 'PPLL1',V_PLL(1) + ! write(output_unit_write,*) 'MU1',V_MU(1) + + call check_if_in_fields_domain_2D_ACC(F,Y_R,Y_PHI,Y_Z,flagCon) + ! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + flagCon,PSIp) + + !write(6,*) 'advance2',pp,Y_R,Y_Z,B_R + !flush(6) + + + call GCEoM1_p_ACC(pp,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU,B_R,B_PHI, & + B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & + gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp,flagCon) + + + k2_R=dt*RHS_R + k2_PHI=dt*RHS_PHI + k2_Z=dt*RHS_Z + k2_PLL=dt*RHS_PLL + k2_MU=dt*RHS_MU + + Y_R=Y0_R+a21*k1_R+a22*k2_R + Y_PHI=Y0_PHI+a21*k1_PHI+a22*k2_PHI + Y_Z=Y0_Z+a21*k1_Z+a22*k2_Z + V_PLL=V0_PLL +a21*k1_PLL+a22*k2_PLL + V_MU=V0_MU +a21*k1_MU+a22*k2_MU + + + + call check_if_in_fields_domain_2D_ACC(F,Y_R,Y_PHI,Y_Z,flagCon) + ! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + flagCon,PSIp) + + + call GCEoM1_p_ACC(pp,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU,B_R,B_PHI, & + B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & + gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp,flagCon) + + + k3_R=dt*RHS_R + k3_PHI=dt*RHS_PHI + k3_Z=dt*RHS_Z + k3_PLL=dt*RHS_PLL + k3_MU=dt*RHS_MU + + Y_R=Y0_R+a31*k1_R+a32*k2_R+a33*k3_R + Y_PHI=Y0_PHI+a31*k1_PHI+a32*k2_PHI+ & + a33*k3_PHI + Y_Z=Y0_Z+a31*k1_Z+a32*k2_Z+a33*k3_Z + + V_PLL=V0_PLL +a31*k1_PLL+a32*k2_PLL+a33*k3_PLL + V_MU=V0_MU +a31*k1_MU+a32*k2_MU+a33*k3_MU + + call check_if_in_fields_domain_2D_ACC(F,Y_R,Y_PHI,Y_Z,flagCon) + ! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + flagCon,PSIp) + + call GCEoM1_p_ACC(pp,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU,B_R,B_PHI, & + B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & + gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp,flagCon) + + k4_R=dt*RHS_R + k4_PHI=dt*RHS_PHI + k4_Z=dt*RHS_Z + k4_PLL=dt*RHS_PLL + k4_MU=dt*RHS_MU + + Y_R=Y0_R+a41*k1_R+a42*k2_R+a43*k3_R+ & + a44*k4_R + Y_PHI=Y0_PHI+a41*k1_PHI+a42*k2_PHI+ & + a43*k3_PHI+a44*k4_PHI + Y_Z=Y0_Z+a41*k1_Z+a42*k2_Z+a43*k3_Z+ & + a44*k4_Z + V_PLL=V0_PLL +a41*k1_PLL+a42*k2_PLL+ & + a43*k3_PLL+a44*k4_PLL + V_MU=V0_MU +a41*k1_MU+a42*k2_MU+ & + a43*k3_MU+a44*k4_MU + + + call check_if_in_fields_domain_2D_ACC(F,Y_R,Y_PHI,Y_Z,flagCon) + ! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + flagCon,PSIp) + + call GCEoM1_p_ACC(pp,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU,B_R,B_PHI, & + B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & + gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp,flagCon) + + k5_R=dt*RHS_R + k5_PHI=dt*RHS_PHI + k5_Z=dt*RHS_Z + k5_PLL=dt*RHS_PLL + k5_MU=dt*RHS_MU + + Y_R=Y0_R+a51*k1_R+a52*k2_R+a53*k3_R+ & + a54*k4_R+a55*k5_R + Y_PHI=Y0_PHI+a51*k1_PHI+a52*k2_PHI+ & + a53*k3_PHI+a54*k4_PHI+a55*k5_PHI + Y_Z=Y0_Z+a51*k1_Z+a52*k2_Z+a53*k3_Z+ & + a54*k4_Z+a55*k5_Z + V_PLL=V0_PLL +a51*k1_PLL+a52*k2_PLL+ & + a53*k3_PLL+a54*k4_PLL+a55*k5_PLL + V_MU=V0_MU +a51*k1_MU+a52*k2_MU+ & + a53*k3_MU+a54*k4_MU+a55*k5_MU + + call check_if_in_fields_domain_2D_ACC(F,Y_R,Y_PHI,Y_Z,flagCon) + ! call interp_fields_p(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + flagCon,PSIp) + + call GCEoM1_p_ACC(pp,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU,B_R,B_PHI, & + B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & + gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp,flagCon) + + + k6_R=dt*RHS_R + k6_PHI=dt*RHS_PHI + k6_Z=dt*RHS_Z + k6_PLL=dt*RHS_PLL + k6_MU=dt*RHS_MU + + Y_R=Y0_R+(b1*k1_R+b2*k2_R+ & + b3*k3_R+b4*k4_R+b5*k5_R+b6*k6_R)* & + REAL(flagCol)*REAL(flagCon0) + Y_PHI=Y0_PHI+(b1*k1_PHI+b2*k2_PHI+ & + b3*k3_PHI+b4*k4_PHI+b5*k5_PHI+b6*k6_PHI)* & + REAL(flagCol)*REAL(flagCon0) + Y_Z=Y0_Z+(b1*k1_Z+b2*k2_Z+ & + b3*k3_Z+b4*k4_Z+b5*k5_Z+b6*k6_Z)* & + REAL(flagCol)*REAL(flagCon0) + V_PLL=V0_PLL+(b1*k1_PLL+b2*k2_PLL+ & + b3*k3_PLL+b4*k4_PLL+b5*k5_PLL+b6*k6_PLL)* & + REAL(flagCol)*REAL(flagCon0) + V_MU=V0_MU+(b1*k1_MU+b2*k2_MU+ & + b3*k3_MU+b4*k4_MU+b5*k5_MU+b6*k6_MU)* & + REAL(flagCol)*REAL(flagCon0) + + Y_R1=Y_R1+(Y_R0-Y_R1)*REAL(flagCon0) + Y_PHI1=Y_PHI1+(Y_PHI0-Y_PHI1)*REAL(flagCon0) + Y_Z1=Y_Z1+(Y_Z0-Y_Z1)*REAL(flagCon0) + + Y_R0=Y_R0+(Y0_R-Y_R0)*REAL(flagCon0) + Y_PHI0=Y_PHI0+(Y0_PHI-Y_PHI0)*REAL(flagCon0) + Y_Z0=Y_Z0+(Y0_Z-Y_Z0)*REAL(flagCon0) + + call check_if_in_fields_domain_2D_ACC(F,Y_R,Y_PHI,Y_Z,flagCon) + call calculate_GCfieldswE_p_ACC(F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI, & + E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z, & + flagCon,PSIp) + +end subroutine advance_GCinterp_psiwE_vars_ACC + FUNCTION fRE_BMC(Nr_a,r_a,nRE,rm) REAL(rp), INTENT(IN) :: rm INTEGER :: Nr_a @@ -7794,6 +8479,7 @@ subroutine advance_GCinterp_psi2x1t_vars(vars,pp,tt,params,random,Y_R,Y_PHI,Y_Z, REAL(rp),DIMENSION(params%pchunk) :: RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU REAL(rp),DIMENSION(params%pchunk) :: V0_PLL,V0_MU REAL(rp),DIMENSION(params%pchunk) :: Te,Zeff + REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species) :: nimp INTEGER(is),DIMENSION(params%pchunk),intent(INOUT) :: flagCon,flagCol REAL(rp),intent(IN) :: q_cache,m_cache @@ -8122,7 +8808,7 @@ subroutine advance_GCinterp_psi2x1t_vars(vars,pp,tt,params,random,Y_R,Y_PHI,Y_Z, if (params%collisions) then call include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & - V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,PSIp) + V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) end if @@ -8145,12 +8831,14 @@ subroutine advance_FPinterp_vars(params,random,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & REAL(rp),intent(in) :: m_cache INTEGER(is),DIMENSION(params%pchunk),intent(INOUT) :: flagCon,flagCol REAL(rp),DIMENSION(params%pchunk), INTENT(OUT) :: ne + REAL(rp),DIMENSION(params%pchunk) :: Te,Zeff + REAL(rp), DIMENSION(params%pchunk,params%num_impurity_species) :: nimp ! write(output_unit_write,'("E_PHI_FP: ",E17.10)') E_PHI do tt=1_ip,params%t_skip call include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & - V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,PSIp) + V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,Te,Zeff,nimp,PSIp) ! write(output_unit_write,'("Collision Loop in FP")') @@ -8483,6 +9171,147 @@ subroutine GCEoM1_p(pchunk,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, & end subroutine GCEoM1_p +subroutine GCEoM1_p_ACC(pp,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, & + B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,V_PLL,V_MU,Y_R,Y_PHI,Y_Z,q_cache,m_cache,PSIp,flag_cache) +!$acc routine seq + IMPLICIT NONE + + TYPE(KORC_PARAMS), INTENT(INOUT) :: params + !! Core KORC simulation parameters. + TYPE(FIELDS), INTENT(IN) :: F + TYPE(PROFILES), INTENT(IN) :: P + REAL(rp) :: Bmag,bhat_R,bhat_PHI,bhat_Z,Bst_R,Bst_PHI + REAL(rp) :: BstdotE,BstdotgradB,EcrossB_R,EcrossB_PHI,bdotBst + REAL(rp) :: bcrossgradB_R,bcrossgradB_PHI,bcrossgradB_Z,gamgc + REAL(rp) :: EcrossB_Z,Bst_Z + REAL(rp) :: pm,xi,tau_R + REAL(rp) :: SR_PLL,SR_MU,BREM_PLL,BREM_MU,BREM_P + REAL(rp),INTENT(in) :: gradB_R,gradB_PHI,gradB_Z,curlb_R + REAL(rp),INTENT(in) :: curlb_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z + REAL(rp),INTENT(OUT) :: RHS_R,RHS_PHI,RHS_Z + REAL(rp),INTENT(OUT) :: RHS_PLL,RHS_MU + REAL(rp),INTENT(IN) :: V_PLL,V_MU,Y_R,Y_PHI,Y_Z,curlb_PHI + REAL(rp),INTENT(IN) :: PSIp + REAL(rp),INTENT(in) :: q_cache,m_cache + INTEGER :: cc + INTEGER(ip),INTENT(IN) :: tt + INTEGER,INTENT(IN) :: pp + INTEGER(is),INTENT(OUT) :: flag_cache + INTEGER(is) :: flagCon + REAL(rp) :: time,re_cache,alpha_cache + INTEGER :: thread_num + REAL(rp) :: Zeff,Te,ne,nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1 + + !$acc routine (analytical_profiles_p_ACC) seq + !$acc routine (analytical_profiles_p_ACC_1) seq + !$acc routine (interp_Hcollision_p_ACC) seq + + ne=-1._rp + Te=-1._rp + Zeff=-1._rp + + Bmag = SQRT(B_R*B_R+B_PHI*B_PHI+B_Z*B_Z) + + bhat_R = B_R/Bmag + bhat_PHI = B_PHI/Bmag + bhat_Z = B_Z/Bmag + + Bst_R=q_cache*B_R+V_PLL*curlb_R + Bst_PHI=q_cache*B_PHI+V_PLL*curlb_PHI + Bst_Z=q_cache*B_Z+V_PLL*curlb_Z + + !write(6,*) pp,'bmag',Bmag,'bhat',bhat_R,bhat_PHI,bhat_Z,'Bst',Bst_R,Bst_PHI,Bst_Z + !flush(6) + + bdotBst=bhat_R*Bst_R+bhat_PHI*Bst_PHI+ & + bhat_Z*Bst_Z + BstdotE=Bst_R*E_R+Bst_PHI*E_PHI+Bst_Z*E_Z + BstdotgradB=Bst_R*gradB_R+Bst_PHI*gradB_PHI+ & + Bst_Z*gradB_Z + + !write(output_unit_write,*) 'bdotBst',bdotBst,BstdotE,BstdotgradB + + Ecrossb_R=E_PHI*bhat_Z-E_Z*bhat_PHI + Ecrossb_PHI=E_Z*bhat_R-E_R*bhat_Z + Ecrossb_Z=E_R*bhat_PHI-E_PHI*bhat_R + + !write(output_unit_write,*) 'Ecrossb',Ecrossb_R,Ecrossb_PHI,Ecrossb_Z + + bcrossgradB_R=bhat_PHI*gradB_Z-bhat_Z*gradB_PHI + bcrossgradB_PHI=bhat_Z*gradB_R-bhat_R*gradB_Z + bcrossgradB_Z=bhat_R*gradB_PHI-bhat_PHI*gradB_R + + ! write(output_unit_write,*) 'bcrossgradB',bcrossgradB_R,bcrossgradB_PHI,bcrossgradB_Z + + gamgc=sqrt(1+V_PLL*V_PLL+2*V_MU*Bmag) + + pm=sqrt(gamgc**2-1) + xi=V_PLL/pm + + RHS_R=(q_cache*Ecrossb_R+(m_cache*V_MU* & + bcrossgradB_R+V_PLL*Bst_R)/(m_cache*gamgc))/ & + bdotBst + RHS_PHI=(q_cache*Ecrossb_PHI+(m_cache*V_MU* & + bcrossgradB_PHI+V_PLL*Bst_PHI)/(m_cache*gamgc))/ & + (Y_R*bdotBst) + RHS_Z=(q_cache*Ecrossb_Z+(m_cache*V_MU* & + bcrossgradB_Z+V_PLL*Bst_Z)/(m_cache*gamgc))/ & + bdotBst + RHS_PLL=(q_cache*BstdotE-V_MU*BstdotgradB/gamgc)/ & + bdotBst + RHS_MU=0._rp + + if (params%radiation.and.(params%GC_rad_model.eq.'SDE')) then + + ! write(output_unit_write,*) 'RHS_PLL',RHS_PLL(1) + + re_cache=C_RE/params%cpp%length + alpha_cache=C_a + + if (.not.params%LargeCollisions) then + time=params%init_time+(params%it-1+tt)*params%dt + else + time=params%init_time+params%it*params%dt+ & + tt*params%coll_per_dump_dt + end if + + if (params%profile_model(1:10).eq.'ANALYTICAL') then + call analytical_profiles_p_ACC_1(time,params,Y_R,Y_Z, & + P,F,ne,Te,Zeff,PSIp) +#ifdef PSPLINE + else if (params%profile_model(10:10).eq.'H') then + call interp_Hcollision_p_ACC(params,Y_R,Y_PHI,Y_Z, & + ne,Te,Zeff,nAr0,nAr1,nAr2,nAr3,nAr4,nD,nD1, & + flagCon) +#endif + endif + + !write(6,*) 'Y: ',Y_R,Y_PHI,Y_Z + !write(6,*) 'Profs: ',ne,Te,Zeff + + tau_R=6*C_PI*E0/(Bmag*Bmag) + + SR_PLL=V_PLL*(1._rp-xi*xi)/tau_R* & + (1._rp/gamgc-gamgc) + SR_MU=-2._rp*V_MU/tau_R* & + (gamgc*(1-xi*xi)+xi*xi/gamgc) + + !Normalizations done here + BREM_P=-4._rp*re_cache**2*ne* & + Zeff*(Zeff+1._rp)*alpha_cache* & + (gamgc-1._rp)*(log(2._rp*gamgc)-1._rp/3._rp) + BREM_PLL=xi*BREM_P + BREM_MU=(1._rp-xi*xi)*V_PLL/ & + (Bmag*xi)*BREM_P + + RHS_PLL=RHS_PLL+SR_PLL+BREM_PLL + RHS_MU=SR_MU+BREM_MU + + end if + +end subroutine GCEoM1_p_ACC + #ifdef FIO subroutine GCEoM1_fio_p(tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, & diff --git a/src/korc_profiles.f90 b/src/korc_profiles.f90 index 0d3104d4..c5cab692 100755 --- a/src/korc_profiles.f90 +++ b/src/korc_profiles.f90 @@ -255,8 +255,15 @@ subroutine initialize_profiles(params,P,F) P%filename = TRIM(filename) P%axisymmetric = axisymmetric + P%ReInterp_2x1t = F%ReInterp_2x1t + P%ind_2x1t=F%ind_2x1t call load_profiles_data_from_hdf5(params,P) + + !write(6,*) 'profiles ne',P%ne_3D(13,1,16) + !write(6,*) 'Te(:,1,:)',P%Te_3D(:,1,:)/C_E + !write(6,*) 'Te(:,1,:)',P%Zeff_3D(:,1,:) + else if (params%profile_model.eq.'UNIFORM') then !open(unit=default_unit_open,file=TRIM(params%path_to_inputs), & ! status='OLD',form='formatted') @@ -313,6 +320,7 @@ subroutine uniform_profiles(vars,P) vars%Zeff = P%Zeffo end subroutine uniform_profiles + subroutine analytical_profiles_p(pchunk,time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp) !! @note Subroutine that calculates the analytical plasma profiles at !! the particles' position. @endnote @@ -337,7 +345,7 @@ subroutine analytical_profiles_p(pchunk,time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp) REAL(rp) :: R0_RE,Z0_RE,sigmaR_RE,sigmaZ_RE,psimax_RE REAL(rp) :: n_REr0,n_tauion,n_lamfront,n_lamback,n_lamshelf REAL(rp) :: n_psifront,n_psiback,n_psishelf - REAL(rp) :: n_tauin,n_tauout,n_shelfdelay,n_shelf + REAL(rp) :: n_tauin,n_tauout,n_shelfdshelf REAL(rp) :: n0t,n_taut REAL(rp) :: PSIp0,PSIp_lim,psiN_0 REAL(rp), DIMENSION(params%pchunk) :: r_a,rm,rm_RE,PSIpN,PSIp_temp @@ -603,6 +611,264 @@ subroutine analytical_profiles_p(pchunk,time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp) end subroutine analytical_profiles_p + subroutine analytical_profiles_p_ACC(time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp) + !$acc routine seq + !! @note Subroutine that calculates the analytical plasma profiles at + !! the particles' position. @endnote + TYPE(KORC_PARAMS), INTENT(IN) :: params + REAL(rp), INTENT(IN) :: Y_R,Y_Z,PSIp + REAL(rp), INTENT(IN) :: time + TYPE(PROFILES), INTENT(IN) :: P + !! An instance of KORC's derived type PROFILES containing all the + !! information about the plasma profiles used in the simulation. + !! See [[korc_types]] and [[korc_profiles]]. + TYPE(FIELDS), INTENT(IN) :: F + REAL(rp),INTENT(OUT) :: ne + !! Background electron density seen by simulated particles. + REAL(rp),INTENT(OUT) :: Te + !! Backgroun temperature density seen by simulated particles. + REAL(rp),INTENT(OUT) :: Zeff + !! Effective atomic charge seen by simulated particles. + INTEGER(ip) :: cc + !! Particle iterator. + REAL(rp) :: R0,Z0,a,ne0,n_ne,Te0,n_Te,Zeff0,R0a + REAL(rp) :: R0_RE,Z0_RE,sigmaR_RE,sigmaZ_RE,psimax_RE + REAL(rp) :: n_REr0,n_tauion,n_lamfront,n_lamback,n_lamshelf + REAL(rp) :: n_psifront,n_psiback,n_psishelf + REAL(rp) :: n_tauin,n_tauout,n_shelfdshelf + REAL(rp) :: n0t,n_taut + REAL(rp) :: PSIp0,PSIp_lim,psiN_0 + REAL(rp) :: r_a,rm,rm_RE,PSIpN,PSIp_temp + + R0=P%R0 + Z0=P%Z0 + a=P%a + R0a=F%AB%Ro + + ne0=P%neo + n_ne=P%n_ne + + Te0=P%Teo + n_Te=P%n_Te + + Zeff0=P%Zeffo + + R0_RE=P%R0_RE + Z0_RE=P%Z0_RE + n_REr0=P%n_REr0 + n_tauion=P%n_tauion + n_tauin=P%n_tauin + n_tauout=P%n_tauout + n_shelfdelay=P%n_shelfdelay + n_lamfront=P%n_lamfront + n_lamback=P%n_lamback + n_lamshelf=P%n_lamshelf + n_psifront=P%n_lamfront*params%cpp%length + n_psiback=P%n_lamback*params%cpp%length + n_psishelf=P%n_lamshelf*params%cpp%length + n_shelf=P%n_shelf + + PSIp_lim=F%PSIp_lim + PSIp0=F%PSIP_min + psiN_0=P%psiN_0 + + +! write(output_unit_write,*) 'PSIp',PSIp(1)*(params%cpp%Bo*params%cpp%length**2) +! write(output_unit_write,*) 'PSIp_lim',PSIp_lim*(params%cpp%Bo*params%cpp%length**2) +! write(output_unit_write,*) 'PSIp0',PSIp0*(params%cpp%Bo*params%cpp%length**2) + +! write(output_unit_write,'("R0_RE: "E17.10)') R0_RE +! write(output_unit_write,'("Z0_RE: "E17.10)') Z0_RE +! write(output_unit_write,'("n_REr0: "E17.10)') n_REr0 + + + SELECT CASE (TRIM(P%ne_profile)) + CASE('FLAT') + ne = ne0 + CASE('FLAT-RAMP') + ne = n_ne+(ne0-n_ne)*time/n_tauion + CASE('TANH-RAMP') + ne = n_ne+(ne0-n_ne)/2*(tanh((time-n_shelfdelay)/n_tauion)+1._rp) + CASE('SINE') + ne = n_ne+(ne0-n_ne)*sin(time/n_tauion) + CASE('SPONG') + rm=sqrt((Y_R-R0)**2+(Y_Z-Z0)**2) + r_a=rm/a + ne = ne0*(1._rp-0.2*r_a**8)+n_ne + CASE('MST_FSA') + rm=sqrt((Y_R-R0a)**2+(Y_Z-Z0)**2) + r_a=rm/a + ne = (ne0-n_ne)*(1._rp-r_a**4._rp)**4._rp+n_ne + CASE('RE-EVO') + rm_RE=sqrt((Y_R-R0_RE)**2+(Y_Z-Z0_RE)**2) + ne = (ne0-n_ne)/4._rp*(1+tanh((rm_RE+ & + n_REr0*(time/n_tauion-1))/n_lamfront))* & + (1+tanh(-(rm_RE-n_REr0)/n_lamback))+n_ne + CASE('RE-EVO1') + rm_RE=sqrt((Y_R-R0_RE)**2+(Y_Z-Z0_RE)**2) + ne = (ne0-n_ne)/8._rp*(1+tanh((rm_RE+ & + n_REr0*(time/n_tauion-1))/n_lamfront))* & + (1+tanh(-(rm_RE-n_REr0)/n_lamback))* & + (2*(n_shelf-n_ne)/(ne0-n_ne)+(ne0-n_shelf)/(ne0-n_ne)* & + (1-tanh((rm_RE+n_REr0*((time-n_shelfdelay)/n_tauin-1))/ & + n_lamshelf)))+n_ne + CASE('RE-EVO-PSI') + PSIpN=(PSIp-PSIp0)/(PSIp_lim-PSIp0) + ne = (ne0-n_ne)/8._rp*(1+tanh((sqrt(abs(PSIpN))+ & + sqrt(abs(psiN_0))*(time/n_tauion-1))/n_psifront))* & + (1+tanh(-(sqrt(abs(PSIpN))-sqrt(abs(psiN_0)))/n_psiback))* & + (2*(n_shelf-n_ne)/(ne0-n_ne)+(ne0-n_shelf)/(ne0-n_ne)* & + (1-tanh((sqrt(abs(PSIpN))+ sqrt(abs(psiN_0))* & + ((time-n_shelfdelay)/n_tauin-1))/n_psishelf)))+n_ne + CASE('RE-EVO-PSIN-SG') + + n0t=(ne0-n_ne)/2._rp*(tanh(time/n_tauin)- & + tanh((time-n_shelfdelay)/n_tauin)) + n_taut=n_psishelf*erf((time+params%dt/100._rp)/n_tauion) + + PSIpN=(PSIp-PSIp0)/(PSIp_lim-PSIp0) + ne = n0t*exp(-(sqrt(abs(PSIpN))-sqrt(abs(psiN_0)))**2._rp/ & + (2._rp*n_taut**2._rp))*(1._rp+erf(-10._rp* & + (sqrt(abs(PSIpN))-sqrt(abs(psiN_0)))/ & + (sqrt(2._rp)*n_taut)))/2._rp+n_ne + + CASE('RE-EVO-PSIP-G') + + n0t=(ne0-n_ne)/2._rp*(tanh((time-n_tauin)/n_tauin)- & + tanh((time-n_shelfdelay)/n_tauout)) + n_taut=n_psishelf*erf((time+params%dt/100._rp)/n_tauion) + + PSIp_temp=PSIp*(params%cpp%Bo*params%cpp%length**2) + ne = n0t*exp(-(sqrt(abs(PSIp_temp))-sqrt(abs(psiN_0)))**2._rp/ & + (2._rp*n_taut**2._rp))+n_ne + + CASE('RE-EVO-PSIP-G1') + +! write(output_unit_write,*) 'time: ',time*params%cpp%time + + n0t=(ne0-n_ne)/2._rp*(tanh((time-1.75*n_tauin)/n_tauin)- & + tanh((time-n_shelfdelay)/n_tauout)) + n_taut=n_psishelf*erf((time+params%dt/100._rp)/n_tauion) + + PSIp_temp=PSIp*(params%cpp%Bo*params%cpp%length**2) + ne = n0t*exp(-(sqrt(abs(PSIp_temp))-sqrt(abs(psiN_0)))**2._rp/ & + (2._rp*n_taut**2._rp))+n_ne + + CASE DEFAULT + ne = ne0 + END SELECT + + SELECT CASE (TRIM(P%Te_profile)) + CASE('FLAT') + Te = Te0 + CASE('SPONG') + rm=sqrt((Y_R-R0)**2+(Y_Z-Z0)**2) + r_a=rm/a + Te = Te0*(1._rp-0.6*r_a**2)**2+Te0*n_Te + CASE('MST_FSA') + + rm=sqrt((Y_R-R0a)**2+(Y_Z-Z0)**2) + r_a=rm/a + Te = (Te0-n_Te)*(1._rp-r_a**8._rp)**4._rp+n_Te + + CASE DEFAULT + + Te = P%Teo + + END SELECT + + SELECT CASE (TRIM(P%Zeff_profile)) + CASE('FLAT') + + Zeff = P%Zeffo + + CASE('SPONG') + + Zeff = P%Zeffo + + CASE DEFAULT + + Zeff = P%Zeffo + + END SELECT + + end subroutine analytical_profiles_p_ACC + + subroutine analytical_profiles_p_ACC_1(time,params,Y_R,Y_Z,P,F,ne,Te,Zeff,PSIp) + !$acc routine seq + !! @note Subroutine that calculates the analytical plasma profiles at + !! the particles' position. @endnote + TYPE(KORC_PARAMS), INTENT(IN) :: params + REAL(rp), INTENT(IN) :: Y_R,Y_Z,PSIp + REAL(rp), INTENT(IN) :: time + TYPE(PROFILES), INTENT(IN) :: P + !! An instance of KORC's derived type PROFILES containing all the + !! information about the plasma profiles used in the simulation. + !! See [[korc_types]] and [[korc_profiles]]. + TYPE(FIELDS), INTENT(IN) :: F + REAL(rp),INTENT(OUT) :: ne + !! Background electron density seen by simulated particles. + REAL(rp),INTENT(OUT) :: Te + !! Backgroun temperature density seen by simulated particles. + REAL(rp),INTENT(OUT) :: Zeff + !! Effective atomic charge seen by simulated particles. + INTEGER(ip) :: cc + !! Particle iterator. + REAL(rp) :: R0,Z0,a,ne0,n_ne,Te0,n_Te,Zeff0,R0a + REAL(rp) :: R0_RE,Z0_RE,sigmaR_RE,sigmaZ_RE,psimax_RE + REAL(rp) :: n_REr0,n_tauion,n_lamfront,n_lamback,n_lamshelf + REAL(rp) :: n_psifront,n_psiback,n_psishelf + REAL(rp) :: n_tauin,n_tauout,n_shelfdshelf + REAL(rp) :: n0t,n_taut + REAL(rp) :: PSIp0,PSIp_lim,psiN_0 + REAL(rp) :: r_a,rm,rm_RE,PSIpN,PSIp_temp + + R0=P%R0 + Z0=P%Z0 + a=P%a + R0a=F%AB%Ro + + ne0=P%neo + n_ne=P%n_ne + + Te0=P%Teo + n_Te=P%n_Te + + Zeff0=P%Zeffo + + R0_RE=P%R0_RE + Z0_RE=P%Z0_RE + n_REr0=P%n_REr0 + n_tauion=P%n_tauion + n_tauin=P%n_tauin + n_tauout=P%n_tauout + n_shelfdelay=P%n_shelfdelay + n_lamfront=P%n_lamfront + n_lamback=P%n_lamback + n_lamshelf=P%n_lamshelf + n_psifront=P%n_lamfront*params%cpp%length + n_psiback=P%n_lamback*params%cpp%length + n_psishelf=P%n_lamshelf*params%cpp%length + n_shelf=P%n_shelf + + PSIp_lim=F%PSIp_lim + PSIp0=F%PSIP_min + psiN_0=P%psiN_0 + + n0t=(ne0-n_ne)/2._rp*(tanh((time-1.75*n_tauin)/n_tauin)- & + tanh((time-n_shelfdelay)/n_tauout)) + n_taut=n_psishelf*erf((time+params%dt/100._rp)/n_tauion) + + PSIp_temp=PSIp*(params%cpp%Bo*params%cpp%length**2) + ne = n0t*exp(-(sqrt(abs(PSIp_temp))-sqrt(abs(psiN_0)))**2._rp/ & + (2._rp*n_taut**2._rp))+n_ne + + Te = Te0 + + Zeff = P%Zeffo + + end subroutine analytical_profiles_p_ACC_1 + subroutine get_analytical_profiles(P,Y,ne,Te,Zeff,flag) !! @note Subroutine that calculates the analytical plasma profiles at !! the particles' position. @endnote @@ -790,7 +1056,7 @@ subroutine load_profiles_data_from_hdf5(params,P) call load_from_hdf5(h5file_id,dset,rdatum) P%dims(1) = INT(rdatum) - if (P%axisymmetric) then + if (P%axisymmetric.and.(.not.P%ReInterp_2x1t)) then P%dims(2) = 0 else dset = "/NPHI" @@ -802,16 +1068,16 @@ subroutine load_profiles_data_from_hdf5(params,P) call load_from_hdf5(h5file_id,dset,rdatum) P%dims(3) = INT(rdatum) - if (P%axisymmetric) then + if (P%axisymmetric.and.(.not.P%ReInterp_2x1t)) then call ALLOCATE_2D_PROFILES_ARRAYS(params,P) else - call ALLOCATE_3D_PROFILES_ARRAYS(P) + call ALLOCATE_3D_PROFILES_ARRAYS(params,P) end if dset = "/R" call load_array_from_hdf5(h5file_id,dset,P%X%R) - if (.NOT.P%axisymmetric) then + if (.NOT.P%axisymmetric.and.(.not.P%ReInterp_2x1t)) then dset = "/PHI" call load_array_from_hdf5(h5file_id,dset,P%X%PHI) end if @@ -820,21 +1086,21 @@ subroutine load_profiles_data_from_hdf5(params,P) call load_array_from_hdf5(h5file_id,dset,P%X%Z) dset = "/FLAG" - if (P%axisymmetric) then + if (P%axisymmetric.and.(.not.P%ReInterp_2x1t)) then call load_array_from_hdf5(h5file_id,dset,P%FLAG2D) else call load_array_from_hdf5(h5file_id,dset,P%FLAG3D) end if dset = "/ne" - if (P%axisymmetric) then + if (P%axisymmetric.and.(.not.P%ReInterp_2x1t)) then call load_array_from_hdf5(h5file_id,dset,P%ne_2D) else call load_array_from_hdf5(h5file_id,dset,P%ne_3D) end if dset = "/Te" - if (P%axisymmetric) then + if (P%axisymmetric.and.(.not.P%ReInterp_2x1t)) then call load_array_from_hdf5(h5file_id,dset,P%Te_2D) P%Te_2D = P%Te_2D*C_E else @@ -845,31 +1111,52 @@ subroutine load_profiles_data_from_hdf5(params,P) !write(output_unit_write,'("Te: ",E17.10)') P%Te_2D(1,1) dset = "/Zeff" - if (P%axisymmetric) then + if (P%axisymmetric.and.(.not.P%ReInterp_2x1t)) then call load_array_from_hdf5(h5file_id,dset,P%Zeff_2D) else call load_array_from_hdf5(h5file_id,dset,P%Zeff_3D) end if if (params%profile_model(10:10).eq.'H') then - - dset = "/RHON" - call load_array_from_hdf5(h5file_id,dset,P%RHON) - dset = "/nRE" - call load_array_from_hdf5(h5file_id,dset,P%nRE_2D) - dset = "/nAr0" - call load_array_from_hdf5(h5file_id,dset,P%nAr0_2D) - dset = "/nAr1" - call load_array_from_hdf5(h5file_id,dset,P%nAr1_2D) - dset = "/nAr2" - call load_array_from_hdf5(h5file_id,dset,P%nAr2_2D) - dset = "/nAr3" - call load_array_from_hdf5(h5file_id,dset,P%nAr3_2D) - dset = "/nD" - call load_array_from_hdf5(h5file_id,dset,P%nD_2D) - dset = "/nD1" - call load_array_from_hdf5(h5file_id,dset,P%nD1_2D) - + if (P%axisymmetric.and.(.not.P%ReInterp_2x1t)) then + dset = "/RHON" + call load_array_from_hdf5(h5file_id,dset,P%RHON_2D) + dset = "/nRE" + call load_array_from_hdf5(h5file_id,dset,P%nRE_2D) + dset = "/nAr0" + call load_array_from_hdf5(h5file_id,dset,P%nAr0_2D) + dset = "/nAr1" + call load_array_from_hdf5(h5file_id,dset,P%nAr1_2D) + dset = "/nAr2" + call load_array_from_hdf5(h5file_id,dset,P%nAr2_2D) + dset = "/nAr3" + call load_array_from_hdf5(h5file_id,dset,P%nAr3_2D) + dset = "/nAr4" + call load_array_from_hdf5(h5file_id,dset,P%nAr4_2D) + dset = "/nD" + call load_array_from_hdf5(h5file_id,dset,P%nD_2D) + dset = "/nD1" + call load_array_from_hdf5(h5file_id,dset,P%nD1_2D) + else + dset = "/RHON" + call load_array_from_hdf5(h5file_id,dset,P%RHON_3D) + dset = "/nRE" + call load_array_from_hdf5(h5file_id,dset,P%nRE_2D) + dset = "/nAr0" + call load_array_from_hdf5(h5file_id,dset,P%nAr0_3D) + dset = "/nAr1" + call load_array_from_hdf5(h5file_id,dset,P%nAr1_3D) + dset = "/nAr2" + call load_array_from_hdf5(h5file_id,dset,P%nAr2_3D) + dset = "/nAr3" + call load_array_from_hdf5(h5file_id,dset,P%nAr3_3D) + dset = "/nAr4" + call load_array_from_hdf5(h5file_id,dset,P%nAr4_3D) + dset = "/nD" + call load_array_from_hdf5(h5file_id,dset,P%nD_3D) + dset = "/nD1" + call load_array_from_hdf5(h5file_id,dset,P%nD1_3D) + endif end if call h5fclose_f(h5file_id, h5error) @@ -895,21 +1182,23 @@ subroutine ALLOCATE_2D_PROFILES_ARRAYS(params,P) ALLOCATE(P%Zeff_2D(P%dims(1),P%dims(3))) if (params%profile_model(10:10).eq.'H') then - ALLOCATE(P%RHON(P%dims(1),P%dims(3))) + ALLOCATE(P%RHON_2D(P%dims(1),P%dims(3))) ALLOCATE(P%nRE_2D(P%dims(1),P%dims(3))) ALLOCATE(P%nAr0_2D(P%dims(1),P%dims(3))) ALLOCATE(P%nAr1_2D(P%dims(1),P%dims(3))) ALLOCATE(P%nAr2_2D(P%dims(1),P%dims(3))) ALLOCATE(P%nAr3_2D(P%dims(1),P%dims(3))) + ALLOCATE(P%nAr4_2D(P%dims(1),P%dims(3))) ALLOCATE(P%nD_2D(P%dims(1),P%dims(3))) ALLOCATE(P%nD1_2D(P%dims(1),P%dims(3))) end if end subroutine ALLOCATE_2D_PROFILES_ARRAYS - subroutine ALLOCATE_3D_PROFILES_ARRAYS(P) + subroutine ALLOCATE_3D_PROFILES_ARRAYS(params,P) !! @note Subroutine that allocates the mesh information and 3-D arrays !! for keeping the data of pre-computed plasma profiles. @endnote + TYPE(KORC_PARAMS), INTENT(IN) :: params TYPE(PROFILES), INTENT(INOUT) :: P !! @param[out] P An instance of KORC's derived type PROFILES containing !! all the information about the plasma profiles used in the @@ -921,6 +1210,18 @@ subroutine ALLOCATE_3D_PROFILES_ARRAYS(P) ALLOCATE(P%ne_3D(P%dims(1),P%dims(2),P%dims(3))) ALLOCATE(P%Te_3D(P%dims(1),P%dims(2),P%dims(3))) ALLOCATE(P%Zeff_3D(P%dims(1),P%dims(2),P%dims(3))) + + if (params%profile_model(10:10).eq.'H') then + ALLOCATE(P%RHON_3D(P%dims(1),P%dims(2),P%dims(3))) + ALLOCATE(P%nRE_2D(P%dims(1),P%dims(3))) + ALLOCATE(P%nAr0_3D(P%dims(1),P%dims(2),P%dims(3))) + ALLOCATE(P%nAr1_3D(P%dims(1),P%dims(2),P%dims(3))) + ALLOCATE(P%nAr2_3D(P%dims(1),P%dims(2),P%dims(3))) + ALLOCATE(P%nAr3_3D(P%dims(1),P%dims(2),P%dims(3))) + ALLOCATE(P%nAr4_3D(P%dims(1),P%dims(2),P%dims(3))) + ALLOCATE(P%nD_3D(P%dims(1),P%dims(2),P%dims(3))) + ALLOCATE(P%nD1_3D(P%dims(1),P%dims(2),P%dims(3))) + end if end subroutine ALLOCATE_3D_PROFILES_ARRAYS subroutine DEALLOCATE_PROFILES_ARRAYS(P) @@ -940,14 +1241,24 @@ subroutine DEALLOCATE_PROFILES_ARRAYS(P) if (ALLOCATED(P%Te_3D)) DEALLOCATE(P%Te_3D) if (ALLOCATED(P%Zeff_3D)) DEALLOCATE(P%Zeff_3D) - if (ALLOCATED(P%RHON)) DEALLOCATE(P%RHON) + if (ALLOCATED(P%RHON_2D)) DEALLOCATE(P%RHON_2D) if (ALLOCATED(P%nRE_2D)) DEALLOCATE(P%nRE_2D) if (ALLOCATED(P%nAr0_2D)) DEALLOCATE(P%nAr0_2D) if (ALLOCATED(P%nAr1_2D)) DEALLOCATE(P%nAr1_2D) if (ALLOCATED(P%nAr2_2D)) DEALLOCATE(P%nAr2_2D) if (ALLOCATED(P%nAr3_2D)) DEALLOCATE(P%nAr3_2D) + if (ALLOCATED(P%nAr4_2D)) DEALLOCATE(P%nAr4_2D) if (ALLOCATED(P%nD_2D)) DEALLOCATE(P%nD_2D) if (ALLOCATED(P%nD1_2D)) DEALLOCATE(P%nD1_2D) + + if (ALLOCATED(P%RHON_3D)) DEALLOCATE(P%RHON_3D) + if (ALLOCATED(P%nAr0_3D)) DEALLOCATE(P%nAr0_3D) + if (ALLOCATED(P%nAr1_3D)) DEALLOCATE(P%nAr1_3D) + if (ALLOCATED(P%nAr2_3D)) DEALLOCATE(P%nAr2_3D) + if (ALLOCATED(P%nAr3_3D)) DEALLOCATE(P%nAr3_3D) + if (ALLOCATED(P%nAr4_3D)) DEALLOCATE(P%nAr4_3D) + if (ALLOCATED(P%nD_3D)) DEALLOCATE(P%nD_3D) + if (ALLOCATED(P%nD1_3D)) DEALLOCATE(P%nD1_3D) #ifdef FIO if (ALLOCATED(P%FIO_nimp)) DEALLOCATE(P%FIO_nimp) diff --git a/src/korc_pspline.f90 b/src/korc_pspline.f90 index 027d8d47..31c27baa 100644 --- a/src/korc_pspline.f90 +++ b/src/korc_pspline.f90 @@ -2800,6 +2800,49 @@ subroutine EZspline_interp2_collision(spline_oBR, spline_oBPHI, & end subroutine EZspline_interp2_collision +subroutine EZspline_interp2_Hcollision(spline_oBR, spline_oBPHI, & + spline_oBZ, spline_oER, spline_oEPHI, spline_oEZ,spline_oEZ1, & + spline_one, spline_oTe, spline_oZeff, p1, p2, fBR, & + fBPHI, fBZ, fER, fEPHI, fEZ, fEZ1, fne, fTe, fZeff, ier) + !$acc routine seq + type(EZspline2) spline_oBR,spline_oBPHI,spline_oBZ + type(EZspline2) spline_oER,spline_oEPHI,spline_oEZ,spline_oEZ1 + type(EZspline2) spline_one,spline_oTe,spline_oZeff + real(fp), intent(in) :: p1, p2 + real(fp), intent(out):: fBR, fBPHI, fBZ + real(fp), intent(out):: fER, fEPHI, fEZ, fEZ1 + real(fp), intent(out):: fne, fTe, fZeff + integer, intent(out) :: ier + integer :: ifail + integer:: iwarn = 0 + + !$acc routine (EZspline_allocated2) seq + !$acc routine (evbicub_Hcollision) seq + + ier = 0 + ifail = 0 + + if( .not.EZspline_allocated2(spline_oBR) .or. spline_oBR%isReady /= 1) then + ier = 94 + return + endif + + call evbicub_Hcollision(p1, p2, & + spline_oBR%x1(1), spline_oBR%n1, & + spline_oBR%x2(1), spline_oBR%n2, & + spline_oBR%ilin1, spline_oBR%ilin2, & + spline_oBR%fspl(1,1,1), spline_oBPHI%fspl(1,1,1), & + spline_oBZ%fspl(1,1,1), spline_oER%fspl(1,1,1), & + spline_oEPHI%fspl(1,1,1), spline_oEZ%fspl(1,1,1), spline_oEZ1%fspl(1,1,1), & + spline_one%fspl(1,1,1), & + spline_oTe%fspl(1,1,1), spline_oZeff%fspl(1,1,1), & + spline_oBR%n1, & + fBR, fBPHI, fBZ, fER, fEPHI, fEZ, fEZ1, fne, fTe, fZeff, ifail) + + if(ifail /= 0) ier = 97 + +end subroutine EZspline_interp2_Hcollision + subroutine EZspline_interp2_FOmars(spline_oA, spline_oBR, spline_oBPHI, & spline_oBZ, spline_oER, spline_oEPHI, spline_oEZ, p1, p2, fA, fBR, & fBPHI, fBZ, fER, fEPHI, fEZ, ier) @@ -3619,6 +3662,166 @@ subroutine evbicub_collision(xget,yget,x,nx,y,ny,ilinx,iliny, & return end subroutine evbicub_collision +subroutine evbicub_Hcollision(xget,yget,x,nx,y,ny,ilinx,iliny, & + fBR,fBPHI,fBZ,fER,fEPHI,fEZ,fEZ1,fne,fTe,fZeff,inf2,fvalBR,fvalBPHI,fvalBZ, & + fvalER,fvalEPHI,fvalEZ,fvalEZ1,fvalne,fvalTe,fvalZeff,ier) + !$acc routine seq + ! + ! evaluate a 2d cubic Spline interpolant on a rectilinear + ! grid -- this is C2 in both directions. + ! + ! this subroutine calls two subroutines: + ! herm2xy -- find cell containing (xget,yget) + ! fvbicub -- evaluate interpolant function and (optionally) derivatives + ! + ! input arguments: + ! ================ + ! + !============ + implicit none + integer inf2 + !============ + integer,intent(in) :: nx,ny ! grid sizes + real(fp) :: xget,yget ! target of this interpolation + real(fp) :: x(nx) ! ordered x grid + real(fp) :: y(ny) ! ordered y grid + integer ilinx ! ilinx=1 => assume x evenly spaced + integer iliny ! iliny=1 => assume y evenly spaced + ! + real(fp) :: fBR(0:3,inf2,ny) ! function data + real(fp) :: fBPHI(0:3,inf2,ny) + real(fp) :: fBZ(0:3,inf2,ny) + real(fp) :: fER(0:3,inf2,ny) + real(fp) :: fEPHI(0:3,inf2,ny) + real(fp) :: fEZ(0:3,inf2,ny) + real(fp) :: fEZ1(0:3,inf2,ny) + real(fp) :: fne(0:3,inf2,ny) + real(fp) :: fTe(0:3,inf2,ny) + real(fp) :: fZeff(0:3,inf2,ny) + ! + ! f 2nd dimension inf2 must be .ge. nx + ! contents of f: + ! + ! f(0,i,j) = f @ x(i),y(j) + ! f(1,i,j) = d2f/dx2 @ x(i),y(j) + ! f(2,i,j) = d2f/dy2 @ x(i),y(j) + ! f(3,i,j) = d4f/dx2dy2 @ x(i),y(j) + ! + ! (these are spline coefficients selected for continuous 2- + ! diffentiability, see mkbicub[w].f90) + ! + ! + ! ict(1)=1 -- return f (0, don't) + ! ict(2)=1 -- return df/dx (0, don't) + ! ict(3)=1 -- return df/dy (0, don't) + ! ict(4)=1 -- return d2f/dx2 (0, don't) + ! ict(5)=1 -- return d2f/dy2 (0, don't) + ! ict(6)=1 -- return d2f/dxdy (0, don't) + ! the number of non zero values ict(1:6) + ! determines the number of outputs... + ! + ! new dmc December 2005 -- access to higher derivatives (even if not + ! continuous-- but can only go up to 3rd derivatives on any one coordinate. + ! if ict(1)=3 -- want 3rd derivatives + ! ict(2)=1 for d3f/dx3 + ! ict(3)=1 for d3f/dx2dy + ! ict(4)=1 for d3f/dxdy2 + ! ict(5)=1 for d3f/dy3 + ! number of non-zero values ict(2:5) gives no. of outputs + ! if ict(1)=4 -- want 4th derivatives + ! ict(2)=1 for d4f/dx3dy + ! ict(3)=1 for d4f/dx2dy2 + ! ict(4)=1 for d4f/dxdy3 + ! number of non-zero values ict(2:4) gives no. of outputs + ! if ict(1)=5 -- want 5th derivatives + ! ict(2)=1 for d5f/dx3dy2 + ! ict(3)=1 for d5f/dx2dy3 + ! number of non-zero values ict(2:3) gives no. of outputs + ! if ict(1)=6 -- want 6th derivatives + ! d6f/dx3dy3 -- one value is returned. + ! + ! output arguments: + ! ================= + ! + real(fp) :: fvalBR ! output data + real(fp) :: fvalBPHI + real(fp) :: fvalBZ + real(fp) :: fvalER + real(fp) :: fvalEPHI + real(fp) :: fvalEZ + real(fp) :: fvalEZ1 + real(fp) :: fvalne + real(fp) :: fvalTe + real(fp) :: fvalZeff + + integer ier ! error code =0 ==> no error + ! + ! fval(1) receives the first output (depends on ict(...) spec) + ! fval(2) receives the second output (depends on ict(...) spec) + ! fval(3) receives the third output (depends on ict(...) spec) + ! fval(4) receives the fourth output (depends on ict(...) spec) + ! fval(5) receives the fourth output (depends on ict(...) spec) + ! fval(6) receives the fourth output (depends on ict(...) spec) + ! + ! examples: + ! on input ict = [1,1,1,0,0,1] + ! on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced. + ! + ! on input ict = [1,0,0,0,0,0] + ! on output fval= [f] ... elements 2 -- 6 never referenced. + ! + ! on input ict = [0,0,0,1,1,0] + ! on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced. + ! + ! on input ict = [0,0,1,0,0,0] + ! on output fval= [df/dy] ... elements 2 -- 6 never referenced. + ! + ! ier -- completion code: 0 means OK + !------------------- + ! local: + ! + integer i,j ! cell indices + ! + ! normalized displacement from (x(i),y(j)) corner of cell. + ! xparam=0 @x(i) xparam=1 @x(i+1) + ! yparam=0 @y(j) yparam=1 @y(j+1) + ! + real(fp) :: xparam,yparam + ! + ! cell dimensions and + ! inverse cell dimensions hxi = 1/(x(i+1)-x(i)), hyi = 1/(y(j+1)-y(j)) + ! + real(fp) :: hx,hy + real(fp) :: hxi,hyi + ! + ! 0 .le. xparam .le. 1 + ! 0 .le. yparam .le. 1 + ! + ! ** the interface is very similar to herm2ev.f90; can use herm2xy ** + !--------------------------------------------------------------------- + !$acc routine (herm2xy) seq + !$acc routine (fvbicub) seq + ! + i=0 + j=0 + call herm2xy(xget,yget,x,nx,y,ny,ilinx,iliny, & + i,j,xparam,yparam,hx,hxi,hy,hyi,ier) + if(ier.ne.0) return + ! + call fvbicub(fvalBR,i,j,xparam,yparam,hx,hxi,hy,hyi,fBR,inf2,ny) + call fvbicub(fvalBPHI,i,j,xparam,yparam,hx,hxi,hy,hyi,fBPHI,inf2,ny) + call fvbicub(fvalBZ,i,j,xparam,yparam,hx,hxi,hy,hyi,fBZ,inf2,ny) + call fvbicub(fvalER,i,j,xparam,yparam,hx,hxi,hy,hyi,fER,inf2,ny) + call fvbicub(fvalEPHI,i,j,xparam,yparam,hx,hxi,hy,hyi,fEPHI,inf2,ny) + call fvbicub(fvalEZ,i,j,xparam,yparam,hx,hxi,hy,hyi,fEZ,inf2,ny) + call fvbicub(fvalEZ1,i,j,xparam,yparam,hx,hxi,hy,hyi,fEZ1,inf2,ny) + call fvbicub(fvalne,i,j,xparam,yparam,hx,hxi,hy,hyi,fne,inf2,ny) + call fvbicub(fvalTe,i,j,xparam,yparam,hx,hxi,hy,hyi,fTe,inf2,ny) + call fvbicub(fvalZeff,i,j,xparam,yparam,hx,hxi,hy,hyi,fZeff,inf2,ny) + ! + return +end subroutine evbicub_Hcollision + subroutine evbicub_FOmars(xget,yget,x,nx,y,ny,ilinx,iliny, & fA,fBR,fBPHI,fBZ,fER,fEPHI,fEZ,inf2,fvalA,fvalBR,fvalBPHI,fvalBZ, & fvalER,fvalEPHI,fvalEZ,ier) diff --git a/src/korc_spatial_distribution.f90 b/src/korc_spatial_distribution.f90 index 4b935385..1d109353 100755 --- a/src/korc_spatial_distribution.f90 +++ b/src/korc_spatial_distribution.f90 @@ -1517,11 +1517,10 @@ subroutine MH_psi(params,random,spp,F) ! write(output_unit_write,'("sample:",I15)') ii -#if DBG_CHECK if (modulo(ii,nsamples/10).eq.0) then write(output_unit_write,'("Sample: ",I10)') ii end if -#endif + PHI_test = 2.0_rp*C_PI*random%uniform%get() @@ -2559,6 +2558,7 @@ subroutine intitial_spatial_distribution(params,random,spp,P,F) INTEGER :: ss !! Species iterator. INTEGER :: mpierr + INTEGER,DIMENSION(34) :: seed=(/1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1/) do ss=1_idef,params%num_species SELECT CASE (TRIM(spp(ss)%spatial_distribution)) @@ -2589,6 +2589,7 @@ subroutine intitial_spatial_distribution(params,random,spp,P,F) spp(ss)%vars%X(:,1)=spp(ss)%Xtrace(1) spp(ss)%vars%X(:,2)=spp(ss)%Xtrace(2) spp(ss)%vars%X(:,3)=spp(ss)%Xtrace(3) + CASE ('SPONG-3D') call Spong_3D(params,random,spp(ss)) CASE ('HOLLMANN-3D') diff --git a/src/korc_types.f90 b/src/korc_types.f90 index 17ae17d6..ee0b4fbf 100755 --- a/src/korc_types.f90 +++ b/src/korc_types.f90 @@ -811,14 +811,23 @@ module korc_types !! Full path to the HDF5 file containing the pre-computed plasma profiles. LOGICAL :: axisymmetric !! Flag to indicate if the plasma profiles are axisymmetric. - REAL(rp), DIMENSION(:,:), ALLOCATABLE :: RHON + REAL(rp), DIMENSION(:,:), ALLOCATABLE :: RHON_2D REAL(rp), DIMENSION(:,:), ALLOCATABLE :: nRE_2D REAL(rp), DIMENSION(:,:), ALLOCATABLE :: nAr0_2D REAL(rp), DIMENSION(:,:), ALLOCATABLE :: nAr1_2D REAL(rp), DIMENSION(:,:), ALLOCATABLE :: nAr2_2D REAL(rp), DIMENSION(:,:), ALLOCATABLE :: nAr3_2D + REAL(rp), DIMENSION(:,:), ALLOCATABLE :: nAr4_2D REAL(rp), DIMENSION(:,:), ALLOCATABLE :: nD_2D REAL(rp), DIMENSION(:,:), ALLOCATABLE :: nD1_2D + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: RHON_3D + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: nAr0_3D + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: nAr1_3D + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: nAr2_3D + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: nAr3_3D + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: nAr4_3D + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: nD_3D + REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: nD1_3D #ifdef FIO INTEGER (C_INT) :: FIO_ne @@ -827,6 +836,9 @@ module korc_types INTEGER (C_INT), DIMENSION(:), ALLOCATABLE :: FIO_nimp INTEGER (C_INT) :: FIO_zeff #endif + + LOGICAL :: ReInterp_2x1t + INTEGER :: ind_2x1t END TYPE PROFILES end module korc_types diff --git a/src/korc_units.f90 b/src/korc_units.f90 index 896403b2..45c61792 100755 --- a/src/korc_units.f90 +++ b/src/korc_units.f90 @@ -210,16 +210,28 @@ subroutine normalize_variables(params,spp,F,P) if (ALLOCATED(P%ne_3D)) P%ne_3D = P%ne_3D/params%cpp%density if (ALLOCATED(P%Te_3D)) P%Te_3D = P%Te_3D/params%cpp%temperature - if (params%profile_model(10:10).eq.'H') then + if (params%profile_model(10:11).eq.'H0') then if (ALLOCATED(P%nRE_2D)) P%nRE_2D = P%nRE_2D/params%cpp%density if (ALLOCATED(P%nAr0_2D)) P%nAr0_2D = P%nAr0_2D/params%cpp%density if (ALLOCATED(P%nAr1_2D)) P%nAr1_2D = P%nAr1_2D/params%cpp%density if (ALLOCATED(P%nAr2_2D)) P%nAr2_2D = P%nAr2_2D/params%cpp%density if (ALLOCATED(P%nAr3_2D)) P%nAr3_2D = P%nAr3_2D/params%cpp%density + if (ALLOCATED(P%nAr3_2D)) P%nAr4_2D = P%nAr4_2D/params%cpp%density if (ALLOCATED(P%nD_2D)) P%nD_2D = P%nD_2D/params%cpp%density if (ALLOCATED(P%nD1_2D)) P%nD1_2D = P%nD1_2D/params%cpp%density endif + if (params%profile_model(10:12).eq.'Hdt') then + if (ALLOCATED(P%nRE_2D)) P%nRE_2D = P%nRE_2D/params%cpp%density + if (ALLOCATED(P%nAr0_3D)) P%nAr0_3D = P%nAr0_3D/params%cpp%density + if (ALLOCATED(P%nAr1_3D)) P%nAr1_3D = P%nAr1_3D/params%cpp%density + if (ALLOCATED(P%nAr2_3D)) P%nAr2_3D = P%nAr2_3D/params%cpp%density + if (ALLOCATED(P%nAr3_3D)) P%nAr3_3D = P%nAr3_3D/params%cpp%density + if (ALLOCATED(P%nAr3_3D)) P%nAr4_3D = P%nAr4_3D/params%cpp%density + if (ALLOCATED(P%nD_3D)) P%nD_3D = P%nD_3D/params%cpp%density + if (ALLOCATED(P%nD1_3D)) P%nD1_3D = P%nD1_3D/params%cpp%density + endif + end if if (params%field_model(1:10) .EQ. 'ANALYTICAL') then diff --git a/src/main.f90 b/src/main.f90 index f9a61542..017030fa 100755 --- a/src/main.f90 +++ b/src/main.f90 @@ -334,7 +334,12 @@ program main else if (params%orbit_model(1:2).eq.'GC') then if (.NOT.(params%restart.OR.params%proceed.or.params%reinit)) then - call GC_init(params,F,spp) + +#ifdef ACC + call GC_init_ACC(params,F,P,spp) +#else + call GC_init(params,randoms,F,P,spp) +#endif else call get_fields(params,spp(1)%vars,F) @@ -650,8 +655,11 @@ program main end if do it=params%ito,params%t_steps,params%t_skip +#ifdef ACC + call adv_GCinterp_psiwE_top_ACC(params,randoms,spp,P,F) +#else call adv_GCinterp_psiwE_top(params,randoms,spp,P,F) - +#endif if (.not.params%LargeCollisions) then params%time = params%init_time & +REAL(it-1_ip+params%t_skip,rp)*params%dt @@ -666,11 +674,13 @@ program main call save_simulation_outputs(params,spp,F) F%ind_2x1t=F%ind_2x1t+1_ip + P%ind_2x1t=F%ind_2x1t if (params%mpi_params%rank .EQ. 0) then write(output_unit_write,*) 'KORC time',params%time*params%cpp%time write(output_unit_write,*) '2x1t_ind time',F%X%PHI(F%ind_2x1t)*params%cpp%time end if call initialize_fields_interpolant(params,F) + call initialize_profiles_interpolant(params,P) if (params%LargeCollisions) then call initialize_collision_params(params,spp,P,F,.false.) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 873d5e40..9ffdaa8c 100755 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -23,6 +23,7 @@ set_property(TARGET xtest PROPERTY LINKER_LANGUAGE Fortran) configure_file(${CMAKE_SOURCE_DIR}/test/egyro/korc_egyro.sh.in ${CMAKE_BINARY_DIR}/egyro_test/korc_egyro.sh) configure_file(${CMAKE_SOURCE_DIR}/test/mars/korc_mars.sh.in ${CMAKE_BINARY_DIR}/mars_test/korc_mars.sh) configure_file(${CMAKE_SOURCE_DIR}/test/aorsa/korc_aorsa.sh.in ${CMAKE_BINARY_DIR}/aorsa_test/korc_aorsa.sh) +configure_file(${CMAKE_SOURCE_DIR}/test/hollmann/korc_hollmann.sh.in ${CMAKE_BINARY_DIR}/hollmann_test/korc_hollmann.sh) #if(USE_PSPLINE) # add_test (NAME mars_test_1 @@ -34,6 +35,7 @@ configure_file(${CMAKE_SOURCE_DIR}/test/aorsa/korc_aorsa.sh.in ${CMAKE_BINARY_DI file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/egyro_test) file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/mars_test) file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/aorsa_test) +file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/hollmann_test) foreach(RANK IN ITEMS 1 2 4 8) #16) if(${RANK} LESS_EQUAL ${MPIEXEC_MAX_NUMPROCS}) @@ -52,6 +54,11 @@ foreach(RANK IN ITEMS 1 2 4 8) #16) COMMAND ${CMAKE_BINARY_DIR}/aorsa_test/korc_aorsa.sh ${RANK} WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bin) set_tests_properties (aorsa_test_${RANK} PROPERTIES PROCESSORS ${RANK} ENVIRONMENT OMP_NUM_THREADS=1) + + add_test (NAME hollmann_test_${RANK} + COMMAND ${CMAKE_BINARY_DIR}/hollmann_test/korc_hollmann.sh ${RANK} + WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bin) + set_tests_properties (hollmann_test_${RANK} PROPERTIES PROCESSORS ${RANK} ENVIRONMENT OMP_NUM_THREADS=1) endif() add_test (NAME unit_testing_${RANK} diff --git a/test/hollmann/EFIT_JET_95128_Hratio.h5 b/test/hollmann/EFIT_JET_95128_Hratio.h5 new file mode 100644 index 00000000..c6f45c2f Binary files /dev/null and b/test/hollmann/EFIT_JET_95128_Hratio.h5 differ diff --git a/test/hollmann/Hollmann_PDF_JET_95128_3.h5 b/test/hollmann/Hollmann_PDF_JET_95128_3.h5 new file mode 100644 index 00000000..33cfff7a Binary files /dev/null and b/test/hollmann/Hollmann_PDF_JET_95128_3.h5 differ diff --git a/test/hollmann/input_file_JET_95128_TEST15b.korc b/test/hollmann/input_file_JET_95128_TEST15b.korc new file mode 100644 index 00000000..b9dacd0d --- /dev/null +++ b/test/hollmann/input_file_JET_95128_TEST15b.korc @@ -0,0 +1,170 @@ +&input_parameters + restart = F + ! Restart simulation that exited before simulation_time reached + proceed = F + ! Append simulation results after previous simulation_time reached + reinit = F + ! Begin a new simulation, reinitializing from restart file state + simulation_time = 1.E-8 + ! Total aimed simulation time in seconds + ! Run 10 mu s If transients exist put 5 extra mu s. + snapshot_frequency = 1.E-9 + ! Time between snapshots in seconds + restart_overwrite_frequency = 1.E-1 + ! Time between overwritting of restart file in seconds + dt = 1.58E0 + ! Time step as fraction of relativistic gyro-period + num_species = 1 + minimum_particle_energy = 1.0E2 + ! Minimum energy of simulated particles in eV + radiation = T + GC_rad_model='SDE' + collisions = T + LargeCollisions = F + collisions_model = 'SINGLE_SPECIES' + ! Options are: 'NONE','SINGLE_SPECIES' and 'MULTIPLE_SPECIES' + bound_electron_model = 'HESSLOW' + ! Options are: 'NO_BOUND', 'HESSLOW', and 'ROSENBLUTH' + field_model = 'EXTERNAL-PSI' + profile_model = 'EXTERNAL_Hdt' + ! The two options for this parameter are 'ANALYTICAL' or 'EXTERNAL'. + ! For 'ANALYTICAL', the magnetic field is calculated based on + ! the parameters given in the "analytic_mag_field_params" section. + ! For 'EXTERNAL', the magnetic field is loaded from the file + ! specified in "filename". + ! 'UNIFORM' A uniform magnetic field used to advance only electrons' + ! velocity. + magnetic_field_filename = 'EFIT_JET_95128_Hratio.h5' + outputs_list = '{X,Y,V,B,E,g,eta,flagCon,flagCol,PSIp}' + ! List of outputs + !'{X,Y,V,E,B,g,mu,eta,Prad,Pin,flagCon,flagCol,gradB, + ! curlb,ne,Te,Zeff,PSIp,nimp}' + HDF5_error_handling = T + orbit_model = 'GC' + ! 'FO' for full orbit, 'GCpre' for guiding center with pre-computed + ! auxiliary fields, 'GCgrad' for guiding center with auxiliary + ! fields computed with PSPLINE. + field_eval = 'interp' + ! Set for plasma_model='ANALYTICAL'. Can be 'interp' or 'eqn', + ! where 'eqn' evaluates particle fields at particle positions and + ! 'interp' interpolates precomputed fields. + SameRandSeed = T + !pchunk = 1 +/ + +&plasma_species + runaway = T + ppp = 100 + ! Number of particles per process (mpi) + q = -1.0 + ! Electric charge + m = 1.0 + ! In units of electron mass + spatial_distribution = 'HOLLMANN-1DTRANSPORT' + ! Options are: 'UNIFORM', 'DISK', 'TORUS', 'EXPONENTIAL-TORUS', + ! 'GAUSSIAN-TORUS', 'ELLIPTIC-TORUS', 'EXPONENTIAL-ELLIPTIC-TORUS', + ! 'GAUSSIAN-ELLIPTICAL-TORUS', '2D-GAUSSIAN-ELLIPTIC-TORUS-MH', + ! 'AVALANCHE-4D','TRACER','SPONG-3D','HOLLMANN-3D' + Ro = 2.6 + PHIo = 0.0 + ! In degrees + Zo = 0.0 + sigmaR = 1.e6 + sigmaZ = 0.2 + theta_gauss = 0.0 + psi_max=1. + ! goes as R^2 for HOLLMANN-3D, is psiN_max for HOLLMANN-3D-PSI + energy_distribution = 'HOLLMANN-1DTRANSPORT' + ! Options are: 'MONOENERGETIC', 'THERMAL', 'AVALANCHE', + ! 'EXPERIMENTAL', and 'UNIFORM' + pitch_distribution = 'HOLLMANN-1DTRANSPORT' + ! Options are: 'MONOPITCH', 'AVALANCHE', 'EXPERIMENTAL', and 'UNIFORM'. + Eno = 10.0E6 + ! Initial energy in eV + etao = 170.0 + ! Initial pitch angle + Eo_lims = 1.0E6,50.0E6 + ! Lower and upper limit of simulated energy range, in eV. + etao_lims = 0.0,20.0 + ! Lower and upper limit of simulated pitch-angle range, in degrees. + Xtrace = 2.6,0.0,0.0 + ! Initial position of tracer particle for debugging with + ! spatial_distribution='TRACER' + dth = 3. + ! Variance of sampling normal variate for pitch angle + dgam = 3. + ! Variance of sampling normal variate for pitch angle + dR = 0.1 + ! Variance of sampling normal variate for R location + dZ = 0.1 + ! Variance of sampling normal variate for Z location +/ + +&analytical_fields_params +/ + +&externalPlasmaModel + Bfield = F + B1field = F + E1field = F + dBfield = F + axisymmetric_fields = T + Bflux = F + Bflux3D = T + Efield = T + Dim2x1t =T + ReInterp_2x1t = T + PSIp_lim= 1.046608805604306e+00 + PSIp_0= 8.838935617719927e-01 + psip_conv=-1. + ! sign appended to \nabla\phi\times\nabla\psi_p for definition of + ! poloidal magnetic field. +1 for DIII-D, -1 for JET +/ + +&plasmaProfiles + axisymmetric = T + filename = 'EFIT_JET_95128_Hratio.h5' +/ + +&CollisionParamsSingleSpecies + Te_sing = 3.e0 + ! Background electron temperature in eV + Ti_sing = 3.e0 + ! Background ion temperature in eV + ne_sing = 9.63E19 + ! Background electron density in m^-3 + Zeff_sing = 3.84 + ! Effective atomic number + dTau_sing = 5.e-2 + ! Subcycling time step in collisional time units (Tau) + !p_therm = 1. + !min_secRE = 'CRIT' + !ConserveLA = T + !Clog_model = 'MCDEVITT' + !sample_test = F +/ + +&CollisionParamsMultipleSpecies + num_impurity_species = 1 + Te_mult = 3.e0 + ! Background electron temperature in eV + ne_mult = 9.36E19 + ! Background electron density in 1/m^3 +/ + +&HollmannPDF + filename_Hollmann = 'Hollmann_PDF_JET_95128_3.h5' + rho_ind = 7 + max_pitch_angle_Hollmann = 90. + ! In degrees + min_pitch_angle_Hollmann = 0. + ! In degrees + min_energy_Hollmann = 2.E3 + ! In eV + max_energy_Hollmann = 8.4E6 + ! For Hollmann_PDF_HR.h5, needs to be less than 80 MeV + ! In eV + current_direction_Hollmann = 'PARALLEL' + +/ + diff --git a/test/hollmann/korc_hollmann.sh.in b/test/hollmann/korc_hollmann.sh.in new file mode 100755 index 00000000..909adf7c --- /dev/null +++ b/test/hollmann/korc_hollmann.sh.in @@ -0,0 +1,26 @@ +#!/bin/bash + +set -ex + +#define input file +INPUT_FILE="${CMAKE_SOURCE_DIR}/test/hollmann/input_file_JET_95128_TEST15b.korc" +#define output directory +OUT_DIR="hollmann_test/rank_$1" + +#check that output directory doesn't exist so bash doesn't complain +if [ ! -d $OUT_DIR ]; then + mkdir -p $OUT_DIR +fi +if [ ! -f EFIT_JET_95128_Hratio.h5 ]; then + ln -s ${CMAKE_SOURCE_DIR}/test/hollmann/EFIT_JET_95128_Hratio.h5 EFIT_JET_95128_Hratio.h5 +fi +if [ ! -f Hollmann_PDF_JET_95128_3.h5 ]; then + ln -s ${CMAKE_SOURCE_DIR}/test/hollmann/Hollmann_PDF_JET_95128_3.h5 Hollmann_PDF_JET_95128_3.h5 +fi + +${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} $1 ./xkorc $INPUT_FILE $OUT_DIR/ + +for i in $(seq 0 $(($1-1))); +do + h5diff -r -p 0.000008 $OUT_DIR/file_$i.h5 ${CMAKE_SOURCE_DIR}/test/hollmann/rank_$1/file_$i.h5 +done \ No newline at end of file diff --git a/test/hollmann/rank_1/file_0.h5 b/test/hollmann/rank_1/file_0.h5 new file mode 100644 index 00000000..d264a668 Binary files /dev/null and b/test/hollmann/rank_1/file_0.h5 differ diff --git a/test/hollmann/rank_2/file_0.h5 b/test/hollmann/rank_2/file_0.h5 new file mode 100644 index 00000000..f1ebba91 Binary files /dev/null and b/test/hollmann/rank_2/file_0.h5 differ diff --git a/test/hollmann/rank_2/file_1.h5 b/test/hollmann/rank_2/file_1.h5 new file mode 100644 index 00000000..d6380b6e Binary files /dev/null and b/test/hollmann/rank_2/file_1.h5 differ diff --git a/test/hollmann/rank_4/file_0.h5 b/test/hollmann/rank_4/file_0.h5 new file mode 100644 index 00000000..312207c5 Binary files /dev/null and b/test/hollmann/rank_4/file_0.h5 differ diff --git a/test/hollmann/rank_4/file_1.h5 b/test/hollmann/rank_4/file_1.h5 new file mode 100644 index 00000000..da515aa3 Binary files /dev/null and b/test/hollmann/rank_4/file_1.h5 differ diff --git a/test/hollmann/rank_4/file_2.h5 b/test/hollmann/rank_4/file_2.h5 new file mode 100644 index 00000000..f1091eca Binary files /dev/null and b/test/hollmann/rank_4/file_2.h5 differ diff --git a/test/hollmann/rank_4/file_3.h5 b/test/hollmann/rank_4/file_3.h5 new file mode 100644 index 00000000..a51b1ec6 Binary files /dev/null and b/test/hollmann/rank_4/file_3.h5 differ diff --git a/test/hollmann/rank_8/file_0.h5 b/test/hollmann/rank_8/file_0.h5 new file mode 100644 index 00000000..9a1e0010 Binary files /dev/null and b/test/hollmann/rank_8/file_0.h5 differ diff --git a/test/hollmann/rank_8/file_1.h5 b/test/hollmann/rank_8/file_1.h5 new file mode 100644 index 00000000..b8e62007 Binary files /dev/null and b/test/hollmann/rank_8/file_1.h5 differ diff --git a/test/hollmann/rank_8/file_2.h5 b/test/hollmann/rank_8/file_2.h5 new file mode 100644 index 00000000..859f58bc Binary files /dev/null and b/test/hollmann/rank_8/file_2.h5 differ diff --git a/test/hollmann/rank_8/file_3.h5 b/test/hollmann/rank_8/file_3.h5 new file mode 100644 index 00000000..2f01dca8 Binary files /dev/null and b/test/hollmann/rank_8/file_3.h5 differ diff --git a/test/hollmann/rank_8/file_4.h5 b/test/hollmann/rank_8/file_4.h5 new file mode 100644 index 00000000..fa1b2266 Binary files /dev/null and b/test/hollmann/rank_8/file_4.h5 differ diff --git a/test/hollmann/rank_8/file_5.h5 b/test/hollmann/rank_8/file_5.h5 new file mode 100644 index 00000000..06b277eb Binary files /dev/null and b/test/hollmann/rank_8/file_5.h5 differ diff --git a/test/hollmann/rank_8/file_6.h5 b/test/hollmann/rank_8/file_6.h5 new file mode 100644 index 00000000..1a5ff19f Binary files /dev/null and b/test/hollmann/rank_8/file_6.h5 differ diff --git a/test/hollmann/rank_8/file_7.h5 b/test/hollmann/rank_8/file_7.h5 new file mode 100644 index 00000000..0dbce7a3 Binary files /dev/null and b/test/hollmann/rank_8/file_7.h5 differ