Merge remote-tracking branch 'origin/Filter_KBM' into SAT3_mtmfix_kbmfix

cgyro/bin/cgyro_plot

@@ -424,7 +424,7 @@ case "$VIS_TYPE" in
 	export OMP_NUM_THREADS=$NOMP ; python -m pygacode.cgyro.vis_wheel $EXT $MOMENT $SPECIES $NX $NY $NZ $TIME $FMIN $FMAX $CMAP $FONTSIZE ;;
 
    supertorus)
-	export OMP_NUM_THREADS=$NOMP ; python -m pygacode.cgyro.vis_supertorus $EXT $MOMENT $SPECIES $NX $NZ $NPHI $PHIMAX $TIME $FMIN $FMAX $CMAP $FONTSIZE $LEGACY_FLAG $DN $MAG $NOZONAL $ONLYZONAL ;;
+	export OMP_NUM_THREADS=$NOMP ; python -m pygacode.cgyro.vis_supertorus $EXT $MOMENT $SPECIES $NX $NZ $NPHI $PHIMAX $TIME $FMIN $FMAX $CMAP $FONTSIZE $LEGACY_FLAG $DN $MAG $NOZONAL $ONLYZONAL $PX $PY ;;
 
     *)
         echo "Unrecognized argument to -vis" ; exit 1 ;;

cgyro/src/Makefile

@@ -21,6 +21,7 @@ OBJECTS = cgyro_globals.o \
           cgyro_parallel_lib.o \
           cgyro_step.o \
           cgyro_io.o \
+          cgyro_restart.o \
           cgyro_advect_wavenumber.o \
 	  cgyro_init_rotation.o \
           cgyro_equilibrium.o \
@@ -58,8 +59,6 @@ OBJECTS = cgyro_globals.o \
 	  cgyro_step_gk_bs5.o \
 	  cgyro_step_gk_v76.o \
           cgyro_write_initdata.o \
-          cgyro_read_restart.o \
-          cgyro_write_restart.o \
           cgyro_write_timedata.o \
           cgyro_write_hosts.o
 

cgyro/src/cgyro_cleanup.F90

@@ -95,7 +95,6 @@ subroutine cgyro_cleanup
   if(allocated(omega_rot_star))      deallocate(omega_rot_star)
   if(allocated(gtime))               deallocate(gtime)
   if(allocated(freq))                deallocate(freq)
-  if(allocated(freq_err))            deallocate(freq_err)
   if(allocated(fcoef))  then
      ccl_del_device(fcoef)     
      deallocate(fcoef)

cgyro/src/cgyro_freq.f90 → cgyro/src/cgyro_freq.F90

@@ -11,6 +11,7 @@
 subroutine cgyro_freq
 
   use cgyro_globals
+  use cgyro_io
 
   implicit none
 
@@ -20,30 +21,41 @@ subroutine cgyro_freq
   integer :: itor
   complex, dimension(nc) :: freq_loc
   complex :: fl,myfr,df,total_weighted_freq
+  complex :: fo1,fo2
+  complex :: icdt
 
   if (i_time == 0) then
 
     freq(:) = 0.0
-    freq_err(:) = 0.0
+    freq_err = 0.0
 
   else
 
     !--------------------------------------------------
     ! Standard method: sum all wavenumbers at a given n
     !--------------------------------------------------
 
+    icdt = i_c/delta_t
+#ifdef __INTEL_COMPILER
+!$omp do ordered
+#endif
     do itor=nt1,nt2
 
      total_weight = 0.0
      total_weighted_freq = (0.0,0.0)
+#ifdef __INTEL_COMPILER
+!$omp do ordered
+#endif
      do ic=1,nc
         ! Use potential to compute frequency
-        mw = abs(field_old(1,ic,itor))
+        fo1 = field_old(1,ic,itor)
+        fo2 = field_old2(1,ic,itor)
+        mw = abs(fo1)
         mode_weight(ic) =  mw
         total_weight = total_weight + mw
         ! Define local frequencies
-        if (abs(field_old(1,ic,itor)) > 1e-12 .and. abs(field_old2(1,ic,itor)) > 1e-12) then
-           fl = (i_c/delta_t)*log(field_old(1,ic,itor)/field_old2(1,ic,itor))
+        if ( (mw > 1e-12) .and. (abs(fo2) > 1e-12) ) then
+           fl = icdt*log(fo1/fo2)
         else
            fl = 0.0
         endif
