Added delta_t_gk to out.cgyro.tag to ensure deterministic restart

gafusion · Sep 17, 2024 · 4c56cd0 · 4c56cd0
1 parent c2d5e58
commit 4c56cd0
Show file tree

Hide file tree

Showing 3 changed files with 25 additions and 14 deletions.
diff --git a/cgyro/src/cgyro_init_kernel.F90 b/cgyro/src/cgyro_init_kernel.F90
@@ -120,6 +120,12 @@ subroutine cgyro_init_kernel
   call timer_lib_out('str_mem')
 
   ! Initialize adaptive time-stepping parameter
-  delta_t_gk = delta_t
+  if (delta_t_last > 0.0) then
+     ! Saved from restart data
+     delta_t_gk = delta_t_last
+  else
+     ! New run
+     delta_t_gk = delta_t
+  endif
 
 end subroutine cgyro_init_kernel
diff --git a/cgyro/src/cgyro_restart.F90 b/cgyro/src/cgyro_restart.F90
@@ -25,6 +25,7 @@ subroutine cgyro_write_restart
   use mpi
   use cgyro_globals
   use cgyro_io
+  use cgyro_step
 
   implicit none
 
@@ -53,6 +54,7 @@ subroutine cgyro_write_restart
      open(unit=io,file=trim(path)//runfile_restart_tag,status='replace')
      write(io,*) i_current
      write(io,fmtstr) t_current
+     write(io,*) delta_t_gk
      close(io)
   endif
 
@@ -350,31 +352,34 @@ subroutine cgyro_read_restart
   use mpi
   use cgyro_globals
   use cgyro_io
+  use cgyro_step
 
   implicit none
 
+  integer :: igk
+
   !---------------------------------------------------------
-  ! Read restart parameters from ASCII file.
+  ! Read restart parameters from ASCII tag file.
   !
   if (restart_flag == 1) then
+     delta_t_last = 0.0
      if (i_proc == 0) then
 
         open(unit=io,file=trim(path)//runfile_restart_tag,status='old')
 
         read(io,*) i_current
         read(io,fmtstr) t_current
+        read(io,*,iostat=igk) delta_t_last
         close(io)
-
+        
      endif
 
-     ! Broadcast to all cores.
-
-     call MPI_BCAST(i_current,&
-          1,MPI_INTEGER,0,CGYRO_COMM_WORLD,i_err)
-
-     call MPI_BCAST(t_current,&
-          1,MPI_DOUBLE_PRECISION,0,CGYRO_COMM_WORLD,i_err)
+     ! Broadcast to all cores
 
+     call MPI_BCAST(i_current,1,MPI_INTEGER,0,CGYRO_COMM_WORLD,i_err)
+     call MPI_BCAST(t_current,1,MPI_DOUBLE_PRECISION,0,CGYRO_COMM_WORLD,i_err)
+     call MPI_BCAST(delta_t_last,1,MPI_DOUBLE_PRECISION,0,CGYRO_COMM_WORLD,i_err)
+
   endif
 
   if (i_proc == 0) then

diff --git a/cgyro/src/cgyro_write_timedata.f90 b/cgyro/src/cgyro_write_timedata.f90
@@ -135,7 +135,7 @@ subroutine cgyro_write_timedata
   endif
   !---------------------------------------------------------------
 
-  ! Linear frequency diagnostics for every value of n
+  ! out.cgyro.freq: linear frequency diagnostics for every value of n
   call cgyro_freq
   fvec(1,:) = real(freq(:)) ; fvec(2,:) = aimag(freq(:))
   if (n_toroidal > 1) then
@@ -148,7 +148,7 @@ subroutine cgyro_write_timedata
   ! Output to screen
   if (printout) call print_scrdata()
 
-  ! Output to files
+  ! out.cgyro.time: time,total_error,rk_error,delta_t
   vec(1) = t_current ; vec(2:3) = integration_error(:) ; vec(4) = delta_t_gk
   call write_ascii(trim(path)//runfile_time,4,vec(1:4))
 
@@ -567,11 +567,11 @@ subroutine write_ascii(datafile,n_fn,fn)
 
      open(unit=io,file=datafile,status='old')
      do i=1,i_current
-        read(io,fmtstr) fn(1)
+        read(io,fmtstrn) fn(:)
      enddo
      endfile(io)
      close(io)
-
+     
   end select
 
 end subroutine write_ascii