Skip to content
Snippets Groups Projects
Commit a4cfb827 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann
Browse files

refresh

parent 3cf05fae
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
wk/benchmark_kperp_scan.m
+ 7
− 7
View file @ a4cfb827
...@@ -12,18 +12,18 @@ OUTPUTS.write_phi = '.true.'; ...@@ -12,18 +12,18 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.'; OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results'''; OUTPUTS.resfile0 = '''results''';
%% Grid parameters %% Grid parameters
GRID.pmaxe = 12; GRID.pmaxe = 15;
GRID.jmaxe = 7; GRID.jmaxe = 6;
GRID.pmaxi = 12; GRID.pmaxi = 15;
GRID.jmaxi = 7; GRID.jmaxi = 6;
GRID.nkr = 1; GRID.nkr = 1;
GRID.krmin = 0.; GRID.krmin = 0.;
GRID.krmax = 0.; GRID.krmax = 0.;
GRID.nkz = 20; GRID.nkz = 20;
GRID.kzmin = 0.1; GRID.kzmin = 0.1;
GRID.kzmax = 2.0; GRID.kzmax = 1.5;
%% Model parameters %% Model parameters
MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0 MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0
% temperature ratio T_a/T_e % temperature ratio T_a/T_e
MODEL.tau_e = 1.0; MODEL.tau_e = 1.0;
...@@ -44,7 +44,7 @@ MODEL.lambdaD = 0.0; ...@@ -44,7 +44,7 @@ MODEL.lambdaD = 0.0;
TIME_INTEGRATION.numerical_scheme = '''RK4'''; TIME_INTEGRATION.numerical_scheme = '''RK4''';
BASIC.nrun = 100000; BASIC.nrun = 100000;
BASIC.dt = 0.05; BASIC.dt = 0.05;
BASIC.tmax = 100.0; BASIC.tmax = 50.0;
INITIAL.initback_moments = 0.01; INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.; INITIAL.initnoise_moments = 0.;
INITIAL.iseed = 42; INITIAL.iseed = 42;
......
wk/benchmark_time_solver.m
+ 48
− 57
View file @ a4cfb827
clear all; close all; clear all; close all;
SIMID = 'benchmark_time_solver'; % Name of the simulations SIMID = 'test_full_coulomb'; % Name of the simulations
addpath(genpath('../matlab')) % ... add addpath(genpath('../matlab')) % ... add
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% outputs options %% outputs options
OUTPUTS.nsave_0d = 0; OUTPUTS.nsave_0d = 0;
...@@ -13,19 +13,19 @@ OUTPUTS.write_phi = '.true.'; ...@@ -13,19 +13,19 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.'; OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results'''; OUTPUTS.resfile0 = '''results''';
%% Grid parameters %% Grid parameters
GRID.pmaxe = 6; GRID.pmaxe = 15;
GRID.jmaxe = 6; GRID.jmaxe = 6;
GRID.pmaxi = 6; GRID.pmaxi = 15;
GRID.jmaxi = 6; GRID.jmaxi = 6;
GRID.nkr = 1; GRID.nkr = 1;
GRID.krmin = 0.; GRID.krmin = 0.;
GRID.krmax = 0.; GRID.krmax = 0.;
GRID.nkz = 1; GRID.nkz = 1;
GRID.kzmin = 0.67; GRID.kzmin = 1.0;
GRID.kzmax = 0.67; GRID.kzmax = 1.0;
%% Model parameters %% Model parameters
MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty) MODEL.CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0 MODEL.nu = 1.0; % collisionality nu*L_perp/Cs0
% temperature ratio T_a/T_e % temperature ratio T_a/T_e
MODEL.tau_e = 1.0; MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0; MODEL.tau_i = 1.0;
...@@ -36,16 +36,16 @@ MODEL.sigma_i = 1.0; ...@@ -36,16 +36,16 @@ MODEL.sigma_i = 1.0;
MODEL.q_e =-1.0; MODEL.q_e =-1.0;
MODEL.q_i = 1.0; MODEL.q_i = 1.0;
% gradients L_perp/L_x % gradients L_perp/L_x
MODEL.eta_n = 1.0; % Density MODEL.eta_n = 0.0; % Density
MODEL.eta_T = 0.0; % Temperature MODEL.eta_T = 0.0; % Temperature
MODEL.eta_B = 0.5; % Magnetic MODEL.eta_B = 0.0; % Magnetic
% Debye length % Debye length
MODEL.lambdaD = 0.0; MODEL.lambdaD = 0.0;
%% Time integration and intialization parameters %% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4'''; TIME_INTEGRATION.numerical_scheme = '''RK4''';
BASIC.nrun = 100000; BASIC.nrun = 100000;
BASIC.dt = 0.05; BASIC.dt = 0.05;
BASIC.tmax = 100.0; BASIC.tmax = 10.0;
INITIAL.initback_moments = 0.01; INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.; INITIAL.initnoise_moments = 0.;
INITIAL.iseed = 42; INITIAL.iseed = 42;
...@@ -86,63 +86,54 @@ default_plots_options ...@@ -86,63 +86,54 @@ default_plots_options
SAVEFIG = 1; SAVEFIG = 1;
filename = 'results_00.h5'; filename = 'results_00.h5';
default_plots_options default_plots_options
TITLE = [', $k_z=',num2str(GRID.kzmin)]; TITLE = [TITLE,', $k_z=',num2str(GRID.kzmin),'$'];
%% Nipj bare = @(p_,j_) (GRID.jmaxe+1)*p_ + j_ + 1;
bari = @(p_,j_) bare(GRID.pmaxe, GRID.jmaxe) + (GRID.jmaxi+1)*p_ + j_ + 1;
%% Load moments
% load HeLaZ data moment = 'moments_i';
moment = 'Ni00';
kr = h5read(filename,['/data/var2d/' moment '/coordkr']);
kz = h5read(filename,['/data/var2d/' moment '/coordkz']);
timeNi = squeeze(h5read(filename,'/data/var2d/time'));
Ni00 = zeros(numel(kr),numel(kz),numel(timeNi));
for it = 1:numel(timeNi) kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
tmp = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]); kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
Ni00(:,:,it) = tmp.real + 1i * tmp.imaginary; time = h5read(filename,'/data/var5d/time');
Nipj = zeros(GRID.pmaxi+1, GRID.jmaxi+1,numel(kr),numel(kz),numel(time));
for it = 1:numel(time)
tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
Nipj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
end end
moment = 'moments_i'; moment = 'moments_e';
kr = h5read(filename,['/data/var5d/' moment '/coordkr']); kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
kz = h5read(filename,['/data/var5d/' moment '/coordkz']); kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
time = h5read(filename,'/data/var5d/time'); time = h5read(filename,'/data/var5d/time');
Nipj = zeros(GRID.pmaxi+1, GRID.jmaxi+1,numel(kr),numel(kz),numel(time)); Nepj = zeros(GRID.pmaxe+1, GRID.jmaxe+1,numel(kr),numel(kz),numel(time));
for it = 1:numel(time) for it = 1:numel(time)
tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]); tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
Nipj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary; Nepj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
end end
err_ = mean(abs(squeeze(Nipj(1,1,1,1,:))-squeeze((Ni00(1,1,:)))));
disp(['2D 5D error :',num2str(err_)])
% growth rate :
gammaHeLaZ = LinearFit_s(squeeze(time),...
squeeze(abs(Nipj(1,1,1,1,:)))) ;
gammaMOLI = LinearFit(results.Napj,results.time,params,options);
disp(['Growth rate HeLaZ : ',num2str(gammaHeLaZ)])
disp(['Growth rate MOLI : ',num2str( gammaMOLI)])
disp(['Error : ',num2str( abs(gammaMOLI-gammaHeLaZ))])
% Plot moments %% Plot moments
fig = figure; fig = figure;
x1 = time; x1 = time;
x2 = results.time; x2 = results.time;
ic = 1; ic = 1;
% Real part % Electrons
subplot(321) subplot(321)
for ip = 0:1 for ip = 0:1
for ij = 0:1 for ij = 0:1
y1 = squeeze(real(Nipj(ip+1,ij+1,1,1,:))); y1 = squeeze(real(Nepj(ip+1,ij+1,1,1,:)));
plot(x1,y1,'-','DisplayName',... plot(x1,y1,'-','DisplayName',...
['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],... ['HeLaZ $N_e^{',num2str(ip),num2str(ij),'}$'],...
'color', line_colors(ic,:)) 'color', line_colors(ic,:))
hold on hold on
y2 = squeeze(real(results.Napj(:,bari(ip,ij)))); y2 = squeeze(real(results.Napj(:,bare(ip,ij))));
plot(x2,y2,'--','DisplayName',... plot(x2,y2,'--','DisplayName',...
['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],... ['MOLI $N_e^{',num2str(ip),num2str(ij),'}$'],...
'color', line_colors(ic,:)) 'color', line_colors(ic,:))
hold on hold on
ic = ic + 1; ic = ic + 1;
...@@ -150,19 +141,19 @@ for ip = 0:1 ...@@ -150,19 +141,19 @@ for ip = 0:1
end end
grid on grid on
xlabel('$t$') xlabel('$t$')
ylabel(['Re$(N_i^{pj})$']) ylabel(['Re$(N_e^{pj})$'])
% Im part % Ions
ic = 1; ic = 1;
subplot(322) subplot(322)
for ip = 0:1 for ip = 0:1
for ij = 0:1 for ij = 0:1
y1 = squeeze(imag(Nipj(ip+1,ij+1,1,1,:))); y1 = squeeze(real(Nipj(ip+1,ij+1,1,1,:)));
plot(x1,y1,'-','DisplayName',... plot(x1,y1,'-','DisplayName',...
['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],... ['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],...
'color', line_colors(ic,:)) 'color', line_colors(ic,:))
hold on hold on
y2 = squeeze(imag(results.Napj(:,bari(ip,ij)))); y2 = squeeze(real(results.Napj(:,bari(ip,ij))));
plot(x2,y2,'--','DisplayName',... plot(x2,y2,'--','DisplayName',...
['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],... ['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],...
'color', line_colors(ic,:)) 'color', line_colors(ic,:))
...@@ -172,11 +163,11 @@ for ip = 0:1 ...@@ -172,11 +163,11 @@ for ip = 0:1
end end
grid on grid on
xlabel('$t$') xlabel('$t$')
ylabel(['Im$(N_i^{pj})$']) ylabel(['Re$(N_i^{pj})$'])
%suptitle(TITLE); %suptitle(TITLE);
%if SAVEFIG; FIGNAME = ['Nipj_kz_',num2str(GRID.kzmin)]; save_figure; end; %if SAVEFIG; FIGNAME = ['Nipj_kz_',num2str(GRID.kzmin)]; save_figure; end;
%% phi % phi
timephi = h5read(filename,'/data/var2d/time'); timephi = h5read(filename,'/data/var2d/time');
kr = h5read(filename,'/data/var2d/phi/coordkr'); kr = h5read(filename,'/data/var2d/phi/coordkr');
kz = h5read(filename,'/data/var2d/phi/coordkz'); kz = h5read(filename,'/data/var2d/phi/coordkz');
...@@ -227,14 +218,14 @@ for it = 1:numel(timephiMOLI) ...@@ -227,14 +218,14 @@ for it = 1:numel(timephiMOLI)
phiMOLI(it) = get_phi(results.Napj(it,:),params,options); phiMOLI(it) = get_phi(results.Napj(it,:),params,options);
end end
ERR1 = abs(real(phiMOLI - phiHeLaZ)); ERR1 = abs(real(phiMOLI(1:numel(timephi)) - phiHeLaZ));
ERR2 = abs(imag(phiMOLI - phiHeLaZ)); ERR2 = abs(imag(phiMOLI(1:numel(timephi)) - phiHeLaZ));
%fig = figure; %fig = figure;
subplot(325); subplot(325);
plot(timephiMOLI,ERR1,'-','DisplayName','Real') plot(timephi,ERR1,'-','DisplayName','Real')
hold on hold on
plot(timephiMOLI,ERR2,'--','DisplayName','Imag') plot(timephi,ERR2,'--','DisplayName','Imag')
%title(TITLE); %title(TITLE);
grid on grid on
xlabel('$t$') xlabel('$t$')
...@@ -242,12 +233,12 @@ ylabel('$e(\phi)$') ...@@ -242,12 +233,12 @@ ylabel('$e(\phi)$')
legend('show') legend('show')
%if SAVEFIG; FIGNAME = ['err_phi_kz_',num2str(GRID.kzmin)]; save_figure; end; %if SAVEFIG; FIGNAME = ['err_phi_kz_',num2str(GRID.kzmin)]; save_figure; end;
%% Growth rate fit quantity (Ni00) % Growth rate fit quantity (Ni00)
subplot(326); subplot(326);
y1 = squeeze(abs(Nipj(1,1,1,1,:))); % HeLaZ y1 = squeeze(abs(Nipj(1,1,1,1,:))); % HeLaZ
y2 = squeeze(abs(results.Napj(:,bari(0,0)))); y2 = squeeze(abs(results.Napj(:,bari(0,0))));
ERR3 = abs(y2-y1); ERR3 = abs(y2(1:numel(timephi))-y1);
yyaxis left yyaxis left
semilogy(x1,y1,'-','DisplayName','HeLaZ',... semilogy(x1,y1,'-','DisplayName','HeLaZ',...
...@@ -265,7 +256,7 @@ grid on ...@@ -265,7 +256,7 @@ grid on
xlabel('$t$') xlabel('$t$')
ylabel('$e(N_i^{00})$') ylabel('$e(N_i^{00})$')
%% Finish and save % Finish and save
suptitle(TITLE); suptitle(TITLE);
if SAVEFIG; FIGNAME = ['kz_',num2str(GRID.kzmin)]; save_figure; end; if SAVEFIG; FIGNAME = ['kz_',num2str(GRID.kzmin)]; save_figure; end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
wk/fort.90
+ 12
− 12
View file @ a4cfb827
&BASIC &BASIC
nrun=100000 nrun=100000
dt=0.05 dt=0.01
tmax=100 ! time normalized to 1/omega_pe tmax=10 ! time normalized to 1/omega_pe
/ /
&GRID &GRID
pmaxe =12 ! Electron Hermite moments pmaxe =15 ! Electron Hermite moments
jmaxe = 7 ! Electron Laguerre moments jmaxe = 6 ! Electron Laguerre moments
pmaxi = 12 ! Ion Hermite moments pmaxi = 15 ! Ion Hermite moments
jmaxi = 7 ! Ion Laguerre moments jmaxi = 6 ! Ion Laguerre moments
nkr = 1 nkr = 1
krmin = 0 krmin = 0
krmax = 0 krmax = 0
nkz = 20 nkz = 20
kzmin = 0.1 kzmin = 1.5
kzmax = 2 kzmax = 2.5
/ /
&OUTPUT_PAR &OUTPUT_PAR
nsave_0d = 0 nsave_0d = 0
...@@ -28,7 +28,7 @@ ...@@ -28,7 +28,7 @@
/ /
&MODEL_PAR &MODEL_PAR
! Collisionality ! Collisionality
CO = -2 ! Collision operator (-1:Full Coulomb, 0: Dougherty) CO = -1 ! Collision operator (-1:Full Coulomb, 0: Dougherty)
nu = 0.01 ! Normalized collision frequency normalized to plasma frequency nu = 0.01 ! Normalized collision frequency normalized to plasma frequency
tau_e = 1 ! T_e/T_e tau_e = 1 ! T_e/T_e
tau_i = 1 ! T_i/T_e temperature ratio tau_i = 1 ! T_i/T_e temperature ratio
...@@ -47,9 +47,9 @@ ...@@ -47,9 +47,9 @@
! Noise amplitude ! Noise amplitude
initnoise_moments=0 initnoise_moments=0
iseed=42 iseed=42
selfmat_file='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_12_Jmaxe_7_Pmaxi_12_Jmaxi_7_pamaxx_10.h5' selfmat_file='../iCa/self_Coll_GKE_0_GKI_0_ESELF_1_ISELF_1_Pmaxe_15_Jmaxe_6_Pmaxi_15_Jmaxi_6_pamaxx_10.h5'
eimat_file='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_12_Jmaxe_7_Pmaxi_12_Jmaxi_7_pamaxx_10_tau_1.0000_mu_0.0233.h5' eimat_file='../iCa/ei_Coll_GKE_0_GKI_0_ETEST_1_EBACK_1_Pmaxe_15_Jmaxe_6_Pmaxi_15_Jmaxi_6_pamaxx_10_tau_1.0000_mu_0.0233.h5'
iemat_file='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_12_Jmaxe_7_Pmaxi_12_Jmaxi_7_pamaxx_10_tau_1.0000_mu_0.0233.h5' iemat_file='../iCa/ie_Coll_GKE_0_GKI_0_ITEST_1_IBACK_1_Pmaxe_15_Jmaxe_6_Pmaxi_15_Jmaxi_6_pamaxx_10_tau_1.0000_mu_0.0233.h5'
/ /
&TIME_INTEGRATION_PAR &TIME_INTEGRATION_PAR
numerical_scheme='RK4' numerical_scheme='RK4'
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment