New selection logic for global mode

cgyro/bin/cgyro_json

@@ -1,5 +1,6 @@
 #!/usr/bin/env python
 
+import sys
 import json
 import numpy as np
 
@@ -12,13 +13,26 @@ try:
       version = file.read().rstrip()
 except:
    version = 'unavailable'
-
 
-#--------------------------------------------------------------
+if sim.n_species == 1:
+   print('Conversion not implemented for adiabatic electrons.')
+   sys.exit()
+
+# Electron index
+for i in range(sim.n_species):
+   if sim.z[i] < 0.0:
+      ielec = i
+
+#==================================================================
 # json.cgyro.localdump
 
 d = {}
 
+# Code-specific parameters
+d['N_THETA'] = sim.n_theta
+d['N_RADIAL'] = sim.n_theta
+
+# Physics parameters
 d['RMIN'] = sim.rmin
 d['RMAJ'] = sim.rmaj
 d['Q'] = sim.q
@@ -52,8 +66,11 @@ d['SHAPE_COS'] = sim.shape_cos[:].tolist()
 # JSON output
 with open('json.cgyro.localdump','w') as f:
    json.dump(d,f,indent=2,sort_keys=True)
+#==================================================================
 
+#==================================================================
 # standard localdump output
+#
 nmax = 0
 with open('out.cgyro.localdump','w') as f:
    for key in d:
@@ -70,17 +87,24 @@ with open('out.cgyro.localdump','w') as f:
       else:
          f.write(key+'='+str(x)+'\n')
 
-#--------------------------------------------------------------
+#==================================================================
+
+#==================================================================
 # json.cgyro.imas
 
 nmax = nmax+1
 
+# Reference dimensions in CGYRO units
 lref=sim.rmaj
 vthref = np.sqrt(2.0)
 bref = sim.b_gs2
 
+# Initialize dictionary
 d = {}
 
+#--------------------------------------------------------------------
+# Species and profiles
+
 d['name'] = 'CGYRO'
 d['version'] = version
 d['r_minor_norm'] = sim.rmin/lref
@@ -92,15 +116,23 @@ 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()
+#--------------------------------------------------------------------
 
-# Rotation parameters
+#--------------------------------------------------------------------
+# Rotation, collisons, beta
 d['velocity_tor_norm'] = sim.mach/vthref
 # -L/v 1/B (r/q Bu) (-dw/dr) = L/v Bu/B (r/q dw/dr) = L/v Bu/B (-GAMMA_E)
 d['shearing_rate_norm'] = -(lref/vthref)/bref*sim.gamma_e
 # -L^2/v dw/dr = -L/v (L dw/dr) = -L/v (-GAMMA_P)
 d['velocity_tor_gradient_norm'] = lref/vthref*sim.gamma_p
-
+# nu_ee
+d['collisionality_norm'] = (lref/vthref)*sim.nu[ielec]
+# electron beta
 d['beta_reference'] = sim.betae_unit/bref**2
+#-----------------------------------------------------------------------
+
+#-----------------------------------------------------------------------
+# Geometry
 
 d['dgeometric_axis_r_dr_minor'] = sim.shift
 d['dgeometric_axis_z_dr_minor'] = sim.dzmag
@@ -119,7 +151,9 @@ d['dc_dr_minor_norm'] = (sim.shape_s_cos[:nmax]*lref/sim.rmin).tolist()
 d['ds_dr_minor_norm'] = (sim.shape_s_sin[:nmax]*lref/sim.rmin).tolist()
 d['ds_dr_minor_norm'][1] = sim.s_delta*lref/sim.rmin/np.cos(np.arcsin(sim.delta))
 d['ds_dr_minor_norm'][2] = -sim.s_zeta*lref/sim.rmin 
+#-----------------------------------------------------------------------
 
 with open('json.cgyro.imas','w') as f:
    json.dump(d,f,indent=2,sort_keys=True)
+#==================================================================
 

cgyro/bin/cgyro_parse.py

@@ -79,9 +79,9 @@
 x.add('HIPREC_FLAG','0')
 x.add('UDSYMMETRY_FLAG','0')
 x.add('SHEAR_METHOD','2')
+x.add('GLOBAL_FLAG','0')
 x.add('N_GLOBAL','4')
 x.add('NU_GLOBAL','15.0')
-x.add('PROFILE_SHEAR_FLAG','0')
 x.add('THETA_PLOT','1')
 x.add('GPU_BIGMEM_FLAG','1')
 x.add('UPWIND_SINGLE_FLAG','0')
@@ -139,7 +139,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('SBETA','0.0',n=n)
 x.add('DLNNDR_SCALE','1.0',n=n)
 x.add('DLNTDR_SCALE','1.0',n=n)
 

cgyro/src/cgyro_advect_wavenumber.F90

@@ -16,10 +16,9 @@ subroutine cgyro_advect_wavenumber(ij)
   integer :: ir,l,ll,j,iccj,ivc,itor,llnt
   complex :: rl,he1,he2
 
-  if (nonlinear_flag == 0) return
+  if (nonlinear_flag == 0 .or. source_flag == 0) return
 