@@ -53,22 +65,28 @@ subroutine cgyro_freq
 
      myfr = total_weighted_freq/total_weight
      freq(itor) = myfr
-
-     ! Fractional Frequency Error
-     dfr = 0.0
-     dfi = 0.0
-     do ic=1,nc
-        mw = mode_weight(ic)
-        df = freq_loc(ic)-myfr
-        dfr = dfr + abs(real(df))*mw
-        dfi = dfi + abs(aimag(df))*mw
-     enddo
-     freq_err(itor) = (dfr+i_c*dfi)/total_weight/abs(myfr)
-
+
+     if (n_toroidal == 1) then 
+        ! Fractional Frequency Error
+        dfr = 0.0
+        dfi = 0.0
+        do ic=1,nc
+           mw = mode_weight(ic)
+           df = freq_loc(ic)-myfr
+           dfr = dfr + abs(real(df))*mw
+           dfi = dfi + abs(aimag(df))*mw
+        enddo
+
+        if (abs(myfr) > 1e-12) then
+           freq_err = (dfr+i_c*dfi)/total_weight/abs(myfr)
+        else
+           ! Trap a division-by-zero error and halt
+           call cgyro_error('Underflow in calculation of frequency error')
+        endif
+        if (abs(freq_err) < freq_tol) signal = 1
+     endif
     enddo
 
-    if (n_toroidal == 1 .and. abs(freq_err(nt1)) < freq_tol) signal=1
-
   endif
 
 

cgyro/src/cgyro_globals.F90

@@ -251,8 +251,6 @@ module cgyro_globals
   ! Restart tags
   character(len=8) :: fmt='(I2.2)'
   character(len=6), dimension(100) :: rtag
-  integer, parameter :: restart_header_size = 1024
-  integer  :: restart_magic
   !
   ! error checking
   integer :: error_status = 0
@@ -284,7 +282,7 @@ module cgyro_globals
   real    :: t_current
   real, dimension(:), allocatable    :: gtime
   complex, dimension(:), allocatable :: freq
-  complex, dimension(:), allocatable :: freq_err
+  complex :: freq_err
   integer(KIND=8) :: kernel_start_time, kernel_exit_time, kernel_count_rate, kernel_count_max
   !---------------------------------------------------------------
 

cgyro/src/cgyro_init_h.f90

@@ -3,11 +3,12 @@ subroutine cgyro_init_h
   use mpi
   use cgyro_globals
   use cgyro_io
+  use cgyro_restart
 
   implicit none
 
   integer :: ir,it,is,ie,ix,itor
-  real :: arg, ang
+  real :: arg,ang
 
   !---------------------------------------------------------------------------
   ! Check to see if we have restart data available
@@ -49,24 +50,16 @@ subroutine cgyro_init_h
 
      call cgyro_info('Restart data found.')
      call cgyro_read_restart
-     if (error_status /=0 ) return
+     if (error_status > 0) return
+
      gtime(:) = 0.0
 
   case (2)
 
      call cgyro_info('Initializing with restart data.')
      call cgyro_read_restart
-
-     !call MPI_ALLREDUCE(sum(abs(h_x)), &
-     !     arg, &
-     !     1, &
-     !     MPI_DOUBLE_PRECISION, &
-     !     MPI_SUM, &
-     !     NEW_COMM_1, &
-     !     i_err)
-     !h_x = h_x/arg
+     if (error_status > 0) return
 
-     if (error_status /=0 ) return
      i_current = 0
      t_current = 0.0
      gtime(:) = 0.0

cgyro/src/cgyro_init_kernel.F90

@@ -58,22 +58,20 @@ subroutine cgyro_init_kernel
      mpi_dt = (aftermpi_time-kernel_start_time+kernel_count_max)/real(kernel_count_rate)
   endif
 
-    ! 2. Profile setup
+  ! 2. Profile setup
   call cgyro_make_profiles
-  !if (error_status > 0) call cgyro_final_kernel
 
   ! 3. Parameter consistency checks
   call cgyro_check
   if (error_status > 0) then
-          call cgyro_final_kernel
-          return
+     call cgyro_final_kernel
+     return
   endif
 
   ! 4. Array initialization and construction
   !    NOTE: On exit, field_old = field 
-
   call cgyro_init_manager
-  if (error_status /=0 ) then
+  if (error_status > 0) then
      ! something went terribly wrong, hard abort, as things may be
      ! in weird state
      call abort

cgyro/src/cgyro_init_manager.F90

