Moved norm to cgyro_plot

gafusion · Nov 28, 2024 · e1f2958 · e1f2958
1 parent 9d562fb
commit e1f2958
Show file tree

Hide file tree

Showing 5 changed files with 45 additions and 26 deletions.
diff --git a/cgyro/bin/cgyro_plot b/cgyro/bin/cgyro_plot
@@ -260,6 +260,7 @@ EXT=screen
 
 SPECIES=0
 FIELD=0
+FNORM=0
 FC=0
 TIME=-1
 TMAX=-1.0
@@ -324,6 +325,7 @@ while [ $# -gt 0 ] ; do
   -species) shift ; SPECIES=$1 ;;
 
   -field) shift ; FIELD=$1 ;;
+  -fnorm) shift ; FNORM=$1 ;;
 
   -fc) FC=1 ;;
 
@@ -411,7 +413,7 @@ PYROOT=$GACODE_ROOT/f2py/pygacode/cgyro
 
 PLT="-m pygacode.cgyro.data_plot_single $TEX $FONTSIZE"
 
-PFLAGS="$PLOT_TYPE $LX $LY $W $NORM $EXT $TIME $FIELD $MOMENT $TMAX $THETA $YMIN $YMAX $KXMIN $KXMAX $N $ABS $FC $LOC $NSCALE $CFLUX $SPECIES $BAR $IE $MESH"
+PFLAGS="$PLOT_TYPE $LX $LY $W $NORM $EXT $TIME $FIELD $FNORM $MOMENT $TMAX $THETA $YMIN $YMAX $KXMIN $KXMAX $N $ABS $FC $LOC $NSCALE $CFLUX $SPECIES $BAR $IE $MESH"
 
 case "$VIS_TYPE" in
 

diff --git a/cgyro/src/cgyro_init_h.f90 b/cgyro/src/cgyro_init_h.f90
@@ -129,9 +129,11 @@ subroutine cgyro_init_h
                  ir = ir_c(ic) 
                  it = it_c(ic)
                  ang = theta(it)+2*pi*px(ir)
-                 if (amp >  0.0) then
+                 if (amp > 0.0) then
+                    ! symmetric
                     h_x(ic,iv_loc,nt1) = rho/(1.0+ang**4)
                  else
+                    ! anti-symmetric
                     h_x(ic,iv_loc,nt1) = rho*ang/(1.0+ang**4)
                  endif
               enddo

diff --git a/cgyro/src/cgyro_write_timedata.f90 b/cgyro/src/cgyro_write_timedata.f90
@@ -120,10 +120,10 @@ subroutine cgyro_write_timedata
         enddo
 
         if (has_balloon) then
-           if (i_field == 1) then
-              it = maxloc(abs(ftemp(:,n_radial/2+1)),dim=1)
-              a_norm = ftemp(it,n_radial/2+1)
-           endif
+           !if (i_field == 1) then
+           !   it = maxloc(abs(ftemp(:,n_radial/2+1)),dim=1)
+           !   a_norm = ftemp(it,n_radial/2+1)
+           !endif
            call extended_ang(ftemp)
         else
            a_norm = 1.0

diff --git a/f2py/pygacode/cgyro/data_plot.py b/f2py/pygacode/cgyro/data_plot.py
@@ -44,6 +44,7 @@ def plot_dictinit(self):
       xin['bar']    = 0
       xin['ie']     = 0
       xin['mesh']   = 0
+      xin['fnorm']  = 0
 
       return xin
 
@@ -241,6 +242,7 @@ def plot_ball(self,xin):
       itime = xin['itime']
       field = xin['field']
       tmax  = xin['tmax']
+      fnorm = xin['fnorm']
 
       if xin['fig'] is None:
          fig = plt.figure(MYDIR,figsize=(xin['lx'],xin['ly']))
@@ -277,11 +279,23 @@ def plot_ball(self,xin):
          else:
             ax.set_xlim([-tmax,tmax])
 
-      y1 = np.real(f)
-      y2 = np.imag(f)
+
+      # theta_* = 0
+      n0 = self.n_radial//2*self.n_theta+self.n_theta//2
+
+      # Normalized real and imag parts
+      if fnorm == 0:
+         f_norm = self.phib[0,n0,itime]+1j*self.phib[1,n0,itime]
+      elif fnorm == 1:
+         f_norm = self.aparb[0,n0,itime]+1j*self.aparb[1,n0,itime]
+      else:
+         f_norm = self.bparb[0,n0,itime]+1j*self.bparb[1,n0,itime]
+
+      y1 = np.real(f/f_norm)
+      y2 = np.imag(f/f_norm)
 
       ax.plot(x,y1,'-o',color='black',markersize=2,label=r'$\mathrm{Re}$')
-      ax.plot(x,y2,'-o',color='red',markersize=2,label=r'$\mathrm{Im}$')
+      ax.plot(x,y2,'-o',color='red'  ,markersize=2,label=r'$\mathrm{Im}$')
 
       ax.legend()
 

diff --git a/f2py/pygacode/cgyro/data_plot_single.py b/f2py/pygacode/cgyro/data_plot_single.py
@@ -33,23 +33,24 @@
 xin['ftype']  = sys.argv[8]
 xin['itime']  = int(sys.argv[9])
 xin['field']  = int(sys.argv[10])
-xin['moment'] = sys.argv[11]
-xin['tmax']   = float(sys.argv[12])
-xin['theta']  = int(sys.argv[13])
-xin['ymin']   = sys.argv[14]
-xin['ymax']   = sys.argv[15]
-xin['kxmin']  = sys.argv[16]
-xin['kxmax']  = sys.argv[17]
-xin['nstr']   = sys.argv[18]
-xin['abs']    = int(sys.argv[19])
-xin['fc']     = int(sys.argv[20])
-xin['loc']    = int(sys.argv[21])
-xin['nscale'] = int(sys.argv[22])
-xin['cflux']  = sys.argv[23]
-xin['spec']   = int(sys.argv[24])
-xin['bar']    = bool(int(sys.argv[25]))
-xin['ie']     = int(sys.argv[26])
-xin['mesh']   = int(sys.argv[27])
+xin['fnorm']  = int(sys.argv[11])
+xin['moment'] = sys.argv[12]
+xin['tmax']   = float(sys.argv[13])
+xin['theta']  = int(sys.argv[14])
+xin['ymin']   = sys.argv[15]
+xin['ymax']   = sys.argv[16]
+xin['kxmin']  = sys.argv[17]
+xin['kxmax']  = sys.argv[18]
+xin['nstr']   = sys.argv[19]
+xin['abs']    = int(sys.argv[20])
+xin['fc']     = int(sys.argv[21])
+xin['loc']    = int(sys.argv[22])
+xin['nscale'] = int(sys.argv[23])
+xin['cflux']  = sys.argv[24]
+xin['spec']   = int(sys.argv[25])
+xin['bar']    = bool(int(sys.argv[26]))
+xin['ie']     = int(sys.argv[27])
+xin['mesh']   = int(sys.argv[28])
 
 # plotting control
 ftype   = xin['ftype']