diff --git a/cgyro/src/Makefile b/cgyro/src/Makefile index ad5293929..df99f755c 100644 --- a/cgyro/src/Makefile +++ b/cgyro/src/Makefile @@ -22,7 +22,6 @@ OBJECTS = cgyro_globals.o \ cgyro_step.o \ cgyro_io.o \ cgyro_restart.o \ - cgyro_advect_wavenumber.o \ cgyro_init_rotation.o \ cgyro_equilibrium.o \ cgyro_check.o \ diff --git a/cgyro/src/cgyro_advect_wavenumber.F90 b/cgyro/src/cgyro_advect_wavenumber.F90 deleted file mode 100644 index a2329824a..000000000 --- a/cgyro/src/cgyro_advect_wavenumber.F90 +++ /dev/null @@ -1,105 +0,0 @@ -!--------------------------------------------------------- -! cgyro_advect_wavenumber.f90 -! -! PURPOSE: -! Manage shearing by wavenumber advection (legacy method) -!--------------------------------------------------------- - -subroutine cgyro_advect_wavenumber(ij) - - use cgyro_globals - use timer_lib - - implicit none - - integer, intent(in) :: ij - integer :: ir,l,ll,j,iccj,ivc,itor,llnt - complex :: rh,rl,he1,he2 - - if (nonlinear_flag == 0) return - - if (source_flag == 1) then - call timer_lib_in('shear') - -#if defined(OMPGPU) -!$omp target teams distribute parallel do simd collapse(4) & -!$omp& private(ivc,ir,l,iccj,j,ll,rl,rh,llnt,he1,he2) -#elif defined(_OPENACC) -!$acc parallel loop collapse(4) gang vector & -!$acc& private(ivc,ir,l,iccj,j,ll,rl,rh,llnt,he1,he2) & -!$acc& present(rhs(:,:,:,ij),omega_ss,field,h_x,c_wave) -#else -!$omp parallel do collapse(3) private(ivc,ir,l,iccj,j,ll,rl,rh,llnt,he1,he2) & -!$omp& firstprivate(shear_method,global_flag) -#endif - do itor=nt1,nt2 - do ivc=1,nv_loc - do ir=1,n_radial - do j=1,n_theta - iccj = (ir-1)*n_theta+j - rh = rhs(iccj,ivc,itor,ij) - - ! Wavenumber advection ExB shear - if (shear_method == 2) then - rl = 0.0 -#if (!defined(OMPGPU)) && defined(_OPENACC) -!$acc loop seq -#endif - do l=1,n_wave - ll = (2*l-1) - llnt = ll*n_theta - ! was he(j,ir+ll) - if ( (ir+ll) <= n_radial ) then - he1 = h_x(iccj+llnt,ivc,itor) - else - he1 = 0.0 - endif - ! was he(j,ir-ll) - if ( (ir-ll) >= 1 ) then - he2 = h_x(iccj-llnt,ivc,itor) - else - he2 = 0.0 - endif - ! Sign throughout paper is incorrect (or gamma -> - gamma) - ! Thus sign below has been checked and is correct - rl = rl+c_wave(l)*(he1-he2) - enddo - rh = rh + omega_eb_base*itor*rl - endif - - ! Wavenumber advection profile shear - if (global_flag == 1) then - iccj = (ir-1)*n_theta+j -#if (!defined(OMPGPU)) && defined(_OPENACC) -!$acc loop seq -#endif - do l=1,n_wave - ll = 2*l-1 - llnt = ll*n_theta - ! was he(j,ir+ll) - if ( (ir+ll) <= n_radial ) then - he1 = sum(omega_ss(:,iccj+llnt,ivc,itor)*field(:,iccj+llnt,itor)) - else - he1 = 0.0 - endif - ! was he(j,ir-ll) - if ( (ir-ll) >= 1 ) then - he2 = sum(omega_ss(:,iccj-llnt,ivc,itor)*field(:,iccj-llnt,itor)) - else - he2 = 0.0 - endif - ! Note opposite sign to ExB shear - rh = rh-c_wave(l)*(he1-he2) - enddo - endif - rhs(iccj,ivc,itor,ij) = rh - enddo - enddo - enddo - enddo - - call timer_lib_out('shear') - - endif - -end subroutine cgyro_advect_wavenumber diff --git a/cgyro/src/cgyro_globalshear.F90 b/cgyro/src/cgyro_globalshear.F90 index 4ace2a0a1..a5437b25d 100644 --- a/cgyro/src/cgyro_globalshear.F90 +++ b/cgyro/src/cgyro_globalshear.F90 @@ -19,14 +19,14 @@ subroutine cgyro_globalshear(ij) if (nonlinear_flag == 0 .or. source_flag == 0) return call timer_lib_in('shear') - + #if defined(OMPGPU) !$omp target teams distribute parallel do simd collapse(4) & !$omp& private(ivc,ir,l,iccj,j,ll,rl,llnt,h1,h2) #elif defined(_OPENACC) !$acc parallel loop collapse(4) gang vector & !$acc& private(ivc,ir,l,iccj,j,ll,rl,llnt,h1,h2) & -!$acc& present(rhs(:,:,:,ij),omega_ss,omega_sbeta,field,cap_h_c,c_wave) +!$acc& present(rhs(:,:,:,ij),omega_ss,omega_sbeta,field,cap_h_c,h_x,c_wave) #else !$omp parallel do collapse(4) private(ivc,ir,l,iccj,j,ll,rl,llnt,h1,h2) #endif @@ -48,7 +48,11 @@ subroutine cgyro_globalshear(ij) if ( (ir+ll) <= n_radial ) then ! ExB shear - h1 = omega_eb_base*itor*cap_h_c(iccj+llnt,ivc,itor) + if (shear_method == 3) then + h1 = omega_eb_base*itor*cap_h_c(iccj+llnt,ivc,itor) + else + h1 = omega_eb_base*itor*h_x(iccj+llnt,ivc,itor) + endif ! beta_star shear h1 = h1-omega_sbeta(iccj+llnt,ivc,itor)*cap_h_c(iccj+llnt,ivc,itor) ! omega_star shear @@ -59,7 +63,11 @@ subroutine cgyro_globalshear(ij) if ( (ir-ll) >= 1 ) then ! ExB shear - h2 = omega_eb_base*itor*cap_h_c(iccj-llnt,ivc,itor) + if (shear_method == 3) then + h2 = omega_eb_base*itor*cap_h_c(iccj-llnt,ivc,itor) + else + h2 = omega_eb_base*itor*h_x(iccj-llnt,ivc,itor) + endif ! beta_star shear h2 = h2-omega_sbeta(iccj-llnt,ivc,itor)*cap_h_c(iccj-llnt,ivc,itor) ! omega_star shear diff --git a/cgyro/src/cgyro_rhs.F90 b/cgyro/src/cgyro_rhs.F90 index 9eb18a87d..f42e6d12b 100644 --- a/cgyro/src/cgyro_rhs.F90 +++ b/cgyro/src/cgyro_rhs.F90 @@ -335,12 +335,8 @@ subroutine cgyro_rhs(ij,update_cap) call cgyro_rhs_comm_test - ! Wavenumber advection shear terms - if (shear_method == 2) then - ! s* Phi + gamma_E h - call cgyro_advect_wavenumber(ij) - else if (shear_method == 3) then - ! s* Phi + (gamma_E + sbeta) H + ! Wavenumber advection (ExB shear and/or global) terms + if (source_flag == 1) then call cgyro_globalshear(ij) endif