From 455cb16657039bd8a8672be8b405110b6de26e1a Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Tue, 18 Jun 2024 16:31:35 -0700 Subject: [PATCH 1/7] First round of global code --- cgyro/src/Makefile | 2 +- cgyro/src/cgyro_cleanup.F90 | 4 + cgyro/src/cgyro_global.F90 | 104 ------------------------ cgyro/src/cgyro_globals.F90 | 1 + cgyro/src/cgyro_globalshear.F90 | 78 ++++++++++++++++++ cgyro/src/cgyro_init_arrays.F90 | 137 ++++++++++++++++++-------------- 6 files changed, 162 insertions(+), 164 deletions(-) delete mode 100644 cgyro/src/cgyro_global.F90 create mode 100644 cgyro/src/cgyro_globalshear.F90 diff --git a/cgyro/src/Makefile b/cgyro/src/Makefile index f446a6025..1c627eafb 100644 --- a/cgyro/src/Makefile +++ b/cgyro/src/Makefile @@ -32,7 +32,7 @@ OBJECTS = cgyro_globals.o \ cgyro_field_coefficients.o \ cgyro_flux.o \ cgyro_freq.o \ - cgyro_global.o \ + cgyro_globalshear.o \ cgyro_upwind.o \ cgyro_init.o \ cgyro_init_h.o \ diff --git a/cgyro/src/cgyro_cleanup.F90 b/cgyro/src/cgyro_cleanup.F90 index fa65a764b..0b4d7cbb4 100644 --- a/cgyro/src/cgyro_cleanup.F90 +++ b/cgyro/src/cgyro_cleanup.F90 @@ -197,6 +197,10 @@ subroutine cgyro_cleanup ccl_del_device(omega_ss) deallocate(omega_ss) endif + if(allocated(omega_beta)) then + ccl_del_device(omega_beta) + deallocate(omega_beta) + endif if(allocated(jvec_c)) then ccl_del_device(jvec_c) deallocate(jvec_c) diff --git a/cgyro/src/cgyro_global.F90 b/cgyro/src/cgyro_global.F90 deleted file mode 100644 index 80797d839..000000000 --- a/cgyro/src/cgyro_global.F90 +++ /dev/null @@ -1,104 +0,0 @@ -!--------------------------------------------------------- -! cgyro_advect_wavenumber.f90 -! -! PURPOSE: -! Manage shearing by wavenumber advection. -!--------------------------------------------------------- - -subroutine cgyro_global(ij) - - use cgyro_globals - use timer_lib - - implicit none - - integer, intent(in) :: ij - integer :: ir,l,ll,j,iccj,ivc,itor,llnt - complex :: 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,llnt,he1,he2) -#elif defined(_OPENACC) -!$acc parallel loop collapse(4) gang vector & -!$acc& private(ivc,ir,l,iccj,j,ll,rl,llnt,he1,he2) & -!$acc& present(rhs(:,:,:,ij),omega_ss,field,h_x,c_wave) -#else -!$omp parallel do collapse(4) private(ivc,ir,l,iccj,j,ll,rl,llnt,he1,he2) -#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 - - ! 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 - rhs(iccj,ivc,itor,ij) = rhs(iccj,ivc,itor,ij) + omega_eb_base*itor*rl - endif - - ! Wavenumber advection profile shear - if (profile_shear_flag == 1) then - iccj = (ir-1)*n_theta+j - rl = rhs(iccj,ivc,itor,ij) -#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 - rl = rl-c_wave(l)*(he1-he2) - enddo - rhs(iccj,ivc,itor,ij) = rl - endif - enddo - enddo - enddo - enddo - - call timer_lib_out('shear') - - endif - -end subroutine cgyro_global diff --git a/cgyro/src/cgyro_globals.F90 b/cgyro/src/cgyro_globals.F90 index 1c6efec0c..542475d77 100644 --- a/cgyro/src/cgyro_globals.F90 +++ b/cgyro/src/cgyro_globals.F90 @@ -348,6 +348,7 @@ module cgyro_globals complex, dimension(:,:,:), allocatable :: omega_cap_h complex, dimension(:,:,:), allocatable :: omega_h complex, dimension(:,:,:,:), allocatable :: omega_s,omega_ss + complex, dimension(:,:), allocatable :: omega_beta complex, dimension(:,:,:), allocatable :: cap_h_c complex, dimension(:,:,:), allocatable :: cap_h_c_dot complex, dimension(:,:,:), allocatable :: cap_h_c_old diff --git a/cgyro/src/cgyro_globalshear.F90 b/cgyro/src/cgyro_globalshear.F90 new file mode 100644 index 000000000..f6a10e8e4 --- /dev/null +++ b/cgyro/src/cgyro_globalshear.F90 @@ -0,0 +1,78 @@ +!--------------------------------------------------------- +! cgyro_advect_wavenumber.f90 +! +! PURPOSE: +! Manage shearing by wavenumber advection. +!--------------------------------------------------------- + +subroutine cgyro_global(ij) + + use cgyro_globals + use timer_lib + + implicit none + + integer, intent(in) :: ij + integer :: ir,l,ll,j,iccj,ivc,itor,llnt + complex :: rl,h1,h2 + + if (nonlinear_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,field,h_x,cap_h_c,c_wave) +#else +!$omp parallel do collapse(4) private(ivc,ir,l,iccj,j,ll,rl,llnt,h1,h2) +#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 + 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 + + if ( (ir+ll) <= n_radial ) then + h1 = omega_eb_base*itor*h_x(iccj+llnt,ivc,itor) + h1 = h1-sum(omega_ss(:,iccj+llnt,ivc,itor)*field(:,iccj+llnt,itor)) + h1 = h1-0*cap_h_c(iccj+llnt,ivc,itor) + else + h1 = 0.0 + endif + + if ( (ir-ll) >= 1 ) then + h2 = omega_eb_base*itor*h_x(iccj-llnt,ivc,itor) + h2 = h2-sum(omega_ss(:,iccj-llnt,ivc,itor)*field(:,iccj-llnt,itor)) + h2 = h2-0*cap_h_c(iccj-llnt,ivc,itor) + else + h2 = 0.0 + endif + + rl = rl+c_wave(l)*(h1-h2) + + enddo + + rhs(iccj,ivc,itor,ij) = rl + + enddo + enddo + enddo + enddo + + call timer_lib_out('shear') + +end subroutine cgyro_global diff --git a/cgyro/src/cgyro_init_arrays.F90 b/cgyro/src/cgyro_init_arrays.F90 index a6884486c..1f01d77af 100644 --- a/cgyro/src/cgyro_init_arrays.F90 +++ b/cgyro/src/cgyro_init_arrays.F90 @@ -8,7 +8,7 @@ subroutine cgyro_init_arrays real :: arg real :: efac - real :: u + real :: u,sm,sb real :: fac integer :: ir,it,is,ie,ix integer :: itm,itl,itor,mytor,itf @@ -25,6 +25,14 @@ subroutine cgyro_init_arrays real, dimension(:,:), allocatable :: res_weight real, external :: spectraldiss + real, dimension(:), allocatable :: gdlnndr,gdlntdr + + allocate(gdlnndr(n_species)) + allocate(gdlntdr(n_species)) + + gdlntdr = 0 + gdlnndr = 0 + !------------------------------------------------------------------------- ! Distributed Bessel-function Gyroaverages @@ -182,7 +190,7 @@ subroutine cgyro_init_arrays enddo enddo endif - + res_loc(:,:,:) = 0.0 !$omp parallel private(ic,iv_loc,is,ix,ie) @@ -390,80 +398,91 @@ subroutine cgyro_init_arrays ! Streaming coefficients (for speed optimization) !$omp parallel do collapse(3) & -!$omp& private(iv,ic,iv_loc,is,ix,ie,ir,it,carg,u) +!$omp& private(iv,ic,iv_loc,is,ix,ie,ir,it,arg,carg,u,sm,sb) do itor=nt1,nt2 - do iv=nv1,nv2 - do ic=1,nc - - iv_loc = iv-nv1+1 - is = is_v(iv) - ix = ix_v(iv) - ie = ie_v(iv) - ir = ir_c(ic) - it = it_c(ic) - - u = (pi/n_toroidal)*itor - - ! omega_dalpha - omega_cap_h(ic,iv_loc,itor) = & - -omega_adrift(it,is)*energy(ie)*(1.0+xi(ix)**2)*& - (n_toroidal*q/pi/rmin)*(i_c*u) - - ! omega_dalpha [UPWIND: iu -> spectraldiss] - omega_h(ic,iv_loc,itor) = & - -abs(omega_adrift(it,is))*energy(ie)*(1.0+xi(ix)**2)*& - (n_toroidal*q/pi/rmin)*spectraldiss(u,nup_alpha)*up_alpha - - ! (i ktheta) components from drifts - omega_cap_h(ic,iv_loc,itor) = omega_cap_h(ic,iv_loc,itor) & - - i_c*k_theta_base*itor*(omega_aprdrift(it,is)*energy(ie)*xi(ix)**2 & - + omega_cdrift(it,is)*vel(ie)*xi(ix) + omega_rot_drift(it,is) & - + omega_rot_edrift(it)) + do iv=nv1,nv2 + do ic=1,nc + + iv_loc = iv-nv1+1 + is = is_v(iv) + ix = ix_v(iv) + ie = ie_v(iv) + ir = ir_c(ic) + it = it_c(ic) + + u = (pi/n_toroidal)*itor + + ! omega_dalpha + omega_cap_h(ic,iv_loc,itor) = & + -omega_adrift(it,is)*energy(ie)*(1.0+xi(ix)**2)*& + (n_toroidal*q/pi/rmin)*(i_c*u) + + ! omega_dalpha [UPWIND: iu -> spectraldiss] + omega_h(ic,iv_loc,itor) = & + -abs(omega_adrift(it,is))*energy(ie)*(1.0+xi(ix)**2)*& + (n_toroidal*q/pi/rmin)*spectraldiss(u,nup_alpha)*up_alpha + + ! (i ktheta) components from drifts + omega_cap_h(ic,iv_loc,itor) = omega_cap_h(ic,iv_loc,itor) & + - i_c*k_theta_base*itor*(omega_aprdrift(it,is)*energy(ie)*xi(ix)**2 & + + omega_cdrift(it,is)*vel(ie)*xi(ix) + omega_rot_drift(it,is) & + + omega_rot_edrift(it)) - ! Note that we shift the dissipation with px0 (ballooning angle linear mode) - u = (2.0*pi/n_radial)*(px(ir)+px0) + ! Note that we shift the dissipation with px0 (ballooning angle linear mode) + u = (2.0*pi/n_radial)*(px(ir)+px0) - ! (d/dr) components from drifts + ! (d/dr) components from drifts - omega_cap_h(ic,iv_loc,itor) = omega_cap_h(ic,iv_loc,itor) & - - (n_radial/length)*i_c*u & - * (omega_rdrift(it,is)*energy(ie)*(1.0+xi(ix)**2) & - + omega_cdrift_r(it,is)*vel(ie)*xi(ix) & - + omega_rot_drift_r(it,is) & - + omega_rot_edrift_r(it)) + omega_cap_h(ic,iv_loc,itor) = omega_cap_h(ic,iv_loc,itor) & + - (n_radial/length)*i_c*u & + * (omega_rdrift(it,is)*energy(ie)*(1.0+xi(ix)**2) & + + omega_cdrift_r(it,is)*vel(ie)*xi(ix) & + + omega_rot_drift_r(it,is) & + + omega_rot_edrift_r(it)) - ! (d/dr) upwind components from drifts [UPWIND: iu -> spectraldiss] - omega_h(ic,iv_loc,itor) = omega_h(ic,iv_loc,itor) & - - (n_radial/length)*spectraldiss(u,nup_radial)*up_radial & - * (abs(omega_rdrift(it,is))*energy(ie)*(1.0+xi(ix)**2) & - + abs(omega_cdrift_r(it,is)*xi(ix))*vel(ie) & - + abs(omega_rot_drift_r(it,is)) & - + abs(omega_rot_edrift_r(it))) + ! (d/dr) upwind components from drifts [UPWIND: iu -> spectraldiss] + omega_h(ic,iv_loc,itor) = omega_h(ic,iv_loc,itor) & + - (n_radial/length)*spectraldiss(u,nup_radial)*up_radial & + * (abs(omega_rdrift(it,is))*energy(ie)*(1.0+xi(ix)**2) & + + abs(omega_cdrift_r(it,is)*xi(ix))*vel(ie) & + + abs(omega_rot_drift_r(it,is)) & + + abs(omega_rot_edrift_r(it))) - ! omega_star - carg = -i_c*k_theta_base*itor*rho*(dlnndr(is)+dlntdr(is)*(energy(ie)-1.5)) & - -i_c*k_theta_base*itor*rho*(vel2(ie)*xi(ix)/vth(is) & - *omega_gammap(it)) -i_c*k_theta_base*itor*rho*omega_rot_star(it,is) + ! omega_star + carg = -i_c*k_theta_base*itor*rho*(dlnndr(is)+dlntdr(is)*(energy(ie)-1.5)) & + -i_c*k_theta_base*itor*rho*(vel2(ie)*xi(ix)/vth(is) & + *omega_gammap(it)) -i_c*k_theta_base*itor*rho*omega_rot_star(it,is) + + omega_s(:,ic,iv_loc,itor) = carg*jvec_c(:,ic,iv_loc,itor) + + ! global-Taylor corrections via wavenumber advection (ix -> d/dp) + ! JC: Re-checked sign and normalization (Oct 2019) - omega_s(:,ic,iv_loc,itor) = carg*jvec_c(:,ic,iv_loc,itor) + ! generalized profile curvature (phi shearing) + sm = sdlnndr(is)+sdlntdr(is)*(energy(ie)-1.5) & + +energy(ie)*gdlntdr(is)**2+gdlnndr(is)**2-gdlntdr(is)*gdlnndr(is) - ! Profile curvature via wavenumber advection (ix -> d/dp) - ! See whiteboard notes. - ! JC: Re-checked sign and normalization (Oct 2019) - carg = -k_theta_base*itor*length*(sdlnndr(is)+sdlntdr(is)*(energy(ie)-1.5))/(2*pi) + ! beta_star (H shearing) + sb = beta_star(0)*energy(ie)*(sdlntdr(is)+sdlnndr(is)-2*gdlntdr(is)*gdlnndr(is)) & + /(dlnndr(is)+dlntdr(is)) - omega_ss(:,ic,iv_loc,itor) = carg*jvec_c(:,ic,iv_loc,itor) + arg = -k_theta_base*itor*length/(2*pi) + omega_ss(:,ic,iv_loc,itor) = arg*sm*jvec_c(:,ic,iv_loc,itor) + omega_beta(iv_loc,itor) = arg*sb + + enddo enddo - enddo enddo #if defined(OMPGPU) -!$omp target enter data map(to:omega_cap_h,omega_h,omega_s,omega_ss) +!$omp target enter data map(to:omega_cap_h,omega_h,omega_s,omega_ss,omega_beta) #elif defined(_OPENACC) -!$acc enter data copyin(omega_cap_h,omega_h,omega_s,omega_ss) +!$acc enter data copyin(omega_cap_h,omega_h,omega_s,omega_ss,omega_beta) #endif !------------------------------------------------------------------------- + deallocate(gdlnndr,gdlntdr) + end subroutine cgyro_init_arrays ! Spectral dissipation function From cf3814f23bf66120b6d2635109b408a0765885f9 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Tue, 18 Jun 2024 16:36:15 -0700 Subject: [PATCH 2/7] Delete .gitmodules --- .gitmodules | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 .gitmodules diff --git a/.gitmodules b/.gitmodules deleted file mode 100644 index e69de29bb..000000000 From 07d0a0dddc2693e74a00c63d2724ba90e4915eb3 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Tue, 18 Jun 2024 17:14:07 -0700 Subject: [PATCH 3/7] Added data template --- f2py/pygacode/cgyro/data_template.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) create mode 100644 f2py/pygacode/cgyro/data_template.py diff --git a/f2py/pygacode/cgyro/data_template.py b/f2py/pygacode/cgyro/data_template.py new file mode 100644 index 000000000..cf0e61d6a --- /dev/null +++ b/f2py/pygacode/cgyro/data_template.py @@ -0,0 +1,14 @@ +import os +import numpy as np +from gacodefuncs import * +from cgyro.data import cgyrodata + +sim = cgyrodata('reg01/') + +# all relevant quantities available in gacode/f2py/pygacode/cgyro/data.py + +print('kappa',sim.kappa) +print('ky*rho',sim.ky0) + +print('omega',sim.freq[0,0,-1]) +print('gamma',sim.freq[1,0,-1]) From be1fe1a16ba7ba5385de46b08cca188b96eb2910 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Tue, 30 Jul 2024 18:31:20 -0700 Subject: [PATCH 4/7] Makefile fixes --- cgyro/src/Makefile | 2 +- cgyro/src/cgyro_globalshear.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cgyro/src/Makefile b/cgyro/src/Makefile index 1c627eafb..cd741f853 100644 --- a/cgyro/src/Makefile +++ b/cgyro/src/Makefile @@ -15,9 +15,9 @@ FC += ${FOMP} ${FACC} EXEC = cgyro OBJECTS = cgyro_globals.o \ - cgyro_timer_lib.o \ cgyro_math.o \ cgyro_globals_math.o \ + cgyro_timer_lib.o \ cgyro_parallel_lib.o \ cgyro_step.o \ cgyro_io.o \ diff --git a/cgyro/src/cgyro_globalshear.F90 b/cgyro/src/cgyro_globalshear.F90 index f6a10e8e4..450fbb8b0 100644 --- a/cgyro/src/cgyro_globalshear.F90 +++ b/cgyro/src/cgyro_globalshear.F90 @@ -5,7 +5,7 @@ ! Manage shearing by wavenumber advection. !--------------------------------------------------------- -subroutine cgyro_global(ij) +subroutine cgyro_globalshear(ij) use cgyro_globals use timer_lib @@ -75,4 +75,4 @@ subroutine cgyro_global(ij) call timer_lib_out('shear') -end subroutine cgyro_global +end subroutine cgyro_globalshear From 6d33e78efe3d79c61c915e9a8a2e41af0a69a53e Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Wed, 31 Jul 2024 18:28:51 -0700 Subject: [PATCH 5/7] Changes to little-h equation done in terms of sdln and sbeta_star parameters --- cgyro/bin/cgyro_parse.py | 1 + cgyro/src/cgyro_cleanup.F90 | 6 +++--- cgyro/src/cgyro_globals.F90 | 3 ++- cgyro/src/cgyro_globalshear.F90 | 8 ++++---- cgyro/src/cgyro_init_arrays.F90 | 8 +++----- cgyro/src/cgyro_init_manager.F90 | 3 ++- cgyro/src/cgyro_read_input.f90 | 1 + cgyro/src/cgyro_rhs.F90 | 3 ++- 8 files changed, 18 insertions(+), 15 deletions(-) diff --git a/cgyro/bin/cgyro_parse.py b/cgyro/bin/cgyro_parse.py index dcd597138..de430e15b 100644 --- a/cgyro/bin/cgyro_parse.py +++ b/cgyro/bin/cgyro_parse.py @@ -136,6 +136,7 @@ x.add('DLNTDR','1.0',n=n) x.add('SDLNNDR','0.0',n=n) x.add('SDLNTDR','0.0',n=n) +x.add('SBETA_STAR','0.0',n=n) x.add('DLNNDR_SCALE','1.0',n=n) x.add('DLNTDR_SCALE','1.0',n=n) diff --git a/cgyro/src/cgyro_cleanup.F90 b/cgyro/src/cgyro_cleanup.F90 index 0b4d7cbb4..d868bd20a 100644 --- a/cgyro/src/cgyro_cleanup.F90 +++ b/cgyro/src/cgyro_cleanup.F90 @@ -197,9 +197,9 @@ subroutine cgyro_cleanup ccl_del_device(omega_ss) deallocate(omega_ss) endif - if(allocated(omega_beta)) then - ccl_del_device(omega_beta) - deallocate(omega_beta) + if(allocated(omega_sbeta)) then + ccl_del_device(omega_sbeta) + deallocate(omega_sbeta) endif if(allocated(jvec_c)) then ccl_del_device(jvec_c) diff --git a/cgyro/src/cgyro_globals.F90 b/cgyro/src/cgyro_globals.F90 index 542475d77..bbf87930b 100644 --- a/cgyro/src/cgyro_globals.F90 +++ b/cgyro/src/cgyro_globals.F90 @@ -135,6 +135,7 @@ module cgyro_globals real, dimension(11) :: dlntdr real, dimension(11) :: sdlnndr real, dimension(11) :: sdlntdr + real, dimension(11) :: sbeta_star integer :: subroutine_flag ! only used for cgyro_read_input @@ -348,7 +349,7 @@ module cgyro_globals complex, dimension(:,:,:), allocatable :: omega_cap_h complex, dimension(:,:,:), allocatable :: omega_h complex, dimension(:,:,:,:), allocatable :: omega_s,omega_ss - complex, dimension(:,:), allocatable :: omega_beta + complex, dimension(:,:,:), allocatable :: omega_sbeta complex, dimension(:,:,:), allocatable :: cap_h_c complex, dimension(:,:,:), allocatable :: cap_h_c_dot complex, dimension(:,:,:), allocatable :: cap_h_c_old diff --git a/cgyro/src/cgyro_globalshear.F90 b/cgyro/src/cgyro_globalshear.F90 index 450fbb8b0..b8698bcad 100644 --- a/cgyro/src/cgyro_globalshear.F90 +++ b/cgyro/src/cgyro_globalshear.F90 @@ -26,7 +26,7 @@ subroutine cgyro_globalshear(ij) #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,field,h_x,cap_h_c,c_wave) +!$acc& present(rhs(:,:,:,ij),omega_ss,omega_sbeta,field,h_x,cap_h_c,c_wave) #else !$omp parallel do collapse(4) private(ivc,ir,l,iccj,j,ll,rl,llnt,h1,h2) #endif @@ -36,7 +36,7 @@ subroutine cgyro_globalshear(ij) do j=1,n_theta iccj = (ir-1)*n_theta+j - rl = 0.0 + rl = rhs(iccj,ivc,itor,ij) #if (!defined(OMPGPU)) && defined(_OPENACC) !$acc loop seq @@ -49,7 +49,7 @@ subroutine cgyro_globalshear(ij) if ( (ir+ll) <= n_radial ) then h1 = omega_eb_base*itor*h_x(iccj+llnt,ivc,itor) h1 = h1-sum(omega_ss(:,iccj+llnt,ivc,itor)*field(:,iccj+llnt,itor)) - h1 = h1-0*cap_h_c(iccj+llnt,ivc,itor) + h1 = h1-omega_sbeta(iccj+llnt,ivc,itor)*cap_h_c(iccj+llnt,ivc,itor) else h1 = 0.0 endif @@ -57,7 +57,7 @@ subroutine cgyro_globalshear(ij) if ( (ir-ll) >= 1 ) then h2 = omega_eb_base*itor*h_x(iccj-llnt,ivc,itor) h2 = h2-sum(omega_ss(:,iccj-llnt,ivc,itor)*field(:,iccj-llnt,itor)) - h2 = h2-0*cap_h_c(iccj-llnt,ivc,itor) + h2 = h2-omega_sbeta(iccj-llnt,ivc,itor)*cap_h_c(iccj-llnt,ivc,itor) else h2 = 0.0 endif diff --git a/cgyro/src/cgyro_init_arrays.F90 b/cgyro/src/cgyro_init_arrays.F90 index 1f01d77af..124026d47 100644 --- a/cgyro/src/cgyro_init_arrays.F90 +++ b/cgyro/src/cgyro_init_arrays.F90 @@ -459,17 +459,15 @@ subroutine cgyro_init_arrays ! JC: Re-checked sign and normalization (Oct 2019) ! generalized profile curvature (phi shearing) - sm = sdlnndr(is)+sdlntdr(is)*(energy(ie)-1.5) & - +energy(ie)*gdlntdr(is)**2+gdlnndr(is)**2-gdlntdr(is)*gdlnndr(is) + sm = sdlnndr(is)+sdlntdr(is)*(energy(ie)-1.5) ! beta_star (H shearing) - sb = beta_star(0)*energy(ie)*(sdlntdr(is)+sdlnndr(is)-2*gdlntdr(is)*gdlnndr(is)) & - /(dlnndr(is)+dlntdr(is)) + sb = sbeta_star(is)*energy(ie)/bmag(it)**2 arg = -k_theta_base*itor*length/(2*pi) omega_ss(:,ic,iv_loc,itor) = arg*sm*jvec_c(:,ic,iv_loc,itor) - omega_beta(iv_loc,itor) = arg*sb + omega_sbeta(ic,iv_loc,itor) = arg*sb enddo enddo diff --git a/cgyro/src/cgyro_init_manager.F90 b/cgyro/src/cgyro_init_manager.F90 index 1876eacdd..aec0a8149 100644 --- a/cgyro/src/cgyro_init_manager.F90 +++ b/cgyro/src/cgyro_init_manager.F90 @@ -250,6 +250,7 @@ subroutine cgyro_init_manager allocate(omega_h(nc,nv_loc,nt1:nt2)) allocate(omega_s(n_field,nc,nv_loc,nt1:nt2)) allocate(omega_ss(n_field,nc,nv_loc,nt1:nt2)) + allocate(omega_sbeta(nc,nv_loc,nt1:nt2)) allocate(jvec_c(n_field,nc,nv_loc,nt1:nt2)) allocate(jvec_v(n_field,nc_loc,nt1:nt2,nv)) allocate(dvjvec_c(n_field,nc,nv_loc,nt1:nt2)) @@ -257,7 +258,7 @@ subroutine cgyro_init_manager allocate(jxvec_c(n_field,nc,nv_loc,nt1:nt2)) allocate(upfac1(nc,nv_loc,nt1:nt2)) allocate(upfac2(nc,nv_loc,nt1:nt2)) - + #if defined(OMPGPU) !$omp target enter data map(alloc:cap_h_c,cap_h_ct,cap_h_c_dot,cap_h_c_old,cap_h_c_old2) !$omp target enter data map(alloc:cap_h_v,dvjvec_c,dvjvec_v) diff --git a/cgyro/src/cgyro_read_input.f90 b/cgyro/src/cgyro_read_input.f90 index faab599b5..480c59b28 100644 --- a/cgyro/src/cgyro_read_input.f90 +++ b/cgyro/src/cgyro_read_input.f90 @@ -126,6 +126,7 @@ subroutine cgyro_read_input call cgyro_readbc_realv(dlntdr,is) call cgyro_readbc_realv(sdlnndr,is) call cgyro_readbc_realv(sdlntdr,is) + call cgyro_readbc_realv(sbeta_star,is) call cgyro_readbc_realv(dlnndr_scale,is) call cgyro_readbc_realv(dlntdr_scale,is) diff --git a/cgyro/src/cgyro_rhs.F90 b/cgyro/src/cgyro_rhs.F90 index 8955746f9..1a9a15028 100644 --- a/cgyro/src/cgyro_rhs.F90 +++ b/cgyro/src/cgyro_rhs.F90 @@ -337,7 +337,8 @@ subroutine cgyro_rhs(ij,update_cap) ! Wavenumber advection shear terms ! depends on h_x and field only, updates rhs - call cgyro_advect_wavenumber(ij) + !call cgyro_advect_wavenumber(ij) + call cgyro_globalshear(ij) call cgyro_rhs_comm_test From de5c320069e2e7228fda90beda4536250df0001d Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Tue, 3 Sep 2024 13:54:48 -0700 Subject: [PATCH 6/7] commented out global mode for merge back into master --- cgyro/src/cgyro_rhs.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cgyro/src/cgyro_rhs.F90 b/cgyro/src/cgyro_rhs.F90 index 1a9a15028..c47a6038f 100644 --- a/cgyro/src/cgyro_rhs.F90 +++ b/cgyro/src/cgyro_rhs.F90 @@ -337,8 +337,8 @@ subroutine cgyro_rhs(ij,update_cap) ! Wavenumber advection shear terms ! depends on h_x and field only, updates rhs - !call cgyro_advect_wavenumber(ij) - call cgyro_globalshear(ij) + call cgyro_advect_wavenumber(ij) + !call cgyro_globalshear(ij) call cgyro_rhs_comm_test From e274a743814f40177726e3a5e33ade86bef1ec28 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Tue, 3 Sep 2024 18:28:07 -0700 Subject: [PATCH 7/7] Cleaned out some old files. Added new JSON input converter. --- cgyro/bin/README_omp.txt | 22 +++---- cgyro/bin/cgyro_b2b.py | 26 -------- cgyro/bin/cgyro_compress | 43 ------------ cgyro/bin/cgyro_json | 42 ++++++++++++ cgyro/tools/zf/getdata.py | 38 ----------- cgyro/tools/zf/timeaverage_fields.py | 98 ---------------------------- 6 files changed, 51 insertions(+), 218 deletions(-) delete mode 100644 cgyro/bin/cgyro_b2b.py delete mode 100755 cgyro/bin/cgyro_compress create mode 100755 cgyro/bin/cgyro_json delete mode 100644 cgyro/tools/zf/getdata.py delete mode 100644 cgyro/tools/zf/timeaverage_fields.py diff --git a/cgyro/bin/README_omp.txt b/cgyro/bin/README_omp.txt index c115d4eb9..326eace84 100644 --- a/cgyro/bin/README_omp.txt +++ b/cgyro/bin/README_omp.txt @@ -1,22 +1,18 @@ Running CGYRO with high nomp ============================ -CGYRO allocates some of its data on the stack, -and the amount is proportional with the -number of omp threads you use. -(omp==openMP) +CGYRO allocates some of its data on the stack, and the amount is +proportional to the number of omp (OpenMP) threads you use. The +default stack limit of 8M is only sufficient for up to about 4 omp +threads. Thus, when using CGYRO with openMP, it is recommended +that you put the following line into your .bashrc -The default stack limit of 8M is only sufficient -for up to about 4 omp threads. +ulimit -s 32768 -Thus, when using CGYRO with openMP, -it is recommended that you put the following line -into your .bashrc -ulimit -s 32768 -(must be set on all the nodes running CGYRO) +(This must be set on all the nodes running CGYRO). Moreover, for +platforms that do not use the system stack limit for OpenMP +threads, also set -Moreover, for platform that do not use the -system stack limit for OpenMP threads, also set export OMP_STACKSIZE=32M This should be sufficient for at least 128 omp threads. diff --git a/cgyro/bin/cgyro_b2b.py b/cgyro/bin/cgyro_b2b.py deleted file mode 100644 index c16ab66b8..000000000 --- a/cgyro/bin/cgyro_b2b.py +++ /dev/null @@ -1,26 +0,0 @@ -# file processed by 2to3 -from __future__ import print_function, absolute_import -from builtins import map, filter, range -import os -import sys -import numpy as np - -if len(sys.argv) == 1: - print('python cgyro_b2b.py ') - sys.exit() - -afile = sys.argv[1] -bfile = sys.argv[2] - -if not os.path.isfile(afile): - print("INFO: (cgyro_b2b) "+afile+' not found') - sys.exit() - -arr = np.fromfile(afile,dtype='float64') -brr = arr.astype('float32') - -print("INFO: (cgyro_b2b) Read float64 data in "+afile) - -brr.tofile(bfile) - -print("INFO: (cgyro_t2b) Wrote float32 data to "+bfile) diff --git a/cgyro/bin/cgyro_compress b/cgyro/bin/cgyro_compress deleted file mode 100755 index 16ccb9901..000000000 --- a/cgyro/bin/cgyro_compress +++ /dev/null @@ -1,43 +0,0 @@ -#!/bin/bash -#============================================================= -# cgyro_compress -# -# PURPOSE: -# Create "small" version of CGYRO directory. -#============================================================= - -if [ -d $1-s ] -then - echo "WARNING: (cgyro_compress) Removed old directory $1-s." - rm -r $1-s -fi - -SDIR=$1-s -echo "INFO: (cgyro) Compressing $2 -> $2-s" - -# Basic content -mkdir $SDIR -cp -p $1/input.cgyro $SDIR -cp -p $1/input.cgyro.gen $SDIR/input.cgyro.gen.original -cp -p $1/out.cgyro.grids $SDIR -cp -p $1/out.cgyro.equilibrium $SDIR -cp -p $1/out.cgyro.time $SDIR -cp -p $1/out.cgyro.timing $SDIR -cp -p $1/out.cgyro.info $SDIR -cp -p $1/*.cgyro.geo $SDIR - -# Linear data -if ! grep -q NONLINEAR_FLAG=1 $SDIR/input.cgyro -then - cp -p $1/*.cgyro.freq $SDIR - cp -p $1/*.cgyro.phib $SDIR -fi - -# Large files (but no restart) -cp -p $1/bin.cgyro.*flux* $SDIR - -# We will have a version if NOT run through TGYRO -if [ -f $1/out.cgyro.version ] -then - cp -p $1/out.cgyro.version $SDIR -fi diff --git a/cgyro/bin/cgyro_json b/cgyro/bin/cgyro_json new file mode 100755 index 000000000..27c8ff893 --- /dev/null +++ b/cgyro/bin/cgyro_json @@ -0,0 +1,42 @@ +#!/usr/bin/env python + +import json +from pygacode.cgyro import data + +sim = data.cgyrodata('./') + +d = {} + +d['RMIN'] = sim.rmin +d['RMAJ'] = sim.rmaj +d['SHEAR'] = sim.shear +d['Q'] = sim.q +d['Z'] = sim.z[:].tolist() +d['MASS'] = sim.mass[:].tolist() +d['DENS'] = sim.dens[:].tolist() +d['TEMP'] = sim.temp[:].tolist() +d['DLNNDR'] = sim.dlnndr[:].tolist() + +with open('json.cgyro.localdump','w') as f: + json.dump(d,f,indent=2,sort_keys=True) + +#-------------------------------------------------------------- + +lref=sim.rmaj +vthref = np.sqrt(2.0) + +d = {} +d['r_minor_norm'] = sim.rmin/lref +d['magnetic_shear_r_minor'] = sim.shear +d['q'] = sim.q +d['charge_norm'] = sim.z[:].tolist() +d['mass_norm'] = sim.mass[:].tolist() +d['density_norm'] = sim.dens[:].tolist() +d['temperature_norm'] = sim.temp[:].tolist() +d['density_log_gradient_norm'] = (lref*sim.dlnndr[:]).tolist() +d['temperature_log_gradient_norm'] = (lref*sim.dlntdr[:]).tolist() +d['velocity_tor_norm'] = + +with open('json.cgyro.imas','w') as f: + json.dump(d,f,indent=2,sort_keys=True) + diff --git a/cgyro/tools/zf/getdata.py b/cgyro/tools/zf/getdata.py deleted file mode 100644 index 2a1693b0f..000000000 --- a/cgyro/tools/zf/getdata.py +++ /dev/null @@ -1,38 +0,0 @@ -# file processed by 2to3 -from __future__ import print_function, absolute_import -from builtins import map, filter, range -import sys -import numpy as np -import matplotlib.pyplot as plt -from gacodefuncs import * -from cgyro.data import cgyrodata - -# Get data directory from command line -data_dir=sys.argv[1] - -# Read data using cgyro data class -data = cgyrodata(data_dir+'/') - -print("Time vector:") -print(data.t) - -print() -print("Theta vector:") -print(data.theta) - -data.getgeo() -print() -print("B(theta):") -print(data.geo[:,2]) - -data.getdata() -print() -print("B_unit:") -print(data.b_unit) - -print() -print("Re[phi(theta)] at last time:") -print(data.phib[0,:,-1]) -print() -print("Im[phi(theta)] at last time:") -print(data.phib[1,:,-1]) diff --git a/cgyro/tools/zf/timeaverage_fields.py b/cgyro/tools/zf/timeaverage_fields.py deleted file mode 100644 index 4dd1fc70c..000000000 --- a/cgyro/tools/zf/timeaverage_fields.py +++ /dev/null @@ -1,98 +0,0 @@ -import sys -import numpy as np -import matplotlib.pyplot as plt -import math -from gacodefuncs import * -from cgyro.data import cgyrodata - -# This is the fraction of the simulation length at the end -# that is considered for the average -frac=0.5 - -# Get data directory from command line -data_dir=sys.argv[1] - -# Read data using cgyro data class -data = cgyrodata(data_dir+'/') - -# print "Time vector:" -# print data.t - -# print -# print "Theta vector:" -# print data.theta - -data.getgeo() - -data.getdata() - -ntheta=len(data.theta) -ntime=len(data.t) - -tmin=math.ceil(ntime*(1-frac)) -tlength=len(data.t[tmin:ntime]) - -rephiaver=math.fsum(list(map(math.fsum, data.phib[0,:,tmin:ntime-1])))/(ntheta*tlength) -reapaver=math.fsum(list(map(math.fsum, data.aparb[0,:,tmin:ntime-1])))/(ntheta*tlength) -rebpaver=math.fsum(list(map(math.fsum, data.bparb[0,:,tmin:ntime-1])))/(ntheta*tlength) - -imphiaver=math.fsum(list(map(math.fsum, data.phib[1,:,tmin:ntime-1])))/(ntheta*tlength) -imapaver=math.fsum(list(map(math.fsum, data.aparb[1,:,tmin:ntime-1])))/(ntheta*tlength) -imbpaver=math.fsum(list(map(math.fsum, data.bparb[1,:,tmin:ntime-1])))/(ntheta*tlength) - - -print("Average Re[delta phi] and Im[delta phi] during time interval") -print(rephiaver) -print(imphiaver) -print("Average Re[delta A||] and Im[delta A||] during time interval") -print(reapaver) -print(imapaver) -print("Average Re[delta B||] and Im[delta B||] during time interval") -print(rebpaver) -print(imbpaver) - -# writing time traces of poloidal averages into files -# Time - -with open("PhiAv.txt","w+") as fphi, open("ApaAv.txt","w+") as fapa, open("BpaAv.txt","w+") as fbpa, open("BpaCos.txt","w+") as cosbpa: - rebpcosaver=0.0 - for i in range(ntime): - rephifsa=math.fsum(data.phib[0,:,i])/ntheta - imphifsa=math.fsum(data.phib[1,:,i])/ntheta - - fphi.write("%11.3e"% (data.t[i])) - fphi.write("%11.3e"% (rephifsa)) - fphi.write("%11.3e\n"% (imphifsa)) - - reapafsa=math.fsum(data.aparb[0,:,i])/ntheta - imapafsa=math.fsum(data.aparb[1,:,i])/ntheta - - fapa.write("%11.3e"% (data.t[i])) - fapa.write("%11.3e"% (reapafsa)) - fapa.write("%11.3e\n"% (imapafsa)) - - rebpafsa=math.fsum(data.bparb[0,:,i])/ntheta - imbpafsa=math.fsum(data.bparb[1,:,i])/ntheta - - fbpa.write("%11.3e"% (data.t[i])) - fbpa.write("%11.3e"% (rebpafsa)) - fbpa.write("%11.3e\n"% (imbpafsa)) - - rebpacos=0.0 - imbpacos=0.0 - for j in range(ntheta): - rebpacos=rebpacos+2.0*math.cos(data.theta[j])*data.bparb[0,j,i]/ntheta - imbpacos=imbpacos+2.0*math.cos(data.theta[j])*data.bparb[1,j,i]/ntheta - - cosbpa.write("%11.3e"% (data.t[i])) - cosbpa.write("%11.3e"% (rebpacos)) - cosbpa.write("%11.3e\n"% (imbpacos)) - - if i >= tmin: - rebpcosaver=rebpcosaver+rebpacos/tlength - -print("Average cos component of Re[delta B||] during time interval") -print(rebpcosaver) - -# print "theta dependence of Re[delta B||] at last time point" -# print data.bparb[0,:,ntime-1]