From 6c7cd43ff5a60553fb2fa5c4501f3151935b824c Mon Sep 17 00:00:00 2001 From: "jiyoon2421@gmail.com" Date: Tue, 14 May 2024 16:17:20 +0900 Subject: [PATCH 01/16] Updated platform files for NEURON - Used nvidia_hpc_sdk compiler for NEURON - Found that CUDA 11.0 and cc80 worked for NEURON --- platform/build/make.inc.NEURON | 3 +-- platform/env/env.NEURON | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/platform/build/make.inc.NEURON b/platform/build/make.inc.NEURON index 930d9cbfa..e17bc0c73 100644 --- a/platform/build/make.inc.NEURON +++ b/platform/build/make.inc.NEURON @@ -14,12 +14,11 @@ NUMAS_PER_NODE=1 # Compilers FC = mpif90 -module ${GACODE_ROOT}/modules -Mpreprocess -DUSE_INLINE -# -I/apps/compiler/pgi/linux86-64/19.1/cudampi/10.0/openmpi/2.3/applib2/ F77 = ${FC} # Compiler options / flags -FACC = -acc -Minfo=accel -Mcuda=cuda10.0 -ta=nvidia:cc70 -Mcudalib=cufft +FACC = -acc -Minfo=accel -Mcuda=cuda11.0 -ta=nvidia:cc80 -Mcudalib=cufft FOMP = -mp -Mstack_arrays FMATH = -r8 FOPT = -fast diff --git a/platform/env/env.NEURON b/platform/env/env.NEURON index c4b71c97b..5e03afc67 100644 --- a/platform/env/env.NEURON +++ b/platform/env/env.NEURON @@ -1,2 +1,2 @@ module purge -module load pgi/19.1 netcdf/4.6.1 cuda/10.0 cudampi/openmpi-3.1.0 fftw_mpi/3.3.7 python/3.7.1 +module load nvidia_hpc_sdk/22.7 python/3.7.1 From 3785eee3a9623d831c15d7075ac6c8b5ddcc5ecd Mon Sep 17 00:00:00 2001 From: Igor Sfiligoi Date: Tue, 14 May 2024 16:05:45 -0400 Subject: [PATCH 02/16] Add restart read and write into a common file, and make it a module --- cgyro/src/Makefile | 3 +- cgyro/src/cgyro_init_h.f90 | 1 + cgyro/src/cgyro_kernel.F90 | 1 + cgyro/src/cgyro_read_restart.F90 | 233 ----------------- ...ro_write_restart.F90 => cgyro_restart.F90} | 244 +++++++++++++++++- 5 files changed, 245 insertions(+), 237 deletions(-) delete mode 100644 cgyro/src/cgyro_read_restart.F90 rename cgyro/src/{cgyro_write_restart.F90 => cgyro_restart.F90} (56%) diff --git a/cgyro/src/Makefile b/cgyro/src/Makefile index 01f8edcf5..dea9cdb44 100644 --- a/cgyro/src/Makefile +++ b/cgyro/src/Makefile @@ -19,6 +19,7 @@ OBJECTS = cgyro_globals.o \ cgyro_math.o \ cgyro_globals_math.o \ cgyro_parallel_lib.o \ + cgyro_restart.o \ cgyro_step.o \ cgyro_io.o \ cgyro_advect_wavenumber.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 diff --git a/cgyro/src/cgyro_init_h.f90 b/cgyro/src/cgyro_init_h.f90 index 7df95651b..f3a37395b 100644 --- a/cgyro/src/cgyro_init_h.f90 +++ b/cgyro/src/cgyro_init_h.f90 @@ -3,6 +3,7 @@ subroutine cgyro_init_h use mpi use cgyro_globals use cgyro_io + use cgyro_restart implicit none diff --git a/cgyro/src/cgyro_kernel.F90 b/cgyro/src/cgyro_kernel.F90 index 18eb16eb2..10f2d3c8d 100644 --- a/cgyro/src/cgyro_kernel.F90 +++ b/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 diff --git a/cgyro/src/cgyro_read_restart.F90 b/cgyro/src/cgyro_read_restart.F90 deleted file mode 100644 index cb50fb418..000000000 --- a/cgyro/src/cgyro_read_restart.F90 +++ /dev/null @@ -1,233 +0,0 @@ -!------------------------------------------------------ -! cgyro_read_restart.f90 -! -! PURPOSE: -! This is the master file controlling the restart -! via MPI-IO. -!------------------------------------------------------ - -subroutine cgyro_read_restart - - use mpi - use cgyro_globals - use cgyro_io - - implicit none - - !--------------------------------------------------------- - ! Read restart parameters from ASCII file. - ! - if (restart_flag == 1) then - 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 - 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) - - endif - - if (i_proc == 0) then - open(unit=io,file=trim(path)//runfile_restart,status='old',iostat=i_err) - close(io) - endif - - call cgyro_read_restart_one - -end subroutine cgyro_read_restart - - -subroutine cgyro_read_restart_verify - - use cgyro_globals - use cgyro_io - - !--------------------------------------------------- - implicit none - - integer :: magic, version, recid - integer :: t_n_theta,t_n_radial,t_n_species,t_n_xi,t_n_energy,t_n_toroidal - - ! Different compilers have a different semantics for RECL... must use inquire - integer :: reclen - integer, dimension(1) :: recltest - - inquire(iolength=reclen) recltest - - open(unit=io,& - file=trim(path)//runfile_restart,& - status='old',access='DIRECT',RECL=reclen,ACTION='READ') - - recid = 1 - - read(io,REC=recid) magic - recid = recid + 1 - if ( magic /= restart_magic) then - close(io) - call cgyro_error('Wrong magic number in restart header') - return - endif - - read(io,REC=recid) version - recid = recid + 1 - if ( version /= 2) then - close(io) - call cgyro_error('Wrong version in restart header, only v2 supported') - return - endif - - read(io,REC=recid) t_n_theta - recid = recid + 1 - read(io,REC=recid) t_n_radial - recid = recid + 1 - read(io,REC=recid) t_n_species - recid = recid + 1 - read(io,REC=recid) t_n_xi - recid = recid + 1 - read(io,REC=recid) t_n_energy - recid = recid + 1 - read(io,REC=recid) t_n_toroidal - recid = recid + 1 - if ( (t_n_theta/=n_theta) .or. (t_n_radial/=n_radial) .or. & - (t_n_species/=n_species) .or. (t_n_xi/=n_xi) .or. & - (t_n_energy/=n_energy) .or. (t_n_toroidal/=n_toroidal) ) then - close(io) - call cgyro_error('Wrong geometry in restart header') - return - endif - - ! follow MPI params... will ignore them for v2, as they do not change the file format - - close(io) - -end subroutine cgyro_read_restart_verify - -subroutine cgyro_read_restart_one - - use mpi - use cgyro_globals - use cgyro_io - - !--------------------------------------------------- - implicit none - ! - ! Required for MPI-IO: - ! - integer :: filemode - integer :: finfo - integer :: fhv - integer :: fstatus(MPI_STATUS_SIZE) - integer(kind=MPI_OFFSET_KIND) :: disp - integer(kind=MPI_OFFSET_KIND) :: offset1 - !--------------------------------------------------- - - character(8) :: sdate - character(10) :: stime - character(len=64) :: platform - integer(KIND=8) :: start_time,cp_time - integer(KIND=8) :: count_rate, count_max - real :: cp_dt - integer :: ic0,j,statusfd - - ! use system_clock to be consistent with cgyro_kernel - call system_clock(start_time,count_rate,count_max) - - ! First read the header, and verify that it is compatible with current setup - - if (i_proc == 0) then - call cgyro_read_restart_verify - if (error_status /=0 ) return - endif - - filemode = MPI_MODE_RDONLY - disp = 0 - - offset1 = size(h_x,kind=MPI_OFFSET_KIND)*(i_proc_1+i_proc_2*n_proc_1) + restart_header_size - if (offset1 < restart_header_size) then - call cgyro_error('Overflow detected in cgyro_read_restart_one') - return - endif - - call MPI_INFO_CREATE(finfo,i_err) - - call MPI_FILE_OPEN(CGYRO_COMM_WORLD,& - trim(path)//runfile_restart,& - filemode,& - finfo,& - fhv,& - i_err) - - if (i_err /= 0) then - call cgyro_error('MPI_FILE_OPEN in cgyro_read_restart_one failed') - return - endif - - call MPI_FILE_SET_VIEW(fhv,& - disp,& - MPI_COMPLEX16,& - MPI_COMPLEX16,& - 'native',& - finfo,& - i_err) - - call MPI_FILE_READ_AT(fhv,& - offset1,& - h_x,& - size(h_x),& - MPI_COMPLEX16,& - fstatus,& - i_err) - - if (i_err /= 0) then - call cgyro_error('MPI_FILE_READ_AT in cgyro_read_restart_one failed') - return - endif - - call MPI_FILE_CLOSE(fhv,i_err) - call MPI_INFO_FREE(finfo,i_err) - - ! Unpack h(0,0) into source - if (source_flag == 1 .and. nt1 == 0) then - ic0 = (n_radial/2)*n_theta - do j=1,n_theta - source(j,:,0) = h_x(ic0+j,:,0) - h_x(ic0+j,:,0) = 0.0 - enddo - sa = 0.0 - do j=1,nint(t_current/delta_t) - sa = 1.0+exp(-delta_t/tau_ave)*sa - enddo - endif - - call system_clock(cp_time,count_rate,count_max) - if (cp_time > start_time) then - cp_dt = (cp_time-start_time)/real(count_rate) - else - cp_dt = (cp_time-start_time+count_max)/real(count_rate) - endif - - if (i_proc == 0) then - call date_and_time(sdate,stime); - call get_environment_variable('GACODE_PLATFORM',platform) - open(NEWUNIT=statusfd,FILE=trim(path)//runfile_startups,action='write',status='unknown',position='append') - write(statusfd,'(14(a),f7.3)') & - sdate(1:4),'/',sdate(5:6),'/',sdate(7:8),' ', & - stime(1:2),':',stime(3:4),':',stime(5:6),' ', & - trim(platform),' [ CHECKPOINT READ] Time =',cp_dt - close(statusfd) - endif - -end subroutine cgyro_read_restart_one - - diff --git a/cgyro/src/cgyro_write_restart.F90 b/cgyro/src/cgyro_restart.F90 similarity index 56% rename from cgyro/src/cgyro_write_restart.F90 rename to cgyro/src/cgyro_restart.F90 index b9852235c..abf315a5f 100644 --- a/cgyro/src/cgyro_write_restart.F90 +++ b/cgyro/src/cgyro_restart.F90 @@ -1,11 +1,17 @@ !------------------------------------------------ -! cgyro_write_restart.f90 +! cgyro_restart.f90 ! ! PURPOSE: -! This is the master file controlling output of +! This is the master file controlling input and output of ! restart data using MPI-IO. !------------------------------------------------ +module cgyro_restart + + implicit none + +contains + subroutine cgyro_write_restart use mpi @@ -281,3 +287,237 @@ subroutine cgyro_write_restart_header_part close(io) end subroutine cgyro_write_restart_header_part +!------------------------------------------------------ +! cgyro_read_restart.f90 +! +! PURPOSE: +! This is the master file controlling the restart +! via MPI-IO. +!------------------------------------------------------ + +subroutine cgyro_read_restart + + use mpi + use cgyro_globals + use cgyro_io + + implicit none + + !--------------------------------------------------------- + ! Read restart parameters from ASCII file. + ! + if (restart_flag == 1) then + 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 + 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) + + endif + + if (i_proc == 0) then + open(unit=io,file=trim(path)//runfile_restart,status='old',iostat=i_err) + close(io) + endif + + call cgyro_read_restart_one + +end subroutine cgyro_read_restart + + +subroutine cgyro_read_restart_verify + + use cgyro_globals + use cgyro_io + + !--------------------------------------------------- + implicit none + + integer :: magic, version, recid + integer :: t_n_theta,t_n_radial,t_n_species,t_n_xi,t_n_energy,t_n_toroidal + + ! Different compilers have a different semantics for RECL... must use inquire + integer :: reclen + integer, dimension(1) :: recltest + + inquire(iolength=reclen) recltest + + open(unit=io,& + file=trim(path)//runfile_restart,& + status='old',access='DIRECT',RECL=reclen,ACTION='READ') + + recid = 1 + + read(io,REC=recid) magic + recid = recid + 1 + if ( magic /= restart_magic) then + close(io) + call cgyro_error('Wrong magic number in restart header') + return + endif + + read(io,REC=recid) version + recid = recid + 1 + if ( version /= 2) then + close(io) + call cgyro_error('Wrong version in restart header, only v2 supported') + return + endif + + read(io,REC=recid) t_n_theta + recid = recid + 1 + read(io,REC=recid) t_n_radial + recid = recid + 1 + read(io,REC=recid) t_n_species + recid = recid + 1 + read(io,REC=recid) t_n_xi + recid = recid + 1 + read(io,REC=recid) t_n_energy + recid = recid + 1 + read(io,REC=recid) t_n_toroidal + recid = recid + 1 + if ( (t_n_theta/=n_theta) .or. (t_n_radial/=n_radial) .or. & + (t_n_species/=n_species) .or. (t_n_xi/=n_xi) .or. & + (t_n_energy/=n_energy) .or. (t_n_toroidal/=n_toroidal) ) then + close(io) + call cgyro_error('Wrong geometry in restart header') + return + endif + + ! follow MPI params... will ignore them for v2, as they do not change the file format + + close(io) + +end subroutine cgyro_read_restart_verify + +subroutine cgyro_read_restart_one + + use mpi + use cgyro_globals + use cgyro_io + + !--------------------------------------------------- + implicit none + ! + ! Required for MPI-IO: + ! + integer :: filemode + integer :: finfo + integer :: fhv + integer :: fstatus(MPI_STATUS_SIZE) + integer(kind=MPI_OFFSET_KIND) :: disp + integer(kind=MPI_OFFSET_KIND) :: offset1 + !--------------------------------------------------- + + character(8) :: sdate + character(10) :: stime + character(len=64) :: platform + integer(KIND=8) :: start_time,cp_time + integer(KIND=8) :: count_rate, count_max + real :: cp_dt + integer :: ic0,j,statusfd + + ! use system_clock to be consistent with cgyro_kernel + call system_clock(start_time,count_rate,count_max) + + ! First read the header, and verify that it is compatible with current setup + + if (i_proc == 0) then + call cgyro_read_restart_verify + if (error_status /=0 ) return + endif + + filemode = MPI_MODE_RDONLY + disp = 0 + + offset1 = size(h_x,kind=MPI_OFFSET_KIND)*(i_proc_1+i_proc_2*n_proc_1) + restart_header_size + if (offset1 < restart_header_size) then + call cgyro_error('Overflow detected in cgyro_read_restart_one') + return + endif + + call MPI_INFO_CREATE(finfo,i_err) + + call MPI_FILE_OPEN(CGYRO_COMM_WORLD,& + trim(path)//runfile_restart,& + filemode,& + finfo,& + fhv,& + i_err) + + if (i_err /= 0) then + call cgyro_error('MPI_FILE_OPEN in cgyro_read_restart_one failed') + return + endif + + call MPI_FILE_SET_VIEW(fhv,& + disp,& + MPI_COMPLEX16,& + MPI_COMPLEX16,& + 'native',& + finfo,& + i_err) + + call MPI_FILE_READ_AT(fhv,& + offset1,& + h_x,& + size(h_x),& + MPI_COMPLEX16,& + fstatus,& + i_err) + + if (i_err /= 0) then + call cgyro_error('MPI_FILE_READ_AT in cgyro_read_restart_one failed') + return + endif + + call MPI_FILE_CLOSE(fhv,i_err) + call MPI_INFO_FREE(finfo,i_err) + + ! Unpack h(0,0) into source + if (source_flag == 1 .and. nt1 == 0) then + ic0 = (n_radial/2)*n_theta + do j=1,n_theta + source(j,:,0) = h_x(ic0+j,:,0) + h_x(ic0+j,:,0) = 0.0 + enddo + sa = 0.0 + do j=1,nint(t_current/delta_t) + sa = 1.0+exp(-delta_t/tau_ave)*sa + enddo + endif + + call system_clock(cp_time,count_rate,count_max) + if (cp_time > start_time) then + cp_dt = (cp_time-start_time)/real(count_rate) + else + cp_dt = (cp_time-start_time+count_max)/real(count_rate) + endif + + if (i_proc == 0) then + call date_and_time(sdate,stime); + call get_environment_variable('GACODE_PLATFORM',platform) + open(NEWUNIT=statusfd,FILE=trim(path)//runfile_startups,action='write',status='unknown',position='append') + write(statusfd,'(14(a),f7.3)') & + sdate(1:4),'/',sdate(5:6),'/',sdate(7:8),' ', & + stime(1:2),':',stime(3:4),':',stime(5:6),' ', & + trim(platform),' [ CHECKPOINT READ] Time =',cp_dt + close(statusfd) + endif + +end subroutine cgyro_read_restart_one + +end module cgyro_restart + From bc2b907d25827ed2a421ccc9d74222e3f1ed84be Mon Sep 17 00:00:00 2001 From: Igor Sfiligoi Date: Wed, 15 May 2024 09:36:28 -0700 Subject: [PATCH 03/16] Fix make order --- cgyro/src/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cgyro/src/Makefile b/cgyro/src/Makefile index dea9cdb44..e93f92091 100644 --- a/cgyro/src/Makefile +++ b/cgyro/src/Makefile @@ -19,9 +19,9 @@ OBJECTS = cgyro_globals.o \ cgyro_math.o \ cgyro_globals_math.o \ cgyro_parallel_lib.o \ - cgyro_restart.o \ cgyro_step.o \ cgyro_io.o \ + cgyro_restart.o \ cgyro_advect_wavenumber.o \ cgyro_init_rotation.o \ cgyro_equilibrium.o \ From 4454b22056cf24b05874572d8ee80f926f8e0080 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Thu, 16 May 2024 11:22:23 -0700 Subject: [PATCH 04/16] Working on supertorus improvements --- cgyro/bin/cgyro_plot | 2 +- cgyro/install/make.ext.MINT_PGI | 2 - cgyro/install/make.ext.TITAN_PGI | 2 - f2py/pygacode/cgyro/vis_supertorus.py | 49 +++--- f2py/vis/vis.f90 | 241 ++++++++++---------------- 5 files changed, 121 insertions(+), 175 deletions(-) delete mode 100644 cgyro/install/make.ext.MINT_PGI delete mode 100644 cgyro/install/make.ext.TITAN_PGI diff --git a/cgyro/bin/cgyro_plot b/cgyro/bin/cgyro_plot index 53950fdf8..97cd5d995 100755 --- a/cgyro/bin/cgyro_plot +++ b/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 ;; diff --git a/cgyro/install/make.ext.MINT_PGI b/cgyro/install/make.ext.MINT_PGI deleted file mode 100644 index f0747935a..000000000 --- a/cgyro/install/make.ext.MINT_PGI +++ /dev/null @@ -1,2 +0,0 @@ -cgyro_nl_fftw.o : cgyro_nl_fftw.gpu.F90 - $(FC) $(FMATH) $(FFLAGS) -o cgyro_nl_fftw.o -c cgyro_nl_fftw.gpu.F90 diff --git a/cgyro/install/make.ext.TITAN_PGI b/cgyro/install/make.ext.TITAN_PGI deleted file mode 100644 index f0747935a..000000000 --- a/cgyro/install/make.ext.TITAN_PGI +++ /dev/null @@ -1,2 +0,0 @@ -cgyro_nl_fftw.o : cgyro_nl_fftw.gpu.F90 - $(FC) $(FMATH) $(FFLAGS) -o cgyro_nl_fftw.o -c cgyro_nl_fftw.gpu.F90 diff --git a/f2py/pygacode/cgyro/vis_supertorus.py b/f2py/pygacode/cgyro/vis_supertorus.py index a61a22777..ad2e24116 100644 --- a/f2py/pygacode/cgyro/vis_supertorus.py +++ b/f2py/pygacode/cgyro/vis_supertorus.py @@ -1,19 +1,19 @@ import struct import sys +import time import numpy as np import os -import matplotlib.pyplot as plt -from matplotlib import rc from matplotlib import cm from ..gacodefuncs import * from .data import cgyrodata from mayavi import mlab + try: from pygacode import geo from vis.vis import vis - print("INFO: (vis_supertorus) Successfully imported vis") + print('INFO: (vis_supertorus) Successfully imported vis') except: - print("ERROR: (vis_supertorus) Please type 'pip install pygacode'") + print('ERROR: (vis_supertorus) Please type make in f2py/vis') sys.exit() ext = sys.argv[1] @@ -33,6 +33,8 @@ mag = float(sys.argv[15]) nozonal = bool(int(sys.argv[16])) onlyzonal = bool(int(sys.argv[17])) +px = int(sys.argv[18]) +py = int(sys.argv[19]) sim = cgyrodata('./') nt = sim.n_time @@ -41,11 +43,14 @@ ns = sim.n_species nth = sim.theta_plot +lovera = sim.length*sim.rho/dn*mag + print('HINT: adjust -dn to match experimental dn (rho/a and Lx/a will shrink)') print('Lx/rho = {:.2f}'.format(sim.length)) print('rho/a = {:.4f}'.format(sim.rho/dn)) -lovera = sim.length*sim.rho/dn*mag print('Lx/a = {:.4f}'.format(lovera)) +print('px = {} py = {}'.format(px,py)) +print('nx = {} nz = {} nphi = {}'.format(nx,nz,nphi)) ivec = time_vector(istr,nt) @@ -170,8 +175,9 @@ def frame(): a = np.reshape(aa,(2,nr,nth,nn),order='F') else: a = np.reshape(aa,(2,nr,nth,ns,nn),order='F') - - mlab.figure(size=(900,900),bgcolor=(1,1,1)) + + #mlab.figure(size=(px,py),bgcolor=(1,1,1)) + mlab.figure(bgcolor=(1,1,1)) if isfield: c = a[0,:,:,:]+1j*a[1,:,:,:] else: @@ -214,30 +220,25 @@ def subframe_tp(): f = np.zeros([2,nz],order='F') for j in range(nphi): - # - for k in range(nz): - phi = -j*dphi - g1[k] = g10[k] - phi - xpp[0,k,j] = xp0[0,k] * np.cos(phi) - zpp[0,k,j] = xp0[0,k] * np.sin(phi) - ypp[0,k,j] = yp[0,k] - xpp[1,k,j] = xp0[nx-1,k] * np.cos(phi) - zpp[1,k,j] = xp0[nx-1,k] * np.sin(phi) - ypp[1,k,j] = yp[nx-1,k] + phi = -j*dphi + g1[:] = g10[:] - phi + xpp[0,:,j] = xp0[0,:] * np.cos(phi) + zpp[0,:,j] = xp0[0,:] * np.sin(phi) + ypp[0,:,j] = yp[0,:] + xpp[1,:,j] = xp0[nx-1,:] * np.cos(phi) + zpp[1,:,j] = xp0[nx-1,:] * np.sin(phi) + ypp[1,:,j] = yp[nx-1,:] # f[:,:] = 0.0 vis.torcut(dn,sim.m_box,sim.q,sim.thetap,g1,g2,c,f) - for k in range(nz): - f_all[0,k,j] = f[0,k] - f_all[1,k,j] = f[1,k] + f_all[:,:,j] = f[:,:] for i in range(2): image = mlab.mesh(xpp[i,:,:],ypp[i,:,:],zpp[i,:,:],scalars=f_all[i,:,:],colormap=colormap,vmin=fmin_all,vmax=fmax_all,opacity=1.0) PREC='f' ; BIT=4 -print("ivec") -print(ivec) +print('ivec {}'.format(ivec)) # Open binary file with open(fdata,'rb') as fbin: @@ -273,12 +274,12 @@ def subframe_tp(): subframe_tp() ################## # View from positive z-axis - mlab.view(azimuth=0, elevation=0) + mlab.view(azimuth=1, elevation=1) if ftype == 'screen': mlab.show() else: # Filename uses frame number - mlab.savefig(pre+str(i)+'.'+ftype) + mlab.savefig(pre+str(i)+'.'+ftype,size=(px,py)) # Close each time to prevent memory accumulation mlab.close() diff --git a/f2py/vis/vis.f90 b/f2py/vis/vis.f90 index 0fabbadef..3537798cd 100644 --- a/f2py/vis/vis.f90 +++ b/f2py/vis/vis.f90 @@ -1,9 +1,8 @@ module vis ! torcut + ! torside ! realfluct - ! wheel1 - ! wheel2 contains @@ -133,256 +132,206 @@ subroutine torcut(dn,m,q,nr,nth,nn,nx,nz,thi,g1,g2,c,f) end subroutine torcut - subroutine realfluct(nr,nn,nx,ny,c,f) + subroutine torside(dn,m,q,nr,nth,nn,nx,nz,thi,g1,g2,c,f) implicit none - integer, intent(in) :: nr,nn,nx,ny - double complex, intent(in) :: c(0:nr-1,0:nn-1) - double precision, intent(inout) :: f(0:nx-1,0:ny-1) + integer, intent(in) :: dn,m + double precision, intent(in) :: q + integer, intent(in) :: nr,nth,nn,nx,nz + double precision, intent(in) :: thi(0:nth-1) + double precision, intent(in) :: g1(0:nz-1),g2(0:nz-1) + double complex, intent(in) :: c(0:nr-1,0:nth-1,0:nn-1) + double precision, intent(inout) :: f(0:nx-1,0:nz-1) - double complex, dimension(:,:), allocatable :: epx,eny + double complex, dimension(:,:), allocatable :: epx,eny,c0 + double complex, dimension(:,:,:), allocatable :: cx + double precision, dimension(:), allocatable :: x,z,th - integer :: i,j,n,p - double precision :: pi,xi,yj,fsum + integer :: i,k,kc + integer :: n,p,pp + double precision :: pi,fsum,x0 double complex :: ic + ! f2py intent(in) dn + ! f2py intent(in) m + ! f2py intent(in) q ! f2py intent(in) nr + ! f2py intent(in) nth ! f2py intent(in) nn ! f2py intent(in) nx - ! f2py intent(in) ny + ! f2py intent(in) nz + ! f2py intent(in) thi + ! f2py intent(in) g1 + ! f2py intent(in) g2 ! f2py intent(in) c ! f2py intent(in,out) f allocate(epx(0:nr-1,0:nx-1)) - allocate(eny(0:nn-1,0:ny-1)) + allocate(eny(0:nn-1,0:nz-1)) + allocate(x(0:nx-1)) + allocate(z(0:nz-1)) + allocate(th(0:nth)) + allocate(c0(0:nr-1,0:nn-1)) + allocate(cx(0:nr-1,0:nth,0:nn-1)) ic = (0d0,1d0) pi = atan(1d0)*4d0 + th(0:nth-1) = thi + th(nth) = -th(0) + do i=0,nx-1 - xi = i*2*pi/(nx-1) + x(i) = i*2*pi/(nx-1)/dn do p=0,nr-1 - epx(p,i) = exp(ic*(p-nr/2)*xi) + epx(p,i) = exp(dn*ic*(p-nr/2)*x(i)) enddo enddo if (nn > 1) then - do j=0,ny-1 - yj = j*2*pi/(ny-1) + do k=0,nz-1 + z(k) = k*2*pi/(nz-1)-pi do n=0,nn-1 - eny(n,j) = exp(-ic*n*yj) + eny(n,k) = exp(dn*ic*n*g1(k)) enddo enddo ! factor of 1/2 for n=0 eny(0,:) = eny(0,:)/2 else - do j=0,ny-1 - yj = j*2*pi/(ny-1) - n=1 - eny(0,j) = exp(-ic*n*yj) + do k=0,nz-1 + z(k) = k*2*pi/(nz-1)-pi + eny(0,k) = exp(dn*ic*g1(k)) enddo endif -!$omp parallel do private(j,i,fsum,n,p) - do j=0,ny-1 - do i=0,nx-1 - fsum = 0d0 - if (nn > 1) then - do n=0,nn-1 - do p=0,nr-1 - fsum = fsum+real(c(p,n)*epx(p,i)*eny(n,j)) - enddo - enddo - else - n=1 - do p=0,nr-1 - fsum = fsum+real(c(p,0)*epx(p,i)*eny(0,j)) - enddo - endif - f(i,j) = fsum - enddo - enddo - - deallocate(epx,eny) - - end subroutine realfluct - - subroutine wheel1(nr,nth,nn,nx,nz,c,f) - - implicit none - - integer, intent(in) :: nr,nth,nn,nx,nz - double complex, intent(in) :: c(0:nr-1,0:nth-1,0:nn-1) - double precision, intent(inout) :: f(0:nx-1,0:nz-1) - - double complex, dimension(:,:), allocatable :: epx,c0 - double precision, dimension(:), allocatable :: eny - double precision, dimension(:), allocatable :: z,th - - integer :: i,k,kc - integer :: n,p - double precision :: pi,fsum,xi - double complex :: ic - - ! f2py intent(in) nr - ! f2py intent(in) nth - ! f2py intent(in) nn - ! f2py intent(in) nx - ! f2py intent(in) nz - ! f2py intent(in) c - ! f2py intent(in,out) f - - allocate(epx(0:nr-1,0:nx-1)) - allocate(eny(0:nn-1)) - allocate(z(0:nz-1)) - allocate(th(0:nth)) - allocate(c0(0:nr-1,0:nn-1)) - - ic = (0d0,1d0) - pi = atan(1d0)*4d0 - - do i=0,nx-1 - xi = i*2*pi/(nx-1) - do p=0,nr-1 - epx(p,i) = exp(ic*(p-nr/2)*xi) + ! Extended coefficients + cx(:,0:nth-1,:) = c(:,:,:) + if (nn > 1) then + do n=0,nn-1 + do p=0,nr-1 + pp = p+n*m + if (pp > nr-1) pp = pp-nr + cx(p,nth,n) = c(pp,0,n)*exp(-2*pi*dn*ic*n*q) + enddo enddo - enddo - - do kc=0,nth - th(kc) = kc*2*pi/nth-pi - enddo - - do k=0,nz-1 - z(k) = k*pi/(nz-1)-pi - enddo - - if (nn > 1) then - ! factor of 1/2 for n=0 - eny(:) = 1d0 - eny(0) = 0.5d0 else - eny(0) = 1d0 + do p=0,nr-1 + pp = p+1*m + if (pp > nr-1) pp = pp-nr + cx(p,nth,0) = c(pp,0,0)*exp(-2*pi*dn*ic*q) + enddo endif -!$omp parallel do private(k,kc,c0,i,n,p,fsum) +!$omp parallel do private(k,kc,c0,i,fsum,x0,n,p) do k=0,nz-1 + do kc=0,nth-1 + ! Interpolation of ballooning potential (on coarse, irregular th-grid) + ! onto fine, equally-spaced z-grid if (z(k) >= th(kc) .and. z(k) <= th(kc+1)) then - c0(:,:) = ( c(:,kc+1,:)*(z(k)-th(kc)) & - + c(:,kc,:)*(th(kc+1)-z(k)) )/(th(1)-th(0)) + c0(:,:) = ( cx(:,kc+1,:)*(z(k)-th(kc)) & + + cx(:,kc,:)*(th(kc+1)-z(k)) )/(th(kc+1)-th(kc)) exit endif enddo + do i=0,nx-1 - fsum = 0d0 + fsum = 0.0 + x0 = m*x(i)*g2(k)/(2*pi) if (nn > 1) then do n=0,nn-1 do p=0,nr-1 - fsum = fsum+real(c0(p,n)*epx(p,i)*eny(n)) + fsum = fsum+real( c0(p,n)*epx(p,i)*eny(n,k)*exp(dn*ic*n*x0) ) enddo enddo else - n=1 do p=0,nr-1 - fsum = fsum+real(c0(p,0)*epx(p,i)*eny(0)) + fsum = fsum+real( c0(p,0)*epx(p,i)*eny(0,k)*exp(dn*ic*x0) ) enddo endif f(i,k) = fsum enddo + enddo - deallocate(epx,eny,z,th,c0) + deallocate(epx,eny,x,z,th,c0,cx) - end subroutine wheel1 + end subroutine torside - subroutine wheel2(nr,nth,nn,ny,nz,c,f) + subroutine realfluct(nr,nn,nx,ny,c,f) implicit none - integer, intent(in) :: nr,nth,nn,ny,nz - double complex, intent(in) :: c(0:nr-1,0:nth-1,0:nn-1) - double precision, intent(inout) :: f(0:ny-1,0:nz-1) + integer, intent(in) :: nr,nn,nx,ny + double complex, intent(in) :: c(0:nr-1,0:nn-1) + double precision, intent(inout) :: f(0:nx-1,0:ny-1) - double complex, dimension(:,:), allocatable :: eny,c0 - double precision, dimension(:), allocatable :: y,z,th + double complex, dimension(:,:), allocatable :: epx,eny - integer :: j,k,kc - integer :: n,p - double precision :: pi,fsum + integer :: i,j,n,p + double precision :: pi,xi,yj,fsum double complex :: ic ! f2py intent(in) nr - ! f2py intent(in) nth ! f2py intent(in) nn + ! f2py intent(in) nx ! f2py intent(in) ny - ! f2py intent(in) nz ! f2py intent(in) c ! f2py intent(in,out) f + allocate(epx(0:nr-1,0:nx-1)) allocate(eny(0:nn-1,0:ny-1)) - allocate(c0(0:nr-1,0:nn-1)) - allocate(th(0:nth)) - allocate(y(0:ny-1)) - allocate(z(0:nz-1)) ic = (0d0,1d0) pi = atan(1d0)*4d0 - do kc=0,nth - th(kc) = kc*2*pi/nth-pi + do i=0,nx-1 + xi = i*2*pi/(nx-1) + do p=0,nr-1 + epx(p,i) = exp(ic*(p-nr/2)*xi) + enddo enddo if (nn > 1) then do j=0,ny-1 - y(j) = j*2*pi/(ny-1) + yj = j*2*pi/(ny-1) do n=0,nn-1 - eny(n,j) = exp(-ic*n*y(j)) + eny(n,j) = exp(-ic*n*yj) enddo enddo ! factor of 1/2 for n=0 eny(0,:) = eny(0,:)/2 else do j=0,ny-1 - y(j) = j*2*pi/(ny-1) + yj = j*2*pi/(ny-1) n=1 - eny(0,j) = exp(-ic*n*y(j)) + eny(0,j) = exp(-ic*n*yj) enddo endif - do k=0,nz-1 - z(k) = k*pi/(nz-1)-pi - enddo - -!$omp parallel do private(k,kc,c0,j,n,p,fsum) - do k=0,nz-1 - do kc=0,nth-1 - if (z(k) >= th(kc) .and. z(k) <= th(kc+1)) then - c0(:,:) = ( c(:,kc+1,:)*(z(k)-th(kc)) & - + c(:,kc,:)*(th(kc+1)-z(k)) )/(th(1)-th(0)) - exit - endif - enddo - do j=0,ny-1 +!$omp parallel do private(j,i,fsum,n,p) + do j=0,ny-1 + do i=0,nx-1 fsum = 0d0 if (nn > 1) then do n=0,nn-1 do p=0,nr-1 - fsum = fsum+real(c0(p,n)*eny(n,j)) + fsum = fsum+real(c(p,n)*epx(p,i)*eny(n,j)) enddo enddo else n=1 do p=0,nr-1 - fsum = fsum+real(c0(p,0)*eny(0,j)) + fsum = fsum+real(c(p,0)*epx(p,i)*eny(0,j)) enddo endif - f(j,k) = fsum + f(i,j) = fsum enddo enddo - deallocate(eny,c0,th,y,z) + deallocate(epx,eny) - end subroutine wheel2 + end subroutine realfluct + end module vis From 8f08a8393e919386a506718eecaf8d86ab2a45af Mon Sep 17 00:00:00 2001 From: Igor Sfiligoi Date: Thu, 16 May 2024 14:58:18 -0700 Subject: [PATCH 05/16] Introduce a new restart format that records more details. Provide backwards compatibilty --- cgyro/src/cgyro_globals.F90 | 2 - cgyro/src/cgyro_init_manager.F90 | 8 - cgyro/src/cgyro_restart.F90 | 299 ++++++++++++++++++++++++++++--- 3 files changed, 270 insertions(+), 39 deletions(-) diff --git a/cgyro/src/cgyro_globals.F90 b/cgyro/src/cgyro_globals.F90 index c27635973..06cf73f76 100644 --- a/cgyro/src/cgyro_globals.F90 +++ b/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 diff --git a/cgyro/src/cgyro_init_manager.F90 b/cgyro/src/cgyro_init_manager.F90 index 49ac16e4a..2476803e8 100644 --- a/cgyro/src/cgyro_init_manager.F90 +++ b/cgyro/src/cgyro_init_manager.F90 @@ -381,14 +381,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 diff --git a/cgyro/src/cgyro_restart.F90 b/cgyro/src/cgyro_restart.F90 index abf315a5f..014e51cf7 100644 --- a/cgyro/src/cgyro_restart.F90 +++ b/cgyro/src/cgyro_restart.F90 @@ -10,6 +10,14 @@ module cgyro_restart implicit none + integer, parameter :: restart_header_size = 1024 + integer, parameter :: restart_magic = 140974129 + integer, parameter :: restart_version = 3 + + integer, private :: t_velocity_order + integer, private :: t_nt_loc + integer, private :: t_nv_loc + contains subroutine cgyro_write_restart @@ -20,11 +28,26 @@ subroutine cgyro_write_restart implicit none + integer :: j,ic0 + ! Print this data on restart steps only; otherwise exit now if (mod(i_time,restart_step*print_step) /= 0) return call cgyro_write_restart_one + ! Unpack h(0,0) into source + if (source_flag == 1 .and. nt1 == 0) then + ic0 = (n_radial/2)*n_theta + do j=1,n_theta + source(j,:,0) = h_x(ic0+j,:,0) + h_x(ic0+j,:,0) = 0.0 + enddo + sa = 0.0 + do j=1,nint(t_current/delta_t) + sa = 1.0+exp(-delta_t/tau_ave)*sa + enddo + endif + ! Write restart tag if (i_proc == 0) then open(unit=io,file=trim(path)//runfile_restart_tag,status='replace') @@ -249,7 +272,6 @@ subroutine cgyro_write_restart_header_part !--------------------------------------------------- implicit none - integer, parameter :: version = 2 integer :: recid ! Different compilers have a different semantics for RECL... must use inquire @@ -265,7 +287,7 @@ subroutine cgyro_write_restart_header_part recid = 1 write(io,REC=recid) restart_magic recid = recid + 1 - write(io,REC=recid) version + write(io,REC=recid) restart_version recid = recid + 1 write(io,REC=recid) n_theta recid = recid + 1 @@ -279,10 +301,38 @@ subroutine cgyro_write_restart_header_part recid = recid + 1 write(io,REC=recid) n_toroidal recid = recid + 1 - write(io,REC=recid) mpi_rank_order + write(io,REC=recid) velocity_order + recid = recid + 1 + write(io,REC=recid) nt_loc + recid = recid + 1 + write(io,REC=recid) nv_loc + recid = recid + 1 + write(io,REC=recid) nc_loc + recid = recid + 1 + write(io,REC=recid) n_toroidal_procs recid = recid + 1 write(io,REC=recid) n_proc recid = recid + 1 + write(io,REC=recid) mpi_rank_order + recid = recid + 1 + write(io,REC=recid) delta_t_method + recid = recid + 1 + write(io,REC=recid) nonlinear_flag + recid = recid + 1 + write(io,REC=recid) n_jtheta + recid = recid + 1 + write(io,REC=recid) nsplit + recid = recid + 1 + write(io,REC=recid) nsplitA + recid = recid + 1 + write(io,REC=recid) nsplitB + recid = recid + 1 + write(io,REC=recid) nup_theta + recid = recid + 1 + write(io,REC=recid) nv + recid = recid + 1 + write(io,REC=recid) nc + recid = recid + 1 write(io,REC=recid) restart_magic ! just to have a clean end close(io) @@ -347,6 +397,7 @@ subroutine cgyro_read_restart_verify integer :: magic, version, recid integer :: t_n_theta,t_n_radial,t_n_species,t_n_xi,t_n_energy,t_n_toroidal + integer :: restart_magic_v2 ! Different compilers have a different semantics for RECL... must use inquire integer :: reclen @@ -362,18 +413,36 @@ subroutine cgyro_read_restart_verify read(io,REC=recid) magic recid = recid + 1 - if ( magic /= restart_magic) then - close(io) - call cgyro_error('Wrong magic number in restart header') - return - endif - - read(io,REC=recid) version - recid = recid + 1 - if ( version /= 2) then - close(io) - call cgyro_error('Wrong version in restart header, only v2 supported') - return + if ( magic == restart_magic) then + read(io,REC=recid) version + recid = recid + 1 + if ( version /= restart_version) then + close(io) + call cgyro_error('Wrong version in restart header, expected 3') + return + endif + else + ! older versions had different magic + if (velocity_order == 1) then + ! traditional ordering + restart_magic_v2 = 140906808 + else + ! alternative ordering, needed different magic + restart_magic_v2 = 140916753 + endif + if ( magic /= restart_magic_v2) then + ! we don't support anything else, abort + close(io) + call cgyro_error('Wrong magic number in restart header') + return + endif + read(io,REC=recid) version + recid = recid + 1 + if ( version /= 2) then + close(io) + call cgyro_error('Wrong version in restart header, expected 2') + return + endif endif read(io,REC=recid) t_n_theta @@ -396,7 +465,21 @@ subroutine cgyro_read_restart_verify return endif - ! follow MPI params... will ignore them for v2, as they do not change the file format + if ( magic == restart_magic) then + read(io,REC=recid) t_velocity_order + recid = recid + 1 + read(io,REC=recid) t_nt_loc + recid = recid + 1 + read(io,REC=recid) t_nv_loc + recid = recid + 1 + else + ! not recorded, so assume they are the same as current run + t_velocity_order = velocity_order + t_nt_loc = nt_loc + t_nv_loc = nv_loc + endif + + ! follow other params... will ignore them, as they do not change the file format close(io) @@ -427,6 +510,7 @@ subroutine cgyro_read_restart_one integer(KIND=8) :: start_time,cp_time integer(KIND=8) :: count_rate, count_max real :: cp_dt + integer, dimension(3) :: mpibuf integer :: ic0,j,statusfd ! use system_clock to be consistent with cgyro_kernel @@ -437,6 +521,29 @@ subroutine cgyro_read_restart_one if (i_proc == 0) then call cgyro_read_restart_verify if (error_status /=0 ) return + mpibuf(1) = t_nt_loc + mpibuf(2) = t_nv_loc + mpibuf(3) = t_velocity_order + endif + + call MPI_Bcast(mpibuf, 3, MPI_INTEGER, 0, CGYRO_COMM_WORLD, i_err) + + if (i_err /= 0) then + call cgyro_error('MPI in restart failed') + return + endif + + if (i_proc /= 0) then + t_nt_loc = mpibuf(1) + t_nv_loc = mpibuf(2) + t_velocity_order = mpibuf(3) + endif + + if ( (t_nt_loc /= nt_loc) .or. (t_nv_loc /= nv_loc) .or. (t_velocity_order /= velocity_order) ) then + ! the layout on disk does not align to current process distribution + ! use the more complicated (and slower) logic + call cgyro_read_restart_slow + return endif filemode = MPI_MODE_RDONLY @@ -486,19 +593,6 @@ subroutine cgyro_read_restart_one call MPI_FILE_CLOSE(fhv,i_err) call MPI_INFO_FREE(finfo,i_err) - ! Unpack h(0,0) into source - if (source_flag == 1 .and. nt1 == 0) then - ic0 = (n_radial/2)*n_theta - do j=1,n_theta - source(j,:,0) = h_x(ic0+j,:,0) - h_x(ic0+j,:,0) = 0.0 - enddo - sa = 0.0 - do j=1,nint(t_current/delta_t) - sa = 1.0+exp(-delta_t/tau_ave)*sa - enddo - endif - call system_clock(cp_time,count_rate,count_max) if (cp_time > start_time) then cp_dt = (cp_time-start_time)/real(count_rate) @@ -519,5 +613,152 @@ subroutine cgyro_read_restart_one end subroutine cgyro_read_restart_one +subroutine cgyro_read_restart_slow + + use mpi + use cgyro_globals + use cgyro_io + + !--------------------------------------------------- + implicit none + ! + ! Required for MPI-IO: + ! + integer :: filemode + integer :: finfo + integer :: fhv + integer :: fstatus(MPI_STATUS_SIZE) + integer(kind=MPI_OFFSET_KIND) :: disp + integer(kind=MPI_OFFSET_KIND) :: offset1,off_block,off_row,off_col + !--------------------------------------------------- + + integer, dimension(:,:,:), allocatable :: t_iv_v + character(8) :: sdate + character(10) :: stime + character(len=64) :: platform + integer(KIND=8) :: start_time,cp_time + integer(KIND=8) :: count_rate, count_max + real :: cp_dt + integer :: ic0,j,statusfd + integer :: ie,ix,is,itor + integer :: t_iv + + ! we already checked what the format is + allocate(t_iv_v(n_energy,n_xi,n_species)) + + ! Velocity pointers for the restart velocity_order + t_iv = 0 + if (t_velocity_order==1) then + !original + do ie=1,n_energy + do ix=1,n_xi + do is=1,n_species + t_iv = t_iv+1 + t_iv_v(ie,ix,is) = t_iv + enddo + enddo + enddo + else if (t_velocity_order==2) then + ! optimized for minimizing species + do is=1,n_species + do ie=1,n_energy + do ix=1,n_xi + t_iv = t_iv+1 + t_iv_v(ie,ix,is) = t_iv + enddo + enddo + enddo + else + call cgyro_error('Unknown restart VELOCITY_ORDER.') + return + endif + + + ! use system_clock to be consistent with cgyro_kernel + call system_clock(start_time,count_rate,count_max) + + filemode = MPI_MODE_RDONLY + disp = 0 + + call MPI_INFO_CREATE(finfo,i_err) + + call MPI_FILE_OPEN(CGYRO_COMM_WORLD,& + trim(path)//runfile_restart,& + filemode,& + finfo,& + fhv,& + i_err) + + if (i_err /= 0) then + call cgyro_error('MPI_FILE_OPEN in cgyro_read_restart_one failed') + return + endif + + call MPI_FILE_SET_VIEW(fhv,& + disp,& + MPI_COMPLEX16,& + MPI_COMPLEX16,& + 'native',& + finfo,& + i_err) + + do itor=nt1,nt2 ! itor is 0-based + do iv=nv1,nv2 + t_iv = t_iv_v(ie_v(iv),ix_v(iv),is_v(iv)) + iv_loc = iv-nv1+1 + off_block = modulo(t_iv-1,t_nv_loc) + modulo(itor,t_nt_loc)*t_nv_loc + off_row = ((t_iv-1)/t_nv_loc)*(t_nv_loc*t_nt_loc) + off_col = (itor/t_nt_loc)*(nv*t_nt_loc) + offset1 = size(h_x(:,iv_loc,itor),kind=MPI_OFFSET_KIND)*(off_block+off_row+off_col) + offset1 = offset1 + restart_header_size + if (offset1 < restart_header_size) then + call MPI_FILE_CLOSE(fhv,i_err) + call MPI_INFO_FREE(finfo,i_err) + call cgyro_error('Overflow detected in cgyro_read_restart_one') + return + endif + + call MPI_FILE_READ_AT(fhv,& + offset1,& + h_x(:,iv_loc,itor),& + size(h_x(:,iv_loc,itor)),& + MPI_COMPLEX16,& + fstatus,& + i_err) + + if (i_err /= 0) then + call MPI_FILE_CLOSE(fhv,i_err) + call MPI_INFO_FREE(finfo,i_err) + call cgyro_error('MPI_FILE_READ_AT in cgyro_read_restart_slow failed') + return + endif + enddo + enddo + + call MPI_FILE_CLOSE(fhv,i_err) + call MPI_INFO_FREE(finfo,i_err) + + call system_clock(cp_time,count_rate,count_max) + if (cp_time > start_time) then + cp_dt = (cp_time-start_time)/real(count_rate) + else + cp_dt = (cp_time-start_time+count_max)/real(count_rate) + endif + + if (i_proc == 0) then + call date_and_time(sdate,stime); + call get_environment_variable('GACODE_PLATFORM',platform) + open(NEWUNIT=statusfd,FILE=trim(path)//runfile_startups,action='write',status='unknown',position='append') + write(statusfd,'(14(a),f7.3)') & + sdate(1:4),'/',sdate(5:6),'/',sdate(7:8),' ', & + stime(1:2),':',stime(3:4),':',stime(5:6),' ', & + trim(platform),' [ CHECKPOINT READ (slow)] Time =',cp_dt + close(statusfd) + endif + + deallocate(t_iv_v) + +end subroutine cgyro_read_restart_slow + end module cgyro_restart From 4d93e7fb68c1e6630ef97066aae7681099bba687 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Thu, 16 May 2024 16:40:54 -0700 Subject: [PATCH 06/16] Improved version of vis_supertorus working --- f2py/pygacode/cgyro/vis_supertorus.py | 143 ++++++++++++-------------- f2py/vis/vis.f90 | 50 +++++---- 2 files changed, 95 insertions(+), 98 deletions(-) diff --git a/f2py/pygacode/cgyro/vis_supertorus.py b/f2py/pygacode/cgyro/vis_supertorus.py index ad2e24116..423d66b74 100644 --- a/f2py/pygacode/cgyro/vis_supertorus.py +++ b/f2py/pygacode/cgyro/vis_supertorus.py @@ -24,8 +24,8 @@ nphi = int(sys.argv[6]) phimax = float(sys.argv[7]) istr = sys.argv[8] -fmin = sys.argv[9] -fmax = sys.argv[10] +fmin_in = sys.argv[9] +fmax_in = sys.argv[10] colormap = sys.argv[11] font = int(sys.argv[12]) legacy = bool(int(sys.argv[13])) @@ -65,6 +65,8 @@ #------------------------------------------------------------------------ # (r,theta)=(x,y) mesh setup # +start = time.time() + if nth == 1: print('WARNING: (vis_supertorus) Should use THETA_PLOT > 1 in CGYRO.') @@ -73,31 +75,31 @@ if nz < 0: nz = 128 -x = np.zeros([nx]) -z = np.zeros([nz]) - -for i in range(nx): - x[i] = i*2*np.pi/(nx-1.0)/dn -for k in range(nz): - z[k] = k*2*np.pi/(nz-1.0)-np.pi +x = np.linspace(0,2*np.pi,nx)/dn +z = np.linspace(0,2*np.pi,nz)-np.pi xp = np.zeros([nx,nz]) yp = np.zeros([nx,nz]) zp = np.zeros([nx,nz]) -for i in range(nx): - r = sim.rmin+(dn*x[i]/(2*np.pi)-0.5)*lovera - for k in range(nz): - a = z[k] + (sim.shape_cos[0] - + sim.shape_cos[1]*np.cos(z[k]) - + sim.shape_cos[2]*np.cos(2*z[k]) - + sim.shape_cos[3]*np.cos(3*z[k]) - + np.arcsin(sim.delta)*np.sin(z[k]) - - sim.zeta*np.sin(2*z[k]) - + sim.shape_sin[3]*np.sin(3*z[k])) - xp[i,k] = sim.rmaj+r*np.cos(a) - yp[i,k] = sim.zmag+sim.kappa*r*np.sin(z[k]) - zp[i,k] = 0.0 +# minor radius +r = sim.rmin+(dn*x/(2*np.pi)-0.5)*lovera + +# MXH angle +a = z + (sim.shape_cos[0] + + sim.shape_cos[1]*np.cos(z) + + sim.shape_cos[2]*np.cos(2*z) + + sim.shape_cos[3]*np.cos(3*z) + + np.arcsin(sim.delta)*np.sin(z) + - sim.zeta*np.sin(2*z) + + sim.shape_sin[3]*np.sin(3*z)) + +# MESH +xp[:,:] = sim.rmaj+r[:,None]*np.cos(a[None,:]) +yp[:,:] = sim.zmag+sim.kappa*r[:,None]*np.sin(z[None,:]) +zp[:,:] = 0.0 + +print('MESH TIME = '+'{:.3e}'.format(time.time()-start)+' s.') # Shape functions geo.signb_in=1 # fix @@ -169,6 +171,7 @@ n_chunk = 2*nr*ns*nn*nth def frame(): + global c if isfield: @@ -176,7 +179,6 @@ def frame(): else: a = np.reshape(aa,(2,nr,nth,ns,nn),order='F') - #mlab.figure(size=(px,py),bgcolor=(1,1,1)) mlab.figure(bgcolor=(1,1,1)) if isfield: c = a[0,:,:,:]+1j*a[1,:,:,:] @@ -189,52 +191,53 @@ def frame(): c[:,:,1:] = 0.0 # This is the logic to generate an (r,theta) frame -def subframe_rt(iframe): - global xp,yp,zp,g1,g2,c,fmin_all,fmax_all +def subframe_torcut(iframe): + + global c,fmin,fmax + phi = -2*np.pi*phimax*iframe + xpp = xp*np.cos(phi) + ypp = yp + zpp = xp*np.sin(phi) + f = np.zeros([nx,nz],order='F') - vis.torcut(dn,sim.m_box,sim.q,sim.thetap,g1,g2,c,f) + vis.torcut(dn,sim.m_box,sim.q,sim.thetap,g1-phi/dn,g2,c,f) if iframe == 0: - if fmin == 'auto': - fmin_all=np.min(f) - fmax_all=np.max(f) + if fmin_in == 'auto': + fmin = np.min(f) + fmax = np.max(f) else: - fmin_all=float(fmin) - fmax_all=float(fmax) + fmin = float(fmin_in) + fmax = float(fmax_in) + print('INFO: (vis_supertorus) min={:.3f} | max={:.3f}'.format(fmin,fmax)) - image = mlab.mesh(xp,yp,zp,scalars=f,colormap=colormap,vmin=fmin_all,vmax=fmax_all,opacity=1.0) - print('INFO: (vis_supertorus) min={:.3f} | max={:.3f}'.format(fmin_all,fmax_all)) + image = mlab.mesh(xpp,ypp,zpp,scalars=f, + colormap=colormap,vmin=fmin,vmax=fmax,opacity=1.0) # This is the logic to generate an (theta,phi) frame -def subframe_tp(): - global xp,yp,zp,g1,g2,c,g10,xp0,fmin_all,fmax_all +def subframe_torside(): - dphi = 2*np.pi*phimax/(nphi-1) + global c,fmin,fmax + + phi = -np.linspace(0,2*np.pi*phimax,nphi) + xpp = np.zeros([nz,nphi]) + ypp = np.zeros([nz,nphi]) + zpp = np.zeros([nz,nphi]) - xpp = np.zeros([2,nz,nphi]) - zpp = np.zeros([2,nz,nphi]) - ypp = np.zeros([2,nz,nphi]) - f_all = np.zeros([2,nz,nphi],order='F') - f = np.zeros([2,nz],order='F') + f = np.zeros([2,nz,nphi],order='F') - for j in range(nphi): - phi = -j*dphi - g1[:] = g10[:] - phi - xpp[0,:,j] = xp0[0,:] * np.cos(phi) - zpp[0,:,j] = xp0[0,:] * np.sin(phi) - ypp[0,:,j] = yp[0,:] - xpp[1,:,j] = xp0[nx-1,:] * np.cos(phi) - zpp[1,:,j] = xp0[nx-1,:] * np.sin(phi) - ypp[1,:,j] = yp[nx-1,:] - # - f[:,:] = 0.0 - vis.torcut(dn,sim.m_box,sim.q,sim.thetap,g1,g2,c,f) - f_all[:,:,j] = f[:,:] - - for i in range(2): - image = mlab.mesh(xpp[i,:,:],ypp[i,:,:],zpp[i,:,:],scalars=f_all[i,:,:],colormap=colormap,vmin=fmin_all,vmax=fmax_all,opacity=1.0) + vis.torside(dn,sim.m_box,sim.q,phimax,sim.thetap,g1,g2,c,f) + + for i in [0,1]: + xpp[:,:] = xp[-i,:,None]*np.cos(phi[None,:]) + ypp[:,:] = yp[-i,:,None] + zpp[:,:] = xp[-i,:,None]*np.sin(phi[None,:]) + + image = mlab.mesh(xpp[:,:],ypp[:,:],zpp[:,:],scalars=f[i,:,:], + colormap=colormap,vmin=fmin,vmax=fmax,opacity=1.0) + PREC='f' ; BIT=4 @@ -253,28 +256,16 @@ def subframe_tp(): print('INFO: (vis_supertorus) Time index {:d} '.format(i)) if i in ivec: frame() - g10 = g1 - xp0 = xp - ################## + start = time.time() # cap at phi=0 - phi = 2*np.pi*0.0 - g1 = g10 - phi - xp = xp0 * np.cos(phi) - zp = xp0 * np.sin(phi) - subframe_rt(0) - ################## - # cap at phi=phimax - phi = -2*np.pi*phimax - g1 = g10 - phi - xp = xp0 * np.cos(phi) - zp = xp0 * np.sin(phi) - subframe_rt(1) - ################## + subframe_torcut(0) + # cap at phi=-phimax + subframe_torcut(1) # inner and outer body of torus - subframe_tp() - ################## + subframe_torside() # View from positive z-axis - mlab.view(azimuth=1, elevation=1) + mlab.view(azimuth=0,elevation=0) + print('MAP TIME = '+'{:.3e}'.format(time.time()-start)+' s.') if ftype == 'screen': mlab.show() else: diff --git a/f2py/vis/vis.f90 b/f2py/vis/vis.f90 index 3537798cd..388a82908 100644 --- a/f2py/vis/vis.f90 +++ b/f2py/vis/vis.f90 @@ -132,34 +132,36 @@ subroutine torcut(dn,m,q,nr,nth,nn,nx,nz,thi,g1,g2,c,f) end subroutine torcut - subroutine torside(dn,m,q,nr,nth,nn,nx,nz,thi,g1,g2,c,f) + subroutine torside(dn,m,q,phimax,nr,nth,nn,nphi,nz,thi,g1,g2,c,f) implicit none integer, intent(in) :: dn,m double precision, intent(in) :: q - integer, intent(in) :: nr,nth,nn,nx,nz + double precision, intent(in) :: phimax + integer, intent(in) :: nr,nth,nn,nphi,nz double precision, intent(in) :: thi(0:nth-1) double precision, intent(in) :: g1(0:nz-1),g2(0:nz-1) double complex, intent(in) :: c(0:nr-1,0:nth-1,0:nn-1) - double precision, intent(inout) :: f(0:nx-1,0:nz-1) + double precision, intent(inout) :: f(0:1,0:nz-1,0:nphi-1) - double complex, dimension(:,:), allocatable :: epx,eny,c0 + double complex, dimension(:,:), allocatable :: eny,c0 double complex, dimension(:,:,:), allocatable :: cx - double precision, dimension(:), allocatable :: x,z,th + double precision, dimension(:), allocatable :: phi,z,th integer :: i,k,kc integer :: n,p,pp - double precision :: pi,fsum,x0 + double precision :: pi,fsum1,fsum2,x1,x2,xb1,xb2 double complex :: ic ! f2py intent(in) dn ! f2py intent(in) m ! f2py intent(in) q + ! f2py intent(in) phimax ! f2py intent(in) nr ! f2py intent(in) nth ! f2py intent(in) nn - ! f2py intent(in) nx + ! f2py intent(in) nphi ! f2py intent(in) nz ! f2py intent(in) thi ! f2py intent(in) g1 @@ -167,9 +169,8 @@ subroutine torside(dn,m,q,nr,nth,nn,nx,nz,thi,g1,g2,c,f) ! f2py intent(in) c ! f2py intent(in,out) f - allocate(epx(0:nr-1,0:nx-1)) allocate(eny(0:nn-1,0:nz-1)) - allocate(x(0:nx-1)) + allocate(phi(0:nphi-1)) allocate(z(0:nz-1)) allocate(th(0:nth)) allocate(c0(0:nr-1,0:nn-1)) @@ -181,11 +182,11 @@ subroutine torside(dn,m,q,nr,nth,nn,nx,nz,thi,g1,g2,c,f) th(0:nth-1) = thi th(nth) = -th(0) - do i=0,nx-1 - x(i) = i*2*pi/(nx-1)/dn - do p=0,nr-1 - epx(p,i) = exp(dn*ic*(p-nr/2)*x(i)) - enddo + xb1 = 0.0 + xb2 = 2*pi/dn + + do i=0,nphi-1 + phi(i) = -i*2*pi/(nphi-1)*phimax/dn enddo if (nn > 1) then @@ -222,7 +223,7 @@ subroutine torside(dn,m,q,nr,nth,nn,nx,nz,thi,g1,g2,c,f) enddo endif -!$omp parallel do private(k,kc,c0,i,fsum,x0,n,p) +!$omp parallel do private(k,kc,c0,i,fsum1,fsum2,x1,x2,n,p) do k=0,nz-1 do kc=0,nth-1 @@ -235,26 +236,31 @@ subroutine torside(dn,m,q,nr,nth,nn,nx,nz,thi,g1,g2,c,f) endif enddo - do i=0,nx-1 - fsum = 0.0 - x0 = m*x(i)*g2(k)/(2*pi) + do i=0,nphi-1 + fsum1 = 0.0 + fsum2 = 0.0 + x1 = m*g2(k)*xb1/(2*pi)-phi(i) + x2 = m*g2(k)*xb2/(2*pi)-phi(i) if (nn > 1) then do n=0,nn-1 do p=0,nr-1 - fsum = fsum+real( c0(p,n)*epx(p,i)*eny(n,k)*exp(dn*ic*n*x0) ) + fsum1 = fsum1+real( c0(p,n)*eny(n,k)*exp(dn*ic*n*x1) ) + fsum2 = fsum2+real( c0(p,n)*eny(n,k)*exp(dn*ic*n*x2) ) enddo enddo else do p=0,nr-1 - fsum = fsum+real( c0(p,0)*epx(p,i)*eny(0,k)*exp(dn*ic*x0) ) + fsum1 = fsum1+real( c0(p,0)*eny(0,k)*exp(dn*ic*x1) ) + fsum2 = fsum2+real( c0(p,0)*eny(0,k)*exp(dn*ic*x2) ) enddo endif - f(i,k) = fsum + f(0,k,i) = fsum1 + f(1,k,i) = fsum2 enddo enddo - deallocate(epx,eny,x,z,th,c0,cx) + deallocate(eny,phi,z,th,c0,cx) end subroutine torside From 4aeb5b40b2e6b8e9a42060423406d0145b885359 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Thu, 16 May 2024 19:11:23 -0700 Subject: [PATCH 07/16] Added helper mesh function for torcut and supertorus --- f2py/pygacode/cgyro/vis_mesh.py | 80 ++++++++++++++++++++++ f2py/pygacode/cgyro/vis_supertorus.py | 99 ++++++++------------------- f2py/pygacode/cgyro/vis_torcut.py | 84 +++-------------------- 3 files changed, 117 insertions(+), 146 deletions(-) create mode 100644 f2py/pygacode/cgyro/vis_mesh.py diff --git a/f2py/pygacode/cgyro/vis_mesh.py b/f2py/pygacode/cgyro/vis_mesh.py new file mode 100644 index 000000000..ce5ebda73 --- /dev/null +++ b/f2py/pygacode/cgyro/vis_mesh.py @@ -0,0 +1,80 @@ +import numpy as np + +def vis_mesh(sim,nx,nz,dn,lovera): + + x = np.linspace(0,2*np.pi/dn,nx) # r + z = np.linspace(0,2*np.pi,nz)-np.pi # theta + + xp = np.zeros([nx,nz]) # R + yp = np.zeros([nx,nz]) # Z + zp = np.zeros([nx,nz]) # phi=0 plane + + # minor radius + r = sim.rmin+(dn*x/(2*np.pi)-0.5)*lovera + + # MXH angle + a = z + sim.shape_cos[0] + + sim.shape_cos[1]*np.cos(z) + + sim.shape_cos[2]*np.cos(2*z) + + sim.shape_cos[3]*np.cos(3*z) + + sim.shape_cos[4]*np.cos(4*z) + + sim.shape_cos[5]*np.cos(5*z) + + sim.shape_cos[6]*np.cos(6*z) + + np.arcsin(sim.delta)*np.sin(z) + - sim.zeta*np.sin(2*z) + + sim.shape_sin[3]*np.sin(3*z) + + sim.shape_sin[4]*np.sin(4*z) + + sim.shape_sin[5]*np.sin(5*z) + + sim.shape_sin[6]*np.sin(6*z) + + # MESH + xp[:,:] = sim.rmaj+r[:,None]*np.cos(a[None,:]) + yp[:,:] = sim.zmag+sim.kappa*r[:,None]*np.sin(z[None,:]) + + return x,z,xp,yp,zp + +def vis_geo(sim,geo): + + geo.signb_in=1 # fix + geo.geo_rmin_in=sim.rmin + geo.geo_rmaj_in=sim.rmaj + geo.geo_drmaj_in=sim.shift + geo.geo_zmag_in=sim.zmag + geo.geo_dzmag_in=sim.dzmag + geo.geo_q_in=sim.q + geo.geo_s_in=sim.shear + geo.geo_kappa_in=sim.kappa + geo.geo_s_kappa_in=sim.s_kappa + geo.geo_beta_star_in=sim.beta_star + # Antisymmetric + geo.geo_shape_cos0_in = sim.shape_cos[0] + geo.geo_shape_cos1_in = sim.shape_cos[1] + geo.geo_shape_cos2_in = sim.shape_cos[2] + geo.geo_shape_cos3_in = sim.shape_cos[3] + geo.geo_shape_cos4_in = sim.shape_cos[4] + geo.geo_shape_cos5_in = sim.shape_cos[5] + geo.geo_shape_cos6_in = sim.shape_cos[6] + # Symmetric + geo.geo_delta_in=sim.delta + geo.geo_zeta_in=sim.zeta + geo.geo_shape_sin3_in = sim.shape_sin[3] + geo.geo_shape_sin4_in = sim.shape_sin[4] + geo.geo_shape_sin5_in = sim.shape_sin[5] + geo.geo_shape_sin6_in = sim.shape_sin[6] + # Derivative of antisymmetric + geo.geo_shape_s_cos0_in = sim.shape_s_cos[0] + geo.geo_shape_s_cos1_in = sim.shape_s_cos[1] + geo.geo_shape_s_cos2_in = sim.shape_s_cos[2] + geo.geo_shape_s_cos3_in = sim.shape_s_cos[3] + geo.geo_shape_s_cos4_in = sim.shape_s_cos[4] + geo.geo_shape_s_cos5_in = sim.shape_s_cos[5] + geo.geo_shape_s_cos6_in = sim.shape_s_cos[6] + # Derivative of symmetric + geo.geo_s_delta_in=sim.s_delta + geo.geo_s_zeta_in=sim.s_zeta + geo.geo_shape_s_sin3_in = sim.shape_s_sin[3] + geo.geo_shape_s_sin4_in = sim.shape_s_sin[4] + geo.geo_shape_s_sin5_in = sim.shape_s_sin[5] + geo.geo_shape_s_sin6_in = sim.shape_s_sin[6] + + return geo diff --git a/f2py/pygacode/cgyro/vis_supertorus.py b/f2py/pygacode/cgyro/vis_supertorus.py index 423d66b74..17b247674 100644 --- a/f2py/pygacode/cgyro/vis_supertorus.py +++ b/f2py/pygacode/cgyro/vis_supertorus.py @@ -6,15 +6,7 @@ from matplotlib import cm from ..gacodefuncs import * from .data import cgyrodata -from mayavi import mlab - -try: - from pygacode import geo - from vis.vis import vis - print('INFO: (vis_supertorus) Successfully imported vis') -except: - print('ERROR: (vis_supertorus) Please type make in f2py/vis') - sys.exit() +from .vis_mesh import * ext = sys.argv[1] moment = sys.argv[2] @@ -36,6 +28,27 @@ px = int(sys.argv[18]) py = int(sys.argv[19]) +try: + from vis.vis import vis + print("INFO: (vis_torcut) Successfully imported vis") +except: + print("ERROR: (vis_torcut) Please 'make vismod' in gacode/f2py") + sys.exit() + +try: + from pygacode import geo + print("INFO: (vis_torcut) Successfully imported geo") +except: + print("ERROR: (vis_torcut) Please 'pip install pygacode'") + sys.exit() + +try: + from mayavi import mlab + print("INFO: (vis_torcut) Successfully imported mayavi") +except: + print("ERROR: (vis_torcut) Mayavi module not found.") + sys.exit() + sim = cgyrodata('./') nt = sim.n_time nr = sim.n_radial @@ -63,7 +76,7 @@ ftype = s[0] #------------------------------------------------------------------------ -# (r,theta)=(x,y) mesh setup +# Mesh setup # start = time.time() @@ -75,72 +88,14 @@ if nz < 0: nz = 128 -x = np.linspace(0,2*np.pi,nx)/dn -z = np.linspace(0,2*np.pi,nz)-np.pi - -xp = np.zeros([nx,nz]) -yp = np.zeros([nx,nz]) -zp = np.zeros([nx,nz]) - -# minor radius -r = sim.rmin+(dn*x/(2*np.pi)-0.5)*lovera - -# MXH angle -a = z + (sim.shape_cos[0] - + sim.shape_cos[1]*np.cos(z) - + sim.shape_cos[2]*np.cos(2*z) - + sim.shape_cos[3]*np.cos(3*z) - + np.arcsin(sim.delta)*np.sin(z) - - sim.zeta*np.sin(2*z) - + sim.shape_sin[3]*np.sin(3*z)) - -# MESH -xp[:,:] = sim.rmaj+r[:,None]*np.cos(a[None,:]) -yp[:,:] = sim.zmag+sim.kappa*r[:,None]*np.sin(z[None,:]) -zp[:,:] = 0.0 +x,z,xp,yp,zp = vis_mesh(sim,nx,nz,dn,lovera) print('MESH TIME = '+'{:.3e}'.format(time.time()-start)+' s.') -# Shape functions -geo.signb_in=1 # fix -geo.geo_rmin_in=sim.rmin -geo.geo_rmaj_in=sim.rmaj -geo.geo_drmaj_in=sim.shift -geo.geo_zmag_in=sim.zmag -geo.geo_dzmag_in=sim.dzmag -geo.geo_q_in=sim.q -geo.geo_s_in=sim.shear -geo.geo_kappa_in=sim.kappa -geo.geo_delta_in=sim.delta -geo.geo_zeta_in=sim.zeta -geo.geo_s_kappa_in=sim.s_kappa -geo.geo_s_delta_in=sim.s_delta -geo.geo_s_zeta_in=sim.s_zeta -geo.geo_beta_star_in=sim.beta_star -geo.geo_shape_cos0_in = sim.shape_cos[0] -geo.geo_shape_cos1_in = sim.shape_cos[1] -geo.geo_shape_cos2_in = sim.shape_cos[2] -geo.geo_shape_cos3_in = sim.shape_cos[3] -geo.geo_shape_cos4_in = 0.0 -geo.geo_shape_cos5_in = 0.0 -geo.geo_shape_cos6_in = 0.0 -geo.geo_shape_sin3_in = sim.shape_sin[3] -geo.geo_shape_sin4_in = 0.0 -geo.geo_shape_sin5_in = 0.0 -geo.geo_shape_sin6_in = 0.0 -geo.geo_shape_s_cos0_in = sim.shape_s_cos[0] -geo.geo_shape_s_cos1_in = sim.shape_s_cos[1] -geo.geo_shape_s_cos2_in = sim.shape_s_cos[2] -geo.geo_shape_s_cos3_in = sim.shape_s_cos[3] -geo.geo_shape_s_cos4_in = 0.0 -geo.geo_shape_s_cos5_in = 0.0 -geo.geo_shape_s_cos6_in = 0.0 -geo.geo_shape_s_sin3_in = sim.shape_s_sin[3] -geo.geo_shape_s_sin4_in = 0.0 -geo.geo_shape_s_sin5_in = 0.0 -geo.geo_shape_s_sin6_in = 0.0 - +# Define shape functions for GEO evaluation +geo = vis_geo(sim,geo) geo.geo_interp(z,True) + if legacy: # s-alpha approximate (apparently used in legacy GYRO movies) # g1 -> q*theta diff --git a/f2py/pygacode/cgyro/vis_torcut.py b/f2py/pygacode/cgyro/vis_torcut.py index b5434eafc..df171dffa 100644 --- a/f2py/pygacode/cgyro/vis_torcut.py +++ b/f2py/pygacode/cgyro/vis_torcut.py @@ -3,10 +3,10 @@ import numpy as np import os import matplotlib.pyplot as plt -from matplotlib import rc -from matplotlib import cm +from matplotlib import cm,rc from ..gacodefuncs import * from .data import cgyrodata +from .vis_mesh import * ext = sys.argv[1] moment = sys.argv[2] @@ -55,11 +55,12 @@ ns = sim.n_species nth = sim.theta_plot +lovera = sim.length*sim.rho/dn*mag + print(f'INFO: (vis_torcut) dn = {dn} (cf. dn in out.cgyro.info)') print('HINT: adjust -dn to match experimental dn (rho/a and Lx/a will shrink)') print('Lx/rho = {:.2f}'.format(sim.length)) print('rho/a = {:.4f}'.format(sim.rho/dn)) -lovera = sim.length*sim.rho/dn*mag print('Lx/a = {:.4f}'.format(lovera)) # To calculate dn we could use the following formula given rhosunit, @@ -80,87 +81,22 @@ #------------------------------------------------------------------------ # Mesh setup # +start = time.time() + if nth == 1: print('WARNING: (vis_torcut) Should use THETA_PLOT > 1 in CGYRO.') - if nx < 0: nx = 128 if nz < 0: nz = 128 -x = np.linspace(0,2*np.pi/dn,nx) # r -z = np.linspace(0,2*np.pi,nz)-np.pi # theta +x,z,xp,yp,zp = vis_mesh(sim,nx,nz,dn,lovera) -xp = np.zeros([nx,nz]) # R -yp = np.zeros([nx,nz]) # Z -zp = np.zeros([nx,nz]) # phi=0 plane -a = np.zeros([nz]) - -for i in range(nx): - r = sim.rmin+(dn*x[i]/(2*np.pi)-0.5)*lovera - a[:] = z[:] + sim.shape_cos[0] - + sim.shape_cos[1]*np.cos(1*z[:]) - + sim.shape_cos[2]*np.cos(2*z[:]) - + sim.shape_cos[3]*np.cos(3*z[:]) - + sim.shape_cos[4]*np.cos(4*z[:]) - + sim.shape_cos[5]*np.cos(5*z[:]) - + sim.shape_cos[6]*np.cos(6*z[:]) - + np.arcsin(sim.delta)*np.sin(z[:]) - - sim.zeta*np.sin(2*z[:]) - + sim.shape_sin[3]*np.sin(3*z[:]) - + sim.shape_sin[4]*np.sin(4*z[:]) - + sim.shape_sin[5]*np.sin(5*z[:]) - + sim.shape_sin[6]*np.sin(6*z[:]) - - xp[i,:] = sim.rmaj+r*np.cos(a[:]) - yp[i,:] = sim.zmag+sim.kappa*r*np.sin(z[:]) - zp[i,:] = 0.0 - - +print('MESH TIME = '+'{:.3e}'.format(time.time()-start)+' s.') + # Define shape functions for GEO evaluation -geo.signb_in=1 # fix -geo.geo_rmin_in=sim.rmin -geo.geo_rmaj_in=sim.rmaj -geo.geo_drmaj_in=sim.shift -geo.geo_zmag_in=sim.zmag -geo.geo_dzmag_in=sim.dzmag -geo.geo_q_in=sim.q -geo.geo_s_in=sim.shear -geo.geo_kappa_in=sim.kappa -geo.geo_s_kappa_in=sim.s_kappa -geo.geo_beta_star_in=sim.beta_star -# Antisymmetric -geo.geo_shape_cos0_in = sim.shape_cos[0] -geo.geo_shape_cos1_in = sim.shape_cos[1] -geo.geo_shape_cos2_in = sim.shape_cos[2] -geo.geo_shape_cos3_in = sim.shape_cos[3] -geo.geo_shape_cos4_in = sim.shape_cos[4] -geo.geo_shape_cos5_in = sim.shape_cos[5] -geo.geo_shape_cos6_in = sim.shape_cos[6] -# Symmetric -geo.geo_delta_in=sim.delta -geo.geo_zeta_in=sim.zeta -geo.geo_shape_sin3_in = sim.shape_sin[3] -geo.geo_shape_sin4_in = sim.shape_sin[4] -geo.geo_shape_sin5_in = sim.shape_sin[5] -geo.geo_shape_sin6_in = sim.shape_sin[6] -# Derivative of antisymmetric -geo.geo_shape_s_cos0_in = sim.shape_s_cos[0] -geo.geo_shape_s_cos1_in = sim.shape_s_cos[1] -geo.geo_shape_s_cos2_in = sim.shape_s_cos[2] -geo.geo_shape_s_cos3_in = sim.shape_s_cos[3] -geo.geo_shape_s_cos4_in = sim.shape_s_cos[4] -geo.geo_shape_s_cos5_in = sim.shape_s_cos[5] -geo.geo_shape_s_cos6_in = sim.shape_s_cos[6] -# Derivative of symmetric -geo.geo_s_delta_in=sim.s_delta -geo.geo_s_zeta_in=sim.s_zeta -geo.geo_shape_s_sin3_in = sim.shape_s_sin[3] -geo.geo_shape_s_sin4_in = sim.shape_s_sin[4] -geo.geo_shape_s_sin5_in = sim.shape_s_sin[5] -geo.geo_shape_s_sin6_in = sim.shape_s_sin[6] - +geo = vis_geo(sim,geo) geo.geo_interp(z,True) if legacy: From cd7dd52c68086caaeecbaed736460bbb248ba798 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Fri, 17 May 2024 14:09:06 -0700 Subject: [PATCH 08/16] Added new error trap for underflow --- cgyro/src/cgyro_freq.f90 | 13 ++++++++++--- cgyro/src/cgyro_init_kernel.F90 | 10 ++++------ cgyro/src/cgyro_init_manager.F90 | 9 ++++----- cgyro/src/cgyro_kernel.F90 | 18 +++++++++--------- cgyro/src/cgyro_restart.F90 | 2 +- cgyro/src/cgyro_run.f90 | 6 +++--- 6 files changed, 31 insertions(+), 27 deletions(-) diff --git a/cgyro/src/cgyro_freq.f90 b/cgyro/src/cgyro_freq.f90 index 1be886ffe..0cc6b853b 100644 --- a/cgyro/src/cgyro_freq.f90 +++ b/cgyro/src/cgyro_freq.f90 @@ -11,6 +11,7 @@ subroutine cgyro_freq use cgyro_globals + use cgyro_io implicit none @@ -53,7 +54,7 @@ subroutine cgyro_freq myfr = total_weighted_freq/total_weight freq(itor) = myfr - + ! Fractional Frequency Error dfr = 0.0 dfi = 0.0 @@ -63,11 +64,17 @@ subroutine cgyro_freq 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 (abs(myfr) > 1e-12) then + ! Trap a division-by-zero error and halt + freq_err(itor) = (dfr+i_c*dfi)/total_weight/abs(myfr) + else + call cgyro_error('Underflow in calculation of frequency error') + endif + enddo - if (n_toroidal == 1 .and. abs(freq_err(nt1)) < freq_tol) signal=1 + if (n_toroidal == 1 .and. abs(freq_err(nt1)) < freq_tol) signal = 1 endif diff --git a/cgyro/src/cgyro_init_kernel.F90 b/cgyro/src/cgyro_init_kernel.F90 index 10d4e1033..5bc3c17c1 100644 --- a/cgyro/src/cgyro_init_kernel.F90 +++ b/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 diff --git a/cgyro/src/cgyro_init_manager.F90 b/cgyro/src/cgyro_init_manager.F90 index 2476803e8..087aec964 100644 --- a/cgyro/src/cgyro_init_manager.F90 +++ b/cgyro/src/cgyro_init_manager.F90 @@ -26,7 +26,6 @@ subroutine cgyro_init_manager implicit none - character(len=128) :: msg integer :: ie,ix @@ -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,7 +347,7 @@ 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 @@ -356,7 +355,7 @@ subroutine cgyro_init_manager 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 @@ -394,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 diff --git a/cgyro/src/cgyro_kernel.F90 b/cgyro/src/cgyro_kernel.F90 index 10f2d3c8d..a93717eae 100644 --- a/cgyro/src/cgyro_kernel.F90 +++ b/cgyro/src/cgyro_kernel.F90 @@ -92,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) @@ -110,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 @@ -126,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') @@ -158,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 > 1) exit enddo !--------------------------------------------------------------------------- diff --git a/cgyro/src/cgyro_restart.F90 b/cgyro/src/cgyro_restart.F90 index 014e51cf7..c021dfc36 100644 --- a/cgyro/src/cgyro_restart.F90 +++ b/cgyro/src/cgyro_restart.F90 @@ -212,7 +212,7 @@ subroutine cgyro_write_restart_one call MPI_BARRIER(CGYRO_COMM_WORLD,i_err) if (i_proc == 0) then call cgyro_write_restart_header_part - if (error_status /=0 ) return + if (error_status > 0) return endif ! now that we know things worked well, move the file in its final location diff --git a/cgyro/src/cgyro_run.f90 b/cgyro/src/cgyro_run.f90 index 1bf808d7e..3a4cc4479 100644 --- a/cgyro/src/cgyro_run.f90 +++ b/cgyro/src/cgyro_run.f90 @@ -34,7 +34,7 @@ subroutine cgyro_run(test_flag_in,var_in,n_species_out,flux_tave_out,tave_min_ou test_flag = test_flag_in ! Re-set max time - if(var_in > 0.0) then + if (var_in > 0.0) then max_time = var_in endif @@ -42,8 +42,8 @@ subroutine cgyro_run(test_flag_in,var_in,n_species_out,flux_tave_out,tave_min_ou call cgyro_kernel if (error_status == 0) then - if (nonlinear_flag == 0 .and. signal == 1) then - ! linear converged + if (nonlinear_flag == 0 .and. signal > 1) then + ! linear converged or underflow status_out = 1 endif ! no errors From 0ef97cfffd30c716537dc870f17ebef776d01e03 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Wed, 22 May 2024 11:23:01 -0700 Subject: [PATCH 09/16] Fixed formatting for negative bunit and rho --- cgyro/src/cgyro_write_initdata.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cgyro/src/cgyro_write_initdata.f90 b/cgyro/src/cgyro_write_initdata.f90 index bf08ad24d..cab06fbc2 100644 --- a/cgyro/src/cgyro_write_initdata.f90 +++ b/cgyro/src/cgyro_write_initdata.f90 @@ -146,7 +146,7 @@ subroutine cgyro_write_initdata dn = rho/(rhos/a_meters) kyrat = abs(q/rmin*rhos/a_meters) write(io,*) - write(io,10) ' a[m]:',a_meters,' b_unit[T]:',b_unit, ' rhos/a:',rhos/a_meters,' dn:',dn + write(io,10) ' a[m]:',a_meters, '|b_unit[T]|:',abs(b_unit), ' |rhos/a|:',abs(rhos)/a_meters,' dn:',dn write(io,10) 'n_norm[e19/m^3]:',dens_norm,'v_norm[m/s]:',vth_norm,'T_norm[keV]:',temp_norm write(io,*) write(io,'(t2,a)') ' n = 1 2 3 4 5 6 7 8' From bf6aebea08ce52d9b6f7c1f09f65f53069ec9140 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Thu, 23 May 2024 19:36:49 -0700 Subject: [PATCH 10/16] Fixed recently-introduced frequency tolerance bug --- cgyro/src/cgyro_freq.f90 | 2 +- cgyro/src/cgyro_init_h.f90 | 16 +++-------- cgyro/src/cgyro_kernel.F90 | 2 +- cgyro/src/cgyro_restart.F90 | 56 ++++++++++++++++++------------------- cgyro/src/cgyro_run.f90 | 2 +- 5 files changed, 35 insertions(+), 43 deletions(-) diff --git a/cgyro/src/cgyro_freq.f90 b/cgyro/src/cgyro_freq.f90 index 0cc6b853b..e133516a2 100644 --- a/cgyro/src/cgyro_freq.f90 +++ b/cgyro/src/cgyro_freq.f90 @@ -66,9 +66,9 @@ subroutine cgyro_freq enddo if (abs(myfr) > 1e-12) then - ! Trap a division-by-zero error and halt freq_err(itor) = (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 diff --git a/cgyro/src/cgyro_init_h.f90 b/cgyro/src/cgyro_init_h.f90 index f3a37395b..c545afce0 100644 --- a/cgyro/src/cgyro_init_h.f90 +++ b/cgyro/src/cgyro_init_h.f90 @@ -8,7 +8,7 @@ subroutine cgyro_init_h implicit none integer :: ir,it,is,ie,ix,itor - real :: arg, ang + real :: arg,ang !--------------------------------------------------------------------------- ! Check to see if we have restart data available @@ -50,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 diff --git a/cgyro/src/cgyro_kernel.F90 b/cgyro/src/cgyro_kernel.F90 index a93717eae..95c08273f 100644 --- a/cgyro/src/cgyro_kernel.F90 +++ b/cgyro/src/cgyro_kernel.F90 @@ -159,7 +159,7 @@ subroutine cgyro_kernel if (mod(i_time,print_step) == 0) call write_timers(trim(path)//runfile_timers) ! Exit if converged (1) or underflow (2) - if (signal > 1) exit + if (signal > 0) exit enddo !--------------------------------------------------------------------------- diff --git a/cgyro/src/cgyro_restart.F90 b/cgyro/src/cgyro_restart.F90 index c021dfc36..cbada0800 100644 --- a/cgyro/src/cgyro_restart.F90 +++ b/cgyro/src/cgyro_restart.F90 @@ -217,37 +217,37 @@ subroutine cgyro_write_restart_one ! now that we know things worked well, move the file in its final location if (i_proc == 0) then - if (restart_preservation_mode>2) then - ! restart_preservation_mode == 3 or 4 - ! First try to save any existing restart file as old - i_err = RENAME(trim(path)//runfile_restart, trim(path)//runfile_restart//".old") - ! NOTE: We will not check if it succeeded... not important, may not even exist (yet) - endif + if (restart_preservation_mode>2) then + ! restart_preservation_mode == 3 or 4 + ! First try to save any existing restart file as old + i_err = RENAME(trim(path)//runfile_restart, trim(path)//runfile_restart//".old") + ! NOTE: We will not check if it succeeded... not important, may not even exist (yet) + endif - ! Rename part into the final expected file name - i_err = RENAME(trim(path)//runfile_restart//".part", trim(path)//runfile_restart) - if (i_err /= 0) then - call cgyro_error('Final rename in cgyro_write_restart failed') - return - endif + ! Rename part into the final expected file name + i_err = RENAME(trim(path)//runfile_restart//".part", trim(path)//runfile_restart) + if (i_err /= 0) then + call cgyro_error('Final rename in cgyro_write_restart failed') + return + endif endif call system_clock(cp_time,count_rate,count_max) if (cp_time.gt.start_time) then - cp_dt = (cp_time-start_time)/real(count_rate) + cp_dt = (cp_time-start_time)/real(count_rate) else - cp_dt = (cp_time-start_time+count_max)/real(count_rate) + cp_dt = (cp_time-start_time+count_max)/real(count_rate) endif if (i_proc == 0) then - call date_and_time(sdate,stime); - call get_environment_variable('GACODE_PLATFORM',platform) - open(NEWUNIT=statusfd,FILE=trim(path)//runfile_startups,action='write',status='unknown',position='append') - write(statusfd,'(14(a),f7.3)') & - sdate(1:4),'/',sdate(5:6),'/',sdate(7:8),' ', & - stime(1:2),':',stime(3:4),':',stime(5:6),' ', & - trim(platform),' [CHECKPOINT WRITE] Time =',cp_dt - close(statusfd) + call date_and_time(sdate,stime); + call get_environment_variable('GACODE_PLATFORM',platform) + open(NEWUNIT=statusfd,FILE=trim(path)//runfile_startups,action='write',status='unknown',position='append') + write(statusfd,'(14(a),f7.3)') & + sdate(1:4),'/',sdate(5:6),'/',sdate(7:8),' ', & + stime(1:2),':',stime(3:4),':',stime(5:6),' ', & + trim(platform),' [CHECKPOINT WRITE] Time =',cp_dt + close(statusfd) endif ! Re-set h(0,0)=0 @@ -413,7 +413,7 @@ subroutine cgyro_read_restart_verify read(io,REC=recid) magic recid = recid + 1 - if ( magic == restart_magic) then + if (magic == restart_magic) then read(io,REC=recid) version recid = recid + 1 if ( version /= restart_version) then @@ -430,7 +430,7 @@ subroutine cgyro_read_restart_verify ! alternative ordering, needed different magic restart_magic_v2 = 140916753 endif - if ( magic /= restart_magic_v2) then + if (magic /= restart_magic_v2) then ! we don't support anything else, abort close(io) call cgyro_error('Wrong magic number in restart header') @@ -465,7 +465,7 @@ subroutine cgyro_read_restart_verify return endif - if ( magic == restart_magic) then + if (magic == restart_magic) then read(io,REC=recid) t_velocity_order recid = recid + 1 read(io,REC=recid) t_nt_loc @@ -520,7 +520,7 @@ subroutine cgyro_read_restart_one if (i_proc == 0) then call cgyro_read_restart_verify - if (error_status /=0 ) return + if (error_status > 0) return mpibuf(1) = t_nt_loc mpibuf(2) = t_nv_loc mpibuf(3) = t_velocity_order @@ -648,7 +648,7 @@ subroutine cgyro_read_restart_slow ! Velocity pointers for the restart velocity_order t_iv = 0 - if (t_velocity_order==1) then + if (t_velocity_order == 1) then !original do ie=1,n_energy do ix=1,n_xi @@ -658,7 +658,7 @@ subroutine cgyro_read_restart_slow enddo enddo enddo - else if (t_velocity_order==2) then + else if (t_velocity_order == 2) then ! optimized for minimizing species do is=1,n_species do ie=1,n_energy diff --git a/cgyro/src/cgyro_run.f90 b/cgyro/src/cgyro_run.f90 index 3a4cc4479..67f043b91 100644 --- a/cgyro/src/cgyro_run.f90 +++ b/cgyro/src/cgyro_run.f90 @@ -42,7 +42,7 @@ subroutine cgyro_run(test_flag_in,var_in,n_species_out,flux_tave_out,tave_min_ou call cgyro_kernel if (error_status == 0) then - if (nonlinear_flag == 0 .and. signal > 1) then + if (nonlinear_flag == 0 .and. signal > 0) then ! linear converged or underflow status_out = 1 endif From 37d8cd38e649126056137ece5f63ccc06f1334bd Mon Sep 17 00:00:00 2001 From: Igor Sfiligoi Date: Fri, 24 May 2024 15:10:23 -0500 Subject: [PATCH 11/16] Compute freq_err only for n_toroidal == 1, since it is not used else --- cgyro/src/cgyro_cleanup.F90 | 1 - cgyro/src/cgyro_freq.f90 | 40 +++++++++++++++--------------- cgyro/src/cgyro_globals.F90 | 2 +- cgyro/src/cgyro_init_manager.F90 | 2 +- cgyro/src/cgyro_write_timedata.f90 | 3 +-- 5 files changed, 23 insertions(+), 25 deletions(-) diff --git a/cgyro/src/cgyro_cleanup.F90 b/cgyro/src/cgyro_cleanup.F90 index 6cf05d21c..8ae59dae8 100644 --- a/cgyro/src/cgyro_cleanup.F90 +++ b/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) diff --git a/cgyro/src/cgyro_freq.f90 b/cgyro/src/cgyro_freq.f90 index e133516a2..54ab2eb1e 100644 --- a/cgyro/src/cgyro_freq.f90 +++ b/cgyro/src/cgyro_freq.f90 @@ -25,7 +25,7 @@ subroutine cgyro_freq if (i_time == 0) then freq(:) = 0.0 - freq_err(:) = 0.0 + freq_err = 0.0 else @@ -54,28 +54,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 - - if (abs(myfr) > 1e-12) then - freq_err(itor) = (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') + + 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 diff --git a/cgyro/src/cgyro_globals.F90 b/cgyro/src/cgyro_globals.F90 index 06cf73f76..1c6efec0c 100644 --- a/cgyro/src/cgyro_globals.F90 +++ b/cgyro/src/cgyro_globals.F90 @@ -282,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 !--------------------------------------------------------------- diff --git a/cgyro/src/cgyro_init_manager.F90 b/cgyro/src/cgyro_init_manager.F90 index 087aec964..1876eacdd 100644 --- a/cgyro/src/cgyro_init_manager.F90 +++ b/cgyro/src/cgyro_init_manager.F90 @@ -151,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 diff --git a/cgyro/src/cgyro_write_timedata.f90 b/cgyro/src/cgyro_write_timedata.f90 index 221105eb4..6e56f7777 100644 --- a/cgyro/src/cgyro_write_timedata.f90 +++ b/cgyro/src/cgyro_write_timedata.f90 @@ -811,11 +811,10 @@ subroutine print_scrdata() '[t: ',t_current,& '][e: ',integration_error(:),']' else - ! NOTE: There could be only one my_toroidal when n_toroidal <=1 print '(a,1pe9.3,a,1pe10.3,1x,1pe10.3,a,1pe10.3,a,1pe9.3,1x,1pe9.3,a)',& '[t: ',t_current,& '][w: ',freq(nt1),& - '][dw:',abs(freq_err(nt1)),& + '][dw:',abs(freq_err),& '][e: ',integration_error(:),']' endif From 24e5b3517bbe2171282d90c7ba91c551269c6c49 Mon Sep 17 00:00:00 2001 From: Igor Sfiligoi Date: Fri, 24 May 2024 16:05:55 -0500 Subject: [PATCH 12/16] Force no parallelization, to work around ifx compiler bug. Will not add significantly to total runtime --- cgyro/src/cgyro_freq.f90 | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/cgyro/src/cgyro_freq.f90 b/cgyro/src/cgyro_freq.f90 index 54ab2eb1e..bf90c7fbb 100644 --- a/cgyro/src/cgyro_freq.f90 +++ b/cgyro/src/cgyro_freq.f90 @@ -21,6 +21,8 @@ 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 @@ -33,18 +35,23 @@ subroutine cgyro_freq ! Standard method: sum all wavenumbers at a given n !-------------------------------------------------- + icdt = i_c/delta_t +!$omp do ordered do itor=nt1,nt2 total_weight = 0.0 total_weighted_freq = (0.0,0.0) +!$omp do ordered 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 From e145265d2a09c00462c213f2505ea6c0bb041c05 Mon Sep 17 00:00:00 2001 From: Igor Sfiligoi Date: Fri, 24 May 2024 16:18:58 -0500 Subject: [PATCH 13/16] Use per-sim out file during regression testing --- shared/bin/gacode_reg_do | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/shared/bin/gacode_reg_do b/shared/bin/gacode_reg_do index 5bbdf7451..76a11f125 100755 --- a/shared/bin/gacode_reg_do +++ b/shared/bin/gacode_reg_do @@ -58,12 +58,12 @@ cd $testdir for sim in $list do - $code -g $sim -p $testdir > out + $code -g $sim -p $testdir > out.$sim rm -f $sim/$precfile if [ $n_omp -eq 1 ] && [ $n_proc -eq 1 ] ; then - $code -e $sim -p $testdir > out + $code -e $sim -p $testdir > out.$sim else - $code -e $sim -n $n_proc -nomp $n_omp -p $testdir > out + $code -e $sim -n $n_proc -nomp $n_omp -p $testdir > out.$sim fi gacode_reg $sim $compdir $precfile $tol if [ $reset -eq 1 ] From 50fea94041e3d32b0283ad389f62f9d7f2f511a7 Mon Sep 17 00:00:00 2001 From: Igor Sfiligoi Date: Fri, 24 May 2024 16:27:39 -0500 Subject: [PATCH 14/16] Force no parallelization in ifx, to work around compiler bug. Will not add significantly to total runtime --- cgyro/src/cgyro_freq.f90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/cgyro/src/cgyro_freq.f90 b/cgyro/src/cgyro_freq.f90 index bf90c7fbb..12de3a6c3 100644 --- a/cgyro/src/cgyro_freq.f90 +++ b/cgyro/src/cgyro_freq.f90 @@ -36,12 +36,16 @@ subroutine cgyro_freq !-------------------------------------------------- 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 fo1 = field_old(1,ic,itor) From 0491502e863f9076f723821953def5f69f3dee52 Mon Sep 17 00:00:00 2001 From: Igor Sfiligoi Date: Fri, 24 May 2024 16:28:36 -0500 Subject: [PATCH 15/16] Rename to enable preprocessor --- cgyro/src/{cgyro_freq.f90 => cgyro_freq.F90} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename cgyro/src/{cgyro_freq.f90 => cgyro_freq.F90} (100%) diff --git a/cgyro/src/cgyro_freq.f90 b/cgyro/src/cgyro_freq.F90 similarity index 100% rename from cgyro/src/cgyro_freq.f90 rename to cgyro/src/cgyro_freq.F90 From c8712d9939eaf6e266f691dd982ec866c0a735a5 Mon Sep 17 00:00:00 2001 From: Anass Najlaoui Date: Fri, 7 Jun 2024 14:45:21 +0200 Subject: [PATCH 16/16] Implementation of new filter that will not filter out KBMs --- tglf/src/tglf_LS.f90 | 120 ++++++++++++++++++++++++++++++++++ tglf/src/tglf_eigensolver.f90 | 1 + 2 files changed, 121 insertions(+) diff --git a/tglf/src/tglf_LS.f90 b/tglf/src/tglf_LS.f90 index ada97c0cc..2d4f237a7 100644 --- a/tglf/src/tglf_LS.f90 +++ b/tglf/src/tglf_LS.f90 @@ -34,6 +34,17 @@ SUBROUTINE tglf_LS INTEGER,ALLOCATABLE,DIMENSION(:) :: di,de INTEGER,ALLOCATABLE,DIMENSION(:) :: ipiv COMPLEX,ALLOCATABLE,DIMENSION(:,:) :: zmat +! Filter numerical instabilities + REAL, DIMENSION(max_plot) :: phi_test + REAL :: eps_even + logical :: is_odd + logical :: is_even + logical :: is_sup_to_drift + REAL :: sum_phi, sum_abs_phi + REAL, DIMENSION(max_plot) :: even_function, odd_function + REAL :: max_freq + COMMON /block/ max_freq + ! ! cputime0=MPI_WTIME() ! @@ -260,6 +271,115 @@ SUBROUTINE tglf_LS enddo endif ! +! Filter out numerical instabilities with a mode frequency higher than filter*drift wave freq +! But not modes with an even structure for Re(phi) + if (filter_in.lt.0.0 .and. ibranch_in.eq.-1) then + is_sup_to_drift = .false. + do imax=1,nmodes_out + if(rr(jmax(imax)).gt.0.0.and.ABS(ri(jmax(imax))).gt.ABS(max_freq))then + is_sup_to_drift = .true. + else + is_sup_to_drift = .false. + endif + + if (is_sup_to_drift) then + if(jmax(imax).gt.0)then + do i=1,iur + v(i) = small + do j=1,iur + zmat(i,j) = beta(jmax(imax))*amat(i,j)- & + (small +alpha(jmax(imax)))*bmat(i,j) + enddo + enddo + call zgesv(iur,1,zmat,iur,ipiv,v,iur,info) + if (info /= 0)CALL tglf_error(1,"ZGESV failed in tglf_LS") + eigenvalue = xi*alpha(jmax(imax))/beta(jmax(imax)) + call get_QL_weights + + do i = 1, max_plot + odd_function(i) = COS(pi * REAL(i-1) / REAL(max_plot)) + even_function(i) = SIN(pi * REAL(i-1) / REAL(max_plot)) + enddo + + do i = 1, max_plot + phi_test(i) = 0.0 + do j = 1, nbasis + if (MOD(j, 2) == 0) then + phi_test(i) = phi_test(i) + odd_function(i) * REAL(field_weight_QL_out(1, j)) + else + phi_test(i) = phi_test(i) + even_function(i) * REAL(field_weight_QL_out(1, j)) + endif + enddo + enddo + + eps_even = 0.1 ! 10% threshold + + ! Compute the sum of the function values and the sum of their absolute values + sum_phi = SUM(phi_test(:)) + sum_abs_phi = SUM(ABS(phi_test(:))) + + ! Check if the function is even or odd based on the sum/integral comparison + is_even = .false. + is_odd = .false. + if (ABS(sum_phi) > (1-eps_even) * sum_abs_phi) then + is_even = .true. + endif + if (ABS(sum_phi) < eps_even * sum_abs_phi) then + is_odd = .true. + endif + endif + + ! Check if the mode is a numerical instability + if (.not. (is_even .and. rr(jmax(imax)) > 0.0 .and. ri(jmax(imax)) > 0) ) then + rr(jmax(imax)) = -rr(jmax(imax)) + endif + endif + enddo +! Re-initialize output to zero + do j1=1,maxmodes + jmax(j1) = 0 + gamma_out(j1) = 0.0 + freq_out(j1) = 0.0 + v_QL_out(j1) = 0.0 + a_par_QL_out(j1) = 0.0 + b_par_QL_out(j1) = 0.0 + phi_bar_out(j1) = 0.0 + a_par_bar_out(j1) = 0.0 + b_par_bar_out(j1) = 0.0 + v_bar_out(j1) = 0.0 + do is=1,nsm + do j=1,3 + particle_QL_out(j1,is,j) = 0.0 + energy_QL_out(j1,is,j) = 0.0 + stress_par_QL_out(j1,is,j) = 0.0 + stress_tor_QL_out(j1,is,j) = 0.0 + exchange_QL_out(j1,is,j) = 0.0 + enddo + N_QL_out(j1,is) = 0.0 + T_QL_out(j1,is) = 0.0 + U_QL_out(j1,is) = 0.0 + Q_QL_out(j1,is) = 0.0 + N_bar_out(j1,is) = 0.0 + T_bar_out(j1,is) = 0.0 + U_bar_out(j1,is) = 0.0 + Q_bar_out(j1,is) = 0.0 + Ns_Ts_phase_out(j1,is) = 0.0 + enddo + ne_te_phase_out(j1) = 0.0 + enddo +! Re-find the top nmodes most unstable modes + CALL sort_eigenvalues(nmodes_in,jmax) + nmodes_out = 0 + do j1=1,nmodes_in + if(jmax(j1).ne.0)then + nmodes_out = nmodes_out + 1 + gamma_out(j1)=rr(jmax(j1)) + freq_out(j1)=-ri(jmax(j1)) + endif + enddo + endif +! End of filtering out numerical instabilities + ! get the fluxes for the most unstable modes if(iflux_in)then do imax=1,nmodes_out diff --git a/tglf/src/tglf_eigensolver.f90 b/tglf/src/tglf_eigensolver.f90 index 117975af6..ae6f42bce 100644 --- a/tglf/src/tglf_eigensolver.f90 +++ b/tglf/src/tglf_eigensolver.f90 @@ -91,6 +91,7 @@ SUBROUTINE tglf_eigensolver COMPLEX,ALLOCATABLE,DIMENSION(:,:) :: vleft,vright COMPLEX,ALLOCATABLE,DIMENSION(:) :: work COMPLEX,ALLOCATABLE,DIMENSION(:) :: zomega + COMMON /block/ max_freq ! ifail = 0 lwork = 8*iur