From 2a5f329a8148d9b9451dfd22378ddcfd7e0c3344 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Fri, 19 Feb 2021 11:20:36 +0100
Subject: [PATCH] scripts update
---
matlab/create_gif_5D.m | 25 ++++++-------
matlab/setup.m | 4 +++
matlab/write_fort90.m | 4 +++
wk/analysis_2D.m | 46 ++++++++++++++++++++----
wk/fig_post_processing.m | 78 ++++++++++++++++++++++++++++++++++++++++
wk/linear_study.m | 10 +++---
wk/local_run.m | 25 +++++++------
wk/marconi_run.m | 4 ++-
8 files changed, 160 insertions(+), 36 deletions(-)
create mode 100644 wk/fig_post_processing.m
diff --git a/matlab/create_gif_5D.m b/matlab/create_gif_5D.m
index c8646c7a..d77fd1f6 100644
--- a/matlab/create_gif_5D.m
+++ b/matlab/create_gif_5D.m
@@ -1,9 +1,3 @@
-
-GIFNAME = ['Nipj_kr',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
-plt = @(x) squeeze(max((abs(x)),[],4));
-FIELD = plt(Nipj); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
-FIELDNAME = 'N_i'; XNAME = '$k_{max}\rho_s$'; YNAME = '$P$';
-
title1 = GIFNAME;
FIGDIR = BASIC.RESDIR;
GIFNAME = [FIGDIR, GIFNAME,'.gif'];
@@ -14,34 +8,37 @@ sz = size(FIELD);
fig = figure('Color','white','Position', [100, 100, sz(2)*400, 400]);
for ij_ = 1:sz(2)
subplot(100+sz(2)*10+ij_)
- pclr = imagesc(X,Y,squeeze(FIELD(:,ij_,:,1)));
+% pclr = imagesc(X,Y,squeeze(FIELD(:,ij_,:,FRAMES(1))));
+ pclr = imagesc(X,Y,squeeze(log(FIELD(:,ij_,:,FRAMES(1)))));
xlabel('$k_r$');
if ij_ == 1
ylabel('$P$(max o. $k_z$)');
else
yticks([])
end
- LEGEND = ['$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$']; title(LEGEND);
+ LEGEND = ['ln$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$']; title(LEGEND);
end
- colormap gray
+% colormap gray
axis tight manual % this ensures that getframe() returns a consistent size
in = 1;
nbytes = fprintf(2,'frame %d/%d',in,numel(FIELD(1,1,1,:)));
for n = FRAMES % loop over selected frames
- scale = max(max(max(abs(FIELD(:,:,:,n)))));
+% scale = max(max(max(abs(FIELD(:,:,:,n)))));
for ij_ = 1:sz(2)
subplot(100+sz(2)*10+ij_)
- pclr = imagesc(X,Y,squeeze(FIELD(:,ij_,:,n))/scale);
+ scale = max(max(max(abs(FIELD(:,ij_,:,n)))));
+ pclr = imagesc(X,Y,squeeze(FIELD(:,ij_,:,n))/scale); caxis([0,1]);
+% pclr = imagesc(X,Y,squeeze(log(FIELD(:,ij_,:,n))));
xlabel(XNAME);
if ij_ == 1
ylabel(YNAME);
else
yticks([])
end
- LEGEND = ['$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$']; title(LEGEND);
+ LEGEND = ['ln$|',FIELDNAME,'^{p',num2str(ij_-1),'}|$'];
+ title([LEGEND,', amp = ',sprintf('%.1e',scale)]);
end
- suptitle(['$t \approx$', sprintf('%.3d',ceil(T(n)))...
- ,', scaling = ',sprintf('%.1e',scale)]);
+ suptitle(['$t \approx$', sprintf('%.3d',ceil(T(n)))]);
drawnow
% Capture the plot as an image
frame = getframe(fig);
diff --git a/matlab/setup.m b/matlab/setup.m
index f85f571d..ed77bc65 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -5,6 +5,8 @@ GRID.pmaxe = PMAXE; % Electron Hermite moments
GRID.jmaxe = JMAXE; % Electron Laguerre moments
GRID.pmaxi = PMAXI; % Ion Hermite moments
GRID.jmaxi = JMAXI; % Ion Laguerre moments
+GRID.p_damp = P_DAMP; % Ion Laguerre moments
+GRID.j_damp = J_DAMP; % Ion Laguerre moments
GRID.Nr = N; % r grid resolution
GRID.Lr = L; % r length
GRID.Nz = N * (1-KREQ0) + KREQ0; % z ''
@@ -16,6 +18,8 @@ MODEL.CO = CO; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -
MODEL.CLOS = CLOS;
if NON_LIN; MODEL.NON_LIN = '.true.'; else; MODEL.NON_LIN = '.false.';end;
MODEL.mu = MU;
+MODEL.mu_p = MU_P;
+MODEL.mu_j = MU_J;
MODEL.nu = NU; % hyper diffusive coefficient nu for HW
% temperature ratio T_a/T_e
MODEL.tau_e = TAU;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index d17d1ec0..6785e5ca 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -16,6 +16,8 @@ fprintf(fid,[' pmaxe =', num2str(GRID.pmaxe),'\n']);
fprintf(fid,[' jmaxe = ', num2str(GRID.jmaxe),'\n']);
fprintf(fid,[' pmaxi = ', num2str(GRID.pmaxi),'\n']);
fprintf(fid,[' jmaxi = ', num2str(GRID.jmaxi),'\n']);
+fprintf(fid,[' p_Damp = ', num2str(GRID.p_damp),'\n']);
+fprintf(fid,[' j_Damp = ', num2str(GRID.j_damp),'\n']);
fprintf(fid,[' Nr = ', num2str(GRID.Nr),'\n']);
fprintf(fid,[' Lr = ', num2str(GRID.Lr),'\n']);
fprintf(fid,[' Nz = ', num2str(GRID.Nz),'\n']);
@@ -41,6 +43,8 @@ fprintf(fid,[' CO = ', num2str(MODEL.CO),'\n']);
fprintf(fid,[' CLOS = ', num2str(MODEL.CLOS),'\n']);
fprintf(fid,[' NON_LIN = ', MODEL.NON_LIN,'\n']);
fprintf(fid,[' mu = ', num2str(MODEL.mu),'\n']);
+fprintf(fid,[' mu_p = ', num2str(MODEL.mu_p),'\n']);
+fprintf(fid,[' mu_j = ', num2str(MODEL.mu_j),'\n']);
fprintf(fid,[' nu = ', num2str(MODEL.nu),'\n']);
fprintf(fid,[' tau_e = ', num2str(MODEL.tau_e),'\n']);
fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']);
diff --git a/wk/analysis_2D.m b/wk/analysis_2D.m
index 825ade8c..65ed6764 100644
--- a/wk/analysis_2D.m
+++ b/wk/analysis_2D.m
@@ -8,7 +8,7 @@ if 0
end
%%
% JOBNUM = 0; load_results;
-% JOBNUM = 2; load_results;
+% JOBNUM = 1; load_results;
compile_results
load_params
@@ -272,7 +272,7 @@ end
if 1
%% Ion moments max mode vs pj
% tf = Ts2D(end-3);
-for tf = [0,1]
+for tf = []
[~,it2] = min(abs(Ts2D-tf)); [~,it5] = min(abs(Ts5D-tf));
% it2 = it2 + 1;
fig = figure; FIGNAME = ['kmaxp_Nipj_',sprintf('t=%.2f',Ts2D(it2))];set(gcf, 'Position', [100, 100, 700, 600])
@@ -298,10 +298,12 @@ end
%%
t0 = 0;
+[~, it02D] = min(abs(Ts2D-t0));
+[~, it05D] = min(abs(Ts5D-t0));
skip_ = 1;
-DELAY = 0.01*skip_;
-FRAMES_2D = floor(t0/(Ts2D(2)-Ts2D(1)))+1:skip_:numel(Ts2D);
-FRAMES_5D = floor(t0/(Ts5D(2)-Ts5D(1)))+1:skip_:numel(Ts5D);
+DELAY = 0.02*skip_;
+FRAMES_2D = it02D:skip_:numel(Ts2D);
+FRAMES_5D = it05D:skip_:numel(Ts5D);
%% GIFS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
%% Density ion
@@ -339,7 +341,7 @@ X = (r); T = Ts2D; YMIN = -1.1; YMAX = 1.1; XMIN = min(r); XMAX = max(r);
FIELDNAME = '$\phi(r=0)$'; XNAME = '$r/\rho_s$';
create_gif_1D
end
-if 1
+if 0
%% Density ion frequency
GIFNAME = ['Ni00',sprintf('_%.2d',JOBNUM)]; INTERP = 0; FRAMES = FRAMES_2D;
FIELD =ifftshift((abs(Ni00)),2); X = fftshift(KR,2); Y = fftshift(KZ,2); T = Ts2D;
@@ -362,4 +364,36 @@ FIELD =plt(Sipj(:,1,:,:,:)); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
FIELDNAME = '$N_i^{p0}$'; XNAME = '$k_{max}\rho_s$'; YNAME = '$P$';
create_gif_imagesc
end
+if 1
+%% maxkz, kr vs p, for all Nipj over time
+GIFNAME = ['Nipj_kr',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+plt = @(x) squeeze(max((abs(x)),[],4));
+FIELD = plt(Nipj); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
+FIELDNAME = 'N_i'; XNAME = '$k_r\rho_s$'; YNAME = '$P$, ${k_z}^{max}$';
+create_gif_5D
+end
+if 1
+%% maxkr, kz vs p, for all Nipj over time
+GIFNAME = ['Nipj_kz',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+plt = @(x) fftshift(squeeze(max((abs(x)),[],3)),3);
+FIELD = plt(Nipj); X = sort(kz'); Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
+FIELDNAME = 'N_i'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, ${k_r}^{max}$';
+create_gif_5D
+end
+if 0
+%% maxkz, kr vs p, for all Nepj over time
+GIFNAME = ['Nepj_kr',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+plt = @(x) squeeze(max((abs(x)),[],4));
+FIELD = plt(Nepj); X = kr'; Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
+FIELDNAME = 'N_e'; XNAME = '$k_r\rho_s$'; YNAME = '$P$, ${k_z}^{max}$';
+create_gif_5D
+end
+if 0
+%% maxkz, kz vs p, for all Nepj over time
+GIFNAME = ['Nepj_kz',sprintf('_%.2d',JOBNUM)]; INTERP = 0;
+plt = @(x) fftshift(squeeze(max((abs(x)),[],3)),3);
+FIELD = plt(Nepj); X = sort(kz'); Y = Pi'; T = Ts5D; FRAMES = FRAMES_5D;
+FIELDNAME = 'N_e'; XNAME = '$k_z\rho_s$'; YNAME = '$P$, ${k_r}^{max}$';
+create_gif_5D
+end
%%
\ No newline at end of file
diff --git a/wk/fig_post_processing.m b/wk/fig_post_processing.m
new file mode 100644
index 00000000..1479162b
--- /dev/null
+++ b/wk/fig_post_processing.m
@@ -0,0 +1,78 @@
+%% Load figure
+figpath = 'C:\Users\antoi\Desktop\gamma_eta_05_nu_1e-01_trunc';
+fig = openfig(figpath);
+
+%% Load data
+axObjs = fig.Children;
+dataObjs = findobj(fig,'-property','YData');
+Nlines = numel(dataObjs);
+
+%% Post processing
+sigma = zeros(Nlines,1);
+mu = sigma;
+tmin = 350;
+for i = 1:Nlines
+ x = dataObjs(i).XData;
+
+ [~, itmin] = min(abs(tmin-x));
+
+ y = dataObjs(i).YData;
+
+ sigma(i) = std(y(itmin:end));
+ mu(i) = mean(y(itmin:end));
+end
+mu = flip(mu')
+sigma = flip(sigma')
+
+%% Plot mean with error bar
+
+
+
+
+if 0
+%% Handwritten results for nu = 1.0, 150x75, L=100, DGGK
+Results_150x75.Gamma = [0.3794 0.3194 0.3226, 0.0098 0.0221 0.0321 0.0400 0.2897 0.2886 0.2569 0.0104 0.0086 0.0276 0.0320 0.1375 0.1633 0.0848];
+Results_150x75.error = [0.1909 0.1207 0.1336, 0.0028 0.0038 0.0058 0.0086 0.0832 0.0624 0.0557 0.0021 0.0023 0.0068 0.0088 0.0821 0.0278 0.0083];
+Results_150x75.P = [2, 4, 6, 2, 4, 6 8 2, 4, 6 2, 4, 6 8 2, 4, 10];
+Results_150x75.J = [1, 2, 3 1 2 3 4 1, 2, 3 1, 2, 3 4 1, 2 5];
+Results_150x75.etaB = [0.49, 0.49 0.49 0.59 0.59 0.59 0.59 0.50, 0.50 0.50 0.60, 0.60 0.60 0.60 0.51 0.51 0.51];
+Results_150x75.nu = [1.0, 1.0 1.0 1.0 1.0 1.0 1.0 0.5, 0.5 0.5 0.5, 0.5 0.5 0.5 0.1 0.1 0.1];
+Results_150x75.mrkr = [ '*', '*', '*', '*', '*', '*', '*', 'o', 'o', 'o', 'o', 'o', 'o', 'o' 's' 's' 's'];
+Results_150x75.iclr = [ 1, 2, 3, 1 2 3 4 1, 2, 3 1, 2, 3 4 1 2 5];
+
+% Ricci_Rogers.Gamma = [2 1e-1];
+% Ricci_Rogers.etaB = [0.5 1.0];
+Ricci_Rogers.Gamma = [10 1e-2];
+Ricci_Rogers.etaB = [0.5 1.25];
+
+if 1
+% Fig 3 of Ricci Rogers 2006
+SCALING = 2*sqrt(2);
+fig = figure;
+semilogy(Ricci_Rogers.etaB,Ricci_Rogers.Gamma,'--','color',[0,0,0]+0.6);
+hold on;
+plot(10,10,'color',line_colors(1,:)); plot(10,10,'color',line_colors(2,:));
+plot(10,10,'color',line_colors(3,:)); plot(10,10,'color',line_colors(4,:));
+plot(10,10,'color',line_colors(5,:));
+plot(10,10,'*k','MarkerSize',10, 'LineWidth',1.0);
+plot(10,10,'ok','MarkerSize',10, 'LineWidth',1.0);
+plot(10,10,'sk','MarkerSize',10, 'LineWidth',1.0);
+
+res = Results_150x75;
+for i = 1:numel(res.Gamma)
+ errorbar(res.etaB(i),res.Gamma(i)*SCALING,res.error(i)*SCALING,...
+ res.mrkr(i),'DisplayName','256x128', 'color', line_colors(res.iclr(i),:),...
+ 'MarkerSize',12, 'LineWidth',2.0);
+ hold on;
+end
+
+ xlabel('$L_n/L_B$'); ylabel('$2\sqrt(2)\Gamma^\infty_{part}$')
+end
+grid on; title('$L=100$, $150\times75$, $\nu_{hyp}=0.1$');
+xlim([0,1.75]); ylim([1e-6,100])
+legend('Mix. Length, Ricci 2006','$P=2$, $J=1$','$P=4$, $J=2$','$P=6$, $J=3$','$P=8$, $J=4$','$P=10$, $J=5$',...
+ '$\nu_{DGGK}=1.0$', '$\nu_{DGGK}=0.5$','$\nu_{DGGK}=0.1$');
+plot([0.3 0.3],[1e-6,1e2],'r')
+plot([1.6 1.6],[1e-6,1e2],'r')
+plot([0.5],[0.3965],'--','color',[0,0,0]+0.6)
+end
\ No newline at end of file
diff --git a/wk/linear_study.m b/wk/linear_study.m
index a6a13d6f..6643ab27 100644
--- a/wk/linear_study.m
+++ b/wk/linear_study.m
@@ -5,18 +5,20 @@ default_plots_options
%% Set Up parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.1; % Collision frequency
+NU = 1.0; % Collision frequency
TAU = 1.0; % e/i temperature ratio
-ETAB = 0.5;
+ETAB = 0.6;
ETAN = 1.0; % Density gradient
ETAT = 0.0; % Temperature gradient
-NU_HYP = 0.1; % Hyperdiffusivity coefficient
+NU_HYP = 1.0; % Hyperdiffusivity coefficient
LAMBDAD = 0.0;
NOISE0 = 1.0e-5;
%% GRID PARAMETERS
N = 150; % Frequency gridpoints (Nkr = N/2)
L = 70; % Size of the squared frequency domain
KREQ0 = 1; % put kr = 0
+P_DAMP = 0; % Hermite damping term -(j/Jmax)^{2*r}Napj (0: no damping, 1: r=1, 2: r=2 etc...)
+J_DAMP = 0; % Laguerre damping term -(j/Jmax)^{2*r}Napj (0: no damping, 1: r=1, 2: r=2 etc...)
%% TIME PARMETERS
TMAX = 100; % Maximal time unit
DT = 1e-2; % Time step
@@ -30,7 +32,7 @@ JOB2LOAD= 00;
SIMID = 'linear_study'; % Name of the simulation
NON_LIN = 0 *(1-KREQ0); % activate non-linearity (is cancelled if KREQ0 = 1)
CO = -3; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
-CLOS = 0; % Closure model (0 : =0 truncation, 1 : n+j = min(nmax,n+j), 2: odd/even adapted)
+CLOS = 2; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
KERN = 0; % Kernel model (0 : GK)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/wk/local_run.m b/wk/local_run.m
index 377fd3ec..e64ad61f 100644
--- a/wk/local_run.m
+++ b/wk/local_run.m
@@ -12,28 +12,31 @@ ETAT = 0.0; % Temperature gradient
HD_CO = 0.5; % Hyper diffusivity cutoff ratio
NU_HYP = 0.1;
%% GRID PARAMETERS
-N = 60; % Frequency gridpoints (Nkr = N/2)
+N = 150; % Frequency gridpoints (Nkr = N/2)
L = 70; % Size of the squared frequency domain
-PMAXE = 2; % Highest electron Hermite polynomial degree
-JMAXE = 1; % Highest '' Laguerre ''
-PMAXI = 2; % Highest ion Hermite polynomial degree
-JMAXI = 1; % Highest '' Laguerre ''
+PMAXE = 4; % Highest electron Hermite polynomial degree
+JMAXE = 4; % Highest '' Laguerre ''
+PMAXI = 4; % Highest ion Hermite polynomial degree
+JMAXI = 4; % Highest '' Laguerre ''
+MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
-TMAX = 50; % Maximal time unit
-DT = 1e-2; % Time step
-SPS0D = 1; % Sampling per time unit for profiler
+TMAX = 600; % Maximal time unit
+DT = 2e-2; % Time step
+SPS0D = 1/DT; % Sampling per time unit for profiler
SPS2D = 1; % Sampling per time unit for 2D arrays
SPS5D = 1; % Sampling per time unit for 5D arrays
-SPSCP = 1/10; % Sampling per time unit for checkpoints/10
+SPSCP = 0; % Sampling per time unit for checkpoints/10
RESTART = 0; % To restart from last checkpoint
JOB2LOAD= 0;
%% OPTIONS
% SIMID = 'local_nu_%0.0e'; % Name of the simulation
% SIMID = sprintf(SIMID,NU);
% SIMID = 'Marconi_DGGK_nu_1e+00'; % Name of the simulation
-SIMID = 'debug_cp_load'; % Name of the simulation
+% SIMID = 'test_sapj_cut'; % Name of the simulation
+SIMID = 'test_moments_damping'; % Name of the simulation
CO = -3; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
-CLOS = 0; % Closure model (0 : =0 truncation, 1 : n+j = min(nmax,n+j), 2: odd/even adapted)
+CLOS = 2; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
KERN = 0; % Kernel model (0 : GK)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/wk/marconi_run.m b/wk/marconi_run.m
index 394387bc..0b59434a 100644
--- a/wk/marconi_run.m
+++ b/wk/marconi_run.m
@@ -21,6 +21,8 @@ N = 150; % Frequency gridpoints (Nkr = N/2)
L = 70; % Size of the squared frequency domain
P = 2; % Electron and Ion highest Hermite polynomial degree
J = 1; % Electron and Ion highest Laguerre polynomial degree
+MU_P = 0; % Hermite hyperdiffusivity -mu_p*(d/dvpar)^4 f
+MU_J = 0; % Laguerre hyperdiffusivity -mu_j*(d/dvperp)^4 f
%% TIME PARAMETERS
TMAX = 500; % Maximal time unit
DT = 2e-2; % Time step
@@ -35,7 +37,7 @@ JOB2LOAD= 0;
SIMID = 'Marconi_restart'; % Name of the simulation
SIMID = sprintf(SIMID,NU);
CO = -3; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty, -3 : GK Dougherty)
-CLOS = 0; % Closure model (0 : =0 truncation, 1 : n+j = min(nmax,n+j), 2: odd/even adapted)
+CLOS = 0; % Closure model (0: =0 truncation, 1: semi coll, 2: Copy closure J+1 = J, P+2 = P)
KERN = 0; % Kernel model (0 : GK)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--
GitLab