From 60c4b5a580da9697033b7b03df51733141a69fe0 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Wed, 24 Jun 2020 17:25:40 +0200
Subject: [PATCH] slight fixes

---
 .gitignore    |   1 +
 wk/fort.90    |  14 +++----
 wk/launcher.m | 103 +++++++++++++++++++++++++++++++++++++++++++-------
 3 files changed, 98 insertions(+), 20 deletions(-)

diff --git a/.gitignore b/.gitignore
index 500b87ee..c6daf83b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,3 +28,4 @@ obj/
 bin/
 .gitignore
 wk/fort.90
+.directory
diff --git a/wk/fort.90 b/wk/fort.90
index 73a62560..84051bf6 100644
--- a/wk/fort.90
+++ b/wk/fort.90
@@ -1,18 +1,18 @@
 &BASIC
   nrun=100000
   dt=0.01
-  tmax=200    ! time normalized to 1/omega_pe
+  tmax=50    ! time normalized to 1/omega_pe
 /
 &GRID
-  pmaxe =15    ! number of Hermite moments 
-  jmaxe = 4     ! number of Hermite moments 
-  pmaxi = 15    ! number of Hermite moments 
-  jmaxi = 4     ! number of Hermite moments
+  pmaxe =81    ! number of Hermite moments 
+  jmaxe = 20     ! number of Hermite moments 
+  pmaxi = 81    ! number of Hermite moments 
+  jmaxi = 20     ! number of Hermite moments
   nkr   = 1
   krmin = 0
   krmax = 0     ! Normalized to cs0/Omega_i
-  nkz   = 9
-  kzmin = 0.1
+  nkz   = 1
+  kzmin = 1
   kzmax = 1     ! Normalized to cs0/Omega_i
 /
 &OUTPUT_PAR
diff --git a/wk/launcher.m b/wk/launcher.m
index fd1c83de..9f031c53 100644
--- a/wk/launcher.m
+++ b/wk/launcher.m
@@ -1,4 +1,4 @@
-SIMID = 'Linear_Zpinch'; % Name of the simulations
+SIMID = 'MOLI_Comparison'; % Name of the simulations
 addpath(genpath('../matlab')) % ... add 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -13,18 +13,18 @@ OUTPUTS.write_phi     = '.true.';
 OUTPUTS.write_doubleprecision = '.true.';
 OUTPUTS.resfile0      = '''results''';
 %% Grid parameters
-GRID.pmaxe = 15;
-GRID.jmaxe = 4;
-GRID.pmaxi = 15;
-GRID.jmaxi = 4;
+GRID.pmaxe = 25;
+GRID.jmaxe = 10;
+GRID.pmaxi = 25;
+GRID.jmaxi = 10;
 GRID.nkr   = 1;
 GRID.krmin = 0.;
 GRID.krmax = 0.;
-GRID.nkz   = 9;
+GRID.nkz   = 1;
 GRID.kzmin = 0.1;
-GRID.kzmax = 1.0;
+GRID.kzmax = 0.1;
 %% Model parameters
-MODEL.nu      = 0.001;
+MODEL.nu      = 0.1;
 MODEL.tau_e   = 1.0;
 MODEL.tau_i   = 1.0;
 MODEL.sigma_e = 0.0233380;
@@ -37,9 +37,9 @@ MODEL.eta_B   = 0.5;
 MODEL.lambdaD = 0.0;
 %% Time integration and intialization parameters
 TIME_INTEGRATION.numerical_scheme  = '''RK4''';
-BASIC.nrun              = 100000;
-BASIC.dt                = 0.01;
-BASIC.tmax              = 200;
+BASIC.nrun                = 100000;
+BASIC.dt                  = 0.01;
+BASIC.tmax                = 200;
 INITIAL.initback_moments  = 0.01;
 INITIAL.initnoise_moments = 0.;
 INITIAL.iseed             = 42;
