From 0597cb3c91dc70a5db05ed679bf80a1e56ffc299 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Mon, 17 Oct 2022 17:29:04 +0200
Subject: [PATCH] save
---
matlab/plot/plot_cosol_mat.m | 27 +-
matlab/plot/plot_metric.m | 2 +-
matlab/plot/statistical_transport_averaging.m | 9 +-
matlab/setup.m | 5 +-
src/diagnose.F90 | 2 +
src/geometry_mod.F90 | 250 ++++--------------
src/miller_mod.F90 | 34 ++-
src/model_mod.F90 | 4 +-
src/moments_eq_rhs_mod.F90 | 12 +-
src/nonlinear_mod.F90 | 4 +-
src/processing_mod.F90 | 4 +-
testcases/miller_example/bug/fort_00.90 | 86 ++++++
testcases/miller_example/bug/fort_01.90 | 86 ++++++
testcases/miller_example/fort_00.90 | Bin 1392 -> 1390 bytes
testcases/miller_example_w_salpha/fort_00.90 | 86 ++++++
testcases/miller_example_w_salpha/fort_01.90 | 86 ++++++
wk/Zpinch_coll_scan_kN_1.7.m | 41 ++-
wk/analysis_gene.m | 3 +-
wk/analysis_gyacomo.m | 26 +-
wk/header_2DZP_results.m | 5 +-
wk/header_3D_results.m | 94 +++----
wk/lin_EPY.m | 178 +++++++++++++
wk/lin_ITG.m | 16 +-
23 files changed, 741 insertions(+), 319 deletions(-)
create mode 100644 testcases/miller_example/bug/fort_00.90
create mode 100644 testcases/miller_example/bug/fort_01.90
create mode 100644 testcases/miller_example_w_salpha/fort_00.90
create mode 100644 testcases/miller_example_w_salpha/fort_01.90
create mode 100644 wk/lin_EPY.m
diff --git a/matlab/plot/plot_cosol_mat.m b/matlab/plot/plot_cosol_mat.m
index f10777e3..686686f8 100644
--- a/matlab/plot/plot_cosol_mat.m
+++ b/matlab/plot/plot_cosol_mat.m
@@ -65,11 +65,11 @@ subplot(224)
%% FCGK
P_ = 4; J_ = 2;
-% mat_file_name = '/home/ahoffman/cosolver/gk.coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
mat_file_name = '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5';
+% mat_file_name = '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5';
% mat_file_name = '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5';
-kp = 2.0;
+kp = 0.0;
kp_a = h5read(mat_file_name,'/coordkperp');
[~,matidx] = min(abs(kp_a-kp));
matidx = sprintf('%5.5i',matidx);
@@ -99,10 +99,11 @@ if 0
mfns = {...
'/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
% '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_4.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
-% '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',...
- '/home/ahoffman/gyacomo/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_30.h5',...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_4.h5',...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5',...
+% '/home/ahoffman/gyacomo/iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5',...
+ '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5',...
% '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
'/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
% '/home/ahoffman/gyacomo/iCa/gk.hacked_sugama_P_10_J_5_N_150_kpm_8.0.h5',...
@@ -114,7 +115,7 @@ CONAME_A = {...
% 'FC 10 5 NFLR 4',...
% 'FC 10 5 NFLR 12',...
% 'FC 10 5 NFLR 12 k<2', ...
- 'FC 10 5 NFLR 30', ...
+ 'LD 6 3 NFLR 20', ...
% 'FC 4 2 NFLR 6',...
'FC 4 2 NFLR 12', ...
% 'Hacked SG A',...
@@ -151,13 +152,13 @@ end
if 0
%%
kperp= 1.5;
-mfns = {'/home/ahoffman/HeLaZ/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
- '/home/ahoffman/HeLaZ/iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12.h5',...
+mfns = {'/home/ahoffman/gyacomo/iCa/gk_sugama_P_20_J_10_N_150_kpm_8.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_6_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5',...
+ '/home/ahoffman/gyacomo/iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5',...
};
-CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'FC 10 5 NFLR 12'};
+CONAME_A = {'SG 20 10','PA 20 10', 'FC 4 2 NFLR 6', 'FC 4 2 NFLR 12', 'LD 6 3 NFLR 12'};
grow = {};
puls = {};
for j_ = 1:numel(mfns)
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 7b95c932..d548c762 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -13,7 +13,7 @@ subplot(311)
for i = 1:5
plot(data.z, geo_arrays(1,:,i),'DisplayName',names{i}); hold on;
end
- xlim([min(data.z),max(data.z)]); legend('show'); title('MoNoLiT geometry');
+ xlim([min(data.z),max(data.z)]); legend('show'); title('GYACOMO geometry');
subplot(312)
for i = 6:10
diff --git a/matlab/plot/statistical_transport_averaging.m b/matlab/plot/statistical_transport_averaging.m
index eb284c69..3a7cc924 100644
--- a/matlab/plot/statistical_transport_averaging.m
+++ b/matlab/plot/statistical_transport_averaging.m
@@ -22,7 +22,9 @@ dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
end
time_seg = (data.Ts0D(it0:it1)-data.Ts0D(it0));
-
+
+ fig = 0;
+if options.NPLOTS > 0
fig = figure;
subplot(211)
plot(time_seg,transp_seg_avg,'-'); hold on;
@@ -40,6 +42,9 @@ dt_const = numel(unique(round(diff(data.Ts0D(it0:it1))*100)))==1;
plot(time_seg,transp_seg_avg,'-'); hold on;
xlabel('Averaging time'); ylabel('$\langle Q_x\rangle_{\tau}$');
legend(['$Q_x^\infty=$',sprintf('%2.2e',transp_seg_avg(end))])
-
+end
+Gx_infty_avg = mean(gamma);
+Gx_infty_std = std (gamma);
+disp(['G_x=',sprintf('%2.2e',Gx_infty_avg),'+-',sprintf('%2.2e',Gx_infty_std)]);
end
diff --git a/matlab/setup.m b/matlab/setup.m
index 966b6fbd..26618413 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -64,8 +64,9 @@ switch CO
COLL.mat_file = '''../../../iCa/gk_pitchangle_8_P_20_J_10_N_150_kpm_8.0.h5''';
case 'LD'
% COLL.mat_file = '''../../../iCa/gk_coulomb_NFLR_12_P_4_J_2_N_50_kpm_4.0.h5''';
-% COLL.mat_file = '''../../../iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5''';
- COLL.mat_file = '''../../../iCa/LDGK_P10_J5_dk_5e-2_km_5_NFLR_30.h5''';
+% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_12_k2trunc.h5''';
+% COLL.mat_file = '''../../../iCa/LDGKii_P10_J5_dk_5e-2_km_5_NFLR_30.h5''';
+ COLL.mat_file = '''../../../iCa/LDGK_P6_J3_dk_5e-2_km_2.5_NFLR_20.h5''';
end
COLL.coll_kcut = COLL_KCUT;
% Time integration and intialization parameters
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index bfd9e9bc..07632168 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -150,12 +150,14 @@ SUBROUTINE diagnose_full(kstep)
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/hatB_NL", hatB_NL(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradxB", gradxB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradyB", gradyB(izs:ize,0:1), (/1, 1, 1/))
CALL putarrnd(fidres, "/data/metric/gradzB", gradzB(izs:ize,0:1), (/1, 1, 1/))
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 9840cc0b..4921117b 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -87,9 +87,6 @@ CONTAINS
call eval_1D_geometry
ELSE
SELECT CASE(geom)
- CASE('circular')
- IF( my_id .eq. 0 ) WRITE(*,*) 'circular geometry'
- call eval_circular_geometry
CASE('s-alpha')
IF( my_id .eq. 0 ) WRITE(*,*) 's-alpha-B geometry'
call eval_salphaB_geometry
@@ -102,7 +99,8 @@ CONTAINS
call set_miller_parameters(kappa,s_kappa,delta,s_delta,zeta,s_zeta)
call get_miller(eps,major_R,major_Z,q0,shear,alpha_MHD,edge_opt,&
C_y,C_xy,dpdx_pm_geom,gxx,gyy,gzz,gxy,gxz,gyz,&
- gradxB,gradyB,hatB,jacobian,gradzB,hatR,hatZ,dxdR,dxdZ)
+ gradxB,gradyB,hatB,jacobian,gradzB,hatR,hatZ,dxdR,dxdZ,&
+ Ckxky,gradz_coeff)
CASE DEFAULT
STOP 'geometry not recognized!!'
END SELECT
@@ -112,18 +110,22 @@ CONTAINS
! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
! normalized to rhos_
DO eo = 0,1
- DO iky = ikys, ikye
- ky = kyarray(iky)
- DO ikx = ikxs, ikxe
- kx = kxarray(ikx)
- DO iz = izgs,izge
- kparray(iky, ikx, iz, eo) = &
- SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
- ! there is a factor 1/B from the normalization; important to match GENE
+ DO iz = izgs,izge
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx = ikxs, ikxe
+ kx = kxarray(ikx)
+ kparray(iky, ikx, iz, eo) = &
+ SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
+ ! there is a factor 1/B from the normalization; important to match GENE
ENDDO
- ENDDO
- ENDDO
+ ENDDO
+ ENDDO
ENDDO
+ ! Factor in front of the nonlinear term
+ hatB_NL(izgs:izge,0:1) = Jacobian(izgs:izge,0:1)&
+ *(gxx(izgs:izge,0:1)*gyy(izgs:izge,0:1) - gxy(izgs:izge,0:1)**2)/hatB(izgs:izge,0:1)
+
! set the mapping for parallel boundary conditions
CALL set_ikx_zBC_map
@@ -171,8 +173,7 @@ CONTAINS
Jacobian(iz,eo) = q0*hatR(iz,eo)
! Relative strengh of modulus of B
- hatB (iz,eo) = 1._dp / hatR(iz,eo)
- hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term
+ hatB(iz,eo) = 1._dp / hatR(iz,eo)
! Derivative of the magnetic field strenght
gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this
@@ -186,7 +187,7 @@ CONTAINS
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)* hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - COS(z)*ky -(shear*z-alpha_MHD*SIN(z))*SIN(z)*ky)/ hatB(iz,eo)
! Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
@@ -201,220 +202,59 @@ CONTAINS
!--------------------------------------------------------------------------------
!
-
- SUBROUTINE eval_circular_geometry
- ! evaluate circular geometry model
- ! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0
+ SUBROUTINE eval_zpinch_geometry
implicit none
- REAL(dp) :: chi, kx, ky, Gx, Gy
+ REAL(dp) :: z, kx, ky, alpha_MHD
+ alpha_MHD = 0._dp
- parity: DO eo = 0,1
- zloop: DO iz = izgs,izge
- chi = zarray(iz,eo) ! = chi
+ parity: DO eo = 0,1
+ zloop: DO iz = izgs,izge
+ z = zarray(iz,eo)
- ! metric in x,y,z
+ ! metric
gxx(iz,eo) = 1._dp
- gyy(iz,eo) = 1._dp + (shear*chi)**2 - 2._dp*eps*COS(chi) - 2._dp*shear*chi*eps*SIN(chi)
- gxy(iz,eo) = shear*chi - eps*SIN(chi)
- gxz(iz,eo) = -SIN(chi)
- gyz(iz,eo) = (1._dp - 2._dp*eps*COS(chi) - shear*chi*eps*SIN(chi))/eps
- gzz(iz,eo) = (1._dp - 2._dp*eps*COS(chi))/eps**2
- dxdR(iz,eo)= COS(chi)
- dxdZ(iz,eo)= SIN(chi)
+ gxy(iz,eo) = 0._dp
+ gxz(iz,eo) = 0._dp
+ gyy(iz,eo) = 1._dp
+ gyz(iz,eo) = 0._dp
+ gzz(iz,eo) = 1._dp
+ dxdR(iz,eo)= COS(z)
+ dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
- hatR(iz,eo) = 1._dp + eps*COS(chi)
- hatZ(iz,eo) = 1._dp + eps*SIN(chi)
+ hatR(iz,eo) = 1._dp
+ hatZ(iz,eo) = 1._dp
! toroidal coordinates
Rc (iz,eo) = hatR(iz,eo)
- phic(iz,eo) = chi
+ phic(iz,eo) = z
Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
- Jacobian(iz,eo) = q0*hatR(iz,eo)
+ Jacobian(iz,eo) = 1._dp
! Relative strengh of modulus of B
- hatB (iz,eo) = 1._dp / hatR(iz,eo)
- hatB_NL(iz,eo) = 1._dp ! Factor in front of the nonlinear term
+ hatB (iz,eo) = 1._dp
+ hatB_NL(iz,eo) = 1._dp
! Derivative of the magnetic field strenght
- gradxB(iz,eo) = -COS(chi) ! Gene put a factor hatB^2 or 1/hatR^2 in this
- gradyB(iz,eo) = 0._dp
- gradzB(iz,eo) = eps * SIN(chi) / hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this
+ gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
+ gradyB(iz,eo) = 0._dp
+ gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
! Curvature operator
- Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(chi) ! Kx
- Gy = -COS(chi) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(chi) ! Ky
- DO iky = ikys, ikye
+ DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(iky, ikx, iz,eo) = (Gx*kx + Gy*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(iky, ikx, iz,eo) = -ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
- ENDDO zloop
- ENDDO parity
-
- END SUBROUTINE eval_circular_geometry
- !
- !--------------------------------------------------------------------------------
- !
-
- SUBROUTINE eval_circular_geometry_GENE
- ! evaluate circular geometry model
- ! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0
- implicit none
- REAL(dp) :: qbar, dxdr_, dpsidr, dqdr, dqbardr, Cy, dCydr_Cy, Cxy, fac2, &
- cost, sint, dchidr, dchidt, sign_Ip, B, dBdt, dBdr, &
- g11, g22, g12, g33, dBdchi, dBdr_c !GENE variables
- REAL(dp) :: X, z, kx, ky, Gamma1, Gamma2, Gamma3, feps
-
- sign_Ip = 1._dp
- feps = SQRT(1._dp-eps**2)
- qbar = q0 * feps
- dxdr_ = 1._dp
- dpsidr = eps/qbar
- dqdr = q0/eps*shear
- dqbardr = feps*q0/eps*(shear - eps**2/feps**2)
- Cy = eps/q0 * sign_Ip
- dCydr_Cy = 0._dp
- Cxy = 1._dp/feps
- fac2 = dCydr_Cy + dqdr/q0
-
-
- parity: DO eo = 0,1
- zloop: DO iz = izgs,izge
- z = zarray(iz,eo)
-
- cost = (COS(z)-eps)/(1._dp-eps*COS(z))
- sint = feps*sin(sign_Ip*z)/(1._dp-eps*COS(z))
- dchidr = -SIN(z)/feps**2
- dchidt = sign_Ip*feps/(1._dp+eps*cost)
- B = SQRT(1._dp+(eps/qbar)**2)/(1._dp+eps*cost)
- dBdt = eps*sint*B/(1._dp+eps*cost)
-
- ! metric in r,chi,Phi
- g11 = 1._dp
- g22 = dchidr**2 + dchidt**2/eps**2
- g12 = dchidr
- g33 = 1._dp/(1._dp+eps*cost)**2
-
- ! magnetic field derivatives in
- dBdchi=dBdt/dchidt
- dBdr_c=(dBdr-g12*dBdt/dchidt)/g11
-
- ! metric in x,y,z
- gxx(iz,eo) = dxdr_**2*g11
- gyy(iz,eo) = (Cy*q0)**2 * (fac2*z)**2*g11 &
- + 2._dp*fac2*z*g12 &
- + Cy**2*g33
- gxy(iz,eo) = dxdr_*Cy*sign_Ip*q0*(fac2*z*g11 + g12)
- gxz(iz,eo) = dxdr_*g12
- gyz(iz,eo) = Cy*q0*sign_Ip*(fac2*z*g12 + g22)
- gzz(iz,eo) = g22
-
- ! Jacobian
- Jacobian(iz,eo) = Cxy*abs(q0)*(1._dp+eps*cost)**2
-
- ! Background equilibrium magnetic field
- hatB (iz,eo) = B
- hatB_NL(iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2) ! In front of the NL term
- gradxB(iz,eo) = dBdr_c/dxdr_ ! Gene put a factor hatB^2 or 1/hatR^2 in this
- gradyB(iz,eo) = 0._dp
- gradzB(iz,eo) = dBdchi! Gene put a factor hatB or 1/hatR in this
-
- ! Cylindrical coordinates derivatives
- dxdR(iz,eo) = cost
- dxdZ(iz,eo) = sint
- hatR(iz,eo) = 1._dp + eps*cost
- hatZ(iz,eo) = 1._dp + eps*sint
-
- ! toroidal coordinates
- Rc (iz,eo) = hatR(iz,eo)
- phic(iz,eo) = X
- Zc (iz,eo) = hatZ(iz,eo)
-
- Gamma1 = gxy(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyy(iz,eo)
- Gamma2 = gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo) ! Kx
- Gamma3 = gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo) ! Ky
- ! Curvature operator
- DO iky = ikys, ikye
- ky = kyarray(iky)
- DO ikx= ikxs, ikxe
- kx = kxarray(ikx)
- ! Ckxky(iky, ikx, iz,eo) = (-SIN(z)*kx - (COS(z) + (shear*z - alpha_MHD*SIN(z))* SIN(z))*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
- Ckxky(iky, ikx, iz,eo) = (-sint*kx - cost*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
- ENDDO
- ENDDO
- ! coefficient in the front of parallel derivative
- gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
-
- ENDDO zloop
- ENDDO parity
-
- END SUBROUTINE eval_circular_geometry_GENE
- !
- !--------------------------------------------------------------------------------
- !
-
- SUBROUTINE eval_zpinch_geometry
- ! evaluate s-alpha geometry model
- implicit none
- REAL(dp) :: z, kx, ky, alpha_MHD
- alpha_MHD = 0._dp
-
- parity: DO eo = 0,1
- zloop: DO iz = izgs,izge
- z = zarray(iz,eo)
-
- ! metric
- gxx(iz,eo) = 1._dp
- gxy(iz,eo) = 0._dp
- gxz(iz,eo) = 0._dp
- gyy(iz,eo) = 1._dp
- gyz(iz,eo) = 0._dp
- gzz(iz,eo) = 1._dp
- dxdR(iz,eo)= COS(z)
- dxdZ(iz,eo)= SIN(z)
-
- ! Relative strengh of radius
- hatR(iz,eo) = 1._dp
- hatZ(iz,eo) = 1._dp
-
- ! toroidal coordinates
- Rc (iz,eo) = hatR(iz,eo)
- phic(iz,eo) = z
- Zc (iz,eo) = hatZ(iz,eo)
-
- ! Jacobian
- Jacobian(iz,eo) = 1._dp
-
- ! Relative strengh of modulus of B
- hatB (iz,eo) = 1._dp
- hatB_NL(iz,eo) = 1._dp
-
- ! Derivative of the magnetic field strenght
- gradxB(iz,eo) = 0._dp ! Gene put a factor hatB^2 or 1/hatR^2 in this
- gradyB(iz,eo) = 0._dp
- gradzB(iz,eo) = 0._dp ! Gene put a factor hatB or 1/hatR in this
-
- ! Curvature operator
- DO iky = ikys, ikye
- ky = kyarray(iky)
- DO ikx= ikxs, ikxe
- kx = kxarray(ikx)
- Ckxky(iky, ikx, iz,eo) = -ky * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
- ENDDO
- ENDDO
- ! coefficient in the front of parallel derivative
- gradz_coeff(iz,eo) = 1._dp / Jacobian(iz,eo) / hatB(iz,eo)
-
- ENDDO zloop
- ENDDO parity
+ ENDDO zloop
+ ENDDO parity
END SUBROUTINE eval_zpinch_geometry
!
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index f75a075e..b33e3a57 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -78,7 +78,7 @@ CONTAINS
!>Get Miller metric, magnetic field, jacobian etc.
subroutine get_miller(trpeps,major_R,major_Z,q0,shat,amhd,edge_opt,&
C_y,C_xy,dpdx_pm_geom,gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,dBdx_,dBdy_,&
- Bfield_,jacobian_,dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_)
+ Bfield_,jacobian_,dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_,Ckxky_,gradz_coeff_)
!!!!!!!!!!!!!!!! GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(dp), INTENT(INOUT) :: trpeps ! eps in gyacomo (inverse aspect ratio)
real(dp), INTENT(INOUT) :: major_R ! major radius
@@ -93,13 +93,16 @@ CONTAINS
real(dp), dimension(izgs:izge,0:1), INTENT(INOUT) :: &
gxx_,gyy_,gzz_,gxy_,gxz_,gyz_,&
dBdx_,dBdy_,Bfield_,jacobian_,&
- dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_
+ dBdz_,R_hat_,Z_hat_,dxdR_,dxdZ_, &
+ gradz_coeff_
+ real(dp), dimension(ikys:ikye,ikxs:ikxe,izgs:izge,0:1), INTENT(INOUT) :: Ckxky_
! No parameter in gyacomo yet
integer :: ikxgs =1 ! the left ghost gpdisc%pi1gl
real(dp) :: sign_Ip_CW=1 ! current sign (only normal current)
real(dp) :: sign_Bt_CW=1 ! current sign (only normal current)
!!!!!!!!!!!!!! END GYACOMO INTERFACE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+ ! Auxiliary variables for curvature computation
+ real(dp) :: G1,G2,G3,Cx,Cy,ky,kx
integer:: np, np_s, Npol_ext, Npol_s
@@ -525,17 +528,36 @@ CONTAINS
CALL update_ghosts_z(gxx_(:,eo))
CALL update_ghosts_z(gyy_(:,eo))
CALL update_ghosts_z(gxz_(:,eo))
- CALL update_ghosts_z(dBdx_(:,eo))
- CALL update_ghosts_z(dBdy_(:,eo))
CALL update_ghosts_z(gxy_(:,eo))
CALL update_ghosts_z(gzz_(:,eo))
CALL update_ghosts_z(Bfield_(:,eo))
- CALL update_ghosts_z(jacobian_(:,eo))
+ CALL update_ghosts_z(dBdx_(:,eo))
+ CALL update_ghosts_z(dBdy_(:,eo))
CALL update_ghosts_z(dBdz_(:,eo))
+ CALL update_ghosts_z(jacobian_(:,eo))
CALL update_ghosts_z(R_hat_(:,eo))
CALL update_ghosts_z(Z_hat_(:,eo))
CALL update_ghosts_z(dxdR_(:,eo))
CALL update_ghosts_z(dxdZ_(:,eo))
+
+ ! Curvature operator (Frei et al. 2022 eq 2.15)
+ DO iz = izgs,izge
+ G1 = gxx_(iz,eo)*gyy_(iz,eo)-gxy_(iz,eo)*gxy_(iz,eo)
+ G2 = gxx_(iz,eo)*gyz_(iz,eo)-gxy_(iz,eo)*gxz_(iz,eo)
+ G3 = gxy_(iz,eo)*gyz_(iz,eo)-gyy_(iz,eo)*gxz_(iz,eo)
+ Cx = (G1*dBdy_(iz,eo) + G2*dBdz_(iz,eo))/Bfield_(iz,eo)
+ Cy = (G3*dBdz_(iz,eo) - G1*dBdx_(iz,eo))/Bfield_(iz,eo)
+
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx= ikxs, ikxe
+ kx = kxarray(ikx)
+ Ckxky_(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)/Bfield_(iz,eo)
+ ENDDO
+ ENDDO
+ ENDDO
+ ! coefficient in the front of parallel derivative
+ gradz_coeff_(iz,eo) = 1._dp / Jacobian_(iz,eo) / Bfield_(iz,eo)
enddo
contains
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 75a27687..5d282836 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -52,6 +52,7 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: qe2_taue, qi2_taui ! factor of the gammaD sum
REAL(dp), PUBLIC, PROTECTED :: q_o_sqrt_tau_sigma_e, q_o_sqrt_tau_sigma_i
REAL(dp), PUBLIC, PROTECTED :: sqrt_tau_o_sigma_e, sqrt_tau_o_sigma_i
+ REAL(dp), PUBLIC, PROTECTED :: dpdx = 0 ! radial pressure gradient
PUBLIC :: model_readinputs, model_outputinputs
CONTAINS
@@ -75,8 +76,9 @@ CONTAINS
ENDIF
IF(beta .GT. 0) THEN
- EM = .TRUE.
IF(my_id.EQ.0) print*, 'Electromagnetic effects are included'
+ EM = .TRUE.
+ dpdx = -(tau_i*(k_Ni + k_Ti) + tau_e*(k_Ne + k_Te))
ENDIF
taue_qe = tau_e/q_e ! factor of the magnetic moment coupling
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 8123fe78..948d66c0 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -33,7 +33,7 @@ SUBROUTINE moments_eq_rhs_e
USE model
USE prec_const
USE collision
- USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatR
+ USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB
USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE
@@ -121,7 +121,9 @@ SUBROUTINE moments_eq_rhs_e
!! Sum of all RHS terms
moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- -imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp&
+ -imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp&
+ ! Perpendicular pressure effects
+ -i_ky*beta*dpdx * Tperp&
! Parallel coupling (Landau Damping)
-gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
@@ -179,7 +181,7 @@ SUBROUTINE moments_eq_rhs_i
USE model
USE prec_const
USE collision
- USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatR
+ USE geometry, ONLY: gradz_coeff, gradzB, Ckxky, hatB_NL, hatB
USE calculus, ONLY : interp_z, grad_z, grad_z2
IMPLICIT NONE
@@ -269,7 +271,9 @@ SUBROUTINE moments_eq_rhs_i
!! Sum of all RHS terms
moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- -imagu*Ckxky(iky,ikx,iz,eo)*hatR(iz,eo) * Tperp &
+ -imagu*Ckxky(iky,ikx,iz,eo)*hatB(iz,eo) * Tperp &
+ ! Perpendicular pressure effects
+ -i_ky*beta*dpdx * Tperp&
! Parallel coupling (Landau damping)
-gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
diff --git a/src/nonlinear_mod.F90 b/src/nonlinear_mod.F90
index 8e9a2ce5..61547176 100644
--- a/src/nonlinear_mod.F90
+++ b/src/nonlinear_mod.F90
@@ -141,7 +141,7 @@ ENDIF
eo = MODULO(parray_i(ip),2)
jloopi: DO ij = ijs_i, ije_i ! Loop over Laguerre moments
j_int=jarray_i(ij)
- IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute
+ IF((CLOS .NE. 1) .OR. (p_int+2*j_int .LE. dmaxi)) THEN !compute for every moments except for closure 1
! Set non linear sum truncation
IF (NL_CLOS .EQ. -2) THEN
nmax = Jmaxi
@@ -152,7 +152,7 @@ ENDIF
ENDIF
bracket_sum_r = 0._dp ! initialize sum over real nonlinear term
nloopi: DO in = 1,nmax+1 ! Loop over laguerre for the sum
-!-----------!! ELECTROMAGNETIC CONTRIBUTION -sqrt(tau)/sigma*{Sum_s dnjs [sqrt(p+1)Nap+1s + sqrt(p)Nap-1s], Kernel psi}
+!-----------!! ELECTROSTATIC CONTRIBUTION
! First convolution terms
F_cmpx(ikys:ikye,ikxs:ikxe) = phi(ikys:ikye,ikxs:ikxe,iz) * kernel_i(in, ikys:ikye,ikxs:ikxe, iz, eo)
! Second convolution terms
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 54f99c63..7a4b2783 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -286,8 +286,8 @@ SUBROUTINE compute_radial_ion_heatflux
ENDDO
ENDIF
! Add polarisation contribution
- integrant(izs:ize) = integrant(izs:ize) + tau_i*imagu*ky_&
- *CONJG(phi(iky,ikx,izs:ize))*phi(iky,ikx,izs:ize) * HF_phi_correction_operator(iky,ikx,izs:ize)
+ ! integrant(izgs:izge) = integrant(izgs:izge) + tau_i*imagu*ky_&
+ ! *CONJG(phi(iky,ikx,izgs:izge))*phi(iky,ikx,izgs:izge) * HF_phi_correction_operator(iky,ikx,izgs:izge)
! Integrate over z
integrant(izgs:izge) = Jacobian(izgs:izge,0)*integrant(izgs:izge)
call simpson_rule_z(integrant,integral)
diff --git a/testcases/miller_example/bug/fort_00.90 b/testcases/miller_example/bug/fort_00.90
new file mode 100644
index 00000000..728a3456
--- /dev/null
+++ b/testcases/miller_example/bug/fort_00.90
@@ -0,0 +1,86 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.01
+ tmax = 100
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 2
+ jmaxe = 1
+ pmaxi = 2
+ jmaxi = 1
+ Nx = 128
+ Lx = 100
+ Ny = 64
+ Ly = 120
+ Nz = 16
+ Npol = 1
+ SG = .false.
+/
+&GEOMETRY
+ geom = 'miller'
+ q0 = 1.4
+ shear = 0.0
+ eps = 0.18
+/
+&OUTPUT_PAR
+ nsave_0d = 50
+ nsave_1d = -1
+ nsave_2d = Inf
+ nsave_3d = 100
+ nsave_5d = 500
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .true.
+ write_Sapj = .false.
+ write_dens = .true.
+ write_temp = .true.
+ job2load = -1
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .false.
+ mu_x = 0.1
+ mu_y = 0.1
+ N_HD = 4
+ mu_z = 0.2
+ mu_p = 0
+ mu_j = 0
+ nu = 0.2
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ni = 2.22
+ K_Ti = 6.92
+ K_Ne = 0
+ K_Te = 0
+ GradB = 1
+ CurvB = 1
+ lambdaD = 0
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ gyrokin_CO = .false.
+ interspecies = .true.
+ mat_file = 'null'
+/
+&INITIAL_CON
+ INIT_OPT = 'ppj'
+ ACT_ON_MODES = 'none'
+ init_background = 0
+ init_noiselvl = 1e-3
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/miller_example/bug/fort_01.90 b/testcases/miller_example/bug/fort_01.90
new file mode 100644
index 00000000..7deb16c3
--- /dev/null
+++ b/testcases/miller_example/bug/fort_01.90
@@ -0,0 +1,86 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.01
+ tmax = 400
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 2
+ jmaxe = 1
+ pmaxi = 2
+ jmaxi = 1
+ Nx = 128
+ Lx = 100
+ Ny = 64
+ Ly = 120
+ Nz = 16
+ Npol = 1
+ SG = .false.
+/
+&GEOMETRY
+ geom = 'miller'
+ q0 = 1.4
+ shear = 0.0
+ eps = 0.18
+/
+&OUTPUT_PAR
+ nsave_0d = 50
+ nsave_1d = -1
+ nsave_2d = Inf
+ nsave_3d = 100
+ nsave_5d = 500
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .true.
+ write_Sapj = .false.
+ write_dens = .true.
+ write_temp = .true.
+ job2load = 0
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .false.
+ mu_x = 0.1
+ mu_y = 0.1
+ N_HD = 4
+ mu_z = 0.2
+ mu_p = 0
+ mu_j = 0
+ nu = 0.2
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ni = 2.22
+ K_Ti = 6.92
+ K_Ne = 0
+ K_Te = 0
+ GradB = 1
+ CurvB = 1
+ lambdaD = 0
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ gyrokin_CO = .false.
+ interspecies = .true.
+ mat_file = 'null'
+/
+&INITIAL_CON
+ INIT_OPT = 'ppj'
+ ACT_ON_MODES = 'none'
+ init_background = 0
+ init_noiselvl = 1e-3
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/miller_example/fort_00.90 b/testcases/miller_example/fort_00.90
index 6958e7ad349208316d3f2e9cb71f5bb53d47b7f5..728a3456cca5d895f068ef8ad01cd766d2211cbd 100644
GIT binary patch
delta 58
zcmeys^^R-8a(P1o11<#xze)uVuvIWK0rGrcJVT?6r-fOV3bH0YU=d}>%g;-l{Eejw
E0Lm8-cK`qY
delta 76
zcmaFI^?_@`auqWpJqu$)OD+Wkze)uVuvIWLvH)^@U|b{9jmLypBr-D9fufGiA@Tlx
T@xJ~puE8L+>I{>=uv7s6<8Tu(
diff --git a/testcases/miller_example_w_salpha/fort_00.90 b/testcases/miller_example_w_salpha/fort_00.90
new file mode 100644
index 00000000..2701de9d
--- /dev/null
+++ b/testcases/miller_example_w_salpha/fort_00.90
@@ -0,0 +1,86 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.01
+ tmax = 100
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 2
+ jmaxe = 1
+ pmaxi = 2
+ jmaxi = 1
+ Nx = 128
+ Lx = 100
+ Ny = 64
+ Ly = 120
+ Nz = 16
+ Npol = 1
+ SG = .false.
+/
+&GEOMETRY
+ geom = 's-alpha'
+ q0 = 1.4
+ shear = 0.0
+ eps = 0.18
+/
+&OUTPUT_PAR
+ nsave_0d = 50
+ nsave_1d = -1
+ nsave_2d = Inf
+ nsave_3d = 100
+ nsave_5d = 500
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .true.
+ write_Sapj = .false.
+ write_dens = .true.
+ write_temp = .true.
+ job2load = -1
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .false.
+ mu_x = 0.1
+ mu_y = 0.1
+ N_HD = 4
+ mu_z = 0.2
+ mu_p = 0
+ mu_j = 0
+ nu = 0.2
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ni = 2.22
+ K_Ti = 6.92
+ K_Ne = 0
+ K_Te = 0
+ GradB = 1
+ CurvB = 1
+ lambdaD = 0
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ gyrokin_CO = .false.
+ interspecies = .true.
+ mat_file = 'null'
+/
+&INITIAL_CON
+ INIT_OPT = 'ppj'
+ ACT_ON_MODES = 'none'
+ init_background = 0
+ init_noiselvl = 1e-3
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/testcases/miller_example_w_salpha/fort_01.90 b/testcases/miller_example_w_salpha/fort_01.90
new file mode 100644
index 00000000..7deb16c3
--- /dev/null
+++ b/testcases/miller_example_w_salpha/fort_01.90
@@ -0,0 +1,86 @@
+&BASIC
+ nrun = 100000000
+ dt = 0.01
+ tmax = 400
+ maxruntime = 356400
+/
+&GRID
+ pmaxe = 2
+ jmaxe = 1
+ pmaxi = 2
+ jmaxi = 1
+ Nx = 128
+ Lx = 100
+ Ny = 64
+ Ly = 120
+ Nz = 16
+ Npol = 1
+ SG = .false.
+/
+&GEOMETRY
+ geom = 'miller'
+ q0 = 1.4
+ shear = 0.0
+ eps = 0.18
+/
+&OUTPUT_PAR
+ nsave_0d = 50
+ nsave_1d = -1
+ nsave_2d = Inf
+ nsave_3d = 100
+ nsave_5d = 500
+ write_doubleprecision = .false.
+ write_gamma = .true.
+ write_hf = .true.
+ write_phi = .true.
+ write_Na00 = .true.
+ write_Napj = .true.
+ write_Sapj = .false.
+ write_dens = .true.
+ write_temp = .true.
+ job2load = 0
+/
+&MODEL_PAR
+ ! Collisionality
+ CLOS = 0
+ NL_CLOS = -1
+ LINEARITY = 'nonlinear'
+ KIN_E = .false.
+ mu_x = 0.1
+ mu_y = 0.1
+ N_HD = 4
+ mu_z = 0.2
+ mu_p = 0
+ mu_j = 0
+ nu = 0.2
+ tau_e = 1
+ tau_i = 1
+ sigma_e = 0.023338
+ sigma_i = 1
+ q_e = -1
+ q_i = 1
+ K_Ni = 2.22
+ K_Ti = 6.92
+ K_Ne = 0
+ K_Te = 0
+ GradB = 1
+ CurvB = 1
+ lambdaD = 0
+ beta = 0
+/
+&COLLISION_PAR
+ collision_model = 'DG'
+ gyrokin_CO = .false.
+ interspecies = .true.
+ mat_file = 'null'
+/
+&INITIAL_CON
+ INIT_OPT = 'ppj'
+ ACT_ON_MODES = 'none'
+ init_background = 0
+ init_noiselvl = 1e-3
+ iseed = 42
+/
+&TIME_INTEGRATION_PAR
+ numerical_scheme = 'RK4'
+/
diff --git a/wk/Zpinch_coll_scan_kN_1.7.m b/wk/Zpinch_coll_scan_kN_1.7.m
index e223e0e4..ad00ba14 100644
--- a/wk/Zpinch_coll_scan_kN_1.7.m
+++ b/wk/Zpinch_coll_scan_kN_1.7.m
@@ -1,5 +1,5 @@
-%%
if 0
+%%
figure
Kn = 1.7;
@@ -11,27 +11,42 @@ Kn = 1.7;
% errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama DK (4,2)'); hold on
-% SUGAMA GK 4,2
+% SUGAMA GK 4,2 200x32
nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0];
-Gavg_a = [2.54e-2 4.66e-2 6.96e-2 8.98e-2 1.06e-1 1.24e-1 1.43e-1 1.52e-1 1.69e-1 1.09e-1];
+Gavg_a = [2.54e-2 4.66e-2 6.96e-2 8.98e-2 1.06e-1 1.24e-1 1.43e-1 1.52e-1 1.69e-1 1.79e-1];
Gstd_a = [3.04e-2 1.42e-2 1.56e-2 1.23e-2 1.20e-2 1.57e-2 1.63e-2 2.06e-2 2.14e-02 1.78e-2];
-errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2)'); hold on
+errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2) 200x32'); hold on
+
+% SUGAMA GK 4,2 256x64
+nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0];
+Gavg_a = [1.99e-2 3.03e-2 4.44e-2 5.87e-2 7.87e-2 1.07e-1 1.22e-1 1.38e-1 1.63e-1 1.75e-1];
+Gstd_a = [1.55e-2 2.22e-2 1.12e-2 1.54e-2 2.20e-2 2.76e-2 2.44e-2 3.20e-2 3.12e-2 3.13e-2];
+errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (4,2) 256x64'); hold on
-% FCGK 4,2
+
+% SUGAMA GK 6,3 200x32
+nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0];
+Gavg_a = [2.35e-2 3.57e-2 5.09e-2 6.97e-2 8.65e-2 1.01e-1 1.16e-1 1.33e-1 1.49e-1 1.62e-1];
+Gstd_a = [3.76e-2 3.91e-2 2.96e-2 1.57e-2 1.73e-2 1.94e-2 2.21e-2 2.01e-2 1.93e-2 2.24e-2];
+
+errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Sugama GK (6,3) 200x32'); hold on
+
+
+% FCGK 4,2 200x32
nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 10.0];
Gavg_a = [8.57e-2 1.45e-1 2.25e-1 2.87e-1 3.48e-1 4.06e-1 4.51e-1 3.65e-1*Kn];
Gstd_a = [2.07e-2 2.61e-2 2.40e-2 3.46e-2 4.30e-2 5.00e-2 5.11e-2 0];
-errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Coulomb (4,2)'); hold on
-
-% LDGK ii 6,3
-nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00];
-Gavg_a = [3.86e-2 1.82e-2 3.08e-2 5.24e-2 7.08e-2 8.26e-2 5.78e-2 7.16e-2 7.96e-2];
-Gstd_a = [3.52e-2 1.87e-2 2.86e-2 2.79e-2 1.72e-2 2.40e-2 2.46e-2 1.01e-2 1.21e-2];
+errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Coulomb (4,2) 200x32'); hold on
-errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Landau ii (6,3)'); hold on
+% % LDGK ii 6,3 200x32
+% nu_a = 1e-2*[1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00];
+% Gavg_a = [3.86e-2 1.82e-2 3.08e-2 5.24e-2 7.08e-2 8.26e-2 5.78e-2 7.16e-2 7.96e-2];
+% Gstd_a = [3.52e-2 1.87e-2 2.86e-2 2.79e-2 1.72e-2 2.40e-2 2.46e-2 1.01e-2 1.21e-2];
+%
+% errorbar(nu_a, Gavg_a/Kn, Gstd_a/Kn,'DisplayName','Landau ii (6,3) 200x32'); hold on
% Collisionless
plot([0 1], 0.02343*[1 1],'--k','DisplayName','$\nu=0$');
@@ -41,7 +56,7 @@ plot([0 1], 0.02343*[1 1],'--k','DisplayName','$\nu=0$');
xlim([0 0.1]);
legend('show');
xlabel('$\nu R/c_s$'); ylabel('$\Gamma_x^\infty/\kappa_N$');
-
+set(gca,'Yscale','log');
end
if 0
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index f45fc16d..de09d169 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -28,7 +28,8 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
% folder = '/misc/gene_results/CBC/KT_9_128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/';
% folder = '/misc/gene_results/CBC/Lapillone_Fig6/';
-folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/';
+% folder = '/misc/gene_results/Z-pinch/HP_kN_1.6_adapt_mu_01/';
+folder = '/misc/gene_results/miller/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index eecf7147..65ba29bd 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -28,26 +28,32 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
% data.scale = 1;%/(data.Nx*data.Ny)^2;
-i_ = 11;
+i_ = 1;
disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
options.TAVG_0 = data.TJOB_SE(i_);%0.4*data.Ts3D(end);
options.TAVG_1 = data.TJOB_SE(i_+1);%0.9*data.Ts3D(end); % Averaging times duration
options.NCUT = 4; % Number of cuts for averaging and error estimation
-options.NMVA = 100; % Moving average for time traces
+options.NMVA = 1; % Moving average for time traces
% options.ST_FIELD = '\Gamma_x'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
-% options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
-% options.INTERP = 1;
+options.ST_FIELD = '\phi'; % chose your field to plot in spacetime diag (e.g \phi,v_x,G_x)
+options.INTERP = 0;
+options.RESOLUTION = 256;
fig = plot_radial_transport_and_spacetime(data,options);
save_figure(data,fig,'.png')
end
if 0
%% statistical transport averaging
-options.T = [200 400];
+for i_ = 1:2:21
+% i_ = 3;
+disp([num2str(data.TJOB_SE(i_)),' ',num2str(data.TJOB_SE(i_+1))])
+disp([num2str(data.NU_EVOL(i_)),' ',num2str(data.NU_EVOL(i_+1))])
+options.T = [data.TJOB_SE(i_) data.TJOB_SE(i_+1)];
+options.NPLOTS = 0;
fig = statistical_transport_averaging(data,options);
end
-
+end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
@@ -55,7 +61,7 @@ options.INTERP = 1;
options.POLARPLOT = 0;
options.NAME = '\phi';
% options.NAME = '\omega_z';
-% options.NAME = 'N_i^{00}';
+% options.NAME = 'N_e^{00}';
% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
@@ -63,10 +69,10 @@ options.NAME = '\phi';
options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
-% options.COMP = 'avg';
+options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
% options.TIME = data.Ts3D;
-options.TIME = [000:50:7000];
+options.TIME = [000:0.1:7000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
options.RESOLUTION = 256;
@@ -146,7 +152,7 @@ end
if 0
%% Time averaged spectrum
-options.TIME = [300 600];
+options.TIME = [2000 3000];
options.NORM =1;
options.NAME = '\phi';
% options.NAME = 'N_i^{00}';
diff --git a/wk/header_2DZP_results.m b/wk/header_2DZP_results.m
index b180dd52..10428ea0 100644
--- a/wk/header_2DZP_results.m
+++ b/wk/header_2DZP_results.m
@@ -191,9 +191,12 @@ resdir ='';
% resdir ='Zpinch_rerun/kN_1.7_SGGK_conv_200x32x7x3_nu_0.01';
% resdir ='Zpinch_rerun/kN_1.7_LDGKii_200x32x7x3_nu_scan';
% resdir ='Zpinch_rerun/nu_0.1_LDGKii_200x48x7x4_kN_scan';
-resdir ='Zpinch_rerun/nu_0.1_FCGK_200x48x5x3_kN_scan';
+% resdir ='Zpinch_rerun/nu_0.1_FCGK_200x48x5x3_kN_scan';
+% resdir ='Zpinch_rerun/nu_0.1_SGGK_200x48x5x3_kN_scan';
% resdir = 'Zpinch_rerun/kN_1.7_FCGK_200x32x5x3_nu_scan';
% resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x7x4_nu_scan';
+% resdir = 'Zpinch_rerun/kN_1.7_SGGK_200x32x5x3_nu_scan';
+resdir = 'Zpinch_rerun/kN_1.7_SGGK_256x64x5x3_nu_scan';
%%
JOBNUMMIN = 00; JOBNUMMAX = 10;
resdir = ['results/',resdir];
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 2a56c2df..189cc0e9 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -1,64 +1,64 @@
% Directory of the code "mypathtoHeLaZ/HeLaZ/"
-helazdir = '/home/ahoffman/HeLaZ/';
+gyacomodir = '/home/ahoffman/gyacomo/';
% Directory of the simulation (from results)
% if 1% Local results
-% outfile ='volcokas/64x32x16x5x3_kin_e_npol_1';
+% resdir ='volcokas/64x32x16x5x3_kin_e_npol_1';
%% Dimits
-% outfile ='shearless_cyclone/128x64x16x5x3_Dim_90';
-% outfile ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60';
-% outfile ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70';
-% outfile ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70';
+% resdir ='shearless_cyclone/128x64x16x5x3_Dim_90';
+% resdir ='shearless_cyclone/128x64x16x9x5_Dim_scan/128x64x16x9x5_Dim_60';
+% resdir ='shearless_cyclone/128x64x16x5x3_Dim_scan/128x64x16x5x3_Dim_70';
+% resdir ='shearless_cyclone/64x32x16x5x3_Dim_scan/64x32x16x5x3_Dim_70';
%% AVS
-% outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1';
-% outfile = 'volcokas/64x32x16x5x3_kin_e_npol_1';
-% outfile = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine';
-% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
-% outfile = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine';
-% outfile = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
+% resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1';
+% resdir = 'volcokas/64x32x16x5x3_kin_e_npol_1';
+% resdir = 'shearless_cyclone/64x32x80x5x3_CBC_Npol_5_kine';
+% resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
+% resdir = 'shearless_cyclone/64x32x160x5x3_CBC_Npol_10_kine';
+% resdir = 'shearless_cyclone/96x32x160x5x3_CBC_Npol_10_kine';
%% shearless CBC
-% outfile ='shearless_cyclone/64x32x16x5x3_CBC_080';
-% outfile ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100';
-% outfile ='shearless_cyclone/64x32x16x5x3_CBC_120';
+% resdir ='shearless_cyclone/64x32x16x5x3_CBC_080';
+% resdir ='shearless_cyclone/64x32x16x5x3_CBC_scan/64x32x16x5x3_CBC_100';
+% resdir ='shearless_cyclone/64x32x16x5x3_CBC_120';
-% outfile ='shearless_cyclone/64x32x16x9x5_CBC_080';
-% outfile ='shearless_cyclone/64x32x16x9x5_CBC_100';
-% outfile ='shearless_cyclone/64x32x16x9x5_CBC_120';
+% resdir ='shearless_cyclone/64x32x16x9x5_CBC_080';
+% resdir ='shearless_cyclone/64x32x16x9x5_CBC_100';
+% resdir ='shearless_cyclone/64x32x16x9x5_CBC_120';
-% outfile = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK';
+% resdir = 'shearless_cyclone/64x32x16x5x3_CBC_CO/64x32x16x5x3_CBC_LRGK';
%% CBC
-% outfile = 'CBC/64x32x16x5x3';
-% outfile = 'CBC/64x128x16x5x3';
-% outfile = 'CBC/128x64x16x5x3';
-% outfile = 'CBC/96x96x16x3x2_Nexc_6';
-% outfile = 'CBC/128x96x16x3x2';
-% outfile = 'CBC/192x96x16x3x2';
-% outfile = 'CBC/192x96x24x13x7';
-% outfile = 'CBC/kT_11_128x64x16x5x3';
-% outfile = 'CBC/kT_9_256x128x16x3x2';
-% outfile = 'CBC/kT_4.5_128x64x16x13x3';
-% outfile = 'CBC/kT_4.5_192x96x24x13x7';
-% outfile = 'CBC/kT_4.5_128x64x16x13x7';
-% outfile = 'CBC/kT_4.5_128x96x24x15x5';
-% outfile = 'CBC/kT_5.3_192x96x24x13x7';
-% outfile = 'CBC/kT_13_large_box_128x64x16x5x3';
-% outfile = 'CBC/kT_11_96x64x16x5x3_ky_0.02';
+% resdir = 'CBC/64x32x16x5x3';
+% resdir = 'CBC/64x128x16x5x3';
+% resdir = 'CBC/128x64x16x5x3';
+% resdir = 'CBC/96x96x16x3x2_Nexc_6';
+% resdir = 'CBC/128x96x16x3x2';
+% resdir = 'CBC/192x96x16x3x2';
+% resdir = 'CBC/192x96x24x13x7';
+% resdir = 'CBC/kT_11_128x64x16x5x3';
+% resdir = 'CBC/kT_9_256x128x16x3x2';
+% resdir = 'CBC/kT_4.5_128x64x16x13x3';
+% resdir = 'CBC/kT_4.5_192x96x24x13x7';
+% resdir = 'CBC/kT_4.5_128x64x16x13x7';
+% resdir = 'CBC/kT_4.5_128x96x24x15x5';
+% resdir = 'CBC/kT_5.3_192x96x24x13x7';
+% resdir = 'CBC/kT_13_large_box_128x64x16x5x3';
+% resdir = 'CBC/kT_11_96x64x16x5x3_ky_0.02';
-% outfile = 'CBC/kT_scan_128x64x16x5x3';
-% outfile = 'CBC/kT_scan_192x96x16x3x2';
-% outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
-% outfile = 'dbg/nexc_dbg';
-outfile = 'CBC/NM_F4_kT_4.5_192x64x24x6x4';
+% resdir = 'CBC/kT_scan_128x64x16x5x3';
+% resdir = 'CBC/kT_scan_192x96x16x3x2';
+% resdir = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
+% resdir = 'dbg/nexc_dbg';
+% resdir = 'CBC/NM_F4_kT_4.5_192x64x24x6x4';
-% outfile = 'CBC_Ke_EM/192x96x24x5x3';
-% outfile = 'CBC_Ke_EM/96x48x16x5x3';
-% outfile = 'CBC_Ke_EM/minimal_res';
+% resdir = 'CBC_Ke_EM/192x96x24x5x3';
+% resdir = 'CBC_Ke_EM/96x48x16x5x3';
+% resdir = 'CBC_Ke_EM/minimal_res';
%% KBM
-% outfile = 'NL_KBM/192x64x24x5x3';
+% resdir = 'NL_KBM/192x64x24x5x3';
%% Linear CBC
-% outfile = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
-
+% resdir = 'linear_CBC/20x2x32_21x11_Lx_62.8319_Ly_31.4159_q0_1.4_e_0.18_s_0.8_kN_2.22_kT_5.3_nu_1e-02_DGDK_adiabe';
+resdir = 'testcases/miller_example_bug';
JOBNUMMIN = 00; JOBNUMMAX = 10;
-run analysis_HeLaZ
+run analysis_gyacomo
diff --git a/wk/lin_EPY.m b/wk/lin_EPY.m
new file mode 100644
index 00000000..8a006017
--- /dev/null
+++ b/wk/lin_EPY.m
@@ -0,0 +1,178 @@
+%% QUICK RUN SCRIPT for linear entropy mode in a Zpinch
+% This script create a directory in /results and run a simulation directly
+% from matlab framework. It is meant to run only small problems in linear
+% for benchmark and debugging purpose since it makes matlab "busy"
+%
+SIMID = 'lin_EPY'; % Name of the simulation
+RUN = 1; % To run or just to load
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+PROGDIR = '/home/ahoffman/gyacomo/';
+EXECNAME = 'gyacomo';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Set Up parameters
+CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% PHYSICAL PARAMETERS
+NU = 0.1; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 2.0; % ele Density '''
+K_Te = 0.5; % ele Temperature '''
+K_Ni = 2.0; % ion Density gradient drive
+K_Ti = 0.5; % ion Temperature '''
+SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.0; % electron plasma beta
+%% GRID PARAMETERS
+P = 4;
+J = P/2;
+PMAXE = P; % Hermite basis size of electrons
+JMAXE = J; % Laguerre "
+PMAXI = P; % " ions
+JMAXI = J; % "
+NX = 2; % real space x-gridpoints
+NY = 100; % '' y-gridpoints
+LX = 2*pi/0.8; % Size of the squared frequency domain
+LY = 2*pi/0.05; % Size of the squared frequency domain
+NZ = 1; % number of perpendicular planes (parallel grid)
+NPOL = 1;
+SG = 0; % Staggered z grids option
+%% GEOMETRY
+GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
+Q0 = 1.4; % safety factor
+SHEAR = 0.8; % magnetic shear
+NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+EPS = 0.18; % inverse aspect ratio
+%% TIME PARMETERS
+TMAX = 50; % Maximal time unit
+DT = 1e-2; % Time step
+SPS0D = 1; % Sampling per time unit for 2D arrays
+SPS2D = 0; % Sampling per time unit for 2D arrays
+SPS3D = 2; % Sampling per time unit for 2D arrays
+SPS5D = 1/5; % Sampling per time unit for 5D arrays
+SPSCP = 0; % Sampling per time unit for checkpoints
+JOB2LOAD= -1;
+%% OPTIONS
+LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+% Collision operator
+% (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+CO = 'LD';
+GKCO = 1; % gyrokinetic operator
+ABCO = 1; % interspecies collisions
+INIT_ZF = 0; ZF_AMP = 0.0;
+CLOS = 0; % Closure model (0: =0 truncation, 1: v^Nmax closure (p+2j<=Pmax))s
+NL_CLOS = 0; % nonlinear closure model (-2:nmax=jmax; -1:nmax=jmax-j; >=0:nmax=NL_CLOS)
+KERN = 0; % Kernel model (0 : GK)
+INIT_OPT= 'mom00'; % Start simulation with a noisy mom00/phi/allmom
+%% OUTPUTS
+W_DOUBLE = 1;
+W_GAMMA = 1; W_HF = 1;
+W_PHI = 1; W_NA00 = 1;
+W_DENS = 1; W_TEMP = 1;
+W_NAPJ = 1; W_SAPJ = 0;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% unused
+HD_CO = 0.0; % Hyper diffusivity cutoff ratio
+MU = 0.0; % Hyperdiffusivity coefficient
+INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
+MU_X = MU; %
+MU_Y = MU; %
+N_HD = 4;
+MU_Z = 2.0; %
+MU_P = 0.0; %
+MU_J = 0.0; %
+LAMBDAD = 0.0;
+NOISE0 = 0.0e-5; % Init noise amplitude
+BCKGD0 = 1.0; % Init background
+GRADB = 1.0;
+CURVB = 1.0;
+%%-------------------------------------------------------------------------
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+if RUN
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',PROGDIR,'bin/',EXECNAME,' 1 6 1 0; cd ../../../wk'])
+end
+
+%% Load results
+%%
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR = [PROGDIR,'results/',filename,'/'];
+% Load outputs from jobnummin up to jobnummax
+JOBNUMMIN = 00; JOBNUMMAX = 00;
+data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
+
+%% Short analysis
+if 1
+%% linear growth rate (adapted for 2D zpinch and fluxtube)
+trange = [0.5 1]*data.Ts3D(end);
+nplots = 1; % 1 for only growth rate and error, 2 for omega local evolution, 3 for plot according to z
+lg = compute_fluxtube_growth_rate(data,trange,nplots);
+[gmax, kmax] = max(lg.g_ky(:,end));
+[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
+msg = sprintf('gmax = %2.2f, kmax = %2.2f',gmax,lg.ky(kmax)); disp(msg);
+msg = sprintf('gmax/k = %2.2f, kmax/k = %2.2f',gmaxok,lg.ky(kmaxok)); disp(msg);
+end
+
+if 0
+%% Ballooning plot
+options.time_2_plot = [120];
+options.kymodes = [0.9];
+options.normalized = 1;
+% options.field = 'phi';
+fig = plot_ballooning(data,options);
+end
+
+if 0
+%% Hermite-Laguerre spectrum
+% options.TIME = 'avg';
+options.P2J = 1;
+options.ST = 1;
+options.PLOT_TYPE = 'space-time';
+% options.PLOT_TYPE = 'Tavg-1D';
+% options.PLOT_TYPE = 'Tavg-2D';
+options.NORMALIZED = 0;
+options.JOBNUM = 0;
+options.TIME = [0 50];
+options.specie = 'i';
+options.compz = 'avg';
+fig = show_moments_spectrum(data,options);
+% fig = show_napjz(data,options);
+save_figure(data,fig)
+end
+
+if 0
+%% linear growth rate for 3D Zpinch (kz fourier transform)
+trange = [0.5 1]*data.Ts3D(end);
+options.keq0 = 1; % chose to plot planes at k=0 or max
+options.kxky = 1;
+options.kzkx = 0;
+options.kzky = 0;
+[lg, fig] = compute_3D_zpinch_growth_rate(data,trange,options);
+save_figure(data,fig)
+end
+if 0
+%% Mode evolution
+options.NORMALIZED = 0;
+options.K2PLOT = 1;
+options.TIME = [0:1000];
+options.NMA = 1;
+options.NMODES = 1;
+options.iz = 'avg';
+fig = mode_growth_meter(data,options);
+save_figure(data,fig,'.png')
+end
+
+
+if 0
+%% RH TEST
+ikx = 2; iky = 2; t0 = 0; t1 = data.Ts3D(end);
+[~, it0] = min(abs(t0-data.Ts3D));[~, it1] = min(abs(t1-data.Ts3D));
+plt = @(x) squeeze(mean(real(x(iky,ikx,:,it0:it1)),3))./squeeze(mean(real(x(iky,ikx,:,it0)),3));
+figure
+plot(data.Ts3D(it0:it1), plt(data.PHI));
+xlabel('$t$'); ylabel('$\phi_z(t)/\phi_z(0)$')
+title(sprintf('$k_x=$%2.2f, $k_y=$%2.2f',data.kx(ikx),data.ky(iky)))
+end
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index bdf3ec05..cffb9cc6 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -7,9 +7,8 @@ SIMID = 'lin_ITG'; % Name of the simulation
RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
-HELAZDIR = '/home/ahoffman/HeLaZ/';
-% EXECNAME = 'helaz3';
-EXECNAME = 'helaz3_dbg';
+HELAZDIR = '/home/ahoffman/gyacomo/';
+EXECNAME = 'gyacomo';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
@@ -33,23 +32,22 @@ JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 11; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
+NY = 32; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 2*pi/0.3; % Size of the squared frequency domain
+LY = 2*pi/0.05; % Size of the squared frequency domain
NZ = 16; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
-% GEOMETRY= 'Z-pinch'; % Z-pinch overwrites q0, shear and eps
GEOMETRY= 's-alpha';
-% GEOMETRY= 'circular';
+% GEOMETRY= 'miller';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 25; % Maximal time unit
-DT = 5e-3; % Time step
+TMAX = 50; % Maximal time unit
+DT = 1e-2; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
SPS3D = 5; % Sampling per time unit for 2D arrays
--
GitLab