diff --git a/matlab/extract_fig_data.m b/matlab/extract_fig_data.m
index e64b370b5414d03f9f3bc4a9d21da581b81f0456..e7a98b9d7301ab9939cfabe60f3d35ac0155008a 100644
--- a/matlab/extract_fig_data.m
+++ b/matlab/extract_fig_data.m
@@ -7,7 +7,7 @@ for i = 1:numel(dataObjs)
X_ = [X_ dataObjs(i).XData];
Y_ = [Y_ dataObjs(i).YData];
end
-n0 = 250;
+n0 = 1;
n1 = numel(X_);
figure;
mvm = @(x) movmean(x,1);
@@ -16,7 +16,7 @@ shift = X_(n0);
% plot(X_(n0:end),Y_(n0:end));
plot(mvm(X_(n0:n1)-shift),mvm(Y_(n0:n1))); hold on;
-t0 = ceil(numel(X_)*0.5); t1 = numel(X_);
+t0 = ceil(numel(X_)*0.2); t1 = numel(X_);
avg= mean(Y_(t0:t1)); dev = std(Y_(t0:t1));
disp(['AVG =',sprintf('%2.2f',avg),'+-',sprintf('%2.2f',dev)]);
%
diff --git a/matlab/load/quick_gene_load.m b/matlab/load/quick_gene_load.m
index 3ac095c94f09c8abae00640542a841853413bdd0..8fb67265b2114225e6482cc0fa57574bf97c1c93 100644
--- a/matlab/load/quick_gene_load.m
+++ b/matlab/load/quick_gene_load.m
@@ -41,13 +41,17 @@
% fname ='GENE_LIN_Kn_1.6_KT_0.4_nu_0_32x16.txt';
% fname ='GENE_LIN_Kn_2.5_KT_0.625_nu_0_32x16.txt';
path = '/home/ahoffman/gene/linear_CBC_results/';
+fname = 'CBC_100_20x1x32x30x14_Lv_3_Lw_12_circ.txt';
% fname = 'CBC_100_20x1x32x32x12_Lv_3_Lw_12.txt';
% fname = 'CBC_KT_4_20x1x32x32x12_Lv_3_Lw_12.txt';
% fname = 'CBC_KT_4_20x1x32x64x24_Lv_6_Lw_24.txt';
% fname = 'CBC_KT_5.3_20x1x32x32x12_Lv_3_Lw_12.txt';
% fname = 'CBC_KT_5.3_32x1x48x40x16_Lv_3_Lw_12.txt';
-fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12.txt';
+% fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12.txt';
% fname = 'CBC_ky_0.3_20x1x32x32x12_Lv_3_Lw_12_nuv_1e-3.txt';
+% fname = 'CBC_KT_11_20x1x16x24x10_Lv_3_Lw_12.txt';
+% fname = 'CBC_KT_11_20x1x32x30x14_Lv_3_Lw_12.txt';
+% fname = 'CBC_ky_0.3_20x1x16x24x10_Lv_3_Lw_12_nuv_1e-3.txt';
data_ = load([path,fname]);
figure
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index 22502ecd02273a5d775ddb3aa8c83cc73155ac0e..6e6df0128c30bf11d8da6142d221ca877a42a098 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -14,14 +14,13 @@ function [FIG] = plot_ballooning(data,options)
ncol = 2;
end
% Apply baollooning tranform
+ nexc = round(data.ky(2)*data.SHEAR*2*pi/data.kx(2));
for iky=ikyarray
dims = size(phi_real);
Nkx = dims(2);
is = max(1,iky-1);
- Npi = (Nkx-1)-2*(is-1);
- if(Npi <= 0)
- break
- elseif(Npi == 1)
+ Npi = (Nkx-1)-2*nexc*(is-1);
+ if(Npi <= 1)
ordered_ikx = 1;
else
tmp_ = (Nkx-is+1):-is:(Nkx/2+2);
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 7ad93c0b6bbb408d6f7420c8312e990407d0b3af..de18d2d53a5d0431135166c8b020ea3fa38c76ba 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -17,7 +17,17 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, OPTIONS)
ikzf = min([ikzf,DATA.Nky]);
Ns3D = numel(DATA.Ts3D);
[KX, KY] = meshgrid(DATA.kx, DATA.ky);
-
+ %% error estimation
+ DT_ = (tend-tstart)/OPTIONS.NCUT;
+ Qx_ee = zeros(1,OPTIONS.NCUT);
+ for i = 1:OPTIONS.NCUT
+ [~,its_] = min(abs(DATA.Ts0D - (tstart+(i-1)*DT_)));
+ [~,ite_] = min(abs(DATA.Ts0D - (tstart+ i *DT_)));
+ Qx_ee(i) = mean(DATA.HFLUX_X(its_:ite_))*SCALE;
+ end
+ Qx_avg = mean(Qx_ee);
+ Qx_err = std(Qx_ee);
+ disp(['Q_avg=',sprintf('%2.2f',Qx_avg),'+-',sprintf('%2.2f',Qx_err)]);
%% computations
% Compute Gamma from ifft matlab
@@ -66,7 +76,7 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
% plot(mvm(DATA.Ts3D),mvm(Qx_t_mtlb),'DisplayName','matlab comp.'); hold on;
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
- 'DisplayName',['$Q_x^{\infty} = $',num2str(Qx_infty_avg),'$\pm$',num2str(Qx_infty_std)]);
+ 'DisplayName',['$Q_{avg}=',sprintf('%2.2f',Qx_avg),'\pm',sprintf('%2.2f',Qx_err),'$']); legend('show');
ylabel('$Q_x$')
ylim([0,5*abs(Qx_infty_avg)]);
xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
diff --git a/matlab/setup.m b/matlab/setup.m
index 07f5a65a77b955ee0d263b22a1652e653b9126b2..8600aa5c1d0ff30bd806f35f4a73c62b0bb90f93 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -7,6 +7,7 @@ GRID.pmaxi = PMAXI; % Ion Hermite moments
GRID.jmaxi = JMAXI; % Ion Laguerre moments
GRID.Nx = NX; % x grid resolution
GRID.Lx = LX; % x length
+GRID.Nexc = NEXC; % to extend Lx when s>0
GRID.Ny = NY; % y ''
GRID.Ly = LY; % y ''
GRID.Nz = NZ; % z resolution
@@ -127,6 +128,11 @@ if (NZ > 1) %3D case
geo_ = [geo_,'_s_',num2str(SHEAR)];
end
end
+switch GEOMETRY
+ case 'circular'
+ geo_ = [geo_,'_circ_'];
+end
+
% put everything together in the param character chain
u_ = '_'; % underscore variable
PARAMS = [res_,HLdeg_,geo_,drives_,coll_,lin_,adiabe_];
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 7dd03641a0628fa4312808e049940e525e0c7dc5..69fe008dc0338f5d324827ab59ed1dbbde32874e 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -17,6 +17,7 @@ fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),'\n']);
fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),'\n']);
fprintf(fid,[' Nx = ', num2str(GRID.Nx),'\n']);
fprintf(fid,[' Lx = ', num2str(GRID.Lx),'\n']);
+fprintf(fid,[' Nexc = ', num2str(GRID.Nexc),'\n']);
fprintf(fid,[' Ny = ', num2str(GRID.Ny),'\n']);
fprintf(fid,[' Ly = ', num2str(GRID.Ly),'\n']);
fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']);
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 2d4c4fffa44e3ae5ba0d869fd89d117a43eb317f..ff0e309e0c7b774856e15fc88a79bc7e19fb9dde 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -116,7 +116,7 @@ CONTAINS
SUBROUTINE eval_salphaB_geometry
! evaluate s-alpha geometry model
implicit none
- REAL(dp) :: z, kx, ky, alpha_MHD
+ REAL(dp) :: z, kx, ky, alpha_MHD, Gx, Gy
alpha_MHD = 0._dp
parity: DO eo = 0,1
@@ -128,7 +128,7 @@ CONTAINS
gxy(iz,eo) = shear*z - alpha_MHD*SIN(z)
gxz(iz,eo) = 0._dp
gyy(iz,eo) = 1._dp + (shear*z - alpha_MHD*SIN(z))**2
- gyz(iz,eo) = 1._dp/(eps + EPSILON(eps)) !avoid 1/0 in Zpinch config
+ gyz(iz,eo) = 1._dp/eps
gzz(iz,eo) = 0._dp
dxdR(iz,eo)= COS(z)
dxdZ(iz,eo)= SIN(z)
@@ -154,12 +154,22 @@ CONTAINS
gradyB(iz,eo) = 0._dp
gradzB(iz,eo) = eps * SIN(z) / hatR(iz,eo) ! 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) = (-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
+ ! ENDDO
+ ! ENDDO
! Curvature operator
- DO iky = ikys, ikye
+ Gx = (gxz(iz,eo) * gxy(iz,eo) - gxx(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Kx
+ Gy = -COS(Z) + (gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)) *eps*SIN(Z) ! Ky
+ 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) = (Gx*kx + Gy*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
@@ -176,58 +186,57 @@ CONTAINS
SUBROUTINE eval_circular_geometry
! evaluate circular geometry model
- ! Ref: Lapilonne et al., PoP, 2009
+ ! Ref: Lapilonne et al., PoP, 2009, GENE circular.F90 applied at r=r0
implicit none
- REAL(dp) :: X, kx, ky, Gamma1, Gamma2, Gamma3
+ REAL(dp) :: chi, kx, ky, Gx, Gy
parity: DO eo = 0,1
zloop: DO iz = izgs,izge
- X = zarray(iz,eo) - eps*SIN(zarray(iz,eo)) ! chi = theta - eps sin(theta)
+ chi = zarray(iz,eo) ! = chi
- ! metric
- gxx(iz,eo) = 1._dp
- gxy(iz,eo) = shear*X - eps*SIN(X)
- gxz(iz,eo) = - SIN(X)
- gyy(iz,eo) = 1._dp + (shear*X)**2 - 2._dp*eps*COS(X) - 2._dp*shear*X*eps*SIN(X)
- gyz(iz,eo) = 1._dp/(eps + EPSILON(eps)) - 2._dp*COS(X) - shear*X*SIN(X)
- gzz(iz,eo) = 1._dp/eps**2 - 2._dp*COS(X)/eps
- dxdR(iz,eo)= COS(X)
- dxdZ(iz,eo)= SIN(X)
+ ! metric in x,y,z
+ 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)
- ! Relative strengh of radius
- hatR(iz,eo) = 1._dp + eps*COS(X)
- hatZ(iz,eo) = 1._dp + eps*SIN(X)
+ ! Relative strengh of radius
+ hatR(iz,eo) = 1._dp + eps*COS(chi)
+ hatZ(iz,eo) = 1._dp + eps*SIN(chi)
- ! toroidal coordinates
- Rc (iz,eo) = hatR(iz,eo)
- phic(iz,eo) = X
- Zc (iz,eo) = hatZ(iz,eo)
+ ! toroidal coordinates
+ Rc (iz,eo) = hatR(iz,eo)
+ phic(iz,eo) = chi
+ Zc (iz,eo) = hatZ(iz,eo)
- ! Jacobian
- Jacobian(iz,eo) = q0*hatR(iz,eo)**2
+ ! Jacobian
+ Jacobian(iz,eo) = q0*hatR(iz,eo)
- ! Relative strengh of modulus of B
- hatB (iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2)
- hatB_NL(iz,eo) = SQRT(gxx(iz,eo)*gyy(iz,eo) - (gxy(iz,eo))**2) ! In front of the NL term
+ ! 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
- ! Derivative of the magnetic field strenght
- gradxB(iz,eo) = -(COS(X) + eps*SIN(X)**2)/hatB(iz,eo)**2 ! Gene put a factor hatB^2 or 1/hatR^2 in this
- gradyB(iz,eo) = 0._dp
- gradzB(iz,eo) = eps * SIN(X) * (1._dp - eps*COS(X)) / hatB(iz,eo)**2 ! Gene put a factor hatB or 1/hatR in this
+ ! 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
- 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)
- Gamma3 = gxz(iz,eo) * gyy(iz,eo) - gxy(iz,eo) * gyz(iz,eo)
- ! Curvature operator
- DO iky = ikys, ikye
- ky = kyarray(iky)
- DO ikx= ikxs, ikxe
- kx = kxarray(ikx)
- Ckxky(iky, ikx, iz,eo) = Gamma1*((kx + shear*X*ky)*gradzB(iz,eo)/eps - gradxB(iz,eo)*ky) ! .. 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)
+ ! 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
+ 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
+ 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
@@ -237,6 +246,102 @@ CONTAINS
!--------------------------------------------------------------------------------
!
+ 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
@@ -367,7 +472,7 @@ CONTAINS
kxmax_ = kx_max
IF(contains_zmax) THEN ! Check if the process is at the end of the FT
DO iky = ikys,ikye
- shift = 2._dp*shear*PI*kyarray(iky)*Npol
+ shift = 2._dp*PI*shear*kyarray(iky)*Npol*Nexc
DO ikx = ikxs,ikxe
kx_shift = kxarray(ikx) + shift
IF(kx_shift .GT. kxmax_) THEN ! outside of the frequ domain
@@ -390,7 +495,7 @@ CONTAINS
kxmin_ = -kx_max+deltakx
IF(contains_zmin) THEN ! Check if the process is at the start of the FT
DO iky = ikys,ikye
- shift = 2._dp*shear*PI*kyarray(iky)*Npol
+ shift = 2._dp*PI*shear*kyarray(iky)*Npol*Nexc
DO ikx = ikxs,ikxe
kx_shift = kxarray(ikx) - shift
IF( kx_shift .LT. kxmin_ ) THEN ! outside of the frequ domain
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index a31128d017ac0141b858a8cfd4d47fe05657c722..f6527f2dd7200f33d4a94a23b20c8196a7f66a27 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -16,6 +16,7 @@ MODULE grid
INTEGER, PUBLIC, PROTECTED :: dmaxi = 1 ! The maximal full GF set of i-moments v^dmax
INTEGER, PUBLIC, PROTECTED :: Nx = 16 ! Number of total internal grid points in x
REAL(dp), PUBLIC, PROTECTED :: Lx = 1._dp ! horizontal length of the spatial box
+ INTEGER, PUBLIC, PROTECTED :: Nexc = 1 ! factor to increase Lx when shear>0 (Lx = Nexc/kymin/shear)
INTEGER, PUBLIC, PROTECTED :: Ny = 16 ! Number of total internal grid points in y
REAL(dp), PUBLIC, PROTECTED :: Ly = 1._dp ! vertical length of the spatial box
INTEGER, PUBLIC, PROTECTED :: Nz = 1 ! Number of total perpendicular planes
@@ -123,7 +124,7 @@ CONTAINS
INTEGER :: lu_in = 90 ! File duplicated from STDIN
NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
- Nx, Lx, Ny, Ly, Nz, Npol, SG
+ Nx, Lx, Nexc, Ny, Ly, Nz, Npol, SG
READ(lu_in,grid)
IF(Nz .EQ. 1) & ! overwrite SG option if Nz = 1 for safety of use
@@ -375,7 +376,7 @@ CONTAINS
INTEGER :: i_, counter
IF(shear .GT. 0._dp) THEN
IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..'
- Lx = Ly/(2._dp*pi*shear*Npol)
+ Lx = Ly/(2._dp*pi*shear*Npol)*Nexc
ENDIF
Nkx = Nx;
! Local data
@@ -542,6 +543,7 @@ CONTAINS
CALL attach(fidres, TRIM(str), "jmaxi", jmaxi)
CALL attach(fidres, TRIM(str), "Nx", Nx)
CALL attach(fidres, TRIM(str), "Lx", Lx)
+ CALL attach(fidres, TRIM(str), "Nexc", Nexc)
CALL attach(fidres, TRIM(str), "Ny", Ny)
CALL attach(fidres, TRIM(str), "Ly", Ly)
CALL attach(fidres, TRIM(str), "Nz", Nz)
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index d185f381226c76dfe6be93f802ef95ed6368e21f..a7209fbe369bc221aa0be269fac88eab5af21cc0 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -24,10 +24,11 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp !
REAL(dp), PUBLIC, PROTECTED :: q_e = -1._dp ! Charge
REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp !
- REAL(dp), PUBLIC, PROTECTED :: K_Ne = 1._dp ! Density drive
- REAL(dp), PUBLIC, PROTECTED :: K_ni = 1._dp ! Density drive
- REAL(dp), PUBLIC, PROTECTED :: K_Te = 0._dp ! Temperature drive
- REAL(dp), PUBLIC, PROTECTED :: K_Ti = 0._dp ! Backg. electric field drive
+ REAL(dp), PUBLIC, PROTECTED :: k_N = 1._dp ! Ion density drive
+ REAL(dp), PUBLIC, PROTECTED :: eta_N = 1._dp ! electron-ion density drive ratio (k_Ne/k_Ni)
+ REAL(dp), PUBLIC, PROTECTED :: k_T = 0._dp ! Temperature drive
+ REAL(dp), PUBLIC, PROTECTED :: eta_T = 0._dp ! electron-ion temperature drive ratio (k_Te/k_Ti)
+ REAL(dp), PUBLIC, PROTECTED :: K_E = 0._dp ! Backg. electric field drive
REAL(dp), PUBLIC, PROTECTED :: GradB = 1._dp ! Magnetic gradient
REAL(dp), PUBLIC, PROTECTED :: CurvB = 1._dp ! Magnetic curvature
REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length
@@ -63,7 +64,7 @@ CONTAINS
NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, &
mu_x, mu_y, N_HD, mu_z, mu_p, mu_j, nu,&
tau_e, tau_i, sigma_e, sigma_i, q_e, q_i,&
- K_Ne, K_Ni, K_Te, K_Ti, GradB, CurvB, lambdaD, beta
+ k_N, eta_N, k_T, eta_T, K_E, GradB, CurvB, lambdaD, beta
READ(lu_in,model_par)
@@ -133,10 +134,11 @@ CONTAINS
CALL attach(fidres, TRIM(str), "sigma_i", sigma_i)
CALL attach(fidres, TRIM(str), "q_e", q_e)
CALL attach(fidres, TRIM(str), "q_i", q_i)
- CALL attach(fidres, TRIM(str), "K_Ne", K_Ne)
- CALL attach(fidres, TRIM(str), "K_Ni", K_Ni)
- CALL attach(fidres, TRIM(str), "K_Te", K_Te)
- CALL attach(fidres, TRIM(str), "K_Ti", K_Ti)
+ CALL attach(fidres, TRIM(str), "k_N", k_N)
+ CALL attach(fidres, TRIM(str), "eta_N", eta_N)
+ CALL attach(fidres, TRIM(str), "k_T", k_T)
+ CALL attach(fidres, TRIM(str), "eta_T", eta_T)
+ CALL attach(fidres, TRIM(str), "K_E", K_E)
CALL attach(fidres, TRIM(str), "lambdaD", lambdaD)
CALL attach(fidres, TRIM(str), "beta", beta)
END SUBROUTINE model_outputinputs
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 760f1505d77c6533cb50471d080f0568faf31520..092b390dfb8cad2170c7880bf3f2cab6f25b2bbd 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -300,7 +300,7 @@ SUBROUTINE add_Maxwellian_background_terms
! 40, 01,02, 21 with background gradient dependences.
USE prec_const
USE time_integration, ONLY : updatetlevel
- USE model, ONLY: taue_qe, taui_qi, K_Ne, K_Ni, K_Te, K_Ti, KIN_E
+ USE model, ONLY: taue_qe, taui_qi, k_N, eta_N, k_T, eta_T, KIN_E
USE array, ONLY: moments_rhs_e, moments_rhs_i
USE grid, ONLY: contains_kx0, contains_ky0, ikx_0, iky_0,&
ips_e,ipe_e,ijs_e,ije_e,ips_i,ipe_i,ijs_i,ije_i,&
@@ -320,28 +320,28 @@ SUBROUTINE add_Maxwellian_background_terms
SELECT CASE (ip-1)
CASE(0) ! Na00 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * (1.5_dp*K_Ne - 1.125_dp*K_Te)
+ +taue_qe * sinz(izs:ize) * (1.5_dp*eta_N*k_N - 1.125_dp*eta_N*k_T)
CASE(2) ! Na20 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*K_Ne - 2.75_dp*K_Te)
+ +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*eta_N*k_N - 2.75_dp*eta_T*k_T)
CASE(4) ! Na40 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*K_Te
+ +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*eta_T*k_T
END SELECT
CASE(1) ! j = 1
SELECT CASE (ip-1)
CASE(0) ! Na01 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taue_qe * sinz(izs:ize) * (K_Ne + 3.5_dp*K_Te)
+ -taue_qe * sinz(izs:ize) * (eta_N*k_N + 3.5_dp*eta_T*k_T)
CASE(2) ! Na21 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taue_qe * sinz(izs:ize) * SQRT2*K_Te
+ -taue_qe * sinz(izs:ize) * SQRT2*eta_T*k_T
END SELECT
CASE(2) ! j = 2
SELECT CASE (ip-1)
CASE(0) ! Na02 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * 2._dp*K_Te
+ +taue_qe * sinz(izs:ize) * 2._dp*eta_T*k_T
END SELECT
END SELECT
ENDDO
@@ -355,28 +355,28 @@ SUBROUTINE add_Maxwellian_background_terms
SELECT CASE (ip-1)
CASE(0) ! Na00 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * (1.5_dp*K_Ni - 1.125_dp*K_Ti)
+ +taui_qi * sinz(izs:ize) * (1.5_dp*k_N - 1.125_dp*k_T)
CASE(2) ! Na20 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*K_Ni - 2.75_dp*K_Ti)
+ +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*k_N - 2.75_dp*k_T)
CASE(4) ! Na40 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*K_Ti
+ +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*k_T
END SELECT
CASE(1) ! j = 1
SELECT CASE (ip-1)
CASE(0) ! Na01 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taui_qi * sinz(izs:ize) * (K_Ni + 3.5_dp*K_Ti)
+ -taui_qi * sinz(izs:ize) * (k_N + 3.5_dp*k_T)
CASE(2) ! Na21 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taui_qi * sinz(izs:ize) * SQRT2*K_Ti
+ -taui_qi * sinz(izs:ize) * SQRT2*k_T
END SELECT
CASE(2) ! j = 2
SELECT CASE (ip-1)
CASE(0) ! Na02 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * 2._dp*K_Ti
+ +taui_qi * sinz(izs:ize) * 2._dp*k_T
END SELECT
END SELECT
ENDDO
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 61db0f4f2398592e0367e1df672f372ec2f294d0..2f47015bd1759d7c2dd6aba23a167445aff08bb5 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -208,7 +208,8 @@ END SUBROUTINE evaluate_ampere_op
SUBROUTINE compute_lin_coeff
USE array
USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, &
- K_Te, K_Ti, K_Ne, K_Ni, CurvB, GradB, KIN_E, tau_e, tau_i, sigma_e, sigma_i
+ k_T, eta_T, k_N, eta_N, CurvB, GradB, KIN_E,&
+ tau_e, tau_i, sigma_e, sigma_i
USE prec_const
USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, &
ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i
@@ -304,11 +305,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij_e(ip,ij) =+K_Ne + 2.*j_dp*K_Te
- xphijp1_e(ip,ij) =-K_Te*(j_dp+1._dp)
- xphijm1_e(ip,ij) =-K_Te* j_dp
+ xphij_e(ip,ij) =+eta_N*k_N + 2.*j_dp*eta_T*k_T
+ xphijp1_e(ip,ij) =-eta_T*k_T*(j_dp+1._dp)
+ xphijm1_e(ip,ij) =-eta_T*k_T* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
- xphij_e(ip,ij) =+K_Te/SQRT2
+ xphij_e(ip,ij) =+eta_T*k_T/SQRT2
xphijp1_e(ip,ij) = 0._dp; xphijm1_e(ip,ij) = 0._dp;
ELSE
xphij_e(ip,ij) = 0._dp; xphijp1_e(ip,ij) = 0._dp
@@ -324,11 +325,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij_i(ip,ij) =+K_Ni + 2._dp*j_dp*K_Ti
- xphijp1_i(ip,ij) =-K_Ti*(j_dp+1._dp)
- xphijm1_i(ip,ij) =-K_Ti* j_dp
+ xphij_i(ip,ij) =+k_N + 2._dp*j_dp*k_T
+ xphijp1_i(ip,ij) =-k_T*(j_dp+1._dp)
+ xphijm1_i(ip,ij) =-k_T* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
- xphij_i(ip,ij) =+K_Ti/SQRT2
+ xphij_i(ip,ij) =+k_T/SQRT2
xphijp1_i(ip,ij) = 0._dp; xphijm1_i(ip,ij) = 0._dp;
ELSE
xphij_i(ip,ij) = 0._dp; xphijp1_i(ip,ij) = 0._dp
@@ -346,11 +347,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 1) THEN ! kronecker p1
- xpsij_e(ip,ij) =+(K_Ne + (2._dp*j_dp+1._dp)*K_Te)* SQRT(tau_e)/sigma_e
- xpsijp1_e(ip,ij) =- K_Te*(j_dp+1._dp) * SQRT(tau_e)/sigma_e
- xpsijm1_e(ip,ij) =- K_Te* j_dp * SQRT(tau_e)/sigma_e
+ xpsij_e(ip,ij) =+(eta_N*k_N + (2._dp*j_dp+1._dp)*eta_T*k_T) * SQRT(tau_e)/sigma_e
+ xpsijp1_e(ip,ij) =- eta_T*k_T*(j_dp+1._dp) * SQRT(tau_e)/sigma_e
+ xpsijm1_e(ip,ij) =- eta_T*k_T* j_dp * SQRT(tau_e)/sigma_e
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
- xpsij_e(ip,ij) =+ K_Te*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e
+ xpsij_e(ip,ij) =+ eta_T*k_T*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e
xpsijp1_e(ip,ij) = 0._dp; xpsijm1_e(ip,ij) = 0._dp;
ELSE
xpsij_e(ip,ij) = 0._dp; xpsijp1_e(ip,ij) = 0._dp
@@ -366,11 +367,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 1) THEN ! kronecker p1
- xpsij_i(ip,ij) =+(K_Ni + (2._dp*j_dp+1._dp)*K_Ti)* SQRT(tau_i)/sigma_i
- xpsijp1_i(ip,ij) =- K_Ti*(j_dp+1._dp) * SQRT(tau_i)/sigma_i
- xpsijm1_i(ip,ij) =- K_Ti* j_dp * SQRT(tau_i)/sigma_i
+ xpsij_i(ip,ij) =+(k_N + (2._dp*j_dp+1._dp)*k_T) * SQRT(tau_i)/sigma_i
+ xpsijp1_i(ip,ij) =- k_T*(j_dp+1._dp) * SQRT(tau_i)/sigma_i
+ xpsijm1_i(ip,ij) =- k_T* j_dp * SQRT(tau_i)/sigma_i
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
- xpsij_i(ip,ij) =+ K_Ti*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i
+ xpsij_i(ip,ij) =+ k_T*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i
xpsijp1_i(ip,ij) = 0._dp; xpsijm1_i(ip,ij) = 0._dp;
ELSE
xpsij_i(ip,ij) = 0._dp; xpsijp1_i(ip,ij) = 0._dp
diff --git a/src/restarts_mod.F90 b/src/restarts_mod.F90
index 581b5d271e938fd3a20c7e2d8c475ef9337c8ff0..0840222c54d7e81a8ea074edc181efa9acff9efd 100644
--- a/src/restarts_mod.F90
+++ b/src/restarts_mod.F90
@@ -12,7 +12,7 @@ IMPLICIT NONE
INTEGER :: rank, sz_, n_
INTEGER :: dims(1) = (/0/)
CHARACTER(LEN=50) :: dset_name
-INTEGER :: pmaxe_cp, jmaxe_cp, pmaxi_cp, jmaxi_cp, n0
+INTEGER :: pmaxe_cp, jmaxe_cp, pmaxi_cp, jmaxi_cp, n0, Nkx_cp, Nky_cp, Nz_cp
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_e_cp
COMPLEX(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: moments_i_cp
@@ -25,6 +25,9 @@ CONTAINS
!******************************************************************************!
SUBROUTINE load_moments
IMPLICIT NONE
+ REAL :: timer_tot_1, timer_find_CP_1, timer_load_mom_1
+ REAL :: timer_tot_2, timer_find_CP_2, timer_load_mom_2
+ CALL cpu_time(timer_tot_1)
! Checkpoint filename
WRITE(rstfile,'(a,a1,i2.2,a3)') TRIM(resfile0),'_',job2load,'.h5'
@@ -33,6 +36,9 @@ CONTAINS
! Open file
CALL openf(rstfile, fidrst,mpicomm=comm0)
! Get the checkpoint moments degrees to allocate memory
+ CALL getatt(fidrst,"/data/input/" , "Nkx", Nkx_cp)
+ CALL getatt(fidrst,"/data/input/" , "Nky", Nky_cp)
+ CALL getatt(fidrst,"/data/input/" , "Nz", Nz_cp)
IF (KIN_E) THEN
CALL getatt(fidrst,"/data/input/" , "pmaxe", pmaxe_cp)
CALL getatt(fidrst,"/data/input/" , "jmaxe", jmaxe_cp)
@@ -47,6 +53,8 @@ CONTAINS
CALL load_output_adapt_pj
ELSE
+ CALL cpu_time(timer_find_CP_1)
+
! Find the last results of the checkpoint file by iteration
n_ = n0+1
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_ ! start with moments_e/000001
@@ -67,18 +75,47 @@ CONTAINS
iframe2d = iframe2d-1; iframe5d = iframe5d-1
IF(my_id.EQ.0) WRITE(*,*) '.. restart from t = ', time
+ CALL cpu_time(timer_find_CP_2)
+ IF(my_id.EQ.0) WRITE(*,*) '** Time find CP : ', timer_find_CP_2 - timer_find_CP_1, ' **'
+
+ CALL cpu_time(timer_load_mom_1)
! Read state of system from checkpoint file
+
+ ! Super slow futils routine in Marconi.... but spare RAM
+ ! IF (KIN_E) THEN
+ ! WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
+ ! CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikys:ikye, ikxs:ikxe, izs:ize, 1),(/1,3,5/))
+ ! ENDIF
+ ! WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
+ ! CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikys:ikye, ikxs:ikxe, izs:ize, 1),(/1,3,5/))
+
+ ! Brute force loading: load the full moments and take what is needed (RAM dangerous...)
IF (KIN_E) THEN
- WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
- CALL getarrnd(fidrst, dset_name, moments_e(ips_e:ipe_e, ijs_e:ije_e, ikys:ikye, ikxs:ikxe, izs:ize, 1),(/1,3,5/))
+ CALL allocate_array(moments_e_cp,1,Nky_cp, 1,Nkx_cp, 1,Nz_cp, 1,pmaxe_cp+1, 1,jmaxe_cp+1)
+ WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_e", n_
+ CALL getarr(fidrst, dset_name, moments_e_cp(:,:,:,:,:))
+ moments_e(ips_e:ipe_e, ijs_e:ije_e, ikys:ikye, ikxs:ikxe, izs:ize, 1) &
+ = moments_e_cp(ips_e:ipe_e, ijs_e:ije_e, ikys:ikye, ikxs:ikxe, izs:ize)
+ DEALLOCATE(moments_e_cp)
ENDIF
+ CALL allocate_array(moments_i_cp, 1,pmaxi_cp+1, 1,jmaxi_cp+1, 1,Nky_cp, 1,Nkx_cp, 1,Nz_cp)
WRITE(dset_name, "(A, '/', i6.6)") "/data/var5d/moments_i", n_
- CALL getarrnd(fidrst, dset_name, moments_i(ips_i:ipe_i, ijs_i:ije_i, ikys:ikye, ikxs:ikxe, izs:ize, 1),(/1,3,5/))
+ CALL getarr(fidrst, dset_name, moments_i_cp(:,:,:,:,:))
+ moments_i(ips_i:ipe_i, ijs_i:ije_i, ikys:ikye, ikxs:ikxe, izs:ize, 1) &
+ = moments_i_cp(ips_i:ipe_i, ijs_i:ije_i, ikys:ikye, ikxs:ikxe, izs:ize)
+ DEALLOCATE(moments_i_cp)
CALL closef(fidrst)
IF (my_id .EQ. 0) WRITE(*,'(3x,a)') "Reading from restart file "//TRIM(rstfile)//" completed!"
ENDIF
+ CALL cpu_time(timer_load_mom_2)
+ CALL cpu_time(timer_tot_2)
+
+ CALL mpi_barrier(MPI_COMM_WORLD, ierr)
+
+ IF(my_id.EQ.0) WRITE(*,*) '** Time load mom : ', timer_load_mom_2 - timer_load_mom_1, ' **'
+ IF(my_id.EQ.0) WRITE(*,*) '** Total load time : ', timer_tot_2 - timer_tot_1, ' **'
END SUBROUTINE load_moments
!******************************************************************************!
diff --git a/wk/CBC_kT_PJ_scan.m b/wk/CBC_kT_PJ_scan.m
index 6d45978e37d91baf212cb7494426659e7c3572e1..11d27a379b34180c28edea1b68c4e73af6856aba 100644
--- a/wk/CBC_kT_PJ_scan.m
+++ b/wk/CBC_kT_PJ_scan.m
@@ -3,7 +3,7 @@ default_plots_options
HELAZDIR = '/home/ahoffman/HeLaZ/';
EXECNAME = 'helaz3';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-KT_a = [3:1.0:7];
+KT_a = [9:2:17];
g_max= KT_a*0;
g_avg= KT_a*0;
g_std= KT_a*0;
@@ -11,15 +11,15 @@ k_max= KT_a*0;
CO = 'DG'; GKCO = 0;
NU = 0.01;
-DT = 4e-3;
-TMAX = 50;
+DT = 1e-2;
+TMAX = 25;
ky_ = 0.3;
SIMID = 'linear_CBC_kT_scan_ky_0.3'; % Name of the simulation
RUN = 0;
figure
-P = 12;
-for J = [0 1 2 4 6]
-% J = P/2;
+% P = 12;
+for P = [2 4 6]
+J = P/2;
i=1;
for K_T = KT_a
diff --git a/wk/CBC_ky_PJ_scan.m b/wk/CBC_ky_PJ_scan.m
new file mode 100644
index 0000000000000000000000000000000000000000..d9bc9e0516236adba01b11d213b395df1063d353
--- /dev/null
+++ b/wk/CBC_ky_PJ_scan.m
@@ -0,0 +1,122 @@
+addpath(genpath('../matlab')) % ... add
+default_plots_options
+HELAZDIR = '/home/ahoffman/HeLaZ/';
+EXECNAME = 'helaz3';
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+KY_a = 0.1:0.1:0.8;
+g_max= KY_a*0;
+g_avg= g_max*0;
+g_std= g_max*0;
+k_max= g_max*0;
+
+CO = 'DG'; GKCO = 0;
+NU = 0.05;
+DT = 2e-3;
+TMAX = 25;
+K_T = 6.96;
+SIMID = 'linear_CBC_circ_conv'; % Name of the simulation
+RUN = 0;
+% figure
+% P = 12;
+for P = [12]
+J = P/2;
+
+i=1;
+for ky_ = KY_a
+
+%Set Up parameters
+for j = 1
+ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
+ TAU = 1.0; % e/i temperature ratio
+ K_N = 2.22; K_Ne = K_N;
+ K_Te = K_T; % Temperature '''
+ SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
+ KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+ BETA = 0e-1; % electron plasma beta
+ PMAXE = P; JMAXE = J;
+ PMAXI = P; JMAXI = J;
+ NX = 12; % real space x-gridpoints
+ NY = 2; % '' y-gridpoints
+ LX = 2*pi/0.1; % Size of the squared frequency domain
+ LY = 2*pi/ky_;
+ NZ = 16; % number of perpendicular planes (parallel grid)
+ NPOL = 1; SG = 0;
+% GEOMETRY= 's-alpha';
+ GEOMETRY= 'circular';
+ Q0 = 1.4; % safety factor
+ SHEAR = 0.8; % magnetic shear (Not implemented yet)
+ EPS = 0.18; % inverse aspect ratio
+ SPS0D = 1; SPS2D = 0; SPS3D = 1;SPS5D= 1/5; SPSCP = 0;
+ JOB2LOAD= -1;
+ LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
+ 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= 'phi'; % Start simulation with a noisy mom00/phi/allmom
+ 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;
+ 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;
+end
+%%-------------------------------------------------------------------------
+% RUN
+setup
+if RUN
+ system(['cd ../results/',SIMID,'/',PARAMS,'/; mpirun -np 6 ',HELAZDIR,'bin/',EXECNAME,' 2 1 3 0; cd ../../../wk'])
+end
+
+% Load results
+filename = [SIMID,'/',PARAMS,'/'];
+LOCALDIR = [HELAZDIR,'results/',filename,'/'];
+data = compile_results(LOCALDIR,0,0); %Compile the results from first output found to JOBNUMMAX if existing
+
+%linear growth rate (adapted for 2D zpinch and fluxtube)
+trange = [0.5 1]*data.Ts3D(end);
+nplots = 0;
+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);
+
+ g_max(i) = gmax;
+ k_max(i) = kmax;
+
+ g_avg(i) = lg.avg_g;
+ g_std(i) = lg.std_g;
+
+ i = i + 1;
+end
+%%
+
+% plot(KT_a,max(g_max,0));
+y_ = g_avg;
+e_ = g_std;
+
+y_ = y_.*(y_-e_>0);
+e_ = e_ .* (y_>0);
+% errorbar(KY_a,y_,e_,...
+% 'LineWidth',1.2,...
+% 'DisplayName',['(',num2str(P),',',num2str(J),')']);
+% hold on;
+% title(['Linear CBC $K_T$ threshold $k_y=$',num2str(ky_),' (CLOS = 1)']);
+% legend('show'); xlabel('$k_y$'); ylabel('$\gamma$');
+% drawnow
+
+fig = plot_ballooning(data,options);
+
+end
+
diff --git a/wk/CBC_nu_CO_scan.m b/wk/CBC_nu_CO_scan.m
index 949847b83d6628c21e9e32a3bfd42ca1a093e36a..d14d69ba31a7db9aad1fe96ad5f9474477d7e13d 100644
--- a/wk/CBC_nu_CO_scan.m
+++ b/wk/CBC_nu_CO_scan.m
@@ -1,20 +1,21 @@
% NU_a = [0.05 0.15 0.25 0.35 0.45];
-NU_a = [0:0.05:0.5];
+NU_a = [0:0.1:0.5];
g_max= NU_a*0;
g_avg= NU_a*0;
g_std= NU_a*0;
k_max= NU_a*0;
-CO = 'LD';
+CO = 'DG';
-K_T = 5.3;
-DT = 2e-3;
-TMAX = 30;
+K_T = 7;
+DT = 5e-3;
+TMAX = 20;
ky_ = 0.3;
-SIMID = 'linear_CBC_nu_scan_kT_5.3_ky_0.3_LDGK'; % Name of the simulation
-RUN = 0;
+SIMID = 'linear_CBC_nu_scan_kT_7_ky_0.3_DGGK'; % Name of the simulation
+% SIMID = 'linear_CBC_nu_scan_kT_11_ky_0.3_DGGK'; % Name of the simulation
+RUN = 1;
figure
-for P = [2 4 6]
+for P = [6 8 10]
i=1;
for NU = NU_a
diff --git a/wk/Dimits_2000_fig3.m b/wk/Dimits_2000_fig3.m
index 30103b346b8f1bd8ec959a0d19b3eee2de799060..13501cd86b7b315a4a4fef5585ac43d6347076c8 100644
--- a/wk/Dimits_2000_fig3.m
+++ b/wk/Dimits_2000_fig3.m
@@ -1,58 +1,64 @@
%% Heat flux Qi [R/rhos^2/cs]
kN = 2.22;
%-------------- GM ---------------
-%192x96x16x3x2 kymin=0.05
+%(P,J)=(2,1)
kT_Qi_GM_32 = ...
[...
- 9.0 1.0e+2 4.3e+1;...
- 7.0 3.6e+1 0;...
- 6.0 2.6e+1 0;...
- 5.0 1.1e+1 0;...
- 4.5 7.8e+0 0;...
+ 13. 1.5e+2 1.6e+1;...%192x96x16x3x2 kymin=0.02
+ 11. 5.1e+2 3.5e+2;...%192x96x16x3x2 kymin=0.02
+ 9.0 9.6e+1 3.0e+1;...%192x96x16x3x2 kymin=0.05
+ 7.0 5.0e+1 6.6e+0;...%192x96x16x3x2 kymin=0.05
+ 6.0 3.0e+1 4.8e+0;...%192x96x16x3x2 kymin=0.05
+ 5.0 1.1e+1 9.4e-1;...%192x96x16x3x2 kymin=0.05
+ 4.5 9.2e+0 1.6e+0;...%192x96x16x3x2 kymin=0.05
];
-%128x64x16x5x3 kymin=0.05
+%(P,J)=(4,2)
kT_Qi_GM_53 = ...
[...
- 13. 1.3e+2 3.5e+1;...
- 11. 9.7e+1 2.2e+1;...
- 9.0 9.0e+1 2.8e+1;...
- 7.0 4.6e+1 0;...
- 5.3 1.9e+1 0;...
- 5.0 1.5e+1 0;...
- 4.5 9.7e+0 0;...9_128
+ 13. 2.0e+2 1.2e+1;...%96x64x16x3x2 kymin=0.02 (large box)
+% 13. 1.1e+2 2.0e+1;...%128x64x16x5x3 kymin=0.02 (large box)
+% 13. 1.3e+2 3.5e+1;...%128x64x16x5x3 kymin=0.05
+ 11. 1.2e+2 1.6e+1;...%96x64x16x3x2 kymin=0.02 (large box)
+% 11. 1.6e+2 1.8e+1;...%128x64x16x5x3 kymin=0.02
+% 11. 9.7e+1 2.2e+1;...%128x64x16x5x3 kymin=0.05
+ 9.0 8.3e+1 2.2e+1;...%128x64x16x5x3 kymin=0.05
+% 9.0 7.6e+1 2.3e+1;...%256x128x16x3x2 kymin=0.05 (high res)
+ 7.0 4.6e+1 2.3e+0;...%128x64x16x5x3 kymin=0.05
+ 6.0 3.7e+1 6.9e+0;...%128x64x16x5x3 kymin=0.05
+ 5.3 1.9e+1 2.0e+0;...%128x64x16x5x3 kymin=0.05
+ 5.0 1.3e+1 3.3e+0;...%128x64x16x5x3 kymin=0.05
+ 4.5 9.3e+0 1.0e+0;...%128x64x16x5x3 kymin=0.05
];
-%128x64x16x5x3 kymin=0.02 (large box)
-kT_Qi_GM_53_LB = ...
+%(P,J)=(12,2) or higher
+kT_Qi_GM_122 = ...
[...
- 13. 1.1e+2 2.0e+1;...
+ 7.0 4.1e+1 6.6e+0;...%192x96x24x13x7 kymin=0.05
+ 4.5 9.6e-1 1.5e-1;...%128x64x16x13x2 kymin=0.05
];
%-------------- GENE ---------------
-%128x64x16x24x12 kymin=0.05
kT_Qi_GENE = ...
[...
- 13. 2.0e+2 6.6e+1;...
- 11. 3.3e+2 1.6e+2;...
- 9.0 1.1e+2 0;...
- 7.0 5.0e+1 0;...
- 5.3 2.4e+1 0;...
- 4.5 1.9e-1 0;...
- ];
-%128x64x16x24x12 kymin=0.02 (large box)
-kT_Qi_GM_53_LB = ...
- [...
- 13. 2.7e+2 2.2e+1;...
- 11. 1.9e+2 1.7e+1;...
+ 13. 2.7e+2 2.2e+1;...%128x64x16x24x12 kymin=0.02 (large box)
+% 13. 2.0e+2 6.6e+1;...%128x64x16x24x12 kymin=0.05
+ 11. 1.9e+2 1.7e+1;...%128x64x16x24x12 kymin=0.02 (large box)
+% 11. 3.3e+2 1.6e+2;...%128x64x16x24x12 kymin=0.05
+ 9.0 1.1e+2 4.2e+1;...%128x64x16x24x12 kymin=0.05
+ 7.0 4.1e+1 2.1e+1;...%128x64x16x24x12 kymin=0.05
+ 5.3 1.1e+1 1.8e+1;...%128x64x16x24x12 kymin=0.05
+ 4.5 1.9e-1 3.0e-2;...%128x64x16x24x12 kymin=0.05
];
%% Heat conductivity Xi [Ln/rhoi^2/cs] computed as Xi = Qi/kT/kN
%init
-kT_Xi_GM_32 = kT_Qi_GM_32;
-kT_Xi_GM_53 = kT_Qi_GM_53;
-kT_Xi_GENE = kT_Qi_GENE;
+kT_Xi_GM_32 = kT_Qi_GM_32;
+kT_Xi_GM_53 = kT_Qi_GM_53;
+kT_Xi_GM_122 = kT_Qi_GM_122;
+kT_Xi_GENE = kT_Qi_GENE;
%scale
for i = 2:3
-kT_Xi_GM_32(:,i) = kT_Qi_GM_32(:,i)./kT_Qi_GM_32(:,1)./kN;
-kT_Xi_GM_53(:,i) = kT_Qi_GM_53(:,i)./kT_Qi_GM_53(:,1)./kN;
-kT_Xi_GENE (:,i) = kT_Qi_GENE (:,i)./kT_Qi_GENE (:,1)./kN;
+kT_Xi_GM_32 (:,i) = kT_Qi_GM_32 (:,i)./kT_Qi_GM_32 (:,1)./kN;
+kT_Xi_GM_53 (:,i) = kT_Qi_GM_53 (:,i)./kT_Qi_GM_53 (:,1)./kN;
+kT_Xi_GM_122(:,i) = kT_Qi_GM_122(:,i)./kT_Qi_GM_122(:,1)./kN;
+kT_Xi_GENE (:,i) = kT_Qi_GENE (:,i)./kT_Qi_GENE (:,1)./kN;
end
%% Dimits fig 3 data
KT_DIM = [4.0 4.5 5.0 6.0 7.0 9.0 12. 14. 16. 18.];
@@ -75,14 +81,21 @@ GFL__98_DIM = [5.0 1.5;...
7.0 5.0;...
9.0 7.5];
%% Plot
+msz = 8.0; lwt = 2.0;
figure;
if 1
xye = kT_Xi_GM_32;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'o-','DisplayName','192x96x16x3x2','LineWidth',2.0); hold on
+errorbar(xye(:,1), xye(:,2),xye(:,3),'>-','DisplayName','192x96x16x3x2',...
+ 'MarkerSize',msz,'LineWidth',lwt); hold on
xye = kT_Xi_GM_53;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'o-','DisplayName','128x64x16x5x3','LineWidth',2.0); hold on
+errorbar(xye(:,1), xye(:,2),xye(:,3),'<-','DisplayName','128x64x16x5x3',...
+ 'MarkerSize',msz,'LineWidth',lwt); hold on
+xye = kT_Xi_GM_122;
+errorbar(xye(:,1), xye(:,2),xye(:,3),'^-','DisplayName','128x64x16x13x3',...
+ 'MarkerSize',msz,'LineWidth',lwt); hold on
xye = kT_Xi_GENE;
-errorbar(xye(:,1), xye(:,2),xye(:,3),'v-.','DisplayName','GENE 128x64x16x24x12','LineWidth',2.0);
+errorbar(xye(:,1), xye(:,2),xye(:,3),'+-.k','DisplayName','GENE 128x64x16x24x12',...
+ 'MarkerSize',msz,'LineWidth',lwt); hold on
end
if 1
plot(LLNL_GK_DIM(:,1),LLNL_GK_DIM(:,2),'dk--','DisplayName','Dimits GK, LLNL'); hold on
diff --git a/wk/analysis_HeLaZ.m b/wk/analysis_HeLaZ.m
index e3eb440df533963cc25d5345873d1f737bd9a530..dfa4d3a87c1e4aa5533c8e02f01c94a0edd77d94 100644
--- a/wk/analysis_HeLaZ.m
+++ b/wk/analysis_HeLaZ.m
@@ -22,8 +22,10 @@ FMT = '.fig';
if 1
%% Space time diagramm (fig 11 Ivanov 2020)
-options.TAVG_0 = 0.5*data.Ts3D(end); data.scale = 1;%/(data.Nx*data.Ny)^2;
-options.TAVG_1 = data.Ts3D(end); % Averaging times duration
+% data.scale = 1;%/(data.Nx*data.Ny)^2;
+options.TAVG_0 = 25;%0.4*data.Ts3D(end);
+options.TAVG_1 = 40;%0.9*data.Ts3D(end); % Averaging times duration
+options.NCUT = 4; % Number of cuts for averaging and error estimation
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)
@@ -37,7 +39,7 @@ if 0
options.T = [200 400];
fig = statistical_transport_averaging(data,options);
end
-
+
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
@@ -49,13 +51,13 @@ options.NAME = '\phi';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
% options.TIME = data.Ts5D(end-30:end);
-options.TIME = data.Ts3D(1:2:end);
-% options.TIME = [550:2:750];
+options.TIME = data.Ts3D;
+% options.TIME = [850:0.1:1000];
data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
@@ -74,11 +76,11 @@ options.NAME = '\psi';
% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xz';
+options.PLAN = 'kxky';
% options.NAME 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [550 650];
+options.TIME = [400 440];
data.a = data.EPS * 2e3;
fig = photomaton(data,options);
% save_figure(data,fig)
@@ -98,13 +100,13 @@ end
if 0
%% Kinetic distribution function sqrt(<f_a^2>xy) (GENE vsp)
-options.SPAR = linspace(-3,3,32)+(6/127/2);
-options.XPERP = linspace( 0,6,32);
-% options.SPAR = gene_data.vp';
-% options.XPERP = gene_data.mu';
+% options.SPAR = linspace(-3,3,32)+(6/127/2);
+% options.XPERP = linspace( 0,6,32);
+options.SPAR = gene_data.vp';
+options.XPERP = gene_data.mu';
options.iz = 'avg';
-options.T = [40];
-options.PLT_FCT = 'pcolor';
+options.T = [300 600];
+options.PLT_FCT = 'contour';
options.ONED = 0;
options.non_adiab = 0;
options.SPECIE = 'e';
@@ -131,10 +133,10 @@ end
if 0
%% Time averaged spectrum
-options.TIME = [100 500];
+options.TIME = [300 600];
options.NORM =1;
-% options.NAME = '\phi';
-options.NAME = 'N_i^{00}';
+options.NAME = '\phi';
+% options.NAME = 'N_i^{00}';
% options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
diff --git a/wk/analysis_gene.m b/wk/analysis_gene.m
index b5eb9436cc656d444edfc11847e523c00ca4aa92..62fb9452a8f0e648191e022e0255c4f74dfddab0 100644
--- a/wk/analysis_gene.m
+++ b/wk/analysis_gene.m
@@ -18,14 +18,15 @@ addpath(genpath([helazdir,'matlab/load'])) % ... add
folder = '/misc/gene_results/LD_zpinch_1.6/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
% folder = '/misc/gene_results/ZP_HP_kn_1.6_nuv_3.2/';
+% folder = '/misc/gene_results/Z-pinch/ZP_HP_kn_1.6_HRES/';
% folder = '/misc/gene_results/ZP_kn_2.5_large_box/';
% folder = '/misc/gene_results/CBC/128x64x16x24x12/';
-% folder = '/misc/gene_results/CBC/196x96x20x32x16_00/';
+% folder = '/misc/gene_results/CBC/196x96x20x32x16_01/';
% folder = '/misc/gene_results/CBC/128x64x16x6x4/';
-% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_00/';
-% folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
-% folder = '/misc/gene_results/CBC/KT_13_128x64x16x24x12/';
-% folder = '/misc/gene_results/CBC/KT_13_large_box_128x64x16x24x12/';
+% folder = '/misc/gene_results/CBC/KT_5.3_128x64x16x24x12_01/';
+folder = '/misc/gene_results/CBC/KT_4.5_128x64x16x24x12_01/';
+% folder = '/misc/gene_results/CBC/KT_11_128x64x16x24x12/';
+% folder = '/misc/gene_results/CBC/KT_11_large_box_128x64x16x24x12/';
gene_data = load_gene_data(folder);
gene_data = invert_kxky_to_kykx_gene_results(gene_data);
if 1
@@ -57,11 +58,11 @@ options.NAME = '\phi';
% options.NAME = 'n_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'kxky';
% options.NAME ='f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [20 30 40 50 60];
+options.TIME = [325 400];
gene_data.a = data.EPS * 2000;
fig = photomaton(gene_data,options);
save_figure(gene_data,fig,'.png')
@@ -70,7 +71,7 @@ end
if 0
%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Options
-options.INTERP = 1;
+options.INTERP = 1;
options.POLARPLOT = 0;
options.NAME = '\phi';
% options.NAME = 'v_y';
@@ -114,11 +115,11 @@ end
if 0
%% Show f_i(vpar,mu)
-options.times = 20:80;
+options.times = 250:500;
options.specie = 'i';
options.PLT_FCT = 'contour';
options.folder = folder;
-options.iz = 1;
+options.iz = 'avg';
options.FIELD = '<f_>';
options.ONED = 0;
% options.FIELD = 'Q_es';
@@ -129,9 +130,9 @@ if 0
%% Time averaged spectrum
options.TIME = 300:600;
options.NORM =1;
-% options.NAME = '\phi';
+options.NAME = '\phi';
% options.NAME = 'n_i';
-options.NAME ='\Gamma_x';
+% options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
options.COMPZ = 'avg';
options.OK = 0;
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index d67101a92705bd14f9d2f5fe5f2d9c8a319b25dc..5ccf79448facc6d60186d0f25e437777194955b3 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -41,16 +41,23 @@ helazdir = '/home/ahoffman/HeLaZ/';
% outfile = 'CBC/64x32x16x5x3';
% outfile = 'CBC/64x128x16x5x3';
% outfile = 'CBC/128x64x16x5x3';
+% outfile = 'CBC/128x96x16x3x2_Nexc_6';
% outfile = 'CBC/192x96x16x3x2';
-% outfile = 'CBC/128x64x16x5x3';
-% outfile = 'CBC/kT_9_128x64x16x5x3';
-outfile = 'CBC/kT_4.5_128x64x16x13x7';
+% outfile = 'CBC/192x96x24x13x7';
+% outfile = 'CBC/kT_11_128x64x16x5x3';
+% outfile = 'CBC/kT_9_256x128x16x3x2';
+% outfile = 'CBC/kT_4.5_128x64x16x13x2';
+% outfile = 'CBC/kT_4.5_192x96x24x13x7';
+% outfile = 'CBC/kT_5.3_192x96x24x13x7';
% outfile = 'CBC/kT_13_large_box_128x64x16x5x3';
-JOBNUMMIN = 00; JOBNUMMAX = 06;
+outfile = 'CBC/kT_13_96x96x16x3x2_Nexc_6';
+% outfile = 'CBC/kT_11_96x64x16x5x3_ky_0.02';
% outfile = 'CBC/kT_scan_128x64x16x5x3';
% outfile = 'CBC/kT_scan_192x96x16x3x2';
%% 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';
+
+JOBNUMMIN = 00; JOBNUMMAX = 20;
run analysis_HeLaZ
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 47ad4a5bd89c139519f1a592907a003501c78bcd..23acafa52edc071b2d22bef29e55ef0881c648d0 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -3,7 +3,9 @@
% from matlab framework. It is meant to run only small problems in linear
% for benchmark and debugging purpose since it makes matlab "busy"
%
-SIMID = 'Kinetic_ballooning_mode'; % Name of the simulation
+% SIMID = 'test_circular_geom'; % Name of the simulation
+% SIMID = 'linear_CBC'; % Name of the simulation
+SIMID = 'dbg'; % Name of the simulation
RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
@@ -23,20 +25,20 @@ K_Te = 4.5; % Temperature '''
K_Ti = 8.0;%0.25*K_N; % Temperature '''
SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
% 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.03; % electron plasma beta
+KIN_E = 0; % 1: kinetic electrons, 2: adiabatic electrons
+BETA = 0.0; % electron plasma beta
%% GRID PARAMETERS
-P = 8;
-J = P/2;
+P = 4;
+J = 2;%P/2;
PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
-NX = 20; % real space x-gridpoints
-NY = 2; % '' y-gridpoints
+NX = 32; % real space x-gridpoints
+NY = 16; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
-LY = 2*pi/0.25; % Size of the squared frequency domain
-NZ = 32; % number of perpendicular planes (parallel grid)
+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
@@ -44,11 +46,12 @@ SG = 0; % Staggered z grids option
GEOMETRY= 's-alpha';
% GEOMETRY= 'circular';
Q0 = 1.4; % safety factor
-SHEAR = 0.8; % magnetic shear (Not implemented yet)
-EPS = 0.18; % inverse aspect ratio
+SHEAR = 0.8; % magnetic shear
+NEXC = 6; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
+EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 40; % Maximal time unit
-DT = 2e-3; % Time step
+TMAX = 25; % 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 = 1; % Sampling per time unit for 2D arrays
@@ -126,7 +129,7 @@ end
if 1
%% Ballooning plot
options.time_2_plot = [120];
-options.kymodes = [0.25];
+options.kymodes = [0.1];
options.normalized = 1;
% options.field = 'phi';
fig = plot_ballooning(data,options);
@@ -182,4 +185,4 @@ 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=0.00$',data.kx(ikx)))
-end
\ No newline at end of file
+end