@@ -49,16 +49,93 @@ INITIAL.iseed             = 42;
 INPUT = write_fort90(OUTPUTS,GRID,MODEL,INITIAL,TIME_INTEGRATION,BASIC);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%% Run solver
+%% Run HeLaZ
 nproc = 1;
 EXEC  = ' ../bin/helaz ';
 RUN   = ['mpirun -np ' num2str(nproc)];
 CMD   = [RUN, EXEC, INPUT];
 system(CMD);
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Run MOLI
+MOLI_time_solver
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Analysis and basic figures
 SAVEFIG = 1;
-helaz_analysis
+filename = 'results_00.h5';
+%% Nipj
+
+moment = 'Ni00';
+
+kr       = h5read(filename,['/data/var2d/' moment '/coordkr']);
+kz       = h5read(filename,['/data/var2d/' moment '/coordkz']);
+timeNi     = h5read(filename,'/data/var2d/time');
+Nipj     = zeros(numel(timeNi),numel(kr),numel(kz));
+for it = 1:numel(timeNi)
+    tmp          = h5read(filename,['/data/var2d/', moment,'/', num2str(it,'%06d')]);
+    Nipj(it,:,:) = tmp.real + 1i * tmp.imaginary; 
+end
+
+%% phi
+timephi  = h5read(filename,'/data/var2d/time');
+kr       = h5read(filename,'/data/var2d/phi/coordkr');
+kz       = h5read(filename,'/data/var2d/phi/coordkz');
+phiHeLaZ      = zeros(numel(timephi),numel(kr),numel(kz));
+for it = 1:numel(timephi)
+    tmp         = h5read(filename,['/data/var2d/phi/' num2str(it,'%06d')]);
+    phiHeLaZ(it,:,:) = tmp.real + 1i * tmp.imaginary;
+end
+
+timephiMOLI = results.timeRK4;
+phiMOLI  = zeros(size(timephiMOLI));
+for it = 1:numel(timephiMOLI)
+    phiMOLI(it) = get_phi(results.NapjRK4(it,:),params,options);
+end
+
+
+%% Error
+nsave = OUTPUTS.nsave_2d;
+if(numel(phiMOLI(1:nsave:end)) == numel(phiHeLaZ))
+    errphi  = abs(phiMOLI(1:nsave:end)-phiHeLaZ)./abs(phiMOLI(1:nsave:end));
+    errNipj = abs(results.NapjRK4(1:nsave:end,1)-Nipj)./abs(results.NapjRK4(1:nsave:end,1));
+    figure
+    plot(timephi,errphi*100,'-','DisplayName','$\epsilon(\phi)$')
+    hold on;
+    plot(timephi,errNipj*100,'--','DisplayName','$\epsilon(N_i^{00})$')
+    title(TITLE);
+    xlabel('$t$');
+    ylabel('$\epsilon$ [\%]')
+    grid on
+    legend('show')
+else
+    figure
+    %HeLaZ results
+    plot(timephi,abs(phiHeLaZ),'-','DisplayName','HeLaZ RK4')
+    title(TITLE);
+    hold on
+    %MOLI results
+    plot(timephiMOLI,abs(phiMOLI),'--','DisplayName','MOLI RK4')
+    grid on
+    xlabel('$t$')
+    ylabel('$|\phi|$')
+    legend('show')
+
+    figure
+    %HeLaZ results
+    x1 = timeNi;
+    y1 = abs(Nipj);
+    plot(x1,y1,'-','DisplayName','HeLaZ RK4')
+    hold on
+    %MOLI results
+    x2 = results.timeRK4;
+    y2 = abs(results.NapjRK4(:,1));
+    plot(x2(1:100:end),y2(1:100:end),'--','DisplayName','MOLI RK4');
+    title(TITLE);
+    grid on
+    legend('show')
+    xlabel('$t$')
+    ylabel(['$|',moment,'|$'])
 
+end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
-- 
GitLab