Merge branch 'master' of github.com:gafusion/gacode

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.

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)
+

cgyro/bin/cgyro_parse.py

@@ -139,6 +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('DLNNDR_SCALE','1.0',n=n)
 x.add('DLNTDR_SCALE','1.0',n=n)
 

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 \
@@ -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 \

cgyro/src/cgyro_cleanup.F90

@@ -197,6 +197,10 @@ subroutine cgyro_cleanup
      ccl_del_device(omega_ss)        
      deallocate(omega_ss)
   endif
+  if(allocated(omega_sbeta)) then
+     ccl_del_device(omega_sbeta)        
+     deallocate(omega_sbeta)
+  endif
   if(allocated(jvec_c))  then
      ccl_del_device(jvec_c)     
      deallocate(jvec_c)

cgyro/src/cgyro_globals.F90

@@ -138,6 +138,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
 
@@ -351,6 +352,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_sbeta
   complex, dimension(:,:,:), allocatable :: cap_h_c
   complex, dimension(:,:,:), allocatable :: cap_h_c_dot
   complex, dimension(:,:,:), allocatable :: cap_h_c_old

cgyro/src/cgyro_globalshear.F90

@@ -0,0 +1,78 @@
+!---------------------------------------------------------
+! cgyro_advect_wavenumber.f90
+!
+! PURPOSE:
+!  Manage shearing by wavenumber advection.
+!---------------------------------------------------------
+
+subroutine cgyro_globalshear(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,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
+  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   = 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
+
+                 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-omega_sbeta(iccj+llnt,ivc,itor)*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-omega_sbeta(iccj-llnt,ivc,itor)*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_globalshear