@@ -26,7 +26,6 @@ subroutine cgyro_init_manager
 
   implicit none
 
-
   character(len=128) :: msg
   integer :: ie,ix
 
@@ -152,7 +151,7 @@ subroutine cgyro_init_manager
 
   allocate(gtime(nt1:nt2))
   allocate(freq(nt1:nt2))
-  allocate(freq_err(nt1:nt2))
+  freq_err = 0
 
   if (test_flag == 0) then
 
@@ -335,7 +334,7 @@ subroutine cgyro_init_manager
   endif
 
   call cgyro_equilibrium
-  if (error_status /=0 ) then
+  if (error_status > 0) then
      ! something went terribly wrong
      return
   endif
@@ -348,15 +347,15 @@ subroutine cgyro_init_manager
 
      call cgyro_init_arrays
      call timer_lib_out('str_init')
-     if (error_status /=0 ) then
+     if (error_status > 0) then
         ! something went terribly wrong
         return
      endif
 
      call timer_lib_in('coll_init')
      call cgyro_init_collision
      call timer_lib_out('coll_init')
-     if (error_status /=0 ) then
+     if (error_status > 0) then
         ! something went terribly wrong
         return
      endif
@@ -381,14 +380,6 @@ subroutine cgyro_init_manager
 
   call cgyro_check_memory(trim(path)//runfile_memory)
 
-  if (velocity_order == 1) then
-    ! traditional ordering
-    restart_magic = 140906808
-  else
-    ! alternative ordering, need different magic
-    restart_magic = 140916753
-  endif
-
   call timer_lib_out('str_init')
 
   ! Write initial data
@@ -402,7 +393,7 @@ subroutine cgyro_init_manager
   ! Initialize h (via restart or analytic IC)
   call timer_lib_in('str_init')
   call cgyro_init_h
-  if (error_status /=0 ) return
+  if (error_status > 0) return
   call timer_lib_out('str_init')
 
   ! Initialize nonlinear dimensions and arrays 

cgyro/src/cgyro_kernel.F90

@@ -17,6 +17,7 @@ subroutine cgyro_kernel
   use cgyro_globals
   use cgyro_step
   use cgyro_io
+  use cgyro_restart
 
   implicit none
 
@@ -91,9 +92,9 @@ subroutine cgyro_kernel
 #endif
 
      if (mod(i_time,print_step) == 0) then
-       ! cap_h_c will not be modified in GPU memory for the rest of the loop
+        ! cap_h_c will not be modified in GPU memory for the rest of the loop
 #if defined(OMPGPU)
-     ! no async for OMPGPU for now
+        ! no async for OMPGPU for now
 !$omp target update from(cap_h_c)
 #elif defined(_OPENACC)
 !$acc update host(cap_h_c) async(4)
@@ -109,8 +110,8 @@ subroutine cgyro_kernel
      ! NOTE: Fluxes are calculated in cgyro_write_timedata
 
 #if (!defined(OMPGPU)) && defined(_OPENACC)
-  call timer_lib_in('coll_mem')
-  ! wait for fields to be synched into system memory, used by cgyro_error_estimate
+     call timer_lib_in('coll_mem')
+     ! wait for fields to be synched into system memory, used by cgyro_error_estimate
 !$acc wait(3)
   call timer_lib_out('coll_mem')
 #endif
@@ -125,13 +126,13 @@ subroutine cgyro_kernel
      ! IO
      !
      if (mod(i_time,print_step) == 0) then
-       call timer_lib_in('coll_mem')
+        call timer_lib_in('coll_mem')
 #if defined(OMPGPU)
-     ! no async for OMPGPU for now
+        ! no async for OMPGPU for now
 !$omp target update from(cap_h_c_dot)
 #elif defined(_OPENACC)
 !$acc update host(cap_h_c_dot)
-       ! wait for cap_h_c to be synched into system memory, used by cgyro_write_timedata
+        ! wait for cap_h_c to be synched into system memory, used by cgyro_write_timedata
 !$acc wait(4)
 #endif
        call timer_lib_out('coll_mem')
@@ -157,8 +158,8 @@ subroutine cgyro_kernel
      ! Don't wrap timer output in a timer
      if (mod(i_time,print_step) == 0) call write_timers(trim(path)//runfile_timers)
 
-     ! Exit if convergenced
-     if (signal == 1) exit
+     ! Exit if converged (1) or underflow (2)
+     if (signal > 0) exit
 
   enddo
   !---------------------------------------------------------------------------