From 88c3c967c89b1ec9ab33f82e760d9175f1d697ba Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 20 Mar 2023 13:30:52 +0100
Subject: [PATCH] benchmarked with nonlinear Zpinch case
---
matlab/compile_results.m | 36 +++++--
matlab/extract_fig_data.m | 2 +-
src/fourier_mod.F90 | 24 ++---
src/grid_mod.F90 | 24 +++--
src/inital.F90 | 14 +--
src/nonlinear_mod.F90 | 42 ++++----
testcases/smallest_problem/fort.90 | 8 +-
testcases/smallest_problem/fort_00.90 | 8 +-
testcases/zpinch_example/fort.90 | 101 +++++++++++++++++++
testcases/zpinch_example/gyacomo_1 | 1 +
testcases/zpinch_example/gyacomo_1_debug | 1 +
testcases/zpinch_example_new/fort.90 | 101 +++++++++++++++++++
testcases/zpinch_example_new/fort_00.90 | 83 +++++++++++++++
testcases/zpinch_example_new/gyacomo_1 | 1 +
testcases/zpinch_example_new/gyacomo_1_debug | 1 +
testcases/zpinch_example_new/gyacomo_debug | 1 +
wk/analysis_gene.m | 4 +-
wk/analysis_gyacomo.m | 5 +-
wk/header_3D_results.m | 4 +-
19 files changed, 392 insertions(+), 69 deletions(-)
create mode 100644 testcases/zpinch_example/fort.90
create mode 120000 testcases/zpinch_example/gyacomo_1
create mode 120000 testcases/zpinch_example/gyacomo_1_debug
create mode 100644 testcases/zpinch_example_new/fort.90
create mode 100644 testcases/zpinch_example_new/fort_00.90
create mode 120000 testcases/zpinch_example_new/gyacomo_1
create mode 120000 testcases/zpinch_example_new/gyacomo_1_debug
create mode 120000 testcases/zpinch_example_new/gyacomo_debug
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index 27c14463..31fa7062 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -293,16 +293,40 @@ else
DATA.scale = 1;%(1/Nx/Ny)^2;
% Fields
DATA.GGAMMA_RI = GGAMMAI_; DATA.PGAMMA_RI = PGAMMAI_; DATA.HFLUX_X = HFLUXI_;
- DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.Nipjz = Nipjz_;
- DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_;
+ if W_NAPJ
+ DATA.Nipj = zeros(Pi_new,Ji_new,Nky,Nkx,Nz,numel(Ts5D_)); DATA.Nipj(:,:,:,:,1:Nz,:) = Nipj_;
+ end
+ if W_NA00
+ DATA.Ni00 = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.Ni00(:,:,1:Nz,:) = Ni00_;
+ DATA.Nipjz = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.Nipjz(:,:,1:Nz,:) = Nipjz_;
+ end
+ if W_DENS
+ DATA.DENS_I = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.DENS_I(:,:,1:Nz,:) = DENS_I_;
+ end
+ if W_TEMP
+ DATA.TEMP_I = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.TEMP_I(:,:,1:Nz,:) = TEMP_I_;
+ end
if(KIN_E)
- DATA.Nepj = Nepj_; DATA.Ne00 = Ne00_; DATA.Nepjz = Nepjz_;
- DATA.DENS_E = DENS_E_; DATA.TEMP_E = TEMP_E_;
+ if W_NAPJ
+ DATA.Nepj = zeros(Pe_new,Je_new,Nky,Nkx,Nz,numel(Ts5D_)); DATA.Nepj(:,:,:,:,1:Nz,:) = Nepj_;
+ end
+ if W_NA00
+ DATA.Ne00 = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.Ne00(:,:,1:Nz,:) = Ne00_;
+ DATA.Nepjz = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.Nepjz(:,:,1:Nz,:) = Nepjz_;
+ end
+ if W_DENS
+ DATA.DENS_E = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.DENS_E(:,:,1:Nz,:) = DENS_E_;
+ end
+ if W_TEMP
+ DATA.TEMP_E = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.TEMP_E(:,:,1:Nz,:) = TEMP_E_;
+ end
DATA.HFLUX_XE = HFLUXE_;
end
DATA.Ts5D = Ts5D_; DATA.Ts3D = Ts3D_; DATA.Ts0D = Ts0D_;
- DATA.PHI = PHI_;
- DATA.PSI = PSI_;
+ DATA.PHI = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.PHI(:,:,1:Nz,:) = PHI_;
+ if BETA>0
+ DATA.PSI = zeros(Nky,Nkx,Nz,numel(Ts3D_)); DATA.PSI(:,:,1:Nz,:) = PSI_;
+ end
DATA.KIN_E=KIN_E;
% grids
DATA.Pe = Pe; DATA.Pi = Pi;
diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index 8c8f0d92..825aff05 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -4,7 +4,7 @@
% tw = [3000 4000];
% tw = [4000 4500];
% tw = [4500 5000];
-tw = [100 180];
+tw = [00 1000];
fig = gcf;
axObjs = fig.Children;
diff --git a/src/fourier_mod.F90 b/src/fourier_mod.F90
index 3a23febf..f6e226c2 100644
--- a/src/fourier_mod.F90
+++ b/src/fourier_mod.F90
@@ -68,22 +68,22 @@ MODULE fourier
!!! Compute the poisson bracket of [F,G] to real space
! - Compute the convolution using the convolution theorem
- SUBROUTINE poisson_bracket_and_sum(kx_, ky_, inv_Nx, inv_Ny, AA_x, AA_y,&
+ SUBROUTINE poisson_bracket_and_sum(ky_, kx_, inv_Ny, inv_Nx, AA_y, AA_x,&
local_nky_ptr, local_nkx_ptr, F_, G_, sum_real_)
IMPLICIT NONE
- INTEGER(C_INTPTR_T), INTENT(IN) :: local_nkx_ptr,local_nky_ptr
- REAL(dp), INTENT(IN) :: inv_Nx, inv_Ny
- REAL(dp), DIMENSION(:), INTENT(IN) :: kx_, ky_, AA_x, AA_y
- COMPLEX(C_DOUBLE_COMPLEX),INTENT(IN) :: F_(:,:), G_(:,:)
- real(C_DOUBLE), pointer, INTENT(INOUT) :: sum_real_(:,:)
+ INTEGER(C_INTPTR_T), INTENT(IN) :: local_nkx_ptr,local_nky_ptr
+ REAL(dp), INTENT(IN) :: inv_Nx, inv_Ny
+ REAL(dp), DIMENSION(local_nky_ptr), INTENT(IN) :: ky_, AA_y
+ REAL(dp), DIMENSION(local_nkx_ptr), INTENT(IN) :: kx_, AA_x
+ COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(local_nky_ptr,local_nkx_ptr), &
+ INTENT(IN) :: F_(:,:), G_(:,:)
+ real(C_DOUBLE), pointer, INTENT(INOUT) :: sum_real_(:,:)
INTEGER :: ikx,iky
! First term df/dx x dg/dy
DO ikx = 1,local_nkx_ptr
DO iky = 1,local_nky_ptr
- cmpx_data_f(ikx,iky) = &
- imagu*kx_(ikx)*F_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
- cmpx_data_g(ikx,iky) = &
- imagu*ky_(iky)*G_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ cmpx_data_f(ikx,iky) = imagu*kx_(ikx)*F_(iky,ikx)*AA_x(ikx)*AA_y(iky)
+ cmpx_data_g(ikx,iky) = imagu*ky_(iky)*G_(iky,ikx)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
@@ -93,9 +93,9 @@ MODULE fourier
DO ikx = 1,local_nkx_ptr
DO iky = 1,local_nky_ptr
cmpx_data_f(ikx,iky) = &
- imagu*ky_(iky)*F_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ imagu*ky_(iky)*F_(iky,ikx)*AA_x(ikx)*AA_y(iky)
cmpx_data_g(ikx,iky) = &
- imagu*kx_(ikx)*G_(iky,ikx)*AA_x(ikx)*AA_y(iky) !Anti aliasing filter
+ imagu*kx_(ikx)*G_(iky,ikx)*AA_x(ikx)*AA_y(iky)
ENDDO
ENDDO
call fftw_mpi_execute_dft_c2r(planb, cmpx_data_f, real_data_f)
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 5a1604b6..367f3ea4 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -508,6 +508,9 @@ CONTAINS
zarray_full(izs) = 0
ELSEIF(Nz .GE. 4) THEN
ngz =4
+ IF(mod(Nz,2) .NE. 0 ) THEN
+ ERROR STOP '>> ERROR << Nz must be an even number for Simpson integration rule !!!!'
+ ENDIF
ELSE
ERROR STOP '>> ERROR << Nz is not appropriate!!'
ENDIF
@@ -559,18 +562,19 @@ CONTAINS
IF(abs(local_zmax(eo) - (zmax+REAL(eo-1,dp)*grid_shift)) .LT. EPSILON(zmax)) &
contains_zmax = .TRUE.
ENDDO
- IF(mod(Nz,2) .NE. 0 ) THEN
- ERROR STOP '>> ERROR << Nz must be an even number for Simpson integration rule !!!!'
- ENDIF
! local weights for Simpson rule
ALLOCATE(zweights_SR(local_nz))
- DO iz = 1,local_nz
- IF(MODULO(iz+local_nz_offset,2) .EQ. 1) THEN ! odd iz
- zweights_SR(iz) = onethird*deltaz*2._dp
- ELSE ! even iz
- zweights_SR(iz) = onethird*deltaz*4._dp
- ENDIF
- ENDDO
+ IF(total_nz .EQ. 1) THEN
+ zweights_SR = 1._dp
+ ELSE
+ DO iz = 1,local_nz
+ IF(MODULO(iz+local_nz_offset,2) .EQ. 1) THEN ! odd iz
+ zweights_SR(iz) = onethird*deltaz*2._dp
+ ELSE ! even iz
+ zweights_SR(iz) = onethird*deltaz*4._dp
+ ENDIF
+ ENDDO
+ ENDIF
END SUBROUTINE set_zgrid
SUBROUTINE set_kparray(gxx, gxy, gyy,hatB)
diff --git a/src/inital.F90 b/src/inital.F90
index 04c4d511..9e9abed6 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -327,7 +327,7 @@ END SUBROUTINE init_phi_ppj
!******************************************************************************!
SUBROUTINE initialize_blob
USE grid, ONLY: local_na, local_np, local_nj, total_nkx, local_nky, local_nz,&
- AA_x, AA_y,&
+ AA_x, AA_y, parray, jarray,&
ngp,ngj,ngz, iky0, ieven, kxarray, kyarray, zarray
USE fields, ONLY: moments
USE prec_const, ONLY: dp
@@ -336,20 +336,22 @@ SUBROUTINE initialize_blob
USE geometry, ONLY: Jacobian, iInt_Jacobian, shear
IMPLICIT NONE
REAL(dp) ::kx, ky, z, sigma_x, sigma_y, gain
- INTEGER :: ia,iky,ikx,iz,ip,ij
+ INTEGER :: ia,iky,ikx,iz,ip,ij, p, j
sigma_y = 1.0
sigma_x = sigma_y
gain = 10.0
DO ia=1,local_na
DO iky=1,local_nky
ky = kyarray(iky)
- DO iz=1,local_nz+ngz
+ DO iz=1+ngz/2,local_nz+ngz/2
z = zarray(iz,ieven)
DO ikx=1,total_nkx
kx = kxarray(ikx) + z*shear*ky
- DO ip=1,local_np+ngp
- DO ij=1,local_nj+ngj
- IF( (iky .NE. iky0) .AND. (ip .EQ. 1) .AND. (ij .EQ. 1)) THEN
+ DO ip=1+ngp/2,local_np+ngp/2
+ p = parray(ip)
+ DO ij=1+ngj/2,local_nj+ngj/2
+ j = jarray(ij)
+ IF( (iky .NE. iky0) .AND. (p .EQ. 0) .AND. (j .EQ. 0)) THEN
moments(ia,ip,ij,iky,ikx,iz, :) = moments(ia,ip,ij,iky,ikx,iz, :) &
+ gain * exp(-((kx/sigma_x)**2+(ky/sigma_y)**2)) &
* AA_x(ikx)*AA_y(iky)* &
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 7155fc24..44f48f84 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -49,7 +49,6 @@ SUBROUTINE compute_Sapj
!! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
!! Sapj = Sum_n (ikx Kn phi)#(iky Sum_s d_njs Naps) - (iky Kn phi)#(ikx Sum_s d_njs Naps)
!! where # denotes the convolution.
-
! Execution time start
CALL cpu_time(t0_Sapj)
@@ -61,21 +60,22 @@ SUBROUTINE compute_Sapj
CASE DEFAULT
ERROR STOP '>> ERROR << Linearity not recognized '
END SELECT
-
! Execution time END
CALL cpu_time(t1_Sapj)
tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
-
END SUBROUTINE compute_Sapj
SUBROUTINE compute_nonlinear
IMPLICIT NONE
- INTEGER :: iz,ij,ip,eo,ia,ikx,iky
+ INTEGER :: iz,ij,ip,eo,ia,ikx,iky,izi,ipi,iji,ini,isi
DO iz = 1,local_nz
+ izi = iz + ngz/2
DO ij = 1,local_nj ! Loop over Laguerre moments
- j_int=jarray(ij+ngj/2)
+ iji = ij + ngj/2
+ j_int=jarray(iji)
DO ip = 1,local_np ! Loop over Hermite moments
- p_int = parray(ip+ngp/2)
+ ipi = ip + ngp/2
+ p_int = parray(ipi)
sqrt_p = SQRT(REAL(p_int,dp))
sqrt_pp1 = SQRT(REAL(p_int,dp) + 1._dp)
eo = min(nzgrid,MODULO(parray(ip),2)+1)
@@ -83,38 +83,40 @@ SUBROUTINE compute_nonlinear
IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmax)) THEN !compute for every moments except for closure 1
! Set non linear sum truncation
IF (NL_CLOS .EQ. -2) THEN
- nmax = local_nj
+ nmax = Jmax
ELSEIF (NL_CLOS .EQ. -1) THEN
- nmax = (Jmax-j_int)+1+ngj/2-local_nj_offset
+ nmax = Jmax-j_int
ELSE
- nmax = NL_CLOS+1+ngj/2-local_nj_offset
+ nmax = min(NL_CLOS,Jmax-j_int)
ENDIF
bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
- DO in = 1,nmax ! Loop over laguerre for the sum
- n_int = parray(in+ngp/2)
+ DO in = 1,nmax+1 ! Loop over laguerre for the sum
+ ini = in+ngj/2
!-----------!! ELECTROSTATIC CONTRIBUTION
! First convolution terms
- F_cmpx(:,:) = phi(:,:,iz+ngz/2) * kernel(ia,in+ngj/2,:,:,iz+ngz/2,eo)
+ F_cmpx(:,:) = phi(:,:,izi) * kernel(ia,ini,:,:,izi,eo)
! Second convolution terms
G_cmpx = 0._dp ! initialization of the sum
- smax = (n_int+j_int)+1+ngj/2-local_nj_offset
- DO is = 1, MIN(smax,local_nj) ! sum truncation on number of moments
+ smax = MIN( (in-1)+(ij-1), Jmax );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ isi = is + ngj/2
G_cmpx(:,:) = G_cmpx(:,:) + &
- dnjs(in,ij,is) * moments(ia,ip,is+ngj/2,:,:,iz+ngz/2,updatetlevel)
+ dnjs(in,ij,is) * moments(ia,ipi,isi,:,:,izi,updatetlevel)
ENDDO
! this function add its result to bracket_sum_r
CALL poisson_bracket_and_sum(kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,local_nky_ptr,local_nkx_ptr,F_cmpx,G_cmpx,bracket_sum_r)
!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
IF(EM) THEN
! First convolution terms
- F_cmpx(:,:) = -sqrt_tau_o_sigma(ia) * psi(:,:,iz+ngz/2) * kernel(ia,in+ngj/2,:,:,iz+ngz/2,eo)
+ F_cmpx(:,:) = -sqrt_tau_o_sigma(ia) * psi(:,:,izi) * kernel(ia,ini,:,:,izi,eo)
! Second convolution terms
G_cmpx = 0._dp ! initialization of the sum
- smax = (n_int+j_int)+1+ngj/2-local_nj_offset
- DO is = 1, MIN(smax,local_nj) ! sum truncation on number of moments
+ smax = MIN( (in-1)+(ij-1), Jmax );
+ DO is = 1, smax+1 ! sum truncation on number of moments
+ isi = is + ngj/2
G_cmpx(:,:) = G_cmpx(:,:) + &
- dnjs(in,ij,is) * (sqrt_pp1*moments(ia,ip+1+ngj/2,is,:,:,iz+ngz/2,updatetlevel)&
- +sqrt_p *moments(ia,ip-1+ngj/2,is,:,:,iz+ngz/2,updatetlevel))
+ dnjs(in,ij,is) * (sqrt_pp1*moments(ia,ipi+1,isi,:,:,izi,updatetlevel)&
+ +sqrt_p *moments(ia,ipi-1,isi,:,:,izi,updatetlevel))
ENDDO
! this function add its result to bracket_sum_r
CALL poisson_bracket_and_sum(kyarray,kxarray,inv_Ny,inv_Nx,AA_y,AA_x,local_nky_ptr,local_nkx_ptr,F_cmpx,G_cmpx,bracket_sum_r)
diff --git a/testcases/smallest_problem/fort.90 b/testcases/smallest_problem/fort.90
index f6a8807c..bc7dce3f 100644
--- a/testcases/smallest_problem/fort.90
+++ b/testcases/smallest_problem/fort.90
@@ -8,9 +8,9 @@
&GRID
pmax = 4
jmax = 1
- Nx = 2
+ Nx = 16
Lx = 200
- Ny = 4
+ Ny = 12
Ly = 60
Nz = 6
SG = .f.
@@ -50,7 +50,7 @@
! Collisionality
CLOS = 0
NL_CLOS = -1
- LINEARITY = 'linear'
+ LINEARITY = 'nonlinear'
Na = 2 ! number of species
mu_x = 0.2
mu_y = 0.4
@@ -90,7 +90,7 @@
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
&INITIAL_CON
- INIT_OPT = 'allmom'
+ INIT_OPT = 'blob'
ACT_ON_MODES = 'donothing'
init_background = 1.0
init_noiselvl = 0.0
diff --git a/testcases/smallest_problem/fort_00.90 b/testcases/smallest_problem/fort_00.90
index d4f127f0..d330d8e0 100644
--- a/testcases/smallest_problem/fort_00.90
+++ b/testcases/smallest_problem/fort_00.90
@@ -9,9 +9,9 @@
jmaxe = 1
pmaxi = 4
jmaxi = 1
- Nx = 2
+ Nx = 16
Lx = 200
- Ny = 4
+ Ny = 12
Ly = 60
Nz = 6
SG = .f.
@@ -43,7 +43,7 @@
! Collisionality
CLOS = 0
NL_CLOS = -1
- LINEARITY = 'linear'
+ LINEARITY = 'nonlinear'
KIN_E = .t.
mu_x = 0.2
mu_y = 0.4
@@ -72,7 +72,7 @@
mat_file = 'gk_sugama_P_20_J_10_N_150_kpm_8.0.h5'
/
&INITIAL_CON
- INIT_OPT = 'allmom'
+ INIT_OPT = 'blob'
ACT_ON_MODES = 'donothing'
init_background = 1.0
init_noiselvl = 0.0
diff --git a/testcases/zpinch_example/fort.90 b/testcases/zpinch_example/fort.90
new file mode 100644
index 00000000..5c5bd290
--- /dev/null
+++ b/testcases/zpinch_example/fort.90
@@ -0,0 +1,101 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 50
+ maxruntime = 356400
+ job2load = -1
+/
+&GRID
+ pmax = 4
+ jmax = 2
+ Nx = 128
+ Lx = 200
+ Ny = 48
+ Ly = 60
+ Nz = 1
+ SG = .f.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0.0
+ shear = 0.0
+ eps = 0.0
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'shearless'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 0.1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 1
+ dtsave_5d = 10
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ Na = 2 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.1
+ beta = 0.0
+ ADIAB_E = .f.
+ tau_e = 1.0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+&SPECIES
+ ! electrons
+ name_ = 'electrons'
+ tau_ = 1.0
+ sigma_= 0.023338
+ q_ = -1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .t.
+ INTERSPECIES = .true.
+ 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.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/zpinch_example/gyacomo_1 b/testcases/zpinch_example/gyacomo_1
new file mode 120000
index 00000000..2fd67105
--- /dev/null
+++ b/testcases/zpinch_example/gyacomo_1
@@ -0,0 +1 @@
+../../../gyacomo_1/bin/gyacomo
\ No newline at end of file
diff --git a/testcases/zpinch_example/gyacomo_1_debug b/testcases/zpinch_example/gyacomo_1_debug
new file mode 120000
index 00000000..f0c8bbfd
--- /dev/null
+++ b/testcases/zpinch_example/gyacomo_1_debug
@@ -0,0 +1 @@
+../../../gyacomo_1/bin/gyacomo_debug
\ No newline at end of file
diff --git a/testcases/zpinch_example_new/fort.90 b/testcases/zpinch_example_new/fort.90
new file mode 100644
index 00000000..5c5bd290
--- /dev/null
+++ b/testcases/zpinch_example_new/fort.90
@@ -0,0 +1,101 @@
+&BASIC
+ nrun = 99999999
+ dt = 0.01
+ tmax = 50
+ maxruntime = 356400
+ job2load = -1
+/
+&GRID
+ pmax = 4
+ jmax = 2
+ Nx = 128
+ Lx = 200
+ Ny = 48
+ Ly = 60
+ Nz = 1
+ SG = .f.
+ Nexc = 1
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0.0
+ shear = 0.0
+ eps = 0.0
+ kappa = 1.0
+ s_kappa= 0.0
+ delta = 0.0
+ s_delta= 0.0
+ zeta = 0.0
+ s_zeta = 0.0
+ parallel_bc = 'shearless'
+ shift_y= 0.0
+/
+&OUTPUT_PAR
+ dtsave_0d = 0.1
+ dtsave_1d = -1
+ dtsave_2d = -1
+ dtsave_3d = 1
+ dtsave_5d = 10
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .t.
+ write_Napj = .t.
+ write_dens = .t.
+ write_fvel = .t.
+ write_temp = .t.
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ Na = 2 ! number of species
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ HYP_V = 'hypcoll'
+ mu_p = 0.0
+ mu_j = 0.0
+ nu = 0.1
+ beta = 0.0
+ ADIAB_E = .f.
+ tau_e = 1.0
+/
+&SPECIES
+ ! ions
+ name_ = 'ions'
+ tau_ = 1.0
+ sigma_= 1.0
+ q_ = 1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+&SPECIES
+ ! electrons
+ name_ = 'electrons'
+ tau_ = 1.0
+ sigma_= 0.023338
+ q_ = -1.0
+ k_N_ = 2.0
+ k_T_ = 0.4
+/
+
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .t.
+ INTERSPECIES = .true.
+ 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.0
+ init_noiselvl = 0.005
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/zpinch_example_new/fort_00.90 b/testcases/zpinch_example_new/fort_00.90
new file mode 100644
index 00000000..5864634f
--- /dev/null
+++ b/testcases/zpinch_example_new/fort_00.90
@@ -0,0 +1,83 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.01
+ tmax = 50
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 4
+ jmaxe = 2
+ pmaxi = 4
+ jmaxi = 2
+ Nx = 128
+ Lx = 200
+ Ny = 48
+ Ly = 60
+ Nz = 1
+ SG = .f.
+/
+&GEOMETRY
+ geom = 'Z-pinch'
+ q0 = 0
+ shear = 0
+ eps = 0
+ parallel_bc = 'shearless'
+/
+&OUTPUT_PAR
+ nsave_0d = 10
+ nsave_1d = -1
+ nsave_2d = -1
+ nsave_3d = 100
+ nsave_5d = 1000
+ write_doubleprecision = .f.
+ write_gamma = .t.
+ write_hf = .t.
+ write_phi = .t.
+ write_Na00 = .f.
+ write_Napj = .t.
+ write_Sapj = .f.
+ write_dens = .t.
+ write_temp = .t.
+ job2load = -1
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .t.
+ mu_x = 1.0
+ mu_y = 1.0
+ N_HD = 4
+ mu_z = 0.0
+ mu_p = 0
+ mu_j = 0
+ nu = 0.1
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ne = 2.0
+ K_Te = 0.4
+ K_Ni = 2.0
+ K_Ti = 0.4
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG' !DG/SG/PA/LD (dougherty, sugama, pitch angle, landau)
+ GK_CO = .true.
+ INTERSPECIES = .true.
+ !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
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/zpinch_example_new/gyacomo_1 b/testcases/zpinch_example_new/gyacomo_1
new file mode 120000
index 00000000..2fd67105
--- /dev/null
+++ b/testcases/zpinch_example_new/gyacomo_1
@@ -0,0 +1 @@
+../../../gyacomo_1/bin/gyacomo
\ No newline at end of file
diff --git a/testcases/zpinch_example_new/gyacomo_1_debug b/testcases/zpinch_example_new/gyacomo_1_debug
new file mode 120000
index 00000000..f0c8bbfd
--- /dev/null
+++ b/testcases/zpinch_example_new/gyacomo_1_debug
@@ -0,0 +1 @@
+../../../gyacomo_1/bin/gyacomo_debug
\ No newline at end of file
diff --git a/testcases/zpinch_example_new/gyacomo_debug b/testcases/zpinch_example_new/gyacomo_debug
new file mode 120000
index 00000000..363e139a
--- /dev/null
+++ b/testcases/zpinch_example_new/gyacomo_debug
@@ -0,0 +1 @@
+../../bin/gyacomo_debug
\ No newline at end of file
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index 9ab388ea..b3aa355a 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -47,14 +47,14 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
%Paper 2
% folder = '/misc/gene_results/CBC/KT_6.96_64x32x32x24x12_Nexc_5/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x8x4_Nexc_5_00/';
-folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
+% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x16x8_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_6.96_128x64x24x32x16_Nexc_5_01/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_00/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x32x16_Nexc_5_01/';
% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x16x8_Nexc_5/';
-% folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
+folder = '/misc/gene_results/CBC/new_sim/KT_5.3_128x64x24x8x4_Nexc_5/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x8x4_Nexc_5_smallvbox/';
% folder = '/misc/gene_results/CBC/new_sim/KT_6.96_128x64x24x16x8_Nexc_5_largexbox/';
% folder = '/misc/gene_results/CBC/KT_5.3_128x64x24x16x8_Muz_0.02/';
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index 5ba7a9f1..ce908eef 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -1,7 +1,8 @@
%% UNCOMMENT FOR TUTORIAL
gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
-% resdir = '.'; %Name of the directory where the results are located
-% JOBNUMMIN = 00; JOBNUMMAX = 10;
+resdir = '../testcases/zpinch_example/'; %Name of the directory where the results are located
+PARTITION ='';
+JOBNUMMIN = 00; JOBNUMMAX = 10;
%%
addpath(genpath([gyacomodir,'matlab'])) % ... add
addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 0f8bc4ae..b2c5664e 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -50,7 +50,7 @@ PARTITION = '/misc/gyacomo_outputs/';
% resdir = 'paper_2_nonlinear/kT_5.3/5x3x128x64x24';
% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_sp';
% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_dp';
-% resdir = 'paper_2_nonlinear/kT_5.3/7x4x192x96x32_dp';
+resdir = 'paper_2_nonlinear/kT_5.3/7x4x192x96x32_dp';
% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x64';
% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_MUxy_0';
% resdir = 'paper_2_nonlinear/kT_5.3/7x4x128x64x24_NL_-1';
@@ -63,7 +63,7 @@ PARTITION = '/misc/gyacomo_outputs/';
% resdir = 'paper_2_nonlinear/kT_5.3/11x6x128x64x24';
% resdir = 'paper_2_nonlinear/kT_5.3/11x6x128x64x24_dp';
% resdir = 'paper_2_nonlinear/kT_5.3/11x6x128x64x64';
-resdir = 'paper_2_nonlinear/kT_5.3/21x11x128x64x24_dp';
+% resdir = 'paper_2_nonlinear/kT_5.3/21x11x128x64x24_dp';
%% nu scans kT=5.3
% resdir = 'paper_2_nonlinear/nu_scan_kT_5.3/DGGK_7x4x128x64x24_dp';
--
GitLab