From 5228831eef37cb597e3abfc81bb2f2b1cb95ce7e Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Wed, 8 Mar 2023 11:42:36 +0100
Subject: [PATCH] code compiles, runs, and output correclty but wrong values
---
src/array_mod.F90 | 3 -
src/diagnose.F90 | 287 ++++++++++-----------
src/grid_mod.F90 | 22 +-
src/inital.F90 | 10 +-
src/memory.F90 | 1 -
src/parallel_mod.F90 | 20 +-
src/processing_mod.F90 | 18 +-
testcases/smallest_problem/fort.90 | 35 ++-
testcases/smallest_problem/fort_00.90 | 17 +-
testcases/smallest_problem/gyacomo_1 | 1 +
testcases/smallest_problem/gyacomo_1_debug | 1 +
testcases/smallest_problem/gyacomo_alpha | 1 -
wk/header_3D_results.m | 2 +-
13 files changed, 204 insertions(+), 214 deletions(-)
create mode 120000 testcases/smallest_problem/gyacomo_1
create mode 120000 testcases/smallest_problem/gyacomo_1_debug
delete mode 120000 testcases/smallest_problem/gyacomo_alpha
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index d6f1fdeb..da4d7c9b 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -48,9 +48,6 @@ MODULE array
REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: inv_pol_ion
! REAL(dp), DIMENSION(:,:,:), ALLOCATABLE :: HF_phi_correction_operator
- ! Gyrocenter density for electron and ions (ia,iky,ikx,iz)
- COMPLEX(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Na00
-
! Kinetic spectrum sum_kx,ky(|Napj(z)|^2), (ia,ip,ij,iz) (should be real)
REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Napjz
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 20872259..a98d5981 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -5,10 +5,8 @@ SUBROUTINE diagnose(kstep)
USE processing, ONLY: pflux_x, hflux_x
USE parallel, ONLY: my_id
IMPLICIT NONE
-
INTEGER, INTENT(in) :: kstep
CALL cpu_time(t0_diag) ! Measuring time
-
!! Basic diagnose loop for reading input file, displaying advancement and ending
IF ((kstep .EQ. 0)) THEN
INQUIRE(unit=lu_in, name=input_fname)
@@ -51,7 +49,6 @@ SUBROUTINE init_outfile(comm,file0,file,fid)
INTEGER, INTENT(OUT) :: fid
CHARACTER(len=256) :: str
INCLUDE 'srcinfo.h'
-
! Writing output filename
WRITE(file,'(a,a1,i2.2,a3)') TRIM(file0) ,'_',jobnum,'.h5'
! 1.1 Initial run
@@ -67,15 +64,14 @@ SUBROUTINE init_outfile(comm,file0,file,fid)
! File group
CALL creatg(fid, "/files", "files")
CALL attach(fid, "/files", "jobnum", jobnum)
-
! Add the code info and parameters to the file
CALL creatg(fid, "/data/input", "input")
CALL creatd(fid, 0,(/0/),"/data/input/codeinfo",'Code Information')
- CALL attach(fid, TRIM(str), "version", VERSION) !defined in srcinfo.h
- CALL attach(fid, TRIM(str), "branch", BRANCH) !defined in srcinfo.h
- CALL attach(fid, TRIM(str), "author", AUTHOR) !defined in srcinfo.h
- CALL attach(fid, TRIM(str), "execdate", EXECDATE) !defined in srcinfo.h
- CALL attach(fid, TRIM(str), "host", HOST) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "version", VERSION) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "branch", BRANCH) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "author", AUTHOR) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "execdate", EXECDATE) !defined in srcinfo.h
+ CALL attach(fid, "/data/input/codeinfo", "host", HOST) !defined in srcinfo.h
CALL basic_outputinputs(fid)
CALL grid_outputinputs(fid)
CALL geometry_outputinputs(fid)
@@ -85,7 +81,6 @@ SUBROUTINE init_outfile(comm,file0,file,fid)
CALL coll_outputinputs(fid)
CALL initial_outputinputs(fid)
CALL time_integration_outputinputs(fid)
-
! Save STDIN (input file) of this run
IF(jobnum .LE. 99) THEN
WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum
@@ -97,38 +92,31 @@ END SUBROUTINE init_outfile
SUBROUTINE diagnose_full(kstep)
USE basic, ONLY: speak,&
- cstep,iframe0d,iframe2d,iframe3d,iframe5d,&
+ cstep,iframe0d,iframe3d,iframe5d,&
start,finish,crashed
USE grid, ONLY: &
- ias,iae, &
- parray_full,pmax,&
- ijs,ije,jarray_full,jmax,&
- ikys,ikye,kyarray_full,&
- ikxs,ikxe,kxarray_full,&
- izs,ize,zarray_full
+ local_nj,local_nky,local_nkx,local_nz,ngj,ngz,&
+ parray_full,pmax,jarray_full,jmax,&
+ kyarray_full,kxarray_full,zarray_full, total_na
USE diagnostics_par
USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf, putarrnd
+ USE species, ONLY: name
USE array
- USE model, ONLY:
- USE initial_par, ONLY:
- USE fields, ONLY:
- USE time_integration,ONLY:
+ USE model, ONLY: EM
USE parallel, ONLY: my_id, comm0
- USE prec_const, ONLY:
USE collision, ONLY: coll_outputinputs
USE geometry, ONLY: gxx,gxy,gyy,gxz,gyz,gzz,hatR,hatZ,hatB,dBdx,dBdy,dBdz,Jacobian,gradz_coeff,Ckxky
IMPLICIT NONE
-
+ CHARACTER :: letter_a
INTEGER, INTENT(in) :: kstep
INTEGER, parameter :: BUFSIZE = 2
- INTEGER :: rank = 0, ierr
+ INTEGER :: rank = 0, ierr, ia
INTEGER :: dims(1) = (/0/)
!____________________________________________________________________________
! 1. Initial diagnostics
IF ((kstep .EQ. 0)) THEN
CALL init_outfile(comm0, resfile0,resfile,fidres)
-
! Profiler time measurement
CALL creatg(fidres, "/profiler", "performance analysis")
CALL creatd(fidres, 0, dims, "/profiler/Tc_rhs", "cumulative rhs computation time")
@@ -143,7 +131,6 @@ SUBROUTINE diagnose_full(kstep)
CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative sym computation time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_step", "cumulative total step computation time")
CALL creatd(fidres, 0, dims, "/profiler/time", "current simulation time")
-
! Grid info
CALL creatg(fidres, "/data/grid", "Grid data")
CALL putarr(fidres, "/data/grid/coordkx", kxarray_full, "kx*rho_s0", ionode=0)
@@ -151,106 +138,87 @@ SUBROUTINE diagnose_full(kstep)
CALL putarr(fidres, "/data/grid/coordz", zarray_full, "z/R", ionode=0)
CALL putarr(fidres, "/data/grid/coordp" , parray_full, "p", ionode=0)
CALL putarr(fidres, "/data/grid/coordj" , jarray_full, "j", ionode=0)
-
! Metric info
CALL creatg(fidres, "/data/metric", "Metric data")
- CALL putarrnd(fidres, "/data/metric/gxx", gxx(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gxy", gxy(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gxz", gxz(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gyy", gyy(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gyz", gyz(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gzz", gzz(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/hatR", hatR(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/hatZ", hatZ(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/hatB", hatB(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/dBdx", dBdx(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/dBdy", dBdy(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/dBdz", dBdz(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(izs:ize,0:1), (/1, 1, 1/))
- CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/1, 1, 3/))
- CALL putarrnd(fidres, "/data/metric/kernel", kernel(ias:iae,ijs:ije,ikys:ikye,ikxs:ikxe,izs:ize,0:1), (/1, 1, 2, 4/))
-
+ CALL putarrnd(fidres, "/data/metric/gxx", gxx(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gxy", gxy(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gxz", gxz(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gyy", gyy(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gyz", gyz(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gzz", gzz(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatR", hatR(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatZ", hatZ(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatB", hatB(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/dBdx", dBdx(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/dBdy", dBdy(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/dBdz", dBdz(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(1+ngz/2:local_nz+ngz/2,:), (/1, 1, 1/))
+ CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,:), (/1, 1, 3/))
+ CALL putarrnd(fidres, "/data/metric/kernel", kernel(1,1+ngj/2:local_nj+ngj/2,1:local_nky,1:local_nkx,1+ngz/2:local_nz+ngz/2,:), (/1, 1, 2, 4/))
! var0d group (gyro transport)
IF (nsave_0d .GT. 0) THEN
CALL creatg(fidres, "/data/var0d", "0d profiles")
CALL creatd(fidres, rank, dims, "/data/var0d/time", "Time t*c_s/R")
CALL creatd(fidres, rank, dims, "/data/var0d/cstep", "iteration number")
-
IF (write_gamma) THEN
- CALL creatd(fidres, rank, dims, "/data/var0d/gflux_x", "Radial gyro transport")
- CALL creatd(fidres, rank, dims, "/data/var0d/pflux_x", "Radial part transport")
+ DO ia=1,total_na
+ letter_a = name(ia)(1:1)
+ CALL creatd(fidres, rank, dims, "/data/var0d/gflux_x"//letter_a, "Radial gyro transport")
+ CALL creatd(fidres, rank, dims, "/data/var0d/pflux_x"//letter_a, "Radial part transport")
+ ENDDO
ENDIF
IF (write_hf) THEN
- CALL creatd(fidres, rank, dims, "/data/var0d/hflux_x", "Radial part ion heat flux")
+ DO ia=1,total_na
+ letter_a = name(ia)(1:1)
+ CALL creatd(fidres, rank, dims, "/data/var0d/hflux_x"//letter_a, "Radial part heat flux")
+ ENDDO
ENDIF
IF (cstep==0) THEN
iframe0d=0
ENDIF
CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
END IF
-
-
- ! var2d group (??)
- IF (nsave_2d .GT. 0) THEN
- CALL creatg(fidres, "/data/var2d", "2d profiles")
- CALL creatd(fidres, rank, dims, "/data/var2d/time", "Time t*c_s/R")
- CALL creatd(fidres, rank, dims, "/data/var2d/cstep", "iteration number")
- IF (cstep==0) THEN
- iframe2d=0
- ENDIF
- CALL attach(fidres,"/data/var2d/" , "frames", iframe2d)
- END IF
-
- ! var3d group (electro. pot., Ni00 moment)
+ ! var3d group (phi,psi, fluid moments, Ni00, Napjz)
IF (nsave_3d .GT. 0) THEN
CALL creatg(fidres, "/data/var3d", "3d profiles")
CALL creatd(fidres, rank, dims, "/data/var3d/time", "Time t*c_s/R")
CALL creatd(fidres, rank, dims, "/data/var3d/cstep", "iteration number")
-
IF (write_phi) CALL creatg(fidres, "/data/var3d/phi", "phi")
- IF (write_phi) CALL creatg(fidres, "/data/var3d/psi", "psi")
-
- IF (write_Na00) THEN
- CALL creatg(fidres, "/data/var3d/Na00", "Na00")
- CALL creatg(fidres, "/data/var3d/Nepjz", "Napjz")
- ENDIF
-
- IF (write_dens) THEN
- CALL creatg(fidres, "/data/var3d/dens", "dens_e")
- ENDIF
-
- IF (write_fvel) THEN
- CALL creatg(fidres, "/data/var3d/upar", "upar")
- CALL creatg(fidres, "/data/var3d/uper", "uper")
- ENDIF
-
- IF (write_temp) THEN
- CALL creatg(fidres, "/data/var3d/Tper_e", "Tper")
- CALL creatg(fidres, "/data/var3d/Tpar_e", "Tpar")
- CALL creatg(fidres, "/data/var3d/temp_e", "temp")
- ENDIF
-
+ IF (write_phi.AND.EM) CALL creatg(fidres, "/data/var3d/psi", "psi")
+ ! Loop to create species related data
+ DO ia=1,total_na
+ letter_a = name(ia)(1:1)
+ IF (write_Na00) THEN
+ CALL creatg(fidres, "/data/var3d/N"//letter_a//"00", "gyroceneter density "//letter_a)
+ CALL creatg(fidres, "/data/var3d/N"//letter_a//"pjz", "pj(z) moment spectrum "//letter_a)
+ ENDIF
+ IF (write_dens) THEN
+ CALL creatg(fidres, "/data/var3d/dens_"//letter_a, "density "//letter_a)
+ ENDIF
+ IF (write_fvel) THEN
+ CALL creatg(fidres, "/data/var3d/upar_"//letter_a, "parallel fluid velocity "//letter_a)
+ CALL creatg(fidres, "/data/var3d/uper_"//letter_a, "perpendicular fluid velocity "//letter_a)
+ ENDIF
+ IF (write_temp) THEN
+ CALL creatg(fidres, "/data/var3d/Tper_"//letter_a, "perpendicular temperature "//letter_a)
+ CALL creatg(fidres, "/data/var3d/Tpar_"//letter_a, "parallel temperature "//letter_a)
+ CALL creatg(fidres, "/data/var3d/temp_"//letter_a, "tiotal temperature "//letter_a)
+ ENDIF
+ ENDDO
IF (cstep==0) THEN
iframe3d=0
ENDIF
CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
END IF
-
! var5d group (moments)
IF (nsave_5d .GT. 0) THEN
CALL creatg(fidres, "/data/var5d", "5d profiles")
CALL creatd(fidres, rank, dims, "/data/var5d/time", "Time t*c_s/R")
CALL creatd(fidres, rank, dims, "/data/var5d/cstep", "iteration number")
-
IF (write_Napj) THEN
- CALL creatg(fidres, "/data/var5d/moments", "moments")
- ENDIF
-
- IF (write_Sapj) THEN
- CALL creatg(fidres, "/data/var5d/Sapj", "Sapj")
+ CALL creatg(fidres, "/data/var5d/moments", "full moments array")
ENDIF
-
IF (cstep==0) THEN
iframe5d=0
END IF
@@ -258,6 +226,7 @@ SUBROUTINE diagnose_full(kstep)
END IF
ENDIF
+
!_____________________________________________________________________________
! 2. Periodic diagnostics
!
@@ -269,24 +238,15 @@ SUBROUTINE diagnose_full(kstep)
CALL diagnose_0d
END IF
END IF
-
- ! 2.2 1d profiles
- ! empty in our case
-
- ! 2.3 2d profiles
- ! empty in our case
-
-
! 2.3 3d profiles
IF (nsave_3d .GT. 0) THEN
IF (MOD(cstep, nsave_3d) == 0) THEN
CALL diagnose_3d
- ! Looks at the folder if the file check_phi exists and spits a snapshot
+ ! Looks at the folder if the file check_phi exists and spits a snapshot
! of the current electrostatic potential in a basic text file
CALL spit_snapshot_check
ENDIF
ENDIF
-
! 2.4 5d profiles
IF (nsave_5d .GT. 0 .AND. cstep .GT. 0) THEN
IF (MOD(cstep, nsave_5d) == 0) THEN
@@ -316,14 +276,16 @@ END SUBROUTINE diagnose_full
!!-------------- Auxiliary routines -----------------!!
SUBROUTINE diagnose_0d
-
USE basic
USE futils, ONLY: append, attach, getatt
USE diagnostics_par
USE prec_const
USE processing
-
+ USE model, ONLY: Na
+ USE species, ONLY: name
IMPLICIT NONE
+ CHARACTER :: letter_a
+ INTEGER :: ia
! Time measurement data
CALL append(fidres, "/profiler/Tc_rhs", tc_rhs,ionode=0)
CALL append(fidres, "/profiler/Tc_adv_field", tc_adv_field,ionode=0)
@@ -340,18 +302,24 @@ SUBROUTINE diagnose_0d
! Processing data
CALL append(fidres, "/data/var0d/time", time,ionode=0)
CALL append(fidres, "/data/var0d/cstep", real(cstep,dp),ionode=0)
- CALL getatt(fidres, "/data/var0d/", "frames",iframe2d)
+ CALL getatt(fidres, "/data/var0d/", "frames",iframe0d)
iframe0d=iframe0d+1
CALL attach(fidres,"/data/var0d/" , "frames", iframe0d)
! Ion transport data
IF (write_gamma) THEN
CALL compute_radial_transport
- CALL append(fidres, "/data/var0d/gflux_x",gflux_x(1),ionode=0)
- CALL append(fidres, "/data/var0d/pflux_x",pflux_x(1),ionode=0)
+ DO ia=1,Na
+ letter_a = name(ia)(1:1)
+ CALL append(fidres, "/data/var0d/gflux_x"//letter_a,gflux_x(ia),ionode=0)
+ CALL append(fidres, "/data/var0d/pflux_x"//letter_a,pflux_x(ia),ionode=0)
+ ENDDO
ENDIF
IF (write_hf) THEN
CALL compute_radial_heatflux
- CALL append(fidres, "/data/var0d/hflux_x",hflux_x(1),ionode=0)
+ DO ia=1,Na
+ letter_a = name(ia)(1:1)
+ CALL append(fidres, "/data/var0d/hflux_x"//letter_a,hflux_x(ia),ionode=0)
+ ENDDO
ENDIF
END SUBROUTINE diagnose_0d
@@ -359,8 +327,8 @@ SUBROUTINE diagnose_3d
USE basic
USE futils, ONLY: append, getatt, attach, putarrnd, putarr
USE fields, ONLY: phi, psi, moments
- USE array, ONLY: Na00,Napjz,dens,upar,uper,Tpar,Tper,temp
- USE grid, ONLY: CONTAINSp0, ip0,ij0, &
+ USE array, ONLY: Napjz,dens,upar,uper,Tpar,Tper,temp
+ USE grid, ONLY: CONTAINSp0, ip0,ij0, local_na,&
total_np, total_nj, total_nky, total_nkx, total_nz, &
local_np, local_nj, local_nky, local_nkx, local_nz, &
ngz
@@ -368,45 +336,65 @@ SUBROUTINE diagnose_3d
USE diagnostics_par
USE prec_const
USE processing, ONLY: compute_fluid_moments, compute_Napjz_spectrum
+ USE model, ONLY: EM
+ USE species, ONLY: name
IMPLICIT NONE
-
+ CHARACTER :: letter_a
+ INTEGER :: ia
+ COMPLEX(dp), DIMENSION(local_nky,local_nkx,local_nz) :: Na00_
+ COMPLEX(dp), DIMENSION(local_nky,local_nkx,local_nz) :: fmom
+ COMPLEX(dp), DIMENSION(local_np, local_nj, local_nz) :: Napjz_
+ ! add current time, cstep and frame
CALL append(fidres, "/data/var3d/time", time,ionode=0)
CALL append(fidres, "/data/var3d/cstep", real(cstep,dp),ionode=0)
CALL getatt(fidres, "/data/var3d/", "frames",iframe3d)
iframe3d=iframe3d+1
CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
-
- IF (write_phi) CALL write_field3d_kykxz(phi (:,:,1+ngz/2:local_nz+ngz/2), 'phi')
- IF (write_phi) CALL write_field3d_kykxz(psi (:,:,1+ngz/2:local_nz+ngz/2), 'psi')
-
+ ! Write current EM fields
+ IF (write_phi) CALL write_field3d_kykxz(phi (:,:,1+ngz/2:local_nz+ngz/2), 'phi')
+ IF (write_phi.AND.EM) CALL write_field3d_kykxz(psi (:,:,1+ngz/2:local_nz+ngz/2), 'psi')
IF (write_Na00) THEN
- IF (CONTAINSp0) &
- Na00(1,:,:,:) = moments(1,ip0,ij0,:,:,1+ngz/2:local_nz+ngz/2,updatetlevel)
- CALL write_field3d_kykxz(Na00(1,:,:,:), 'Na00')
-
CALL compute_Napjz_spectrum
- CALL write_field3d_pjz(Napjz(1,:,:,:), 'Napjz')
+ DO ia=1,local_na
+ letter_a = name(ia)(1:1)
+ IF (CONTAINSp0) THEN
+ ! gyrocenter density
+ Na00_ = moments(ia,ip0,ij0,:,:,1+ngz/2:local_nz+ngz/2,updatetlevel)
+ CALL write_field3d_kykxz(Na00_, 'N'//letter_a//'00')
+ ENDIF
+ ! <<Napj>x>y spectrum
+ Napjz_ = Napjz(ia,:,:,:)
+ CALL write_field3d_pjz(Napjz_, 'N'//letter_a//'pjz')
+ ENDDO
ENDIF
!! Fuid moments
IF (write_dens .OR. write_fvel .OR. write_temp) &
CALL compute_fluid_moments
- IF (write_dens) THEN
- CALL write_field3d_kykxz(dens(1,:,:,:), 'dens')
- ENDIF
-
- IF (write_fvel) THEN
- CALL write_field3d_kykxz(upar(1,:,:,:), 'upar')
- CALL write_field3d_kykxz(uper(1,:,:,:), 'uper')
- ENDIF
+ DO ia=1,local_na
+ letter_a = name(ia)(1:1)
+ IF (write_dens) THEN
+ fmom = dens(ia,:,:,:)
+ CALL write_field3d_kykxz(fmom, 'dens_'//letter_a)
+ ENDIF
- IF (write_temp) THEN
- CALL write_field3d_kykxz(Tpar(1,:,:,:), 'Tpar')
- CALL write_field3d_kykxz(Tper(1,:,:,:), 'Tper')
- CALL write_field3d_kykxz(temp(1,:,:,:), 'temp')
- ENDIF
+ IF (write_fvel) THEN
+ fmom = upar(ia,:,:,:)
+ CALL write_field3d_kykxz(fmom, 'upar_'//letter_a)
+ fmom = uper(ia,:,:,:)
+ CALL write_field3d_kykxz(fmom, 'uper_'//letter_a)
+ ENDIF
+ IF (write_temp) THEN
+ fmom = Tpar(ia,:,:,:)
+ CALL write_field3d_kykxz(fmom, 'Tpar_'//letter_a)
+ fmom = Tper(ia,:,:,:)
+ CALL write_field3d_kykxz(fmom, 'Tper_'//letter_a)
+ fmom = temp(ia,:,:,:)
+ CALL write_field3d_kykxz(fmom, 'temp_'//letter_a)
+ ENDIF
+ ENDDO
CONTAINS
SUBROUTINE write_field3d_kykxz(field, text)
USE basic, ONLY : GATHERV_OUTPUT
@@ -434,8 +422,8 @@ SUBROUTINE diagnose_3d
SUBROUTINE write_field3d_pjz(field, text)
USE parallel, ONLY : gather_pjz, num_procs
IMPLICIT NONE
- REAL(dp), DIMENSION(local_np,local_nj,local_nz), INTENT(IN) :: field
- REAL(dp), DIMENSION(total_np,total_nj,total_nz) :: field_full
+ COMPLEX(dp), DIMENSION(local_np,local_nj,local_nz), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(total_np,total_nj,total_nz) :: field_full
CHARACTER(*), INTENT(IN) :: text
CHARACTER(LEN=50) :: dset_name
field_full = 0;
@@ -456,12 +444,9 @@ SUBROUTINE diagnose_5d
USE basic, ONLY: time, iframe5d,cstep
USE futils, ONLY: append, getatt, attach, putarrnd, putarr
USE fields, ONLY: moments
- USE array, ONLY: Sapj
- USE grid, ONLY:ips, ipe, ijs, ije, &
- total_np, total_nj, total_nky, total_nkx, total_nz, &
+ USE grid, ONLY:total_np, total_nj, total_nky, total_nkx, total_nz, &
local_np, local_nj, local_nky, local_nkx, local_nz, &
- ngp, ngj, ngz,&
- ikxs,ikxe,ikys,ikye,izs,ize
+ ngp, ngj, ngz
USE time_integration, ONLY: updatetlevel
USE diagnostics_par
USE prec_const, ONLY: dp
@@ -474,12 +459,7 @@ SUBROUTINE diagnose_5d
CALL attach(fidres,"/data/var5d/" , "frames", iframe5d)
IF (write_Napj) THEN
- CALL write_field5d(moments(1,1+ngp/2:local_np+ngp/2,1+ngj/2:local_nj+ngj/2,&
- :,:,1+ngz/2:local_nz+ngz/2,updatetlevel), 'moments')
- ENDIF
-
- IF (write_Sapj) THEN
- CALL write_field5d(Sapj(1,:,:,:,:,:), 'Sapj')
+ CALL write_field5d(moments, 'moments')
ENDIF
CONTAINS
@@ -490,19 +470,22 @@ SUBROUTINE diagnose_5d
USE parallel, ONLY: gather_pjxyz, num_procs
USE prec_const, ONLY: dp
IMPLICIT NONE
- COMPLEX(dp), DIMENSION(local_np,local_nj,local_nky,local_nkx,local_nz), INTENT(IN) :: field
+ COMPLEX(dp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN) :: field
CHARACTER(*), INTENT(IN) :: text
+ COMPLEX(dp), DIMENSION(local_np,local_nj,local_nky,local_nkx,local_nz) :: field_sub
COMPLEX(dp), DIMENSION(total_np,total_nj,total_nky,total_nkx,total_nz) :: field_full
CHARACTER(LEN=50) :: dset_name
+ field_sub = field(1,(1+ngp/2):(local_np+ngp/2),(1+ngj/2):(local_nj+ngj/2),&
+ 1:local_nky,1:local_nkx, (1+ngz/2):(local_nz+ngz/2),updatetlevel)
field_full = 0;
WRITE(dset_name, "(A, '/', A, '/', i6.6)") "/data/var5d", TRIM(text), iframe5d
IF (num_procs .EQ. 1) THEN
- CALL putarr(fidres, dset_name, field(ips:ipe,ijs:ije,ikys:ikye,ikxs:ikxe,izs:ize), ionode=0)
+ CALL putarr(fidres, dset_name, field_sub, ionode=0)
ELSEIF(GATHERV_OUTPUT) THEN ! output using one node (gatherv)
- CALL gather_pjxyz(field,field_full,local_np,total_np,total_nj,local_nky,total_nky,total_nkx,local_nz,total_nz)
- CALL putarr(fidres, dset_name, field_full(1:total_np,1:total_nj,1:total_nky,1:total_nkx,1:total_nz), ionode=0)
+ CALL gather_pjxyz(field_sub,field_full,local_np,total_np,total_nj,local_nky,total_nky,total_nkx,local_nz,total_nz)
+ CALL putarr(fidres, dset_name, field_full, ionode=0)
ELSE
- CALL putarrnd(fidres, dset_name, field(ips:ipe,ijs:ije,ikys:ikye,ikxs:ikxe,izs:ize), (/1,3,5/))
+ CALL putarrnd(fidres, dset_name, field_sub, (/1,3,5/))
ENDIF
CALL attach(fidres, dset_name, 'cstep', cstep)
CALL attach(fidres, dset_name, 'time', time)
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index a344e079..75a2a524 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -163,6 +163,12 @@ CONTAINS
CALL set_kygrid(LINEARITY,N_HD)
CALL set_kxgrid(shear,Npol,LINEARITY,N_HD)
CALL set_zgrid (Npol)
+
+ print*, 'p:',parray
+ print*, 'j:',jarray
+ print*, 'ky:',kyarray
+ print*, 'kx:',kxarray
+ print*, 'z:',zarray
END SUBROUTINE set_grids
SUBROUTINE init_1Dgrid_distr
@@ -237,8 +243,8 @@ CONTAINS
local_pmin = parray(1+ngp/2)
! Fill the ghosts
DO ig = 1,ngp/2
- parray(ig) = local_pmin-ngp/2+(ig-1)
- parray(local_np+ig) = local_pmax+ig
+ parray(ig) = local_pmin-ngp/2+(ig-1)
+ parray(local_np+ngp/2+ig) = local_pmax+ig
ENDDO
IF(CONTAINSp0) SOLVE_POISSON = .TRUE.
IF(CONTAINSp1) SOLVE_AMPERE = .TRUE.
@@ -276,12 +282,12 @@ CONTAINS
DO ij = 1+ngj/2,local_nj+ngj/2
jarray(ij) = ij-1-ngj/2+local_nj_offset
END DO
- local_jmax = jarray(local_np+ngj/2)
+ local_jmax = jarray(local_nj+ngj/2)
local_jmin = jarray(1+ngj/2)
! Fill the ghosts
DO ig = 1,ngj/2
- jarray(ig) = local_jmin-ngj/2+(ig-1)
- jarray(local_nj+ig) = local_jmax+ig
+ jarray(ig) = local_jmin-ngj/2+(ig-1)
+ jarray(local_nj+ngj/2+ig) = local_jmax+ig
ENDDO
! Precomputations
jmax_dp = real(jmax,dp)
@@ -511,11 +517,13 @@ CONTAINS
CALL speak('--2 staggered z grids--')
grid_shift = deltaz/2._dp
! indices for even p and odd p grids (used in kernel, jacobian, gij etc.)
- ieven = 1; iodd = 2
+ ieven = 1
+ iodd = 2
Nzgrid = 2
ELSE
grid_shift = 0._dp
- ieven = 1; iodd = 1
+ ieven = 1
+ iodd = 1
Nzgrid = 1
ENDIF
! Build the full grids on process 0 to diagnose it without comm
diff --git a/src/inital.F90 b/src/inital.F90
index 11f6433f..f3bbe1b1 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -172,14 +172,14 @@ SUBROUTINE init_gyrodens
! Seed random number generator
iseedarr(:)=iseed
CALL RANDOM_SEED(PUT=iseedarr+my_id)
-
+ moments = 0._dp
!**** Broad noise initialization *******************************************
DO ia=1,local_na
- DO ip=1,local_np+ngp
- DO ij=1,local_nj+ngj
+ DO ip=1+ngp/2,local_np+ngp/2
+ DO ij=1+ngj/2,local_nj+ngj/2
DO ikx=1,total_nkx
DO iky=1,local_nky
- DO iz=1,local_nz+ngz
+ DO iz=1+ngz/2,local_nz+ngz/2
CALL RANDOM_NUMBER(noise)
IF ( (parray(ip) .EQ. 1) .AND. (jarray(ij) .EQ. 1) ) THEN
moments(ia,ip,ij,iky,ikx,iz,:) = (init_background + init_noiselvl*(noise-0.5_dp))
@@ -196,6 +196,8 @@ SUBROUTINE init_gyrodens
ENDIF
END DO
END DO
+ print*, SUM(moments)
+ stop
! Putting to zero modes that are not in the 2/3 Orszag rule
IF (LINEARITY .NE. 'linear') THEN
DO ikx=1,total_nkx
diff --git a/src/memory.F90 b/src/memory.F90
index c9d856d1..97cf7191 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -20,7 +20,6 @@ SUBROUTINE memory
! CALL allocate_array(HF_phi_correction_operator, 1,local_nky, 1,local_nkx, 1,local_Nz)
!Moments related arrays
- CALL allocate_array( Na00, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
CALL allocate_array( dens, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
CALL allocate_array( upar, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
CALL allocate_array( uper, 1,local_Na, 1,local_nky, 1,local_nkx, 1,local_Nz)
diff --git a/src/parallel_mod.F90 b/src/parallel_mod.F90
index 6e4fc9c6..d90918ef 100644
--- a/src/parallel_mod.F90
+++ b/src/parallel_mod.F90
@@ -238,12 +238,12 @@ CONTAINS
SUBROUTINE gather_pjz(field_loc,field_tot,np_loc,np_tot,nj_tot,nz_loc,nz_tot)
IMPLICIT NONE
INTEGER, INTENT(IN) :: np_loc,np_tot, nj_tot, nz_loc,nz_tot
- REAL(dp), DIMENSION(np_loc,nj_tot,nz_loc), INTENT(IN) :: field_loc
- REAL(dp), DIMENSION(np_tot,nj_tot,nz_tot), INTENT(OUT) :: field_tot
- REAL(dp), DIMENSION(np_tot,nz_loc) :: buffer_pt_zl !full p, local z
- REAL(dp), DIMENSION(np_tot,nz_tot) :: buffer_pt_zt !full p, full z
- REAL(dp), DIMENSION(np_loc) :: buffer_pl_zc !local p, constant z
- REAL(dp), DIMENSION(np_tot) :: buffer_pt_zc !full p, constant z
+ COMPLEX(dp), DIMENSION(np_loc,nj_tot,nz_loc), INTENT(IN) :: field_loc
+ COMPLEX(dp), DIMENSION(np_tot,nj_tot,nz_tot), INTENT(OUT) :: field_tot
+ COMPLEX(dp), DIMENSION(np_tot,nz_loc) :: buffer_pt_zl !full p, local z
+ COMPLEX(dp), DIMENSION(np_tot,nz_tot) :: buffer_pt_zt !full p, full z
+ COMPLEX(dp), DIMENSION(np_loc) :: buffer_pl_zc !local p, constant z
+ COMPLEX(dp), DIMENSION(np_tot) :: buffer_pt_zc !full p, constant z
INTEGER :: snd_p, snd_z, root_p, root_z, root_ky, ij, iz
snd_p = np_loc ! Number of points to send along p (per z)
@@ -255,16 +255,16 @@ CONTAINS
DO iz = 1,nz_loc
! fill a buffer to contain a slice of data at constant j and z
buffer_pl_zc(1:np_loc) = field_loc(1:np_loc,ij,iz)
- CALL MPI_GATHERV(buffer_pl_zc, snd_p, MPI_DOUBLE_PRECISION, &
- buffer_pt_zc, rcv_p, dsp_p, MPI_DOUBLE_PRECISION, &
+ CALL MPI_GATHERV(buffer_pl_zc, snd_p, MPI_DOUBLE_COMPLEX, &
+ buffer_pt_zc, rcv_p, dsp_p, MPI_DOUBLE_COMPLEX, &
root_p, comm_p, ierr)
buffer_pt_zl(1:np_tot,iz) = buffer_pt_zc(1:np_tot)
ENDDO
! send the full line on y contained by root_kyas
IF(rank_p .EQ. 0) THEN
- CALL MPI_GATHERV(buffer_pt_zl, snd_z, MPI_DOUBLE_PRECISION, &
- buffer_pt_zt, rcv_zp, dsp_zp, MPI_DOUBLE_PRECISION, &
+ CALL MPI_GATHERV(buffer_pt_zl, snd_z, MPI_DOUBLE_COMPLEX, &
+ buffer_pt_zt, rcv_zp, dsp_zp, MPI_DOUBLE_COMPLEX, &
root_z, comm_z, ierr)
ENDIF
! ID 0 (the one who ouptut) rebuild the whole array
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index ecfe01c7..001993a6 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -47,7 +47,7 @@ SUBROUTINE compute_radial_transport
IMPLICIT NONE
COMPLEX(dp) :: pflux_local, gflux_local, integral
REAL(dp) :: buffer(2)
- INTEGER :: i_, root, iky, ikx, ia, iz, in, iodd, ieven
+ INTEGER :: i_, root, iky, ikx, ia, iz, in
COMPLEX(dp), DIMENSION(local_nz+Ngz) :: integrant
DO ia = 1,local_na
pflux_local = 0._dp ! particle flux
@@ -330,13 +330,15 @@ SUBROUTINE compute_Napjz_spectrum
! z-moment spectrum
! build local sum
local_sum = 0._dp
- DO ikx = 1,local_nkx
- DO iky = 1,local_nky
- DO ij = 1,local_nj
- DO ip = 1,local_np
- local_sum(ip,ij,iz) = local_sum(ip,ij,iz) + &
- (moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel) &
- * CONJG(moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel)))
+ DO iz = 1,local_nz
+ DO ikx = 1,local_nkx
+ DO iky = 1,local_nky
+ DO ij = 1,local_nj
+ DO ip = 1,local_np
+ local_sum(ip,ij,iz) = local_sum(ip,ij,iz) + &
+ (moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel) &
+ * CONJG(moments(ia,ip+Ngp/2,ij+Ngj/2,iky,ikx,iz+Ngz/2,updatetlevel)))
+ ENDDO
ENDDO
ENDDO
ENDDO
diff --git a/testcases/smallest_problem/fort.90 b/testcases/smallest_problem/fort.90
index 7c491e5d..40604ef7 100644
--- a/testcases/smallest_problem/fort.90
+++ b/testcases/smallest_problem/fort.90
@@ -1,5 +1,5 @@
&BASIC
- nrun = 500
+ nrun = 1
dt = 0.01
tmax = 5
maxruntime = 356400
@@ -40,9 +40,8 @@
write_gamma = .t.
write_hf = .t.
write_phi = .t.
- write_Na00 = .f.
+ write_Na00 = .t.
write_Napj = .t.
- write_Sapj = .f.
write_dens = .t.
write_fvel = .t.
write_temp = .t.
@@ -53,13 +52,13 @@
NL_CLOS = -1
LINEARITY = 'linear'
Na = 2 ! number of species
- mu_x = 0.0
- mu_y = 0.0
+ mu_x = 0
+ mu_y = 0
N_HD = 4
- mu_z = 0.1
+ mu_z = 0
mu_p = 0
mu_j = 0
- nu = 1
+ nu = 0
beta = 0
ADIAB_E = .f.
tau_e = 1.0
@@ -70,8 +69,8 @@
tau_ = 1.0
sigma_= 1.0
q_ = 1.0
- k_N_ = 2.22
- k_T_ = 6.96
+ k_N_ = 0.0!2.22
+ k_T_ = 0.0!6.96
/
&SPECIES
! electrons
@@ -79,22 +78,22 @@
tau_ = 1.0
sigma_= 0.023338
q_ = 1.0
- k_N_ = 2.22
- k_T_ = 6.96
+ k_N_ = 0.0!2.22
+ k_T_ = 0.0!6.96
/
&COLLISION_PAR
collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
- GK_CO = .f.
+ GK_CO = .f.
INTERSPECIES = .true.
- !mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
+ !mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
&INITIAL_CON
- INIT_OPT = 'phi'
- ACT_ON_MODES = 'donothing'
- init_background = 0
- init_noiselvl = 0.005
- iseed = 42
+ INIT_OPT = 'mom00'
+ ACT_ON_MODES = 'donothing'
+ init_background = 1.0
+ init_noiselvl = 0.0
+ iseed = 42
/
&TIME_INTEGRATION_PAR
numerical_scheme = 'RK4'
diff --git a/testcases/smallest_problem/fort_00.90 b/testcases/smallest_problem/fort_00.90
index 622af92a..9ea911de 100644
--- a/testcases/smallest_problem/fort_00.90
+++ b/testcases/smallest_problem/fort_00.90
@@ -33,9 +33,8 @@
write_gamma = .t.
write_hf = .t.
write_phi = .t.
- write_Na00 = .f.
+ write_Na00 = .t.
write_Napj = .t.
- write_Sapj = .f.
write_dens = .t.
write_temp = .t.
job2load = -1
@@ -49,20 +48,20 @@
mu_x = 0.0
mu_y = 0.0
N_HD = 4
- mu_z = 0.1
+ mu_z = 0!0.1
mu_p = 0
mu_j = 0
- nu = 1
+ nu = 0!1.0
tau_e = 1
tau_i = 1
sigma_e = 0.023338
sigma_i = 1
q_e = -1
q_i = 1
- K_Ne = 1!6.96
- K_Te = 1!2.22
- K_Ni = 1!6.96
- K_Ti = 1!2.22
+ K_Ne = 0!6.96
+ K_Te = 0!2.22
+ K_Ni = 0!6.96
+ K_Ti = 0!2.22
beta = 0
/
&COLLISION_PAR
@@ -72,7 +71,7 @@
!mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
&INITIAL_CON
- INIT_OPT = 'phi'
+ INIT_OPT = 'mom00'
ACT_ON_MODES = 'donothing'
init_background = 1.0
init_noiselvl = 0.0
diff --git a/testcases/smallest_problem/gyacomo_1 b/testcases/smallest_problem/gyacomo_1
new file mode 120000
index 00000000..d087f907
--- /dev/null
+++ b/testcases/smallest_problem/gyacomo_1
@@ -0,0 +1 @@
+../../bin/gyacomo_1
\ No newline at end of file
diff --git a/testcases/smallest_problem/gyacomo_1_debug b/testcases/smallest_problem/gyacomo_1_debug
new file mode 120000
index 00000000..f5e7015f
--- /dev/null
+++ b/testcases/smallest_problem/gyacomo_1_debug
@@ -0,0 +1 @@
+../../bin/gyacomo_1_debug
\ No newline at end of file
diff --git a/testcases/smallest_problem/gyacomo_alpha b/testcases/smallest_problem/gyacomo_alpha
deleted file mode 120000
index 0663323a..00000000
--- a/testcases/smallest_problem/gyacomo_alpha
+++ /dev/null
@@ -1 +0,0 @@
-../../bin/gyacomo_alpha
\ No newline at end of file
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 58375962..25115ab7 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -56,7 +56,7 @@ PARTITION = '/misc/gyacomo_outputs/';
% resdir = 'paper_2_nonlinear/kT_5.3/9x5x128x64x24';
% resdir = 'paper_2_nonlinear/kT_5.3/9x5x128x64x64';
% resdir = 'paper_2_nonlinear/kT_5.3/11x6x128x64x24';
-resdir = 'paper_2_nonlinear/kT_5.3/11x6x128x64x64';
+% resdir = 'paper_2_nonlinear/kT_5.3/11x6x128x64x64';
%% Old stuff
% resdir = 'CBC/kT_4.5_128x64x16x13x7/';
--
GitLab