diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 6f5e6402419bf5c5293139d0df894f0635ec0ccd..c4687314a67ee23d4956e615b9cdbff69723bb7e 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -83,7 +83,9 @@ mvm = @(x) movmean(x,OPTIONS.NMVA);
else
plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',...
'DisplayName',['$\Gamma_{avg}=',sprintf('%2.2f',Gx_infty_avg),'\pm',sprintf('%2.2f',Gx_infty_std),'$']); legend('show');
- ylim([0,5*abs(Gx_infty_avg)]);
+ try ylim([0,5*abs(Gx_infty_avg)]);
+ catch
+ end
end
xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
grid on; set(gca,'xticklabel',[]);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 29304ef0a05df30e763f385b202b102b8142eaa4..adc09fc54c9eb18de2b0bd80b53d967fc9249305 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -17,12 +17,12 @@ 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']);
fprintf(fid,[' Npol = ', num2str(GRID.Npol),'\n']);
fprintf(fid,[' SG = ', GRID.SG,'\n']);
+fprintf(fid,[' Nexc = ', num2str(GRID.Nexc),'\n']);
fprintf(fid,'/\n');
fprintf(fid,'&GEOMETRY\n');
@@ -73,8 +73,8 @@ fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']);
fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']);
fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']);
fprintf(fid,[' K_Ne = ', num2str(MODEL.K_Ne),'\n']);
-fprintf(fid,[' K_Ni = ', num2str(MODEL.K_Ni),'\n']);
fprintf(fid,[' K_Te = ', num2str(MODEL.K_Te),'\n']);
+fprintf(fid,[' K_Ni = ', num2str(MODEL.K_Ni),'\n']);
fprintf(fid,[' K_Ti = ', num2str(MODEL.K_Ti),'\n']);
fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']);
fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index f1869e91a2a5571063923d604acd3be588437944..912baa9e7d45a1505aa3e5f8370d28f0b83f32a0 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -3,6 +3,7 @@ SUBROUTINE diagnose(kstep)
USE basic
USE diagnostics_par
+ USE processing, ONLY: gflux_ri, hflux_xi
IMPLICIT NONE
INTEGER, INTENT(in) :: kstep
@@ -16,7 +17,8 @@ SUBROUTINE diagnose(kstep)
IF (kstep .GE. 0) THEN
! Terminal info
IF (MOD(cstep, INT(1.0/dt)) == 0 .AND. (my_id .EQ. 0)) THEN
- WRITE(*,"(F6.0,A,F6.0)") time,"/",tmax
+ ! WRITE(*,"(F6.0,A,F6.0)") time,"/",tmax
+ WRITE(*,"(A,F6.0,A1,F6.0,A8,G10.2,A8,G10.2,A)")'|t/tmax = ', time,"/",tmax,'| Gxi = ',gflux_ri,'| Qxi = ',hflux_xi,'|'
ENDIF
ELSEIF (kstep .EQ. -1) THEN
CALL cpu_time(finish)
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 06e47107275bd6b75900abb8f7aa833d8f8b7d39..191b319f872b6c5cffaef65faeaab6dedb804110 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -126,7 +126,7 @@ CONTAINS
INTEGER :: lu_in = 90 ! File duplicated from STDIN
NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
- Nx, Lx, Nexc, Ny, Ly, Nz, Npol, SG
+ Nx, Lx, Ny, Ly, Nz, Npol, SG, Nexc
READ(lu_in,grid)
IF(Nz .EQ. 1) & ! overwrite SG option if Nz = 1 for safety of use
@@ -517,7 +517,8 @@ CONTAINS
IF(zarray_full(iz) .LT. zmin) zmin = zarray_full(iz)
END DO
!! Parallel data distribution
- IF( num_procs_z .GT. 1) stop '>>STOPPED<< Cannot have multiple core in z-direction (Nz = 1)'
+ IF( (Nz .EQ. 1) .AND. (num_procs_z .GT. 1) ) &
+ stop '>>STOPPED<< Cannot have multiple core in z-direction (Nz = 1)'
! Local data distribution
CALL decomp1D(total_nz, num_procs_z, rank_z, izs, ize)
local_nz = ize - izs + 1
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 25e7f573b99eac5ca4692cdf66ab5994e5c887b9..d866c947018af203826be4c286909b1eadfc1259 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -140,12 +140,12 @@ SUBROUTINE moments_eq_rhs_e
! Nonlinear term
-hatB_NL(iz,eo) * Sepj(ip,ij,iky,ikx,iz)
- IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) &
- ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
- ! (not used often so not parallelized)
- moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
- moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) &
- + mu_p * moments_e(ip-4,ij,iky,ikx,iz,updatetlevel)
+ ! IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) &
+ ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
+ ! ! (not used often so not parallelized)
+ ! moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = &
+ ! moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) &
+ ! + mu_p * moments_e(ip-4,ij,iky,ikx,iz,updatetlevel)
ELSE
moments_rhs_e(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
@@ -291,12 +291,12 @@ SUBROUTINE moments_eq_rhs_i
! Nonlinear term with a (gxx*gxy - gxy**2)^1/2 factor
-hatB_NL(iz,eo) * Sipj(ip,ij,iky,ikx,iz)
- IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) &
- ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
- ! (not used often so not parallelized)
- moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
- moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) &
- + mu_p * moments_i(ip-4,ij,iky,ikx,iz,updatetlevel)
+ ! IF( (ip-4 .GT. 0) .AND. (num_procs_p .EQ. 1) ) &
+ ! ! Numerical parallel velocity hyperdiffusion "+ dvpar4 g_a" see Pueschel 2010 (eq 33)
+ ! ! (not used often so not parallelized)
+ ! moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = &
+ ! moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) &
+ ! + mu_p * moments_i(ip-4,ij,iky,ikx,iz,updatetlevel)
ELSE
moments_rhs_i(ip,ij,iky,ikx,iz,updatetlevel) = 0._dp
ENDIF
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index 7203bdf0538ec784de84dbf4cc3ad74f2f32d93f..12b829109a1306f72778851be780902fe3bd667a 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -62,10 +62,12 @@ CONTAINS
SELECT CASE (numerical_scheme)
CASE ('RK4')
CALL RK4
+ CASE ('SSPx_RK3')
+ CALL SSPx_RK3
+ CASE ('SSPx_RK2')
+ CALL SSPx_RK2
CASE ('DOPRI5')
- CALL DOPRI5
- CASE ('DOPRI5_ADAPT')
- CALL DOPRI5_ADAPT
+ CALL DOPRI5
CASE DEFAULT
IF (my_id .EQ. 0) WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
END SELECT
@@ -104,7 +106,73 @@ CONTAINS
END SUBROUTINE RK4
+ SUBROUTINE SSPx_RK3
+ ! Butcher coeff for modified strong stability preserving RK3
+ ! used in GX (Hammett 2022, Mandell et al. 2022)
+ USE basic
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER,PARAMETER :: nbstep = 3
+ REAL(dp) :: a1, a2, a3, w1, w2, w3
+ a1 = (1._dp/6._dp)**(1._dp/3._dp)! (1/6)^(1/3)
+ ! a1 = 0.5503212081491044571635029569733887910843_dp ! (1/6)^(1/3)
+ a2 = a1
+ a3 = a1
+ w1 = 0.5_dp*(-1._dp + SQRT( 9._dp - 2._dp * 6._dp**(2._dp/3._dp))) ! (-1 + sqrt(9-2*6^(2/3)))/2
+ ! w1 = 0.2739744023885328783052273138309828937054_dp ! (sqrt(9-2*6^(2/3))-1)/2
+ w2 = 0.5_dp*(-1._dp + 6._dp**(2._dp/3._dp) - SQRT(9._dp - 2._dp * 6._dp**(2._dp/3._dp))) ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
+ ! w2 = 0.3769892220587804931852815570891834795475_dp ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
+ w3 = 1._dp/a1 - w2 * (1._dp + w1)
+ ! w3 = 1.3368459739528868457369981115334667265415_dp
+
+ CALL allocate_array_dp1(c_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
+
+ ntimelevel = 3
+
+ c_E(1) = 0.0_dp
+ c_E(2) = 1.0_dp/2.0_dp
+ c_E(3) = 1.0_dp/2.0_dp
+
+ b_E(1) = a1 * (w1*w2 + w3)
+ b_E(2) = a2 * w2
+ b_E(3) = a3
+
+ A_E(2,1) = a1
+ A_E(3,1) = a1 * w1
+ A_E(3,2) = a2
+
+ END SUBROUTINE SSPx_RK3
+
+ SUBROUTINE SSPx_RK2
+ ! Butcher coeff for modified strong stability preserving RK2
+ ! used in GX (Hammett 2022, Mandell et al. 2022)
+
+ USE basic
+ USE prec_const
+ IMPLICIT NONE
+ INTEGER,PARAMETER :: nbstep = 2
+ REAL(dp) :: alpha, beta
+ alpha = 1._dp/SQRT(2._dp)
+ beta = SQRT(2._dp) - 1._dp
+
+ CALL allocate_array_dp1(c_E,1,nbstep)
+ CALL allocate_array_dp1(b_E,1,nbstep)
+ CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
+
+ ntimelevel = 2
+
+ c_E(1) = 0.0_dp
+ c_E(2) = 1.0_dp/2.0_dp
+
+ b_E(1) = alpha*beta/2._dp
+ b_E(2) = alpha/2._dp
+
+ A_E(2,1) = alpha
+
+ END SUBROUTINE SSPx_RK2
SUBROUTINE DOPRI5
! Butcher coeff for DOPRI5 --> Stiffness detection
diff --git a/wk/analysis_gyacomo.m b/wk/analysis_gyacomo.m
index 8d494f0aac22af2be8841637e8ef7fbe735d1021..e637d9dd64a137c983b4c1b7b083cd5039c7cd1b 100644
--- a/wk/analysis_gyacomo.m
+++ b/wk/analysis_gyacomo.m
@@ -10,7 +10,7 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add
%% Load the results
LOCALDIR = [gyacomodir,resdir,'/'];
-MISCDIR = ['/misc/gyacomo_outputs/',resdir,'/']; %For long term storage
+MISCDIR = ['/misc/',resdir,'/']; %For long term storage
system(['mkdir -p ',MISCDIR]);
system(['mkdir -p ',LOCALDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
@@ -62,14 +62,14 @@ if 0
% Options
options.INTERP = 1;
options.POLARPLOT = 0;
-options.NAME = '\phi';
+% options.NAME = '\phi';
% options.NAME = '\omega_z';
-% options.NAME = 'N_i^{00}';
+options.NAME = 'N_e^{00}';
% options.NAME = 'v_x';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-options.PLAN = 'kxky';
+options.PLAN = 'xy';
% options.NAME = 'f_i';
% options.PLAN = 'sx';
options.COMP = 'avg';
diff --git a/wk/header_3D_results.m b/wk/header_3D_results.m
index 8a872def7c8ecc1c43bf6be710d72542bbaf9a4c..5e9a28c9d7d6d7dcfb740205eddcb8a67b992117 100644
--- a/wk/header_3D_results.m
+++ b/wk/header_3D_results.m
@@ -35,7 +35,7 @@ gyacomodir = '/home/ahoffman/gyacomo/';
% resdir = 'CBC/96x96x16x3x2_Nexc_6';
% resdir = 'CBC/128x96x16x3x2';
% resdir = 'CBC/192x96x16x3x2';
-resdir = 'CBC/192x96x24x13x7';
+% resdir = 'CBC/192x96x24x13x7';
% resdir = 'CBC/kT_11_128x64x16x5x3';
% resdir = 'CBC/kT_9_256x128x16x3x2';
% resdir = 'CBC/kT_4.5_128x64x16x13x3';
@@ -61,6 +61,11 @@ resdir = 'CBC/192x96x24x13x7';
% 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';
% resdir = 'Miller/128x256x3x2_CBC_dt_5e-3';
+%% CBC Miller
+% resdir = 'GCM_CBC/daint/Miller_GX_comparison';
+%% RK3
+resdir = 'dbg/SSPx_RK3_test';
+resdir = 'dbg/SSPx_RK3_test/RK4';
resdir = ['results/',resdir];
JOBNUMMIN = 00; JOBNUMMAX = 10;
run analysis_gyacomo
diff --git a/wk/lin_ETPY.m b/wk/lin_ETPY.m
index 101f81c3062828803d7266099b64f051ce557c3a..9b6c383a95db4e8c877c6b1f79a8a3ceacc9e7c6 100644
--- a/wk/lin_ETPY.m
+++ b/wk/lin_ETPY.m
@@ -3,9 +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 = 'lin_ETPY'; % Name of the simulation
-% SIMID = 'dbg'; % Name of the simulation
-RUN = 0; % To run or just to load
+% SIMID = 'lin_ETPY'; % 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
PROGDIR = '/home/ahoffman/gyacomo/';
@@ -15,13 +15,14 @@ EXECNAME = 'gyacomo';
%% Set Up parameters
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.01; % Collision frequency
+NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 2.5; % ele Density '''
-K_Te = K_Ne/4; % ele Temperature '''
-K_Ni = K_Ne; % ion Density gradient drive
-K_Ti = K_Ni/4; % ion Temperature '''
+K_Ne = 2.0; % ele Density '''
+K_Te = 0.4; % ele Temperature '''
+K_Ni = 2.0; % ion Density gradient drive
+K_Ti = 0.4; % 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
@@ -33,24 +34,24 @@ JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
NX = 2; % real space x-gridpoints
-NY = 100; % '' y-gridpoints
+NY = 50; % '' y-gridpoints
LX = 2*pi/0.8; % Size of the squared frequency domain
-LY = 200;%2*pi/0.05; % 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
+GEOMETRY= 'Z-pinch';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
-KAPPA = 0.0; % elongation
+KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 500; % Maximal time unit
-DT = 5e-2; % Time step
+TMAX = 100; % 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
@@ -61,14 +62,14 @@ JOB2LOAD= -1;
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
+CO = 'DG';
+GKCO = 0; % 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
+INIT_OPT= 'phi'; % Start simulation with a noisy mom00/phi/allmom
%% OUTPUTS
W_DOUBLE = 1;
W_GAMMA = 1; W_HF = 1;
@@ -84,7 +85,7 @@ INIT_BLOB = 0; WIPE_TURB = 0; ACT_ON_MODES = 0;
MU_X = MU; %
MU_Y = MU; %
N_HD = 4;
-MU_Z = 2.0; %
+MU_Z = 0.0; %
MU_P = 0.0; %
MU_J = 0.0; %
LAMBDAD = 0.0;
diff --git a/wk/lin_ITG.m b/wk/lin_ITG.m
index 6f18bee5beed67944cf3a2eb4e0a483b1dda1a09..ebca810e6abee81ea9bba946a9a3cc6e9aa056b4 100644
--- a/wk/lin_ITG.m
+++ b/wk/lin_ITG.m
@@ -3,8 +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 = 'lin_ITG'; % Name of the simulation
-RUN = 0; % To run or just to load
+% SIMID = 'lin_ITG'; % 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
gyacomodir = '/home/ahoffman/gyacomo/';
@@ -18,9 +19,9 @@ CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%% PHYSICAL PARAMETERS
NU = 0.05; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-K_Ne = 2.22; % ele Density '''
+K_Ne = 0*2.22; % ele Density '''
K_Te = 6.96; % ele Temperature '''
-K_Ni = 2.22; % ion Density gradient drive
+K_Ni = 0*2.22; % ion Density gradient drive
K_Ti = 6.96; % ion 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)
@@ -45,13 +46,13 @@ SG = 0; % Staggered z grids option
GEOMETRY= 'miller';
Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear
-KAPPA = 0.0; % elongation
+KAPPA = 1.0; % elongation
DELTA = 0.0; % triangularity
ZETA = 0.0; % squareness
NEXC = 1; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 1; % Maximal time unit
+TMAX = 25; % Maximal time unit
DT = 3e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays