Qlmaster

Makefile

@@ -10,11 +10,13 @@ all:
 	cd tgyro ; make
 	cd profiles_gen ; make
 	cd f2py ; make
+	cd qlgyro ; make
 	@echo "GACODE build done"
 
 clean:
 	rm -f .VERSION
 	cd shared ; make clean
+	cd shared/hybridtest ; make clean
 	cd neo ; make clean
 	cd vgen ; make clean
 	cd tglf ; make clean
@@ -23,6 +25,7 @@ clean:
 	cd tgyro ; make clean
 	cd profiles_gen ; make clean
 	cd f2py ; make clean
+	cd qlgyro ; make clean
 	rm -f f2py/*/*.pyc
 	rm -f f2py/*.pyc
 	rm -f f2py/pygacode/*/*.pyc
@@ -43,6 +46,7 @@ distclean:
 	cd tgyro ; make clean
 	cd profiles_gen ; make clean
 	cd f2py ; make clean
+	cd qlgyro ; make clean
 	rm -f python/*/*.pyc
 	rm -f python/*.pyc
 	rm -f modules/*genmod*

README.rst

@@ -1,4 +1,4 @@
 GACODE repository
 =================
 
-For documentation, please refer to the `GACODE website <http://gafusion.github.io/doc/>`_.
+For documentation, please refer to the `GACODE website <http://gacode.io>`_.

cgyro/src/cgyro_write_initdata.f90

@@ -13,7 +13,7 @@ subroutine cgyro_write_initdata
   implicit none
 
   integer :: p,in,is,it
-  real :: kymax,kyrat,dn
+  real :: kymax,kyrat,dn,dt
   real, external :: spectraldiss
   character(len=50) :: msg
 
@@ -120,13 +120,25 @@ subroutine cgyro_write_initdata
 
      write(io,*)
      write(io,'(a)') &
-          ' i  z  n/n_norm   T/T_norm   m/m_norm     a/Ln       a/Lt       nu        s(a/Ln)    s(a/Lt)'
+          ' i  z  n/n_norm   T/T_norm   m/m_norm     a/Ln       a/Lt       nu'
      do is=1,n_species
         write(io,'(t1,i2,1x,i2,3(2x,1pe9.3),2(1x,1pe10.3),(2x,1pe9.3),2(1x,1pe10.3))') &
-             is,int(z(is)),dens(is),temp(is),mass(is),dlnndr(is),dlntdr(is),nu(is),&
-             sdlnndr(is),sdlntdr(is)
+             is,int(z(is)),dens(is),temp(is),mass(is),dlnndr(is),dlntdr(is),nu(is)
      enddo
 
+     ! Profile shear
+     if (profile_shear_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 ' 
+        do is=1,n_species
+           dn = sdlnndr(is)*length/rho/4
+           dt = sdlntdr(is)*length/rho/4
+           write(io,'(t1,i2,3(1x,1pe9.3),2x,3(1x,1pe9.3))') &
+             is,sdlnndr(is),dlnndr(is)-dn,dlnndr(is)+dn,sdlntdr(is),dlntdr(is)-dt,dlntdr(is)+dt
+        enddo
+     endif
+
+     ! Running from input.gacode
      if (profile_model == 2) then
         dn = rho/(rhos/a_meters)
         kyrat = abs(q/rmin*rhos/a_meters)

f2py/expro/expro.f90

@@ -180,6 +180,7 @@ module expro
   integer :: expro_ctrl_quasineutral_flag 
   integer :: expro_ctrl_numeq_flag
   integer :: expro_error=0
+  integer :: expro_jerr=0
 
   ! Header information
   character(len=70) :: expro_head_original =  '#  *original : null'

f2py/expro/expro_util.f90

@@ -208,7 +208,7 @@ subroutine expro_compute_derived
 
   ! r_min = a [m]
   r_min = expro_rmin(expro_n_exp)
-
+  
   do i=2,expro_n_exp
 
      ! Parameters to be passed to geo library   
@@ -262,9 +262,11 @@ subroutine expro_compute_derived
         geo_fourier_in(1:4,0:geo_nfourier_in) = expro_geo(:,:,i)/r_min
         geo_fourier_in(5:8,0:geo_nfourier_in) = expro_dgeo(:,:,i)
      endif
+
      call geo_interp(1,theta,.true.)
      if (minval(geov_jac_r) <= 0d0) then
-        print '(a,i3,a)','WARNING: (expro_util) Negative Jacobian for i =',i,' in input.gacode'
+        print '(a,i3,a,f5.2)','WARNING: (expro_util) J < 0 for i=',i,' in input.gacode: ',minval(geov_jac_r)
+        expro_jerr = 1
      endif
 
      ! V, dV/dr and S (note that S=dV/dr only in a circle)
@@ -415,9 +417,9 @@ subroutine expro_compute_derived
         expro_ni_new(:) = expro_ni_new(:)+expro_z(is)*expro_ni(is,:)
         expro_dlnnidr_new(:) = expro_dlnnidr_new(:)+expro_z(is)*expro_ni(is,:)*expro_dlnnidr(is,:)
      enddo
-     
+
      ! New quasineutral ion 1 profiles:
-     
+
      ! density  
      expro_ni_new(:) = (expro_ne(:)-expro_ni_new(:))/expro_z(1)
      ! gradient scale length 1/L_ni = -dln(ni)/dr (1/m)
@@ -427,7 +429,7 @@ subroutine expro_compute_derived
      call bound_deriv(expro_sdlnnidr_new(:),expro_ni_new(:)*expro_dlnnidr_new(:),&
           expro_rmin,expro_n_exp)
      expro_sdlnnidr_new(:) = expro_sdlnnidr_new(:)/expro_ni_new(:)*expro_rhos(:)
-     
+
      if (minval(expro_ni_new(:)) <= 0d0) expro_error = 1
 
   else
@@ -498,7 +500,7 @@ subroutine expro_compute_derived
 
   ! Reset geo variables to standard values to prevent expro from leaving
   ! geo in a "mysterious" state
-  
+
   geo_rmin_in      = 1d0
   geo_rmaj_in      = 0.5d0
   geo_drmaj_in     = 3d0

f2py/pygacode/cgyro/data.py

@@ -201,10 +201,7 @@ def getxflux(self):
 
    def xfluxave(self,w,moment,e=0.2,nscale=0):
 
-      """
-      Do complicated spatial averages for xflux
-      RESULT: self.lky_flux_ave
-      """
+      # Averaging functionfor global fluxes
 
       print('INFO: (xfluxave) Computing partial-domain averages')
 
@@ -221,41 +218,39 @@ def xfluxave(self,w,moment,e=0.2,nscale=0):
       else:
          sc[:] = 1.0
 
+      # Sum moments over fields (3) and toroidal modes (4)
+      # NOTE: lky_flux_n[2,ng,ns,nfield,n_n,nt]
       if moment == 'n':
          z = np.sum(self.lky_flux_n,axis=(3,4))
       elif moment == 'e':
          z = np.sum(self.lky_flux_e,axis=(3,4))
       elif moment == 'v':
          z = np.sum(self.lky_flux_v,axis=(3,4))
-      else:
-         raise ValueError('(xfluxave) Invalid moment.')
-
+
       #--------------------------------------------
-      # Useful arrays required outside this routine
-      self.lky_xr = np.zeros((ns,ng))
-      self.lky_xi = np.zeros((ns,ng))
-      self.lky_flux_ave = np.zeros((ns,2))
+      # Arrays required outside this routine
+      self.lky_xr = np.zeros([ng,ns])
+      self.lky_xi = np.zeros([ng,ns])
+      self.lky_flux_ave = np.zeros([ns,2])
       #--------------------------------------------
 
-      imin,imax=time_index(self.t,w)
-
-      for ispec in range(ns):
-         for l in range(ng):
-            self.lky_xr[ispec,l] = time_average(z[0,l,ispec,:],self.t,imin,imax)*sc[ispec]
-            self.lky_xi[ispec,l] = time_average(z[1,l,ispec,:],self.t,imin,imax)*sc[ispec]
+      imin,imax = time_index(self.t,w)
+
+      # Time averages of real and imaginary parts
+      self.lky_xr[:,:] = time_average(z[0,:,:,:],self.t,imin,imax)
+      self.lky_xi[:,:] = time_average(z[1,:,:,:],self.t,imin,imax)
 
+      l = np.arange(1,ng)
+      u = (2*np.pi*e)*l
+      for ispec in range(ns):
          # Flux partial average over [-e,e]
-         g0 = self.lky_xr[ispec,0]
-         g1 = g0
-         for l in range(1,ng):
-            u = 2*np.pi*l*e
-            g0 = g0+2*np.sin(u)*self.lky_xr[ispec,l]/u
-            g1 = g1+2*np.sin(u)*self.lky_xr[ispec,l]/u*(-1)**l
+         g0 = self.lky_xr[0,ispec]+2*np.sum(np.sin(u)*self.lky_xr[l,ispec]/u)
+         g1 = self.lky_xr[0,ispec]+2*np.sum(np.sin(u)*self.lky_xr[l,ispec]/u*(-1)**l)
 
          # Average over true (positive) interval
-         self.lky_flux_ave[ispec,0] = g0
+         self.lky_flux_ave[ispec,0] = g0*sc[ispec]
          # Average over negative interval
-         self.lky_flux_ave[ispec,1] = g1
+         self.lky_flux_ave[ispec,1] = g1*sc[ispec]
 
    def getbigfield(self):
 

f2py/pygacode/cgyro/data_plot.py

@@ -1006,40 +1006,29 @@ def plot_xflux(self,xin):
       nscale = xin['nscale']
       ymin   = xin['ymin']
       ymax   = xin['ymax']
+      mirror = xin['abs']
 
       if xin['fig'] is None:
          fig = plt.figure(MYDIR,figsize=(xin['lx'],xin['ly']))
 
       self.getxflux()
 
       ns = self.n_species
-      nl = self.n_global+1
+      ng = self.n_global+1
       t  = self.t
 
-      ky  = self.ky
-      ave = np.zeros((self.n_n,ns))
-
       if moment == 'phi':
          moment = 'e'
-
-      # NOTE: lky_flux_* -> [ 2, nl , ns , n_n , nt ]
-      #                       0  1    2     3    4 
-
+
       if moment == 'n':
          ntag = 'Density~flux'
          mtag = '\Gamma'
-         z = np.sum(self.lky_flux_n,axis=3)
-         ftag = 'xflux_n'
       elif moment == 'e':
          ntag = 'Energy~flux'
          mtag = 'Q'
-         z = np.sum(self.lky_flux_e,axis=3)
-         ftag = 'xflux_e'
       elif moment == 'v':
          ntag = 'Momentum~flux'
          mtag = '\Pi'
-         z = np.sum(self.lky_flux_v,axis=3)
-         ftag = 'xflux_v'
       else:
          print('ERROR: (plot_xflux) Invalid moment.')
          sys.exit()
@@ -1055,7 +1044,6 @@ def plot_xflux(self,xin):
       else:
          mnorm = ''
 
-
       #============================================================
       # Otherwise plot
       ax = fig.add_subplot(111)
@@ -1065,51 +1053,53 @@ def plot_xflux(self,xin):
 
       color = ['k','m','b','c','g','r']
 
-      imin,imax=time_index(t,w)
+      imin,imax = time_index(t,w)
       mpre,mwin = wintxt(imin,imax,t)
 
       ax.set_title(r'$\mathrm{'+ntag+'} \quad $'+mwin)
-
-      a = -np.pi+2*np.pi*np.arange(0.0,1.0,0.001)
+
+      na = 128
+      a  = np.linspace(-np.pi,np.pi,na)
+      ah = np.linspace(0,-2*np.pi,na)
 
       for ispec in range(ns):
 
          u = specmap(self.mass[ispec],self.z[ispec])
 
-         # Flux curve
-         g = np.zeros(len(t))
-         g = self.lky_xr[ispec,0] 
-         for l in range(1,nl):
-            g = g+2*(np.cos(l*a)*self.lky_xr[ispec,l]-np.sin(l*a)*self.lky_xi[ispec,l])
+         #---------------------------------
+         # Global flux versus x
+         g = np.zeros(na) ; g[:] = self.lky_xr[0,ispec]          
+         for l in range(1,ng):
+            g[:] = g[:]+2*(np.cos(l*a)*self.lky_xr[l,ispec]-np.sin(l*a)*self.lky_xi[l,ispec])
          ax.plot(a/(2*np.pi),g,color=color[ispec])
-
+         g = np.zeros(na) ; g[:] = self.lky_xr[0,ispec]          
+         for l in range(1,ng):
+            g[:] = g[:]+2*(np.cos(l*ah)*self.lky_xr[l,ispec]-np.sin(l*ah)*self.lky_xi[l,ispec])
+         ax.plot(a/(2*np.pi),g,color=color[ispec],linestyle='--')
          #---------------------------------
-         # Flux partial average over [-e,e]
+
+         #---------------------------------
+         # Flux partial average over "positive" [-e,e] interval
          g0 = self.lky_flux_ave[ispec,0]
          label = r'$'+mtag+mnorm+'_'+u+'/'+mtag+'_\mathrm{GB}: '+str(round(g0,3))+'$'
          ax.plot([-e,e],[g0,g0],'o-',color=color[ispec],alpha=0.2,linewidth=3,label=label)
          #---------------------------------
 
          #---------------------------------
-         # Flux partial average over "negative" interval
+         # Flux partial average over "negative" [-e,e] interval
          g1 = self.lky_flux_ave[ispec,1]
          ax.plot([0.5-e,0.5],[g1,g1],'o--',color=color[ispec],alpha=0.2,linewidth=3)
          ax.plot([-0.5,-0.5+e],[g1,g1],'o--',color=color[ispec],alpha=0.2,linewidth=3)
          #---------------------------------
-
-         #---------------------------------
-         # Flux spectral average
-         gs = self.lky_xr[ispec,0]+2*np.pi/4*self.lky_xr[ispec,1]
-         #---------------------------------
-
+
          #---------------------------------
          # Flux domain average
-         ga = self.lky_xr[ispec,0]
+         ga = self.lky_xr[0,ispec]
          ax.plot([-0.5,0.5],[ga,ga],color=color[ispec],alpha=0.5)
          #---------------------------------
 
-         print('INFO: (plot_xflux) Ave [inner/inner_spec, outer, domain] = '
-               '{:.2f}/{:.2f}, {:.2f}, {:.2f}'.format(g0,gs,g1,ga)) 
+         print('INFO: (plot_xflux) Ave [inner-e, outer-e, domain] = '
+               '{:.2f}, {:.2f}, {:.2f}'.format(g0,g1,ga)) 
 
          if ymax != 'auto':
             ax.set_ylim(top=float(ymax))

f2py/pygacode/profiles_gen/prgen_shapeprofile.py

@@ -8,18 +8,14 @@
 from .prgen_contour import *
 from .prgen_shape import *
 
-if len(sys.argv) > 1:
-   gfile   = sys.argv[1]
-   mag     = sys.argv[2]
-   narc    = int(sys.argv[3])
-   npsi    = int(sys.argv[4])
-   nharm   = int(sys.argv[5])
-   nfourier = int(sys.argv[6])
-   plotpng = bool(int(sys.argv[7]))
-   psinorm = float(sys.argv[8])
-else:
-   print('Usage: python prgen_shapeprofile.py <gfile> <mag> <narc> <npsi> <nfourier>')
-   sys.exit()
+gfile = sys.argv[1]
+mag   = sys.argv[2]
+narc  = int(sys.argv[3])
+npsi  = int(sys.argv[4])
+nharm = int(sys.argv[5])
+nfourier = int(sys.argv[6])
+plotpng  = bool(int(sys.argv[7]))
+psinorm  = float(sys.argv[8])
 
 efit = prgen_geqdsk(gfile)
 psi0 = efit['SIMAG']

neo/src/neo_make_profiles.f90

@@ -441,6 +441,7 @@ subroutine neo_make_profiles
      enddo
      close(io)
 
+     ! JC: New file, introduced in Dec 2023 for neo_nice
      open(unit=io,file=trim(path)//'out.neo.species',status='replace')
      do is=1,n_species
         write (io,'(e16.8)',advance='no') mass(is)

platform/build/make.inc.CORG

@@ -1,5 +1,5 @@
 #---------------------------------------------------
-# Linux Mint 
+# CORG (Ubuntu server)
 #---------------------------------------------------
 
 IDENTITY="corg.gat.com gfortran+mpich+OpenBLAS (/usr)"

platform/build/make.inc.OSX_CATALINA

@@ -17,7 +17,7 @@
 
 # Compilers and flags
 
-FC  = mpif90-mpich-gcc9 -I$(GACODE_ROOT)/modules -J$(GACODE_ROOT)/modules -fPIC -framework Accelerate
+FC  = mpif90-mpich-gcc11 -I$(GACODE_ROOT)/modules -J$(GACODE_ROOT)/modules -fPIC -framework Accelerate
 F77 = mpif77 -std=legacy
 CC  = ${FC}
 

platform/exec/exec.MINT

@@ -9,5 +9,5 @@ numa=${5}
 mpinuma=${6}
 
 cd $simdir
-mpiexec -env OMP_NUM_THREADS $nomp -n $nmpi $exec
+mpiexec -env OMP_NUM_THREADS $nomp -n $nmpi $exec