From a4cfb827a85290fc3fead0e718ea986750f49e51 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Wed, 8 Jul 2020 10:07:11 +0200
Subject: [PATCH] refresh
---
wk/benchmark_kperp_scan.m | 14 ++---
wk/benchmark_time_solver.m | 105 +++++++++++++++++--------------------
wk/fort.90 | 24 ++++-----
3 files changed, 67 insertions(+), 76 deletions(-)
diff --git a/wk/benchmark_kperp_scan.m b/wk/benchmark_kperp_scan.m
index e14e54eb..b5432bf4 100644
--- a/wk/benchmark_kperp_scan.m
+++ b/wk/benchmark_kperp_scan.m
@@ -12,18 +12,18 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
%% Grid parameters
-GRID.pmaxe = 12;
-GRID.jmaxe = 7;
-GRID.pmaxi = 12;
-GRID.jmaxi = 7;
+GRID.pmaxe = 15;
+GRID.jmaxe = 6;
+GRID.pmaxi = 15;
+GRID.jmaxi = 6;
GRID.nkr = 1;
GRID.krmin = 0.;
GRID.krmax = 0.;
GRID.nkz = 20;
GRID.kzmin = 0.1;
-GRID.kzmax = 2.0;
+GRID.kzmax = 1.5;
%% 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
% temperature ratio T_a/T_e
MODEL.tau_e = 1.0;
@@ -44,7 +44,7 @@ MODEL.lambdaD = 0.0;
TIME_INTEGRATION.numerical_scheme = '''RK4''';
BASIC.nrun = 100000;
BASIC.dt = 0.05;
-BASIC.tmax = 100.0;
+BASIC.tmax = 50.0;
INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.;
INITIAL.iseed = 42;
diff --git a/wk/benchmark_time_solver.m b/wk/benchmark_time_solver.m
index 2eaaaa26..641f0f0c 100644
--- a/wk/benchmark_time_solver.m
+++ b/wk/benchmark_time_solver.m
@@ -1,6 +1,6 @@
clear all; close all;
-SIMID = 'benchmark_time_solver'; % Name of the simulations
-addpath(genpath('../matlab')) % ... add
+SIMID = 'test_full_coulomb'; % Name of the simulations
+addpath(genpath('../matlab')) % ... add
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% outputs options
OUTPUTS.nsave_0d = 0;
@@ -13,19 +13,19 @@ OUTPUTS.write_phi = '.true.';
OUTPUTS.write_doubleprecision = '.true.';
OUTPUTS.resfile0 = '''results''';
%% Grid parameters
-GRID.pmaxe = 6;
+GRID.pmaxe = 15;
GRID.jmaxe = 6;
-GRID.pmaxi = 6;
+GRID.pmaxi = 15;
GRID.jmaxi = 6;
GRID.nkr = 1;
GRID.krmin = 0.;
GRID.krmax = 0.;
GRID.nkz = 1;
-GRID.kzmin = 0.67;
-GRID.kzmax = 0.67;
+GRID.kzmin = 1.0;
+GRID.kzmax = 1.0;
%% Model parameters
-MODEL.CO = -2; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
-MODEL.nu = 0.01; % collisionality nu*L_perp/Cs0
+MODEL.CO = -1; % Collision operator (0 : L.Bernstein, -1 : Full Coulomb, -2 : Dougherty)
+MODEL.nu = 1.0; % collisionality nu*L_perp/Cs0
% temperature ratio T_a/T_e
MODEL.tau_e = 1.0;
MODEL.tau_i = 1.0;
@@ -36,16 +36,16 @@ MODEL.sigma_i = 1.0;
MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
% 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_B = 0.5; % Magnetic
+MODEL.eta_B = 0.0; % Magnetic
% Debye length
MODEL.lambdaD = 0.0;
%% Time integration and intialization parameters
TIME_INTEGRATION.numerical_scheme = '''RK4''';
BASIC.nrun = 100000;
BASIC.dt = 0.05;
-BASIC.tmax = 100.0;
+BASIC.tmax = 10.0;
INITIAL.initback_moments = 0.01;
INITIAL.initnoise_moments = 0.;
INITIAL.iseed = 42;
@@ -86,63 +86,54 @@ default_plots_options
SAVEFIG = 1;
filename = 'results_00.h5';
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 = '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));
+moment = 'moments_i';
-for it = 1:numel(timeNi)
- tmp = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]);
- Ni00(:,:,it) = tmp.real + 1i * tmp.imaginary;
+kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
+kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
+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
-moment = 'moments_i';
+moment = 'moments_e';
-kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
-kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
+kr = h5read(filename,['/data/var5d/' moment '/coordkr']);
+kz = h5read(filename,['/data/var5d/' moment '/coordkz']);
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)
tmp = h5read(filename,['/data/var5d/', moment,'/', num2str(it,'%06d')]);
- Nipj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
+ Nepj(:,:,:,:,it) = tmp.real + 1i * tmp.imaginary;
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;
x1 = time;
x2 = results.time;
ic = 1;
-% Real part
+% Electrons
subplot(321)
for ip = 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',...
- ['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ ['HeLaZ $N_e^{',num2str(ip),num2str(ij),'}$'],...
'color', line_colors(ic,:))
hold on
- y2 = squeeze(real(results.Napj(:,bari(ip,ij))));
+ y2 = squeeze(real(results.Napj(:,bare(ip,ij))));
plot(x2,y2,'--','DisplayName',...
- ['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],...
+ ['MOLI $N_e^{',num2str(ip),num2str(ij),'}$'],...
'color', line_colors(ic,:))
hold on
ic = ic + 1;
@@ -150,19 +141,19 @@ for ip = 0:1
end
grid on
xlabel('$t$')
-ylabel(['Re$(N_i^{pj})$'])
+ylabel(['Re$(N_e^{pj})$'])
-% Im part
+% Ions
ic = 1;
subplot(322)
for ip = 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',...
['HeLaZ $N_i^{',num2str(ip),num2str(ij),'}$'],...
'color', line_colors(ic,:))
hold on
- y2 = squeeze(imag(results.Napj(:,bari(ip,ij))));
+ y2 = squeeze(real(results.Napj(:,bari(ip,ij))));
plot(x2,y2,'--','DisplayName',...
['MOLI $N_i^{',num2str(ip),num2str(ij),'}$'],...
'color', line_colors(ic,:))
@@ -172,11 +163,11 @@ for ip = 0:1
end
grid on
xlabel('$t$')
-ylabel(['Im$(N_i^{pj})$'])
+ylabel(['Re$(N_i^{pj})$'])
%suptitle(TITLE);
%if SAVEFIG; FIGNAME = ['Nipj_kz_',num2str(GRID.kzmin)]; save_figure; end;
-%% phi
+% phi
timephi = h5read(filename,'/data/var2d/time');
kr = h5read(filename,'/data/var2d/phi/coordkr');
kz = h5read(filename,'/data/var2d/phi/coordkz');
@@ -227,14 +218,14 @@ for it = 1:numel(timephiMOLI)
phiMOLI(it) = get_phi(results.Napj(it,:),params,options);
end
-ERR1 = abs(real(phiMOLI - phiHeLaZ));
-ERR2 = abs(imag(phiMOLI - phiHeLaZ));
+ERR1 = abs(real(phiMOLI(1:numel(timephi)) - phiHeLaZ));
+ERR2 = abs(imag(phiMOLI(1:numel(timephi)) - phiHeLaZ));
%fig = figure;
subplot(325);
-plot(timephiMOLI,ERR1,'-','DisplayName','Real')
+plot(timephi,ERR1,'-','DisplayName','Real')
hold on
-plot(timephiMOLI,ERR2,'--','DisplayName','Imag')
+plot(timephi,ERR2,'--','DisplayName','Imag')
%title(TITLE);
grid on
xlabel('$t$')
@@ -242,12 +233,12 @@ ylabel('$e(\phi)$')
legend('show')
%if SAVEFIG; FIGNAME = ['err_phi_kz_',num2str(GRID.kzmin)]; save_figure; end;
-%% Growth rate fit quantity (Ni00)
+% Growth rate fit quantity (Ni00)
subplot(326);
y1 = squeeze(abs(Nipj(1,1,1,1,:))); % HeLaZ
y2 = squeeze(abs(results.Napj(:,bari(0,0))));
-ERR3 = abs(y2-y1);
+ERR3 = abs(y2(1:numel(timephi))-y1);
yyaxis left
semilogy(x1,y1,'-','DisplayName','HeLaZ',...
@@ -265,7 +256,7 @@ grid on
xlabel('$t$')
ylabel('$e(N_i^{00})$')
-%% Finish and save
+% Finish and save
suptitle(TITLE);
if SAVEFIG; FIGNAME = ['kz_',num2str(GRID.kzmin)]; save_figure; end;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/wk/fort.90 b/wk/fort.90
index 2017fc90..449f4371 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,19 +1,19 @@
&BASIC
nrun=100000
- dt=0.05
- tmax=100 ! time normalized to 1/omega_pe
+ dt=0.01
+ tmax=10 ! time normalized to 1/omega_pe
/
&GRID
- pmaxe =12 ! Electron Hermite moments
- jmaxe = 7 ! Electron Laguerre moments
- pmaxi = 12 ! Ion Hermite moments
- jmaxi = 7 ! Ion Laguerre moments
+ pmaxe =15 ! Electron Hermite moments
+ jmaxe = 6 ! Electron Laguerre moments
+ pmaxi = 15 ! Ion Hermite moments
+ jmaxi = 6 ! Ion Laguerre moments
nkr = 1
krmin = 0
krmax = 0
nkz = 20
- kzmin = 0.1
- kzmax = 2
+ kzmin = 1.5
+ kzmax = 2.5
/
&OUTPUT_PAR
nsave_0d = 0
@@ -28,7 +28,7 @@
/
&MODEL_PAR
! 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
tau_e = 1 ! T_e/T_e
tau_i = 1 ! T_i/T_e temperature ratio
@@ -47,9 +47,9 @@
! Noise amplitude
initnoise_moments=0
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'
- 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'
- 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'
+ 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_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_15_Jmaxe_6_Pmaxi_15_Jmaxi_6_pamaxx_10_tau_1.0000_mu_0.0233.h5'
/
&TIME_INTEGRATION_PAR
numerical_scheme='RK4'
--
GitLab