-  if (source_flag == 1) then
-     call timer_lib_in('shear')
+  call timer_lib_in('shear')
 
 #if defined(OMPGPU)
 !$omp target teams distribute parallel do simd collapse(4) &
@@ -31,14 +30,14 @@ subroutine cgyro_advect_wavenumber(ij)
 #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
+  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
+              ! Wavenumber advection ExB shear
+              if (shear_method == 2) then
                  rl = 0.0
 #if (!defined(OMPGPU)) && defined(_OPENACC)
 !$acc loop seq
@@ -63,10 +62,10 @@ subroutine cgyro_advect_wavenumber(ij)
                     rl = rl+c_wave(l)*(he1-he2)
                  enddo
                  rhs(iccj,ivc,itor,ij) = rhs(iccj,ivc,itor,ij) + omega_eb_base*itor*rl
-           endif
+              endif
 
-           ! Wavenumber advection profile shear
-           if (profile_shear_flag == 1) then
+              ! Wavenumber advection profile shear
+              if (global_flag == 1) then
                  iccj = (ir-1)*n_theta+j
                  rl = rhs(iccj,ivc,itor,ij)
 #if (!defined(OMPGPU)) && defined(_OPENACC)
@@ -91,14 +90,12 @@ subroutine cgyro_advect_wavenumber(ij)
                     rl = rl-c_wave(l)*(he1-he2)
                  enddo
                  rhs(iccj,ivc,itor,ij) = rl
-           endif
-         enddo
-       enddo
-      enddo
+              endif
+           enddo
+        enddo
      enddo
+  enddo
 
-     call timer_lib_out('shear')
-
-  endif
+  call timer_lib_out('shear')
 
 end subroutine cgyro_advect_wavenumber

cgyro/src/cgyro_globals.F90

@@ -91,10 +91,10 @@ module cgyro_globals
   integer :: hiprec_flag
   integer :: udsymmetry_flag
   integer :: shear_method
+  integer :: global_flag
   integer :: n_global
   real    :: nu_global
   integer :: psym_flag
-  integer :: profile_shear_flag
   integer :: theta_plot
   integer :: gpu_bigmem_flag
   integer :: upwind_single_flag
@@ -138,7 +138,7 @@ module cgyro_globals
   real, dimension(11) :: dlntdr
   real, dimension(11) :: sdlnndr
   real, dimension(11) :: sdlntdr
-  real, dimension(11) :: sbeta_star
+  real, dimension(11) :: sbeta
 
   integer :: subroutine_flag  ! only used for cgyro_read_input
 

cgyro/src/cgyro_globalshear.F90

@@ -16,7 +16,7 @@ subroutine cgyro_globalshear(ij)
   integer :: ir,l,ll,j,iccj,ivc,itor,llnt
   complex :: rl,h1,h2
 
-  if (nonlinear_flag == 0) return
+  if (nonlinear_flag == 0 .or. source_flag == 0) return
 
   call timer_lib_in('shear')
 
@@ -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,omega_sbeta,field,h_x,cap_h_c,c_wave)
+!$acc&         present(rhs(:,:,:,ij),omega_ss,omega_sbeta,field,cap_h_c,c_wave)
 #else
 !$omp parallel do collapse(4) private(ivc,ir,l,iccj,j,ll,rl,llnt,h1,h2)
 #endif
@@ -48,22 +48,22 @@ subroutine cgyro_globalshear(ij)
 
                  if ( (ir+ll) <= n_radial ) then
                     ! ExB shear
-                    h1 = omega_eb_base*itor*h_x(iccj+llnt,ivc,itor)
-                    ! omega_star shear
-                    h1 = h1-sum(omega_ss(:,iccj+llnt,ivc,itor)*field(:,iccj+llnt,itor))
+                    h1 = omega_eb_base*itor*cap_h_c(iccj+llnt,ivc,itor)
                     ! beta_star shear
                     h1 = h1-omega_sbeta(iccj+llnt,ivc,itor)*cap_h_c(iccj+llnt,ivc,itor)
+                    ! omega_star shear
+                    h1 = h1-sum(omega_ss(:,iccj+llnt,ivc,itor)*field(:,iccj+llnt,itor))
                  else
                     h1 = 0.0
                  endif
 
                  if ( (ir-ll) >= 1 ) then
                     ! ExB shear
-                    h2 = omega_eb_base*itor*h_x(iccj-llnt,ivc,itor)
-                    ! omega_star shear
-                    h2 = h2-sum(omega_ss(:,iccj-llnt,ivc,itor)*field(:,iccj-llnt,itor))
+                    h2 = omega_eb_base*itor*cap_h_c(iccj-llnt,ivc,itor)
                     ! beta_star shear
                     h2 = h2-omega_sbeta(iccj-llnt,ivc,itor)*cap_h_c(iccj-llnt,ivc,itor)
+                    ! omega_star shear
+                    h2 = h2-sum(omega_ss(:,iccj-llnt,ivc,itor)*field(:,iccj-llnt,itor))
                  else
                     h2 = 0.0
                  endif

cgyro/src/cgyro_init_arrays.F90

@@ -462,7 +462,7 @@ subroutine cgyro_init_arrays
            sm = sdlnndr(is)+sdlntdr(is)*(energy(ie)-1.5)
 
            ! generalized beta/drift shear (acts on H)
-           sb = sbeta_star(is)*energy(ie)/bmag(it)**2
+           sb = sbeta(is)*energy(ie)/bmag(it)**2
 
            arg = -k_theta_base*itor*length/(2*pi)
 

cgyro/src/cgyro_make_profiles.F90

@@ -128,7 +128,7 @@ subroutine cgyro_make_profiles
      dlntdr(1:n_species)  = dlntdr_loc(1:n_species)
      sdlnndr(1:n_species) = sdlnndr_loc(1:n_species)     
      sdlntdr(1:n_species) = sdlntdr_loc(1:n_species)     
-     sbeta_star(1:n_species) = sbeta_loc(1:n_species)     
+     sbeta(1:n_species)   = sbeta_loc(1:n_species)     
 
      if (ae_flag == 1) then
         is_ele = n_species+1
@@ -418,7 +418,10 @@ subroutine cgyro_make_profiles
      case (1)
         call cgyro_info('ExB shear: Hammett discrete shift (do not use for production simulations)') 
      case (2)
-        call cgyro_info('ExB shear: Wavenumber advection') 
+        call cgyro_info('ExB shear: Wavenumber advection (h)') 
+        source_flag = 1
+     case (3)
+        call cgyro_info('ExB shear: Wavenumber advection (H)') 
         source_flag = 1
      case default
         call cgyro_error('Unknown ExB shear method') 
@@ -429,12 +432,13 @@ subroutine cgyro_make_profiles
      call cgyro_info('ExB shear: OFF') 
   endif
 
-  if (profile_shear_flag == 1) then
-     call cgyro_info('Profile shear: Wavenumber advection') 
+  if (global_flag == 1) then
+     call cgyro_info('Global profile variation: ON') 
      source_flag = 1
   else
      sdlnndr(1:n_species) = 0.0
      sdlntdr(1:n_species) = 0.0
+     sbeta(1:n_species)   = 0.0
   endif
 
 

cgyro/src/cgyro_read_input.f90

@@ -82,9 +82,9 @@ subroutine cgyro_read_input
   call cgyro_readbc_int(hiprec_flag,'HIPREC_FLAG')
   call cgyro_readbc_int(udsymmetry_flag,'UDSYMMETRY_FLAG')
   call cgyro_readbc_int(shear_method,'SHEAR_METHOD')
+  call cgyro_readbc_int(global_flag,'GLOBAL_FLAG')
   call cgyro_readbc_int(n_global,'N_GLOBAL')
   call cgyro_readbc_real(nu_global)
-  call cgyro_readbc_int(profile_shear_flag,'PROFILE_SHEAR_FLAG')
   call cgyro_readbc_int(theta_plot,'THETA_PLOT')
   call cgyro_readbc_int(gpu_bigmem_flag,'GPU_BIGMEM_FLAG')
   call cgyro_readbc_int(upwind_single_flag,'UPWIND_SINGLE_FLAG')
@@ -129,7 +129,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(sbeta,is)
   call cgyro_readbc_realv(dlnndr_scale,is)
   call cgyro_readbc_realv(dlntdr_scale,is)
 

cgyro/src/cgyro_rhs.F90

@@ -336,10 +336,14 @@ subroutine cgyro_rhs(ij,update_cap)
   call cgyro_rhs_comm_test
 
   ! Wavenumber advection shear terms
-  ! depends on h_x and field only, updates rhs
-  call cgyro_advect_wavenumber(ij)
-  !call cgyro_globalshear(ij)
-
+  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     
+     call cgyro_globalshear(ij)
+  endif
+
   call cgyro_rhs_comm_test
 
   ! Nonlinear evaluation [f,g]

cgyro/src/cgyro_write_initdata.f90

@@ -130,14 +130,14 @@ subroutine cgyro_write_initdata
      !  k(x) = k(0) + s * r/rho
      !
      ! where x = r/rho. The half-domain has -L/4 < r < L/4.
-     if (profile_shear_flag == 1) then
+     if (global_flag == 1) then
         write(io,*)
         write(io,'(a)') ' i  s(a/Ln)  (a/Ln)_L  (a/Ln)_R  |  s(a/Lt)  (a/Lt)_L  (a/Lt)_R  |  sbeta_*' 
         do is=1,n_species
            dn = sdlnndr(is)*length/rho/4
            dt = sdlntdr(is)*length/rho/4
            write(io,'(t1,i2,3(1x,1pe9.2),2x,3(1x,1pe9.2),2x,1(1x,1pe9.2))') &
-             is,sdlnndr(is),dlnndr(is)-dn,dlnndr(is)+dn,sdlntdr(is),dlntdr(is)-dt,dlntdr(is)+dt,sbeta_star(is)
+             is,sdlnndr(is),dlnndr(is)-dn,dlnndr(is)+dn,sdlntdr(is),dlntdr(is)-dt,dlntdr(is)+dt,sbeta(is)
         enddo
      endif