From 5e48668823b5e7613da65795d3201b8b48604fb9 Mon Sep 17 00:00:00 2001
From: Antoine Cyril David Hoffmann <ahoffman@spcpc606.epfl.ch>
Date: Tue, 5 Apr 2022 17:27:10 +0200
Subject: [PATCH] first commit of z parallel version, still buggy
---
Makefile | 16 +-
matlab/compile_results.m | 74 ++-
matlab/compute/compute_fluxtube_growth_rate.m | 72 +++
matlab/compute/ifourier_GENE.m | 57 ++-
.../compute/rotate_c_plane_nxnky_to_nkxny.m | 52 ++
matlab/load/load_gene_data.m | 72 +++
matlab/load/read_namelist.m | 219 +++++++++
matlab/plot/gene_transport_spectrum.m | 6 +-
matlab/plot/plot_fluxes.m | 5 +
matlab/plot/plot_metric.m | 26 +
.../plot_radial_transport_and_spacetime.m | 72 ++-
matlab/plot/show_geometry.m | 10 +-
matlab/plot/spectrum_1D.m | 146 +++---
matlab/process_field.m | 40 ++
src/advance_field.F90 | 4 +-
src/array_mod.F90 | 24 +-
src/auxval.F90 | 1 -
src/basic_mod.F90 | 9 +-
src/calculus_mod.F90 | 190 +++++---
src/closure_mod.F90 | 132 ++++--
src/collision_mod.F90 | 13 +-
src/diagnose.F90 | 75 ++-
src/diagnostics_par_mod.F90 | 4 +-
src/geometry_mod.F90 | 43 +-
src/ghosts_mod.F90 | 162 ++++++-
src/grid_mod.F90 | 210 +++++----
src/inital.F90 | 90 ++--
src/memory.F90 | 76 +--
...ents_eq_rhs.F90 => moments_eq_rhs_mod.F90} | 162 +++----
src/{compute_Sapj.F90 => nonlinear_mod.F90} | 39 +-
src/numerics_mod.F90 | 12 +-
src/poisson.F90 | 28 +-
src/ppinit.F90 | 17 +-
src/prec_const_mod.F90 | 1 +
src/processing_mod.F90 | 444 +++++++++++++-----
src/srcinfo.h | 4 +-
src/srcinfo/srcinfo.h | 4 +-
src/stepon.F90 | 42 +-
src/utility_mod.F90 | 8 +-
wk/analysis_3D.m | 87 ++--
wk/analysis_header.m | 164 +------
wk/analysis_header_2DZP.m | 171 +++++++
wk/benchmark_HeLaZ_gene_transport_zpinch.m | 6 +-
wk/gene_analysis_3D.m | 84 ++++
wk/tmp.txt | 0
45 files changed, 2259 insertions(+), 914 deletions(-)
create mode 100644 matlab/compute/compute_fluxtube_growth_rate.m
create mode 100644 matlab/compute/rotate_c_plane_nxnky_to_nkxny.m
create mode 100644 matlab/load/load_gene_data.m
create mode 100644 matlab/load/read_namelist.m
create mode 100644 matlab/plot/plot_fluxes.m
create mode 100644 matlab/plot/plot_metric.m
rename src/{moments_eq_rhs.F90 => moments_eq_rhs_mod.F90} (61%)
rename src/{compute_Sapj.F90 => nonlinear_mod.F90} (96%)
create mode 100644 wk/analysis_header_2DZP.m
create mode 100644 wk/gene_analysis_3D.m
delete mode 100644 wk/tmp.txt
diff --git a/Makefile b/Makefile
index 8901e87b..8066b0e2 100644
--- a/Makefile
+++ b/Makefile
@@ -55,10 +55,10 @@ $(OBJDIR)/diagnose.o : src/srcinfo.h
FOBJ=$(OBJDIR)/advance_field.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o $(OBJDIR)/basic_mod.o \
$(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o $(OBJDIR)/collision_mod.o \
-$(OBJDIR)/compute_Sapj.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
+$(OBJDIR)/nonlinear_mod.o $(OBJDIR)/control.o $(OBJDIR)/fourier_mod.o \
$(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o $(OBJDIR)/fields_mod.o \
$(OBJDIR)/inital.o $(OBJDIR)/initial_par_mod.o $(OBJDIR)/geometry_mod.o \
-$(OBJDIR)/main.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs.o \
+$(OBJDIR)/main.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o $(OBJDIR)/moments_eq_rhs_mod.o \
$(OBJDIR)/numerics_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/prec_const_mod.o \
$(OBJDIR)/processing_mod.o $(OBJDIR)/readinputs.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o \
$(OBJDIR)/restarts_mod.o $(OBJDIR)/stepon.o $(OBJDIR)/tesend.o $(OBJDIR)/time_integration_mod.o \
@@ -94,8 +94,8 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/collision_mod.o : src/collision_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/collision_mod.F90 -o $@
- $(OBJDIR)/compute_Sapj.o : src/compute_Sapj.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/compute_Sapj.F90 -o $@
+ $(OBJDIR)/nonlinear_mod.o : src/nonlinear_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/nonlinear_mod.F90 -o $@
$(OBJDIR)/control.o : src/control.F90 $(OBJDIR)/auxval.o $(OBJDIR)/geometry_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/ppexit.o $(OBJDIR)/ppinit.o $(OBJDIR)/readinputs.o $(OBJDIR)/tesend.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/control.F90 -o $@
@@ -139,13 +139,13 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/model_mod.o : src/model_mod.F90 $(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/model_mod.F90 -o $@
- $(OBJDIR)/moments_eq_rhs.o : src/moments_eq_rhs.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
- $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs.F90 -o $@
+ $(OBJDIR)/moments_eq_rhs_mod.o : src/moments_eq_rhs_mod.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/calculus_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/moments_eq_rhs_mod.F90 -o $@
$(OBJDIR)/numerics_mod.o : src/numerics_mod.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o $(OBJDIR)/utility_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/numerics_mod.F90 -o $@
- $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
+ $(OBJDIR)/poisson.o : src/poisson.F90 $(OBJDIR)/array_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/poisson.F90 -o $@
$(OBJDIR)/ppexit.o : src/ppexit.F90 $(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o
@@ -166,7 +166,7 @@ $(OBJDIR)/utility_mod.o
$(OBJDIR)/restarts_mod.o : src/restarts_mod.F90 $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/restarts_mod.F90 -o $@
- $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
+ $(OBJDIR)/stepon.o : src/stepon.F90 $(OBJDIR)/initial_par_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/advance_field.o $(OBJDIR)/basic_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/numerics_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/ghosts_mod.o $(OBJDIR)/moments_eq_rhs_mod.o $(OBJDIR)/poisson.o $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
$(OBJDIR)/tesend.o : src/tesend.F90 $(OBJDIR)/basic_mod.o $(OBJDIR)/prec_const_mod.o
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index b511eb92..03aa07db 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -19,6 +19,16 @@ PGAMMA_ = [];
PHI_ = [];
DENS_E_ = [];
DENS_I_ = [];
+UPAR_E_ = [];
+UPAR_I_ = [];
+UPER_E_ = [];
+UPER_I_ = [];
+TPAR_E_ = [];
+TPAR_I_ = [];
+TPER_E_ = [];
+TPER_I_ = [];
+TEMP_E_ = [];
+TEMP_I_ = [];
TEMP_E_ = [];
TEMP_I_ = [];
Ts0D_ = [];
@@ -39,15 +49,14 @@ while(CONTINUE)
DT_SIM = h5readatt(filename,'/data/input','dt');
[Pe, Je, Pi, Ji, kx, ky, z] = load_grid_data(filename);
W_GAMMA = strcmp(h5readatt(filename,'/data/input','write_gamma'),'y');
- W_HF = 0;%strcmp(h5readatt(filename,'/data/input','write_hf' ),'y');
+ W_HF = strcmp(h5readatt(filename,'/data/input','write_hf' ),'y');
W_PHI = strcmp(h5readatt(filename,'/data/input','write_phi' ),'y');
W_NA00 = strcmp(h5readatt(filename,'/data/input','write_Na00' ),'y');
W_NAPJ = strcmp(h5readatt(filename,'/data/input','write_Napj' ),'y');
- W_SAPJ = 0;%strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
+ W_SAPJ = strcmp(h5readatt(filename,'/data/input','write_Sapj' ),'y');
W_DENS = strcmp(h5readatt(filename,'/data/input','write_dens' ),'y');
W_TEMP = strcmp(h5readatt(filename,'/data/input','write_temp' ),'y');
KIN_E = strcmp(h5readatt(filename,'/data/input', 'KIN_E' ),'y');
-% KIN_E = 1;
% Check polynomials degrees
Pe_new= numel(Pe); Je_new= numel(Je);
@@ -127,11 +136,31 @@ while(CONTINUE)
[DENS_I, Ts3D, ~] = load_3D_data(filename, 'dens_i');
DENS_I_ = cat(4,DENS_I_,DENS_I); clear DENS_I
end
+ if 0
+ if KIN_E
+ [UPAR_E, Ts3D, ~] = load_3D_data(filename, 'upar_e');
+ UPAR_E_ = cat(4,UPAR_E_,UPAR_E); clear UPAR_E
+ [UPER_E, Ts3D, ~] = load_3D_data(filename, 'uper_e');
+% UPER_E_ = cat(4,UPER_E_,UPER_E); clear UPER_E
+ end
+ [UPAR_I, Ts3D, ~] = load_3D_data(filename, 'upar_i');
+ UPAR_I_ = cat(4,UPAR_I_,UPAR_I); clear UPAR_I
+ [UPER_I, Ts3D, ~] = load_3D_data(filename, 'uper_i');
+ UPER_I_ = cat(4,UPER_I_,UPER_I); clear UPER_I
+ end
if W_TEMP
if KIN_E
+% [TPAR_E, Ts3D, ~] = load_3D_data(filename, 'Tpar_e');
+% TPAR_E_ = cat(4,TPAR_E_,TPAR_E); clear TPAR_E
+% [TPER_E, Ts3D, ~] = load_3D_data(filename, 'Tper_e');
+% TPER_E_ = cat(4,TPER_E_,TPER_E); clear TPER_E
[TEMP_E, Ts3D, ~] = load_3D_data(filename, 'temp_e');
TEMP_E_ = cat(4,TEMP_E_,TEMP_E); clear TEMP_E
end
+% [TPAR_I, Ts3D, ~] = load_3D_data(filename, 'Tpar_i');
+% TPAR_I_ = cat(4,TPAR_I_,TPAR_I); clear TPAR_I
+% [TPER_I, Ts3D, ~] = load_3D_data(filename, 'Tper_i');
+% TPER_I_ = cat(4,TPER_I_,TPER_I); clear TPER_I
[TEMP_I, Ts3D, ~] = load_3D_data(filename, 'temp_i');
TEMP_I_ = cat(4,TEMP_I_,TEMP_I); clear TEMP_I
end
@@ -178,19 +207,38 @@ while(CONTINUE)
end
%% Build grids
-Nkx = numel(kx); Nky = numel(ky);
+
+Nkx = numel(kx);
+if Nkx > 1
+ dkx = kx(2);
+ Lx = 2*pi/dkx;
+else
+ dkx = 0;
+ Lx = 0;
+end
+[~,ikx0] = min(abs(kx));
+Nx = 2*Nkx-1;
+x = linspace(-Lx/2,Lx/2,Nx+1); x = x(1:end-1);
+
+Nky = numel(ky);
+if Nky > 1
+ dky = ky(2);
+ Ly = 2*pi/dky;
+else
+ dky = 0;
+ Ly = 0;
+end
+[~,iky0] = min(abs(ky));
+Ny = Nky;
+y = linspace(-Ly/2,Ly/2,Ny+1); y = y(1:end-1);
+
+Nz = numel(z);
+
[KY,KX] = meshgrid(ky,kx);
-Lkx = max(kx)-min(kx); Lky = max(ky)-min(ky);
-dkx = Lkx/(Nkx-1); dky = Lky/(Nky-1);
KPERP2 = KY.^2+KX.^2;
-[~,ikx0] = min(abs(kx)); [~,iky0] = min(abs(ky));
-
-Nx = 2*(Nkx-1); Ny = Nky; Nz = numel(z);
-Lx = 2*pi/dkx; Ly = 2*pi/dky;
-dx = Lx/Nx; dy = Ly/Ny; dz = 2*pi/Nz;
-x = dx*(-Nx/2:(Nx/2-1)); Lx = max(x)-min(x);
-y = dy*(-Ny/2:(Ny/2-1)); Ly = max(y)-min(y);
%% Add everything in output structure
+% scaling
+DATA.scale = -(1/Nx/Ny)^2;
% Fields
DATA.GGAMMA_RI = GGAMMA_; DATA.PGAMMA_RI = PGAMMA_; DATA.HFLUX_X = HFLUX_;
DATA.Nipj = Nipj_; DATA.Ni00 = Ni00_; DATA.DENS_I = DENS_I_; DATA.TEMP_I = TEMP_I_;
diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
new file mode 100644
index 00000000..60fb9437
--- /dev/null
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -0,0 +1,72 @@
+function [ linear_gr ] = compute_fluxtube_growth_rate(DATA, TRANGE, PLOT)
+
+% Remove AA part
+if DATA.Nx > 1
+ ikxnz = abs(DATA.PHI(:,1,1,1)) > 0;
+else
+ ikxnz = abs(DATA.PHI(:,1,1,1)) > -1;
+end
+ikynz = abs(DATA.PHI(1,:,1,1)) > 0;
+
+phi = DATA.PHI(ikxnz,ikynz,:,:);
+t = DATA.Ts3D;
+
+[~,its] = min(abs(t-TRANGE(1)));
+[~,ite] = min(abs(t-TRANGE(end)));
+
+w_ky = zeros(sum(ikynz),ite-its);
+ce = zeros(sum(ikynz),ite-its);
+
+is = 1;
+for it = its+1:ite
+ phi_n = phi(:,:,:,it);
+ phi_nm1 = phi(:,:,:,it-1);
+ dt = t(it)-t(it-1);
+ ZS = sum(sum(phi_nm1,1),3);
+
+ wl = log(phi_n./phi_nm1)/dt;
+ w_ky(:,is) = sum(sum(wl.*phi_nm1,1),3)./ZS;
+
+ for iky = 1:numel(w_ky(:,is))
+ ce(iky,is) = abs(sum(sum(abs(w_ky(iky,is)-wl(:,iky,:)).^2.*phi_nm1(:,iky,:),1),3)./ZS(:,iky,:));
+ end
+ is = is + 1;
+end
+[kys, Is] = sort(DATA.ky(ikynz));
+
+linear_gr.trange = t(its:ite);
+linear_gr.g_ky = real(w_ky(Is,:));
+linear_gr.w_ky = imag(w_ky(Is,:));
+linear_gr.ce = abs(ce);
+linear_gr.ky = kys;
+if PLOT >0
+ figure
+ plot(linear_gr.ky,linear_gr.g_ky(:,end),'DisplayName','$Re(\omega_{k_y})$'); hold on;
+ plot(linear_gr.ky,linear_gr.w_ky(:,end),'DisplayName','$Im(\omega_{k_y})$'); hold on;
+ plot(linear_gr.ky,linear_gr.ce (:,end),'DisplayName','$\epsilon$'); hold on;
+ xlim([min(linear_gr.ky) max(linear_gr.ky)]);
+ xlabel('$k_y$');
+ legend('show');
+ if PLOT > 1
+ [YY,XX] = meshgrid(kys,t(its+1:ite));
+ figure
+ subplot(311)
+% imagesc(t(its:ite),kys,real(w_ky)); set(gca,'YDir','normal');
+ pclr = pcolor(XX',YY',real(w_ky)); set(pclr, 'edgecolor','none');
+ xlabel('$t$'); ylabel('$k_y$');
+ title('Re($\omega$)')
+
+ subplot(312)
+ pclr = pcolor(XX',YY',imag(w_ky)); set(pclr, 'edgecolor','none');
+ xlabel('$t$'); ylabel('$k_y$');
+ title('Im($\omega$)')
+
+ subplot(313)
+ pclr = pcolor(XX',YY',abs(w_ky)); set(pclr, 'edgecolor','none');
+ xlabel('$t$'); ylabel('$k_y$');
+ title('|$\omega$|')
+ end
+
+end
+
+end
\ No newline at end of file
diff --git a/matlab/compute/ifourier_GENE.m b/matlab/compute/ifourier_GENE.m
index 15a94dab..f5d98f75 100644
--- a/matlab/compute/ifourier_GENE.m
+++ b/matlab/compute/ifourier_GENE.m
@@ -1,26 +1,47 @@
-function [ field_r ] = ifourier_GENE( field_c, dims )
+function [ field_r ] = ifourier_GENE( field_c )
%IFOURIER_FIELD_GENE method of ifourier used in GENE post processing
% This fourier transform back from the half complex plane to a full real
% space, in 2 or 3D (dim). Compied from computeFXYZ of GENE for
% comparison purpose.
-sz = size(field_c);
-nkx = sz(1);
-field_r = zeros(dims);
-if numel(dims) == 2 %2D
- nx = dims(1); ny = dims(2);
- %note, we need one extra point which we set to zero for the ifft
- spectrumKxKyZ = zeros(nx,ny);
- spectrumKxKyZ(1:nkx,:) = field_c;
- spectrumKxKyZ((nkx):(nx),1) = conj(field_c(nkx:-1:2));
- spectrumKxKyZ((nkx):(nx),2:ny) = conj(field_c(nkx:-1:2,ny:-1:2));
-
-% spectrumKxYZ = ny*ifft(spectrumKxKyZ,[],2);
-% field_r = nx*ifft(spectrumKxYZ,[],1,'symmetric');
-
- spectrumKxYZ = ifft(spectrumKxKyZ,[],2);
- field_r = ifft(spectrumKxYZ,[],1,'symmetric');
-end
+%% Original
+% [nx,nky,nz]=size(field_c);
+% %extend to whole ky by imposing reality condition
+% ny=2*nky-1;
+%
+% if ny~=1
+% %note, we need one extra point which we set to zero for the ifft
+% spectrumKxKyZ=zeros(nx,ny,nz);
+% spectrumKxKyZ(:,1:nky,:)=field_c(:,:,:);
+% spectrumKxKyZ(1,(nky+1):(ny),:)=conj(field_c(1,nky:-1:2,:));
+% spectrumKxKyZ(2:nx,(nky+1):(ny),:)=conj(field_c(nx:-1:2,nky:-1:2,:));
+% else
+% %pad with zeros to interpolate on fine scale
+% ny=20;
+% spectrumKxKyZ=zeros(nx,ny,nz);
+% spectrumKxKyZ(:,2,:)=field_c(:,:,:);
+% end
+%
+% %inverse fft, symmetric as we are using real data
+% spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1);
+% field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric');
+% clear spectrumKxKyZ
+
+%% Adapted for HeLaZ representation
+[nkx,ny,nz]=size(field_c);
+%extend to whole ky by imposing reality condition
+nx=2*nkx-1;
+
+%note, we need one extra point which we set to zero for the ifft
+spectrumKxKyZ=zeros(nx,ny,nz);
+spectrumKxKyZ(1:nkx,:,:)=field_c(:,:,:);
+spectrumKxKyZ((nkx+1):(nx),1,:)=conj(field_c(nkx:-1:2,1,:));
+spectrumKxKyZ((nkx+1):(nx),2:ny,:)=conj(field_c(nkx:-1:2,ny:-1:2,:));
+%inverse fft, symmetric as we are using real data
+spectrumXKyZ=nx*ifft(spectrumKxKyZ,[],1);
+field_r=ny*ifft(spectrumXKyZ,[],2,'symmetric');
+clear spectrumKxKyZ
+
end
diff --git a/matlab/compute/rotate_c_plane_nxnky_to_nkxny.m b/matlab/compute/rotate_c_plane_nxnky_to_nkxny.m
new file mode 100644
index 00000000..0bcdb1e6
--- /dev/null
+++ b/matlab/compute/rotate_c_plane_nxnky_to_nkxny.m
@@ -0,0 +1,52 @@
+function [ data ] = rotate_c_plane_nxnky_to_nkxny( data )
+%to go from nx nky Gene representation to nkx ny HeLaZ one
+
+kx_full = data.kx;
+
+ky_full = [data.ky; -data.ky(end-1:-1:2)];
+
+dens = data.DENS_I; temp = data.TEMP_I; phi = data.PHI;
+
+dims = size(dens); nz = dims(3); nt = dims(4);
+
+nkx = numel(data.kx);
+nky = numel(data.ky);
+nx = nkx;
+ny = 2*numel(data.ky)-1;
+
+%note, we need one extra point which we set to zero for the ifft
+dens_full = zeros(nx,ny,nz,nt);
+dens_full(:,1:nky,:,:) = dens(:,:,:,:);
+dens_full(1,(nky+1):(ny),:,:) = conj(dens(1,nky:-1:2,:,:));
+dens_full(2:nx,(nky+1):(ny),:,:) = conj(dens(nx:-1:2,nky:-1:2,:,:));
+dens_full = cat(2,dens(1:data.Nkx,:,:,:),conj(dens(1:data.Nkx,end-1:-1:2,:,:)));
+
+temp_full = zeros(nx,ny,nz,nt);
+temp_full(:,1:nky,:,:) = temp(:,:,:,:);
+temp_full(1,(nky+1):(ny),:,:) = conj(temp(1,nky:-1:2,:,:));
+temp_full(2:nx,(nky+1):(ny),:,:) = conj(temp(nx:-1:2,nky:-1:2,:,:));
+temp_full = cat(2,temp(1:data.Nkx,:,:,:),conj(temp(1:data.Nkx,end-1:-1:2,:,:)));
+
+phi_full = zeros(nx,ny,nz,nt);
+phi_full(:,1:nky,:,:) = phi(:,:,:,:);
+phi_full(1,(nky+1):(ny),:,:) = conj(phi(1,nky:-1:2,:,:));
+phi_full(2:nx,(nky+1):(ny),:,:) = conj(phi(nx:-1:2,nky:-1:2,:,:));
+phi_full = cat(2, phi(1:data.Nkx,:,:,:),conj( phi(1:data.Nkx,end-1:-1:2,:,:)));
+
+data.DENS_I = dens_full(1:data.Nkx/2+1,:,:,:);
+data.TEMP_I = temp_full(1:data.Nkx/2+1,:,:,:);
+data.PHI = phi_full (1:data.Nkx/2+1,:,:,:);
+data.kx = kx_full(1:data.Nkx/2+1);
+data.ky = ky_full;
+data.Nkx = numel(data.kx);
+data.Nky = numel(data.ky);
+data.Nx = nx+1;
+data.Ny = ny-1;
+
+dkx = data.kx(2); dky = data.ky(2);
+Lx = 2*pi/dkx; Ly = 2*pi/dky;
+x = linspace(-Lx/2,Lx/2,data.Nx+1); x = x(1:end-1);
+y = linspace(-Ly/2,Ly/2,data.Ny+1); y = y(1:end-1);
+data.x = x; data.y = y; data.Lx = Lx; data.Ly = Ly;
+end
+
diff --git a/matlab/load/load_gene_data.m b/matlab/load/load_gene_data.m
new file mode 100644
index 00000000..2fdb98a4
--- /dev/null
+++ b/matlab/load/load_gene_data.m
@@ -0,0 +1,72 @@
+function [ DATA ] = load_gene_data( folder )
+%to load gene data as for HeLaZ results
+namelist = read_namelist([folder,'parameters.dat']);
+DATA.namelist = namelist;
+%% Grid
+file = 'coord.dat.h5';
+DATA.vp = h5read([folder,file],'/coord/vp');
+DATA.Nvp = numel(DATA.vp);
+
+DATA.mu = h5read([folder,file],'/coord/mu');
+DATA.Nmu = numel(DATA.mu);
+
+DATA.kx = h5read([folder,file],'/coord/kx');
+DATA.Nkx = numel(DATA.kx);
+DATA.Nx = DATA.Nkx;
+
+DATA.ky = h5read([folder,file],'/coord/ky');
+DATA.Nky = numel(DATA.ky);
+DATA.Ny = DATA.Nky*2-1;
+
+DATA.z = h5read([folder,file],'/coord/z');
+DATA.Nz = numel(DATA.z);
+
+dkx = DATA.kx(2); dky = DATA.ky(2);
+Lx = 2*pi/dkx; Ly = 2*pi/dky;
+x = linspace(-Lx/2,Lx/2,DATA.Nx+1); x = x(1:end-1);
+y = linspace(-Ly/2,Ly/2,DATA.Ny+1); y = y(1:end-1);
+DATA.x = x; DATA.y = y; DATA.Lx = Lx; DATA.Ly = Ly;
+%% Transport
+file = 'nrg.dat.h5';
+DATA.Ts0D = h5read([folder,file],'/nrgions/time');
+DATA.PGAMMA_RI = h5read([folder,file],'/nrgions/Gamma_es');
+DATA.HFLUX_X = h5read([folder,file],'/nrgions/Q_es');
+
+%% fields and moments
+file = 'mom_ions.dat.h5';
+DATA.Ts3D = h5read([folder,file],'/mom_ions/time');
+DATA.DENS_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
+DATA.TPER_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
+DATA.TPAR_I = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
+DATA.PHI = zeros(DATA.Nkx,DATA.Nky,DATA.Nz,numel(DATA.Ts3D));
+
+for jt = 1:numel(DATA.Ts3D)
+ t = DATA.Ts3D(jt);
+ [~, it] = min(abs(DATA.Ts3D-t));
+
+ tmp = h5read([folder,file],['/mom_ions/dens/',sprintf('%10.10d',it-1)]);
+ DATA.DENS_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+
+ tmp = h5read([folder,file],['/mom_ions/T_par/',sprintf('%10.10d',it-1)]);
+ DATA.TPAR_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+
+ tmp = h5read([folder,file],['/mom_ions/T_perp/',sprintf('%10.10d',it-1)]);
+ DATA.TPER_I(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+
+ tmp = h5read([folder,'field.dat.h5'],['/field/phi/',sprintf('%10.10d',it-1)]);
+ DATA.PHI(:,:,:,it) = tmp.real + 1i*tmp.imaginary;
+
+end
+
+DATA.TEMP_I = (DATA.TPAR_I + 2*DATA.TPER_I)/3.0 - DATA.DENS_I;
+DATA.scale = 1;
+DATA.PARAMS = ['GENE'];
+DATA.param_title = 'GENE';
+DATA.localdir = folder;
+%% Geometry
+CMD = ['tail -n ',num2str(namelist.box.nz0),' ',folder,namelist.geometry.magn_geometry{1},'.dat > tmp.txt'];
+system(CMD);
+DATA.geo_arrays = load('tmp.txt');
+system('rm tmp.txt');
+end
+
diff --git a/matlab/load/read_namelist.m b/matlab/load/read_namelist.m
new file mode 100644
index 00000000..74f64c69
--- /dev/null
+++ b/matlab/load/read_namelist.m
@@ -0,0 +1,219 @@
+
+function S = read_namelist(filename)
+% S = READ_NAMELIST(FILENAME) returns the struct S containg namelists and
+% variables in the file FILENAME organised in hierachical way:
+%
+% |--VAR1
+% |--VAR2
+% |-- NMLST_A--|...
+% | |--VARNa
+% |
+% | |--VAR1
+% |-- NMLST_B--|--VAR2
+% | |...
+% S --| ... |--VARNb
+% |
+% | |--VAR1
+% |-- NMLST_M--|--VAR2
+% |...
+% |--VARNm
+%
+% Note: The function can read multidimensioal variables as well. The
+% function assumes that there is no more than one namelist section per
+% line. At this time there is no syntax checking functionality so the
+% function will crash in case of errors.
+%
+% Example:
+% NMLST = read_namelist('OPTIONS.nam');
+% NMLST.NAM_FRAC.XUNIF_NATURE = 0.1;
+% write_namelist(NMlST, 'MOD_OPTIONS.nam');
+%
+% Written by: Darien Pardinas Diaz (darien.pardinas-diaz@monash.edu)
+% Version: 1.0
+% Date: 16 Dec 2011
+S = struct();
+% Open and read the text file containing the namelists
+fid = fopen(filename,'r');
+c = 0;
+lines = cell(1);
+% Read all the text lines in namelist file
+while ~feof(fid)
+ line = fgetl(fid);
+ % Remove comments if any on the line
+ idx = find(line == '!');
+ if ~isempty(idx),
+ line = line(1:idx(1)-1);
+ end
+ if ~isempty(line),
+ c = c + 1;
+ lines{c} = line;
+ end
+end
+fclose(fid);
+i = 0;
+while i < c;
+ % Find a record
+ i = i + 1;
+ line = lines{i};
+ idx = find(line == '&');
+ if ~isempty(idx), % i.e. a namelist start
+ line = line(idx(1) + 1:end);
+ % find next space
+ idx = find(line == ' ');
+ if ~isempty(idx),
+ namelst = line(1:idx(1) - 1);
+ line = line(idx(1) + 1:end);
+ else
+ namelst = line;
+ line = [];
+ end
+ nmlst_bdy = [];
+ idx = strfind(line,'/');
+ % Get the variable specification section
+ while isempty(idx) && i < c,
+ nmlst_bdy = [nmlst_bdy ' ' line];
+ i = i + 1;
+ line = lines{i};
+ idx = strfind(line,'/');
+ end
+ if ~isempty(idx) && idx(1) > 1,
+ nmlst_bdy = [nmlst_bdy ' ' line];
+ end
+ % Parse current namelist (set of variables)
+ switch namelst
+ case 'info'
+ otherwise
+ S.(namelst) = parse_namelist(nmlst_bdy);
+ end
+ end
+end
+function S = parse_namelist(strng)
+% Internal function to parse the body text of a namelist section.
+% Limitations: the following patterns are prohibited inside the literal
+% strings: '.t.' '.f.' '.true.' '.false.' '(:)'
+% Get all .true., .t. and .false., .f. to T and F
+strng = regexprep(strng,'\.true\.' ,'T','ignorecase');
+strng = regexprep(strng,'\.false\.','F','ignorecase');
+strng = regexprep(strng,'\.t\.','T','ignorecase');
+strng = regexprep(strng,'\.f\.','F','ignorecase');
+% Make evaluable the (:) expression in MATLAB if any
+strng = regexprep(strng, '\(:\)', '(1,:)');
+[strng, islit] = parse_literal_strings([strng ' ']);
+% Find the position of all the '='
+eq_idx = find(strng == '=');
+nvars = length(eq_idx);
+arg_start = eq_idx + 1;
+arg_end = zeros(size(eq_idx));
+vars = cell(nvars,1);
+S = struct;
+% Loop through every variable
+for k = 1:nvars,
+ i = eq_idx(k) - 1;
+ % Move to the left and discard blank spaces
+ while strng(i) == ' ', i = i - 1; end
+ % Now we are over the variable name or closing parentesis
+ j = i;
+ if strng(i) == ')',
+ while strng(i) ~= '(', i = i - 1; end
+ i = i - 1;
+ % Move to the left and discard any possible blank spaces
+ while strng(i) == ' ', i = i - 1; end
+ end
+
+ % Now we are over the last character of the variable name
+ while strng(i) ~= ' ', i = i - 1; end
+
+ if k > 1, arg_end(k - 1) = i; end
+ vars{k} = ['S.' strng(i + 1: j)];
+end
+if numel(arg_end) > 0
+arg_end(end) = length(strng);
+end
+% This variables are used in the eval function to evaluate True/False,
+% so don't remove it!
+T = '.true.';
+F = '.false.';
+% Loop through every variable guess varible type
+for k = 1:nvars,
+ arg = strng(arg_start(k):arg_end(k));
+ arglit = islit(arg_start(k):arg_end(k))';
+ % Remove commas in non literal string...
+ commas = ~arglit & arg == ',';
+ if any(commas)
+ arg(commas) = ' ';
+ end
+
+ if any(arglit),
+ % We are parsing a variable that is literal string
+ arg = ['{' arg '};'];
+ elseif ~isempty(find( arg == 'T' | arg == 'F', 1)),
+ % We are parsing a boolean variable
+ arg = ['{' arg '};'];
+ else
+ % We are parsing a numerical array
+ arg = ['[' arg '];'];
+ end
+ % Eval the modified syntax in Matlab
+ eval([vars{k} ' = ' arg]);
+end
+function [strng, is_lit] = parse_literal_strings(strng)
+% Parse the literal declarations of strings and change to Matlab syntax
+len = length(strng);
+add_squote = []; % Positions to add a scape single quote on syntax
+rem_dquote = []; % Positions to remove a double quote scape on syntax
+i = 1;
+while i < len,
+ if strng(i) == '''', % Opening string with single quote...
+ i = i + 1;
+ while i < len && strng(i) ~= '''' || strcmp(strng(i:i+1),'''''') ,
+ i = i + 1;
+ if strcmp(strng(i-1:i),''''''),
+ i = i + 1;
+ end
+ end
+ end
+ if strng(i) == '"', % Opening string with double quote...
+ strng(i) = ''''; % Change to single quote
+ i = i + 1;
+ while strng(i) ~= '"' || strcmp(strng(i:i+1),'""') && i < len,
+ % Check for a possible sigle quote here
+ if strng(i) == '''',
+ add_squote = [add_squote i];
+ end
+ i = i + 1;
+ if strcmp(strng(i-1:i),'""'),
+ rem_dquote = [rem_dquote i-1];
+ i = i + 1;
+ end
+ end
+ strng(i) = ''''; % Change to single quote
+ end
+ i = i + 1;
+end
+for i = 1:length(add_squote);
+ strng = [strng(1:add_squote(i)) strng(add_squote(i):end)];
+end
+for i = 1:length(rem_dquote);
+ strng = [strng(1:rem_dquote(i)-1) strng(rem_squote(i)+1:end)];
+end
+% Now everything should be in Matlab string syntax
+% Classify syntax as literal or regular expression
+i = 1;
+len = length(strng);
+is_lit = zeros(len,1);
+while i < len,
+ if strng(i) == '''', % Opening string with single quote...
+ is_lit(i) = 1;
+ i = i + 1;
+ while i < len && strng(i) ~= '''' || strcmp(strng(i:i+1),''''''),
+ is_lit(i) = 1;
+ i = i + 1;
+ if strcmp(strng(i-1:i),''''''),
+ is_lit(i) = 1;
+ i = i + 1;
+ end
+ end
+ is_lit(i) = 1;
+ end
+ i = i + 1;
+end
diff --git a/matlab/plot/gene_transport_spectrum.m b/matlab/plot/gene_transport_spectrum.m
index 66563af5..e6244070 100644
--- a/matlab/plot/gene_transport_spectrum.m
+++ b/matlab/plot/gene_transport_spectrum.m
@@ -1,9 +1,9 @@
% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK_36x20/';
% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
-folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
+% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
% folder = '/misc/gene_results/HP_fig_2c_gyroLES/';
% folder = '/misc/gene_results/NL_Zpinch_Kn_1.8_eta_0.25_nuSG_5e-2_mu_1e-2_SGDK_128x32x48x32/';
-
+folder = '/misc/gene_results/shearless_cyclone/output/';
file = 'coord.dat.h5';
vp = h5read([folder,file],'/coord/vp');
mu = h5read([folder,file],'/coord/mu');
@@ -17,7 +17,7 @@ file = 'field.dat.h5';
time = h5read([folder,file],'/field/time');
KDIM = 'ky';
-TIMES = 2500:100:3500;
+TIMES = time;
switch KDIM
case 'kx'
sdim = 2;
diff --git a/matlab/plot/plot_fluxes.m b/matlab/plot/plot_fluxes.m
new file mode 100644
index 00000000..02905b11
--- /dev/null
+++ b/matlab/plot/plot_fluxes.m
@@ -0,0 +1,5 @@
+figure
+subplot(211)
+plot(data.Ts0D,data.PGAMMA_RI');
+subplot(212)
+plot(data.Ts0D,-data.HFLUX_X);
\ No newline at end of file
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
new file mode 100644
index 00000000..b4b9e080
--- /dev/null
+++ b/matlab/plot/plot_metric.m
@@ -0,0 +1,26 @@
+function [ fig ] = plot_metric( data )
+
+names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
+ '$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',...
+ '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'};
+data.geo_arrays = load([data.localdir,'geometry.dat']);
+fig = figure;
+subplot(311)
+ for i = 1:6
+ plot(data.z, data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show'); title('MoNoLiT geometry');
+
+subplot(312)
+ for i = 7:10
+ plot(data.z, data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show');
+
+subplot(313)
+ for i = 11:16
+ plot(data.z, data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(data.z),max(data.z)]); legend('show');
+end
+
diff --git a/matlab/plot/plot_radial_transport_and_spacetime.m b/matlab/plot/plot_radial_transport_and_spacetime.m
index 9a3234df..be4680f6 100644
--- a/matlab/plot/plot_radial_transport_and_spacetime.m
+++ b/matlab/plot/plot_radial_transport_and_spacetime.m
@@ -1,14 +1,17 @@
-function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stfname,Nmvm)
+function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stfname,Nmvm, COMPZ)
%Compute steady radial transport
tend = TAVG_1; tstart = TAVG_0;
[~,its0D] = min(abs(DATA.Ts0D-tstart));
[~,ite0D] = min(abs(DATA.Ts0D-tend));
- SCALE = (1/DATA.Nx/DATA.Ny)^2;
- gamma_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
- gamma_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
+ SCALE = DATA.scale;%(1/DATA.Nx/DATA.Ny)^2;
+ Gx_infty_avg = mean(DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
+ Gx_infty_std = std (DATA.PGAMMA_RI(its0D:ite0D))*SCALE;
+ Qx_infty_avg = mean(DATA.HFLUX_X(its0D:ite0D))*SCALE;
+ Qx_infty_std = std (DATA.HFLUX_X(its0D:ite0D))*SCALE;
[~,ikzf] = max(squeeze(mean(abs(squeeze(DATA.PHI(:,1,1,:))),2)));
Ns3D = numel(DATA.Ts3D);
[KY, KX] = meshgrid(DATA.ky, DATA.kx);
+ z = DATA.z;
%% computations
@@ -16,12 +19,22 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stf
Gx = zeros(DATA.Nx,DATA.Ny,numel(DATA.Ts3D));
for it = 1:numel(DATA.Ts3D)
for iz = 1:DATA.Nz
- Gx(:,:,it) = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)),[DATA.Nx,DATA.Ny])...
- .*ifourier_GENE(DATA.DENS_I(:,:,iz,it),[DATA.Nx,DATA.Ny]);
+ Gx(:,:,it) = Gx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
+ .*ifourier_GENE(DATA.DENS_I(:,:,iz,it));
end
Gx(:,:,it) = Gx(:,:,it)/DATA.Nz;
end
- Gamma_t_mtlb = squeeze(mean(mean(Gx,1),2));
+ Gx_t_mtlb = squeeze(mean(mean(Gx,1),2));
+ % Compute Heat flux from ifft matlab
+ Qx = zeros(DATA.Nx,DATA.Ny,numel(DATA.Ts3D));
+ for it = 1:numel(DATA.Ts3D)
+ for iz = 1:DATA.Nz
+ Qx(:,:,it) = Qx(:,:,it) + ifourier_GENE(-1i*KY.*(DATA.PHI(:,:,iz,it)))...
+ .*ifourier_GENE(DATA.TEMP_I(:,:,iz,it));
+ end
+ Qx(:,:,it) = Qx(:,:,it)/DATA.Nz;
+ end
+ Qx_t_mtlb = squeeze(mean(mean(Qx,1),2));
% zonal vs nonzonal energies for phi(t)
E_Zmode_SK = zeros(1,Ns3D);
@@ -35,16 +48,25 @@ function [FIGURE] = plot_radial_transport_and_spacetime(DATA, TAVG_0, TAVG_1,stf
%% Figure
mvm = @(x) movmean(x,Nmvm);
- FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1200, 600])
+ FIGURE.fig = figure; FIGURE.FIGNAME = ['ZF_transport_drphi','_',DATA.PARAMS]; set(gcf, 'Position', [100, 100, 1000, 600])
subplot(311)
-% yyaxis left
+ yyaxis left
plot(mvm(DATA.Ts0D),mvm(DATA.PGAMMA_RI*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
- plot(mvm(DATA.Ts3D),mvm(Gamma_t_mtlb),'DisplayName','matlab comp.'); hold on;
- plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*gamma_infty_avg, '-k',...
- 'DisplayName',['$\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
- grid on; set(gca,'xticklabel',[]); ylabel('$\Gamma_x$')
- ylim([0,5*abs(gamma_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
- title([DATA.param_title,', $\Gamma^{\infty} = $',num2str(gamma_infty_avg),'$\pm$',num2str(gamma_infty_std)]);
+ plot(mvm(DATA.Ts3D),mvm(Gx_t_mtlb),'DisplayName','matlab comp.'); hold on;
+ plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Gx_infty_avg, '-k',...
+ 'DisplayName',['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$\pm$',num2str(Gx_infty_std)]);
+ ylabel('$\Gamma_x$')
+ ylim([0,5*abs(Gx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
+ yyaxis right
+ plot(mvm(DATA.Ts0D),mvm(DATA.HFLUX_X*SCALE),'DisplayName','$\langle n_i \partial_y\phi \rangle_y$'); hold on;
+ plot(mvm(DATA.Ts3D),mvm(Qx_t_mtlb),'DisplayName','matlab comp.'); hold on;
+ plot(DATA.Ts0D(its0D:ite0D),ones(ite0D-its0D+1,1)*Qx_infty_avg, '-k',...
+ 'DisplayName',['$Q_x^{\infty} = $',num2str(Qx_infty_avg),'$\pm$',num2str(Qx_infty_std)]);
+ ylabel('$Q_x$')
+ ylim([0,5*abs(Qx_infty_avg)]); xlim([DATA.Ts0D(1),DATA.Ts0D(end)]);
+ grid on; set(gca,'xticklabel',[]);
+ title({DATA.param_title,...
+ ['$\Gamma^{\infty} = $',num2str(Gx_infty_avg),'$, Q^{\infty} = $',num2str(Qx_infty_avg)]});
%% radial shear radial profile
% computation
Ns3D = numel(DATA.Ts3D);
@@ -53,29 +75,39 @@ mvm = @(x) movmean(x,Nmvm);
kycut = max(DATA.ky);
kxcut = max(DATA.kx);
LP = (abs(KY)<kycut).*(abs(KX)<kxcut); %Low pass filter
+ switch COMPZ
+ case 'avg'
+ compz = @(f) mean(f,3);
+ nameL = '$\langle$ ';
+ nameR = ' $\rangle_{z}$';
+ otherwise
+ compz = @(f) f(:,:,COMPZ);
+ nameL = '';
+ nameR = [' $(z=',num2str(z(COMPZ)),')$'];
+ end
switch stfname
case 'phi'
phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
for it = 1:numel(DATA.Ts3D)
- phi(:,:,1,it) = ifourier_GENE(DATA.PHI(:,:,1,it),[DATA.Nx,DATA.Ny]);
+ phi(:,:,1,it) = ifourier_GENE(compz(DATA.PHI(:,:,:,it)));
end
f2plot = phi; fname = '$\phi$';
case 'v_y'
dxphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
for it = 1:numel(DATA.Ts3D)
- dxphi(:,:,1,it) = ifourier_GENE(-1i*KX.*(DATA.PHI(:,:,1,it)).*LP,[DATA.Nx,DATA.Ny]);
+ dxphi(:,:,1,it) = ifourier_GENE(-1i*KX.*compz(DATA.PHI(:,:,:,it)).*LP);
end
f2plot = dxphi; fname = '$\langle \partial_x\phi\rangle_y$';
case 'v_x'
dyphi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
for it = 1:numel(DATA.Ts3D)
- dyphi(:,:,1,it) = ifourier_GENE(1i*KY.*(DATA.PHI(:,:,1,it)).*LP,[DATA.Nx,DATA.Ny]);
+ dyphi(:,:,1,it) = ifourier_GENE(1i*KY.*compz(DATA.PHI(:,:,:,it)).*LP);
end
f2plot = dyphi; fname = '$\langle \partial_y\phi\rangle_y$';
case 'szf'
dx2phi = zeros(DATA.Nx,DATA.Ny,1,Ns3D);
for it = 1:numel(DATA.Ts3D)
- dx2phi(:,:,1,it) = ifourier_GENE(-KX.^2.*(DATA.PHI(:,:,1,it)).*LP,[DATA.Nx,DATA.Ny]);
+ dx2phi(:,:,1,it) = ifourier_GENE(-KX.^2.*compz(DATA.PHI(:,:,:,it)).*LP);
end
f2plot = dx2phi; fname = '$\langle \partial_x^2\phi\rangle_y$';
end
@@ -84,7 +116,7 @@ mvm = @(x) movmean(x,Nmvm);
[TY,TX] = meshgrid(DATA.x,DATA.Ts3D);
pclr = pcolor(TX,TY,squeeze(plt(f2plot))');
set(pclr, 'edgecolor','none');
- legend(fname) %colorbar;
+ legend([nameL,fname,nameR]) %colorbar;
caxis(clim*[-1 1]);
cmap = bluewhitered(256);
colormap(cmap)
diff --git a/matlab/plot/show_geometry.m b/matlab/plot/show_geometry.m
index 3794d8c8..926eef86 100644
--- a/matlab/plot/show_geometry.m
+++ b/matlab/plot/show_geometry.m
@@ -64,11 +64,13 @@ FIGURE.fig = figure; FIGURE.FIGNAME = ['geometry','_',DATA.PARAMS]; set(gcf, 'Po
for it_ = 1:numel(OPTIONS.TIME)
subplot(1,numel(OPTIONS.TIME),it_)
%plot magnetic geometry
+ if OPTIONS.PLT_MTOPO
magnetic_topo=surf(x_tor, y_tor, z_tor); hold on;alpha 1.0;%light('Position',[-1 1 1],'Style','local')
- set(magnetic_topo,'edgecolor',[1 1 1]*0.8,'facecolor','none')
+ set(magnetic_topo,'edgecolor',[1 1 1]*0.7,'facecolor','none')
+ end
%plot field line
theta = linspace(-Nturns*pi, Nturns*pi, 512) ; % Poloidal angle
- plot3(Xfl(theta),Yfl(theta),Zfl(theta))
+ plot3(Xfl(theta),Yfl(theta),Zfl(theta)); hold on;
%plot vector basis
theta = DATA.z ; % Poloidal angle
plot3(Xfl(theta),Yfl(theta),Zfl(theta),'ok'); hold on;
@@ -100,9 +102,9 @@ subplot(1,numel(OPTIONS.TIME),it_)
ylim([-0.1 0.1]);
zlim([-0.2 0.2]);
end
- end
+end
%%
- % axis equal
+ axis equal
view([1,-2,1])
end
diff --git a/matlab/plot/spectrum_1D.m b/matlab/plot/spectrum_1D.m
index d57f6755..dea4c3f2 100644
--- a/matlab/plot/spectrum_1D.m
+++ b/matlab/plot/spectrum_1D.m
@@ -1,8 +1,10 @@
function [ FIGURE ] = spectrum_1D( data, options )
options.PLAN = 'kxky';
-options.COMP = 1;
+options.COMP = 'avg';
+options.INTERP = 0;
options.POLARPLOT = 0;
+options.AXISEQUAL = 1;
toplot = process_field(data,options);
t = data.Ts3D; frames = toplot.FRAMES;
@@ -23,6 +25,7 @@ switch options.NAME
yname = '$\sum_{k_y}|n_e|$'; fieldname = 'elec. dens.';
end
+PLOT2D = 0;
switch options.COMPXY
case 'avg'
compx = @(x) mean(x,1);
@@ -34,12 +37,14 @@ switch options.COMPXY
compx = @(x) max(x,1);
compy = @(x) max(x,2);
otherwise
- compx = @(x) x(data.kx==0,:);
- compy = @(x) x(:,data.ky==0);
+ compx = @(x) x(:,:);
+ compy = @(x) x(:,:);
+ PLOT2D= 1;
end
-set(gcf, 'Position', [20 50 5000 2000])
-subplot(1,2,1)
+if ~PLOT2D
+ set(gcf, 'Position', [20 50 5000 2000])
+ subplot(1,2,1)
k = data.kx;
xname = '$k_x$';
@@ -48,82 +53,89 @@ subplot(1,2,1)
shiftx = @(x) x;%(1:nmax);
shifty = @(x) x;%(1:nmax);
switch options.COMPT
- case 'avg'
- it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
- Gk = compy(abs(mean(toplot.FIELD(:,:,it0:it1),3)));
- Gk = squeeze(Gk);
- if options.NORM
- Gk = Gk./max(abs(Gk));
- end
- X = shiftx(k);
- if options.OK
- Y = shifty(Gk)./X;
- else
- Y = shifty(Gk);
- end
- plot(X,Y,'DisplayName','t-averaged')
- otherwise
- for it = 1:numel(toplot.FRAMES)
- Gk = compy(abs(toplot.FIELD(:,:,toplot.FRAMES(it))));
- Gk = squeeze(Gk);
- if options.NORM
- Gk = Gk./max(abs(Gk));
- end
- X = shiftx(k);
- if options.OK
- Y = shifty(Gk)./X;
- else
- Y = shifty(Gk);
- end
- plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
- 'Color',colors(it,:)); hold on;
+ case 'avg'
+ it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
+ Gk = compy(abs(mean(toplot.FIELD(:,:,it0:it1),3)));
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(max(abs(Gk)));
+ end
+ X = shiftx(k);
+ if options.OK
+ Y = shifty(Gk)./X;
+ else
+ Y = shifty(Gk);
+ end
+ plot(X,Y,'DisplayName','t-averaged')
+ otherwise
+ for it = 1:numel(toplot.FRAMES)
+ Gk = compy(abs(toplot.FIELD(:,:,toplot.FRAMES(it))));
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(max(abs(Gk)));
end
+ X = shiftx(k);
+ if options.OK
+ Y = shifty(Gk)./X;
+ else
+ Y = shifty(Gk);
+ end
+ plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
end
grid on
title(['HeLaZ $k_x$ ',fieldname,' spectrum']); legend('show','Location','eastoutside')
xlabel(xname); ylabel(yname)
-
-subplot(1,2,2)
+
+ subplot(1,2,2)
k = data.ky;
xname = '$k_y$';
nmax = floor(data.Nky/2*2/3);
- AA = @(x) x(1:nmax);
- shiftx = @(x) x;%AA(x);
- shifty = @(x) x;%AA(ifftshift(x));
switch options.COMPT
- case 'avg'
- it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
- Gk = compx(abs(mean(toplot.FIELD(:,:,it0:it1),3)))';
- Gk = squeeze(Gk);
- if options.NORM
- Gk = Gk./max(abs(Gk));
- end
- X = k(k>0); X = X(1:end-1);
- Yp= Gk(k>0);
- Ym= Gk(k<0);
- Y = Yp(1:end-1) + Ym(end:-1:1);
- Y = Y(end:-1:1);
- plot(X,Y,'DisplayName','t-averaged')
- otherwise
- for it = 1:numel(toplot.FRAMES)
- Gk = compx(abs(toplot.FIELD(:,:,toplot.FRAMES(it))));
- Gk = squeeze(Gk);
- if options.NORM
- Gk = Gk./max(abs(Gk));
- end
- X = k(k>0); X = X(1:end-1);
- Yp= Gk(k>0);
- Ym= Gk(k<0);
- Y = Yp(1:end-1) + Ym(end:-1:1);
- Y = Y(end:-1:1);
- plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
- 'Color',colors(it,:)); hold on;
+ case 'avg'
+ it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
+ Gk = compx(abs(mean(toplot.FIELD(:,:,it0:it1),3)))';
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(max(abs(Gk)));
+ end
+ X = k(k>0); X = X(1:end-1);
+ Yp= Gk(k>0);
+ Ym= Gk(k<0);
+ Y = Yp(1:end-1) + Ym(end:-1:1);
+ Y = Y(end:-1:1);
+ plot(X,Y,'DisplayName','t-averaged')
+ otherwise
+ for it = 1:numel(toplot.FRAMES)
+ Gk = compx(abs(toplot.FIELD(:,:,toplot.FRAMES(it))));
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(max(abs(Gk)));
end
+ X = k(k>0); X = X(1:end-1);
+ Yp= Gk(k>0);
+ Ym= Gk(k<0);
+ Y = Yp(1:end-1) + Ym(end:-1:1);
+ Y = Y(end:-1:1);
+ plot(X,Y,'DisplayName',['$t=',num2str(t(frames(it))),'$'],...
+ 'Color',colors(it,:)); hold on;
+ end
end
grid on
-% title('HeLaZ $k_y$ transport spectrum');
+ % title('HeLaZ $k_y$ transport spectrum');
legend('show','Location','eastoutside');
xlabel(xname); ylabel(yname)
+else
+ it0 = toplot.FRAMES(1); it1 = toplot.FRAMES(end);
+ Gk = mean(abs(toplot.FIELD(:,:,it0:it1)),3);
+ Gk = squeeze(Gk);
+ if options.NORM
+ Gk = Gk./max(max(abs(Gk)));
+ end
+ pclr = pcolor(toplot.X,toplot.Y,Gk);
+ set(pclr, 'edgecolor','none');
+
end
diff --git a/matlab/process_field.m b/matlab/process_field.m
index ff8d00ee..f19f1b2d 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -151,6 +151,26 @@ switch OPTIONS.NAME
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
+ case 'N_i^{00}'
+ NAME = 'Ni00';
+ if COMPDIM == 3
+ for it = FRAMES
+ tmp = squeeze(compr(DATA.Ni00(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(DATA.Ni00(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
case 'n_e'
NAME = 'ne';
if COMPDIM == 3
@@ -234,6 +254,26 @@ switch OPTIONS.NAME
FIELD(:,:,it) = squeeze(compr(tmp));
end
end
+ case 'T_i'
+ NAME = 'Ti';
+ if COMPDIM == 3
+ for it = FRAMES
+ tmp = squeeze(compr(DATA.TEMP_I(:,:,:,it)));
+ FIELD(:,:,it) = squeeze(process(tmp));
+ end
+ else
+ if REALP
+ tmp = zeros(Nx,Ny,Nz);
+ else
+ tmp = zeros(DATA.Nkx,DATA.Nky,Nz);
+ end
+ for it = FRAMES
+ for iz = 1:numel(DATA.z)
+ tmp(:,:,iz) = squeeze(process(DATA.TEMP_I(:,:,iz,it)));
+ end
+ FIELD(:,:,it) = squeeze(compr(tmp));
+ end
+ end
case 'n_i^{NZ}'
NAME = 'niNZ';
if COMPDIM == 3
diff --git a/src/advance_field.F90 b/src/advance_field.F90
index 071d1217..9ebdba45 100644
--- a/src/advance_field.F90
+++ b/src/advance_field.F90
@@ -32,7 +32,7 @@ CONTAINS
p_int = parray_e(ip)
DO ij=ijs_e,ije_e
IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxe))&
- CALL advance_field(moments_e(ip,ij,:,:,:,:), moments_rhs_e(ip,ij,:,:,:,:))
+ CALL advance_field(moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:), moments_rhs_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:))
ENDDO
ENDDO
ENDIF
@@ -41,7 +41,7 @@ CONTAINS
DO ij=ijs_i,ije_i
j_int = jarray_i(ij)
IF((CLOS .NE. 1) .OR. (ip-1+2*(ij-1)+1 .LE. dmaxi))&
- CALL advance_field(moments_i(ip,ij,:,:,:,:), moments_rhs_i(ip,ij,:,:,:,:))
+ CALL advance_field(moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:), moments_rhs_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,:))
ENDDO
ENDDO
! Execution time end
diff --git a/src/array_mod.F90 b/src/array_mod.F90
index d7765312..2193df1d 100644
--- a/src/array_mod.F90
+++ b/src/array_mod.F90
@@ -47,22 +47,22 @@ MODULE array
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: ynipp1j, ynipm1j, ynipp1jm1, ynipm1jm1 ! mirror lin coeff for non adiab mom
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: zNipm1j, zNipm1jp1, zNipm1jm1 ! mirror lin coeff for adiab mom
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: xphij, xphijp1, xphijm1
- ! Geoemtrical operators
+ ! Geometrical operators
! Curvature
REAL(dp), DIMENSION(:,:,:,:), ALLOCATABLE :: Ckxky ! dimensions: kx, ky, z, odd/even p
! Jacobian
REAL(dp), DIMENSION(:,:), ALLOCATABLE :: Jacobian ! dimensions: z, odd/even p
! Metric
- REAL(dp), DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gyy, gxz, gyz ! dimensions: z, odd/even p
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: gxx, gxy, gxz, gyy, gyz, gzz ! dimensions: z, odd/even p
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: dxdr, dxdZ, Rc, phic, Zc
! derivatives of magnetic field strength
- REAL(dp), DIMENSION(:,:), allocatable :: gradzB ! dimensions: z, odd/even p
- REAL(dp), DIMENSION(:,:), allocatable :: gradxB
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: gradxB, gradyB, gradzB
! Relative magnetic field strength
- REAL(dp), DIMENSION(:,:), allocatable :: hatB
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: hatB
! Relative strength of major radius
- REAL(dp), DIMENSION(:,:), allocatable :: hatR
+ REAL(dp), DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients
- REAL(dp), DIMENSION(:,:) , allocatable :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
+ REAL(dp), DIMENSION(:,:) , ALLOCATABLE :: gradz_coeff ! 1 / [ J_{xyz} \hat{B} ]
! Kernel function evaluation (ij,ikx,iky,iz,odd/even p)
REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_e
REAL(dp), DIMENSION(:,:,:,:,:), ALLOCATABLE :: kernel_i
@@ -78,7 +78,17 @@ MODULE array
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dens_e
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: dens_i
+ ! particle fluid velocity for electron and ions (ikx,iky,iz)
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: upar_e
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: upar_i
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: uper_e
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: uper_i
+
! particle temperature for electron and ions (ikx,iky,iz)
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Tpar_e
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Tpar_i
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Tper_e
+ COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: Tper_i
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp_e
COMPLEX(dp), DIMENSION(:,:,:), ALLOCATABLE :: temp_i
diff --git a/src/auxval.F90 b/src/auxval.F90
index 1d978124..e826615f 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -71,7 +71,6 @@ subroutine auxval
' izs = ', izs , ', ize = ', ize
! write(*,*) 'local kx =', kxarray
! write(*,*) 'local ky =', kyarray
- ! write(*,*) 'local iz =', izarray
IF (my_id .NE. num_procs-1) WRITE (*,*) ''
IF (my_id .EQ. num_procs-1) WRITE(*,*) '------------------------------------------'
ENDIF
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index ecda6e4a..75e2a6d4 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -12,7 +12,7 @@ MODULE basic
real(dp) :: time = 0 ! Current simulation time (Init from restart file)
INTEGER :: comm0 ! Default communicator with a topology
- INTEGER :: comm_p, comm_kx ! Communicators for 1-dim cartesian subgrids of comm0
+ INTEGER :: comm_p, comm_kx, comm_z! Communicators for 1-dim cartesian subgrids of comm0
INTEGER :: commr_p0 ! Communicators along kx for only rank 0 on p
INTEGER :: jobnum = 0 ! Job number
@@ -26,9 +26,12 @@ MODULE basic
INTEGER :: num_procs ! number of MPI processes
INTEGER :: num_procs_p ! Number of processes in p
INTEGER :: num_procs_kx ! Number of processes in r
- INTEGER :: rank_0, rank_p, rank_kx! Ranks in comm0, comm_p, comm_kx
+ INTEGER :: num_procs_z ! Number of processes in z
+ INTEGER :: rank_0 ! Ranks in comm0
+ INTEGER :: rank_p, rank_kx, rank_z! Ranks in comm_p, comm_kx, comm_z
INTEGER :: nbr_L, nbr_R ! Left and right neighbours (along p)
INTEGER :: nbr_T, nbr_B ! Top and bottom neighbours (along kx)
+ INTEGER :: nbr_U, nbr_D ! Upstream and downstream neighbours (along z)
INTEGER :: iframe0d ! counting the number of times 0d datasets are outputed (for diagnose)
INTEGER :: iframe1d ! counting the number of times 1d datasets are outputed (for diagnose)
@@ -83,7 +86,7 @@ CONTAINS
INTEGER :: nargs, fileid, l
LOGICAL :: mlexist
nargs = COMMAND_ARGUMENT_COUNT()
- IF((nargs .EQ. 1) .OR. (nargs .EQ. 3)) THEN
+ IF((nargs .EQ. 1) .OR. (nargs .EQ. 4)) THEN
CALL GET_COMMAND_ARGUMENT(nargs, str, l, ierr)
READ(str(1:l),'(i3)') fileid
WRITE(input_file,'(a,a1,i2.2,a3)') 'fort','_',fileid,'.90'
diff --git a/src/calculus_mod.F90 b/src/calculus_mod.F90
index dd2eea85..c06c9ca4 100644
--- a/src/calculus_mod.F90
+++ b/src/calculus_mod.F90
@@ -8,12 +8,14 @@ MODULE calculus
(/ onetwelfth, -twothird, &
0._dp, &
twothird, -onetwelfth /) ! fd4 centered stencil
+ REAL(dp), dimension(-2:2) :: dz4_usu = &
+ (/ 1._dp, -4._dp, 6._dp, -4._dp, 1._dp /)* 0.0625 ! 4th derivative, 2nd order (for parallel hypdiff)
REAL(dp), dimension(-2:1) :: dz_o2e = &
(/ onetwentyfourth,-nineeighths, nineeighths,-onetwentyfourth /) ! fd4 odd to even stencil
REAL(dp), dimension(-1:2) :: dz_e2o = &
(/ onetwentyfourth,-nineeighths, nineeighths,-onetwentyfourth /) ! fd4 odd to even stencil
- PUBLIC :: simpson_rule_z, interp_z, grad_z
+ PUBLIC :: simpson_rule_z, interp_z, grad_z, grad_z4
CONTAINS
@@ -24,8 +26,8 @@ SUBROUTINE grad_z(target,f,ddzf)
! not staggered : the derivative results must be on the same grid as the field
! staggered : the derivative is computed from a grid to the other
INTEGER, INTENT(IN) :: target
- COMPLEX(dp),dimension( izs:ize ), intent(in) :: f
- COMPLEX(dp),dimension( izs:ize ), intent(out) :: ddzf
+ COMPLEX(dp),dimension( izgs:izge ), intent(in) :: f
+ COMPLEX(dp),dimension( izs:ize ), intent(out) :: ddzf
IF(Nz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
IF(SG) THEN
IF(TARGET .EQ. 0) THEN
@@ -34,60 +36,55 @@ SUBROUTINE grad_z(target,f,ddzf)
CALL grad_z_e2o(f,ddzf)
ENDIF
ELSE ! No staggered grid -> usual centered finite differences
- DO iz = 3,Nz-2
+ DO iz = izs,ize
ddzf(iz) = dz_usu(-2)*f(iz-2) + dz_usu(-1)*f(iz-1) &
+dz_usu( 1)*f(iz+1) + dz_usu( 2)*f(iz+2)
ENDDO
- ! Periodic boundary conditions
- ddzf(Nz-1) = dz_usu(-2)*f(Nz-3) + dz_usu(-1)*f(Nz-2) &
- +dz_usu(+1)*f(Nz ) + dz_usu(+2)*f(1 )
- ddzf(Nz ) = dz_usu(-2)*f(Nz-2) + dz_usu(-1)*f(Nz-1) &
- +dz_usu(+1)*f(1 ) + dz_usu(+2)*f(2 )
- ddzf(1 ) = dz_usu(-2)*f(Nz-1) + dz_usu(-1)*f(Nz ) &
- +dz_usu(+1)*f(2 ) + dz_usu(+2)*f(3 )
- ddzf(2 ) = dz_usu(-2)*f(Nz ) + dz_usu(-1)*f(1 ) &
- +dz_usu(+1)*f(3 ) + dz_usu(+2)*f(4)
ENDIF
ELSE
ddzf = 0._dp
ENDIF
-END SUBROUTINE grad_z
+CONTAINS
+ SUBROUTINE grad_z_o2e(fo,ddzfe) ! Paruta 2018 eq (27)
+ ! gives the gradient of a field from the odd grid to the even one
+ implicit none
+ COMPLEX(dp),dimension(izgs:izge), intent(in) :: fo
+ COMPLEX(dp),dimension( izs:ize ), intent(out) :: ddzfe !
+ DO iz = izs,ize
+ ddzfe(iz) = dz_o2e(-2)*fo(iz-2) + dz_o2e(-1)*fo(iz-1) &
+ +dz_o2e( 0)*fo(iz ) + dz_o2e( 1)*fo(iz+1)
+ ENDDO
+ ddzfe(:) = ddzfe(:) * inv_deltaz
+ END SUBROUTINE grad_z_o2e
-SUBROUTINE grad_z_o2e(fo,dzfe) ! Paruta 2018 eq (27)
- ! gives the value of a field from the odd grid to the even one
- implicit none
- COMPLEX(dp),dimension(1:Nz), intent(in) :: fo
- COMPLEX(dp),dimension(1:Nz), intent(out) :: dzfe !
- DO iz = 3,Nz-1
- dzfe(iz) = dz_o2e(-2)*fo(iz-2) + dz_o2e(-1)*fo(iz-1) &
- +dz_o2e( 0)*fo(iz ) + dz_o2e( 1)*fo(iz+1)
- ENDDO
- ! Periodic boundary conditions
- dzfe(Nz) = dz_o2e(-2)*fo(Nz-2) + dz_o2e(-1)*fo(Nz-1) &
- +dz_o2e( 0)*fo(Nz ) + dz_o2e( 1)*fo(1)
- dzfe(1 ) = dz_o2e(-2)*fo(Nz-1) + dz_o2e(-1)*fo(Nz ) &
- +dz_o2e( 0)*fo(1 ) + dz_o2e( 1)*fo(2)
- dzfe(2 ) = dz_o2e(-2)*fo(Nz ) + dz_o2e(-1)*fo(1 ) &
- +dz_o2e( 0)*fo(2 ) + dz_o2e( 1)*fo(3)
-END SUBROUTINE grad_z_o2e
+ SUBROUTINE grad_z_e2o(fe,ddzfo) ! n2v for Paruta 2018 eq (28)
+ ! gives the gradient of a field from the even grid to the odd one
+ implicit none
+ COMPLEX(dp),dimension(izgs:izge), intent(in) :: fe
+ COMPLEX(dp),dimension( izs:ize ), intent(out) :: ddzfo
+ DO iz = izs,ize
+ ddzfo(iz) = dz_e2o(-1)*fe(iz-1) + dz_e2o(0)*fe(iz ) &
+ +dz_e2o( 1)*fe(iz+1) + dz_e2o(2)*fe(iz+2)
+ ENDDO
+ ddzfo(:) = ddzfo(:) * inv_deltaz
+ END SUBROUTINE grad_z_e2o
+END SUBROUTINE grad_z
-SUBROUTINE grad_z_e2o(fe,dzfo)
- ! gives the value of a field from the even grid to the odd one
+SUBROUTINE grad_z4(f,ddz4f)
implicit none
- COMPLEX(dp),dimension(1:Nz), intent(in) :: fe
- COMPLEX(dp),dimension(1:Nz), intent(out) :: dzfo
- DO iz = 2,Nz-2
- dzfo(iz) = dz_e2o(-1)*fe(iz-1) + dz_e2o(0)*fe(iz ) &
- +dz_e2o( 1)*fe(iz+1) + dz_e2o(2)*fe(iz+2)
- ENDDO
- ! Periodic boundary conditions
- dzfo(Nz-1) = dz_e2o(-1)*fe(Nz-2) + dz_e2o(0)*fe(Nz-1) &
- +dz_e2o( 1)*fe(Nz ) + dz_e2o(2)*fe(1 )
- dzfo(Nz ) = dz_e2o(-1)*fe(Nz-1) + dz_e2o(0)*fe(Nz ) &
- +dz_e2o( 1)*fe(1 ) + dz_e2o(2)*fe(2 )
- dzfo(1 ) = dz_e2o(-1)*fe(Nz ) + dz_e2o(0)*fe(1 ) &
- +dz_e2o( 1)*fe(2 ) + dz_e2o(2)*fe(3 )
-END SUBROUTINE grad_z_e2o
+ ! Compute the second order fourth derivative for periodic boundary condition
+ COMPLEX(dp),dimension(izgs:izge), intent(in) :: f
+ COMPLEX(dp),dimension( izs:ize ), intent(out) :: ddz4f
+ IF(Nz .GT. 3) THEN ! Cannot apply four points stencil on less than four points grid
+ DO iz = izs,ize
+ ddz4f(iz) = dz4_usu(-2)*f(iz-2) + dz4_usu(-1)*f(iz-1) &
+ +dz4_usu( 1)*f(iz+1) + dz4_usu( 2)*f(iz+2)
+ ENDDO
+ ELSE
+ ddz4f = 0._dp
+ ENDIF
+ ddz4f = ddz4f * inv_deltaz**4
+END SUBROUTINE grad_z4
SUBROUTINE interp_z(target,f_in,f_out)
@@ -96,10 +93,10 @@ SUBROUTINE interp_z(target,f_in,f_out)
! If Staggered Grid flag (SG) is false, returns identity
implicit none
INTEGER, intent(in) :: target ! target grid : 0 for even grid, 1 for odd
- COMPLEX(dp),dimension(1:Nz), intent(in) :: f_in
- COMPLEX(dp),dimension(1:Nz), intent(out) :: f_out !
+ COMPLEX(dp),dimension(izgs:izge), intent(in) :: f_in
+ COMPLEX(dp),dimension( izs:ize ), intent(out) :: f_out !
IF(SG) THEN
- IF(TARGET .EQ. 0) THEN
+ IF(target .EQ. 0) THEN
CALL interp_o2e_z(f_in,f_out)
ELSE
CALL interp_e2o_z(f_in,f_out)
@@ -107,38 +104,36 @@ SUBROUTINE interp_z(target,f_in,f_out)
ELSE ! No staggered grid -> identity
f_out(:) = f_in(:)
ENDIF
-END SUBROUTINE interp_z
-
-SUBROUTINE interp_o2e_z(fo,fe)
- ! gives the value of a field from the odd grid to the even one
- implicit none
- COMPLEX(dp),dimension(1:Nz), intent(in) :: fo
- COMPLEX(dp),dimension(1:Nz), intent(out) :: fe !
- DO iz = 2,Nz
- fe(iz) = 0.5_dp*(fo(iz)+fo(iz-1))
- ENDDO
- ! Periodic boundary conditions
- fe(1) = 0.5_dp*(fo(1) + fo(Nz))
-END SUBROUTINE interp_o2e_z
+CONTAINS
+ SUBROUTINE interp_o2e_z(fo,fe)
+ ! gives the value of a field from the odd grid to the even one
+ implicit none
+ COMPLEX(dp),dimension(izgs:izge), intent(in) :: fo
+ COMPLEX(dp),dimension( izs:ize ), intent(out) :: fe !
+ ! 4th order interp
+ DO iz = izs,ize
+ fe(iz) = onesixteenth * (-fo(iz-2) + 9._dp*(fo(iz-1) + fo(iz)) - fo(iz+1))
+ ENDDO
+ END SUBROUTINE interp_o2e_z
-SUBROUTINE interp_e2o_z(fe,fo)
- ! gives the value of a field from the even grid to the odd one
- implicit none
- COMPLEX(dp),dimension(1:Nz), intent(in) :: fe
- COMPLEX(dp),dimension(1:Nz), intent(out) :: fo
- DO iz = 1,Nz-1
- fo(iz) = 0.5_dp*(fe(iz+1)+fe(iz))
- ENDDO
- ! Periodic boundary conditions
- fo(Nz) = 0.5_dp*(fe(1) + fe(Nz))
-END SUBROUTINE interp_e2o_z
+ SUBROUTINE interp_e2o_z(fe,fo)
+ ! gives the value of a field from the even grid to the odd one
+ implicit none
+ COMPLEX(dp),dimension(izgs:izge), intent(in) :: fe
+ COMPLEX(dp),dimension( izs:ize ), intent(out) :: fo
+ ! 4th order interp
+ DO iz = 1,Nz
+ fo(iz) = onesixteenth * (-fe(iz-1) + 9._dp*(fe(iz) + fe(iz+1)) - fe(iz+2))
+ ENDDO
+ END SUBROUTINE interp_e2o_z
+END SUBROUTINE interp_z
SUBROUTINE simpson_rule_z(f,intf)
! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
!from molix BJ Frei
implicit none
- complex(dp),dimension(izs:ize), intent(in) :: f
+ complex(dp),dimension(izgs:izge), intent(in) :: f
COMPLEX(dp), intent(out) :: intf
COMPLEX(dp) :: buff_
IF(Nz .GT. 1) THEN
@@ -159,4 +154,49 @@ SUBROUTINE simpson_rule_z(f,intf)
ENDIF
END SUBROUTINE simpson_rule_z
+! SUBROUTINE simpson_rule_z(f,intf)
+! ! integrate f(z) over z using the simpon's rule. Assume periodic boundary conditions (f(ize+1) = f(izs))
+! !from molix BJ Frei
+! implicit none
+! complex(dp),dimension(izs:ize), intent(in) :: f
+! COMPLEX(dp), intent(out) :: intf
+! COMPLEX(dp) :: buffer(1:2), local_int
+! INTEGER :: root, i_
+!
+! IF(Nz .EQ. 1) THEN !2D zpinch simulations
+! intf = f(izs)
+!
+! ELSE !3D fluxtube
+! IF(mod(Nz,2) .ne. 0 ) THEN
+! ERROR STOP 'Simpson rule: Nz must be an even number !!!!'
+! ENDIF
+! ! Buil local sum using the weights of composite Simpson's rule
+! local_int = 0._dp
+! DO iz = izs,ize
+! local_int = local_int + zweights_SR(iz)*f(iz)
+! ENDDO
+!
+! buffer(2) = local_int
+! root = 0
+! !Gather manually among the rank_z=0 processes and perform the sum
+! intf = 0._dp
+! IF (num_procs_z .GT. 1) THEN
+! !! Everyone sends its local_sum to root = 0
+! IF (rank_z .NE. root) THEN
+! CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_COMPLEX, root, 5678, comm_z, ierr)
+! ELSE
+! ! Recieve from all the other processes
+! DO i_ = 0,num_procs_z-1
+! IF (i_ .NE. rank_kx) &
+! CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_COMPLEX, i_, 5678, comm_z, MPI_STATUS_IGNORE, ierr)
+! intf = intf + buffer(2)
+! ENDDO
+! ENDIF
+! ELSE
+! intf = local_int
+! ENDIF
+! intf = onethird*deltaz*intf
+! ENDIF
+! END SUBROUTINE simpson_rule_z
+
END MODULE calculus
diff --git a/src/closure_mod.F90 b/src/closure_mod.F90
index b1a5a7ef..40be81ee 100644
--- a/src/closure_mod.F90
+++ b/src/closure_mod.F90
@@ -8,7 +8,7 @@ USE fields, ONLY: moments_e, moments_i
USE time_integration, ONLY: updatetlevel
IMPLICIT NONE
-PUBLIC :: apply_closure_model, ghosts_truncation
+PUBLIC :: apply_closure_model
CONTAINS
@@ -19,7 +19,7 @@ SUBROUTINE apply_closure_model
CALL cpu_time(t0_clos)
IF (CLOS .EQ. 0) THEN
! zero truncation, An+1=0 for n+1>nmax only
- CALL ghosts_truncation
+ CALL ghosts_upper_truncation
ELSEIF (CLOS .EQ. 1) THEN
@@ -32,15 +32,15 @@ SUBROUTINE apply_closure_model
DO iky = ikys,ikye
DO iz = izs,ize
IF(KIN_E) THEN
- DO ip = ipsg_e,ipeg_e
- DO ij = ijsg_e,ijeg_e
+ DO ip = ipgs_e,ipge_e
+ DO ij = ijgs_e,ijge_e
IF ( parray_e(ip)+2*jarray_e(ip) .GT. dmaxe) &
moments_e(ip,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
ENDDO
ENDIF
- DO ip = ipsg_i,ipeg_i
- DO ij = ijsg_i,ijeg_i
+ DO ip = ipgs_i,ipge_i
+ DO ij = ijgs_i,ijge_i
IF ( parray_i(ip)+2*jarray_i(ip) .GT. dmaxi) &
moments_i(ip,ij,ikx,iky,iz,updatetlevel) = 0._dp
ENDDO
@@ -49,58 +49,102 @@ SUBROUTINE apply_closure_model
ENDDO
ENDDO
! + ghosts truncation
- CALL ghosts_truncation
+ CALL ghosts_upper_truncation
ELSE
- if(my_id .EQ. 0) write(*,*) '! Closure scheme not found !'
+ ERROR STOP '! Closure scheme not found !'
ENDIF
+ CALL ghosts_lower_truncation
+
CALL cpu_time(t1_clos)
tc_clos = tc_clos + (t1_clos - t0_clos)
END SUBROUTINE apply_closure_model
! Positive Oob indices are approximated with a model
-SUBROUTINE ghosts_truncation
+SUBROUTINE ghosts_upper_truncation
IMPLICIT NONE
! zero truncation, An+1=0 for n+1>nmax
- DO ikx = ikxs,ikxe
- DO iky = ikys,ikye
- DO iz = izs,ize
- IF(KIN_E) THEN
- DO ip = ipsg_e,ipeg_e
- moments_e(ip,ijsg_e,ikx,iky,iz,updatetlevel) = 0._dp
- moments_e(ip,ijeg_e,ikx,iky,iz,updatetlevel) = 0._dp
- ENDDO
- DO ij = ijsg_e,ijeg_e
- moments_e(ipsg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_e(ipeg_e ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
- moments_e(ipsg_e+1,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_e(ipeg_e-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
- ENDIF
- ENDDO
- kernel_e(ijsg_e,ikx,iky,iz,:) = 0._dp
- kernel_e(ijeg_e,ikx,iky,iz,:) = 0._dp
- ENDIF
- DO ip = ipsg_i,ipeg_i
- moments_i(ip,ijsg_i,ikx,iky,iz,updatetlevel) = 0._dp
- moments_i(ip,ijeg_i,ikx,iky,iz,updatetlevel) = 0._dp
- ENDDO
- DO ij = ijsg_i,ijeg_i
- moments_i(ipsg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_i(ipeg_i ,ij,ikx,iky,iz,updatetlevel) = 0._dp
- IF(deltapi .EQ. 1) THEN ! Must truncate the second stencil
- moments_i(ipsg_i+1,ij,ikx,iky,iz,updatetlevel) = 0._dp
- moments_i(ipeg_i-1,ij,ikx,iky,iz,updatetlevel) = 0._dp
- ENDIF
- ENDDO
- kernel_i(ijsg_i,ikx,iky,iz,:) = 0._dp
- kernel_i(ijeg_i,ikx,iky,iz,:) = 0._dp
- ENDDO
+ ! Electrons
+ IF(KIN_E) THEN
+ ! applies only for the process that has largest j index
+ IF(ije_e .EQ. Jmaxe+1) THEN
+ DO ip = ipgs_e,ipge_e
+ moments_e(ip,ije_e+1,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ ENDDO
+ ENDIF
+ ! applies only for the process that has largest p index
+ IF(ipe_e .EQ. Pmaxe+1) THEN
+ DO ij = ijgs_e,ijge_e
+ moments_e(ipe_e+1,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
+ moments_e(ipe_e+2,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ ENDIF
ENDDO
+ ENDIF
+ ENDIF
+
+ ! Ions
+ ! applies only for the process that has largest j index
+ IF(ije_i .EQ. Jmaxi+1) THEN
+ DO ip = ipgs_i,ipge_i
+ moments_i(ip,ije_i+1,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
ENDDO
+ ENDIF
+ ! applies only for the process that has largest p index
+ IF(ipe_i .EQ. Pmaxi+1) THEN
+ DO ij = ijgs_i,ijge_i
+ moments_i(ipe_i+1,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
+ moments_i(ipe_i+2,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ ENDIF
+ ENDDO
+ ENDIF
+
+END SUBROUTINE ghosts_upper_truncation
+
+! Negative OoB indices are 0
+SUBROUTINE ghosts_lower_truncation
+ IMPLICIT NONE
+
+! zero truncation, An=0 for n<0
+ ! Electrons
+ IF(KIN_E) THEN
+ ! applies only for the process that has lowest j index
+ IF(ijs_e .EQ. 1) THEN
+ DO ip = ipgs_e,ipge_e
+ moments_e(ip,ijs_e-1,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ ENDDO
+ ENDIF
+ ! applies only for the process that has lowest p index
+ IF(ips_e .EQ. 1) THEN
+ DO ij = ijgs_e,ijge_e
+ moments_e(ips_e-1,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
+ moments_e(ips_e-2,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ ENDIF
+ ENDDO
+ ENDIF
+ ENDIF
+
+ ! Ions
+ IF(ijs_i .EQ. 1) THEN
+ ! applies only for the process that has lowest j index
+ DO ip = ipgs_i,ipge_i
+ moments_i(ip,ije_i-1,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ ENDDO
+ ENDIF
+ ! applies only for the process that has lowest p index
+ IF(ips_i .EQ. 1) THEN
+ DO ij = ijgs_i,ijge_i
+ moments_i(ips_i-1,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ IF(deltape .EQ. 1) THEN ! Must truncate the second stencil
+ moments_i(ips_i-2,ij,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel) = 0._dp
+ ENDIF
+ ENDDO
+ ENDIF
-END SUBROUTINE ghosts_truncation
+END SUBROUTINE ghosts_lower_truncation
END module closure
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 4cd4a596..451787e7 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -79,15 +79,26 @@ CONTAINS
IMPLICIT NONE
! Execution time start
CALL cpu_time(t0_coll)
+
IF (nu .NE. 0) THEN
SELECT CASE(collision_model)
CASE ('LB')
CALL compute_lenard_bernstein
CASE ('DG')
CALL compute_dougherty
- CASE DEFAULT
+ CASE ('SG','LR','LD')
CALL compute_cosolver_coll
+ CASE ('none')
+ IF(KIN_E) &
+ TColl_e = 0._dp
+ TColl_i = 0._dp
+ CASE DEFAULT
+ ERROR STOP 'Error stop: collision operator not recognized!!'
END SELECT
+ ELSE
+ IF(KIN_E) &
+ TColl_e = 0._dp
+ TColl_i = 0._dp
ENDIF
! Execution time end
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 056661f9..84eb57d3 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -5,6 +5,7 @@ SUBROUTINE diagnose(kstep)
USE grid
USE diagnostics_par
USE futils, ONLY: creatf, creatg, creatd, closef, putarr, putfile, attach, openf, putarrnd
+ USE array
USE model
USE initial_par
USE fields
@@ -70,6 +71,24 @@ SUBROUTINE diagnose(kstep)
CALL putarr(fidres, "/data/grid/coordp_i" , parray_i_full, "p_i", ionode=0)
CALL putarr(fidres, "/data/grid/coordj_i" , jarray_i_full, "j_i", ionode=0)
+ ! Metric info
+ CALL creatg(fidres, "/data/metric", "Metric data")
+ CALL putarrnd(fidres, "/data/metric/gxx", gxx(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gxy", gxy(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gyy", gyy(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gyz", gyz(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gzz", gzz(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatR", hatR(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatZ", hatZ(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/hatB", hatB(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gradxB", gradxB(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gradyB", gradyB(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gradzB", gradzB(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/Jacobian", Jacobian(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/gradz_coeff", gradz_coeff(izs:ize,0:1), (/1, 1/))
+ CALL putarrnd(fidres, "/data/metric/Ckxky", Ckxky(ikxs:ikxe,ikys:ikye,izs:ize,0:1), (/1, 3/))
+ CALL putarrnd(fidres, "/data/metric/kernel_i", kernel_i(ijs_i:ije_i,ikxs:ikxe,ikys:ikye,izs:ize,0:1), (/1, 3/))
+
! var0d group (gyro transport)
IF (nsave_0d .GT. 0) THEN
CALL creatg(fidres, "/data/var0d", "0d profiles")
@@ -121,10 +140,24 @@ SUBROUTINE diagnose(kstep)
CALL creatg(fidres, "/data/var3d/dens_i", "dens_i")
ENDIF
+ IF (write_fvel) THEN
+ IF(KIN_E) THEN
+ CALL creatg(fidres, "/data/var3d/upar_e", "upar_e")
+ CALL creatg(fidres, "/data/var3d/uper_e", "uper_e")
+ ENDIF
+ CALL creatg(fidres, "/data/var3d/upar_i", "upar_i")
+ CALL creatg(fidres, "/data/var3d/uper_i", "uper_i")
+ ENDIF
+
IF (write_temp) THEN
- IF(KIN_E)&
- CALL creatg(fidres, "/data/var3d/temp_e", "temp_e")
- CALL creatg(fidres, "/data/var3d/temp_i", "temp_i")
+ IF(KIN_E) THEN
+ CALL creatg(fidres, "/data/var3d/Tper_e", "Tper_e")
+ CALL creatg(fidres, "/data/var3d/Tpar_e", "Tpar_e")
+ CALL creatg(fidres, "/data/var3d/temp_e", "temp_e")
+ ENDIF
+ CALL creatg(fidres, "/data/var3d/Tper_i", "Tper_i")
+ CALL creatg(fidres, "/data/var3d/Tpar_i", "Tpar_i")
+ CALL creatg(fidres, "/data/var3d/temp_i", "temp_i")
ENDIF
IF (cstep==0) THEN
@@ -190,6 +223,7 @@ SUBROUTINE diagnose(kstep)
CALL attach(fidres, TRIM(str), "write_Napj", write_Napj)
CALL attach(fidres, TRIM(str), "write_Sapj", write_Sapj)
CALL attach(fidres, TRIM(str), "write_dens", write_dens)
+ CALL attach(fidres, TRIM(str), "write_fvel", write_fvel)
CALL attach(fidres, TRIM(str), "write_temp", write_temp)
CALL grid_outputinputs(fidres, str)
@@ -338,7 +372,7 @@ SUBROUTINE diagnose_2d
USE basic
USE futils, ONLY: append, getatt, attach, putarrnd
USE fields
- USE array, ONLY: Ne00, Ni00, dens_e, dens_i, temp_e, temp_i
+ USE array
USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx, ikx, iky, ips_e, ips_i
USE time_integration
USE diagnostics_par
@@ -393,7 +427,7 @@ SUBROUTINE diagnose_3d
USE basic
USE futils, ONLY: append, getatt, attach, putarrnd
USE fields
- USE array, ONLY: Ne00, Ni00, dens_e, dens_i, temp_e, temp_i
+ USE array
USE grid, ONLY: ikxs,ikxe, ikys,ikye, Nkx, Nky, local_nkx, ikx, iky, ips_e, ips_i
USE time_integration
USE diagnostics_par
@@ -411,32 +445,49 @@ SUBROUTINE diagnose_3d
iframe3d=iframe3d+1
CALL attach(fidres,"/data/var3d/" , "frames", iframe3d)
- IF (write_phi) CALL write_field3d(phi (:,:,:), 'phi')
+ IF (write_phi) CALL write_field3d(phi (ikxs:ikxe,ikys:ikye,izs:ize), 'phi')
IF (write_Na00) THEN
IF(KIN_E)THEN
IF (ips_e .EQ. 1) THEN
Ne00(ikxs:ikxe,ikys:ikye,izs:ize) = moments_e(ips_e,1,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
ENDIF
- CALL manual_3D_bcast(Ne00(ikxs:ikxe,ikys:ikye,izs:ize))
+ ! CALL manual_3D_bcast(Ne00(ikxs:ikxe,ikys:ikye,izs:ize))
CALL write_field3d(Ne00(ikxs:ikxe,ikys:ikye,izs:ize), 'Ne00')
ENDIF
IF (ips_i .EQ. 1) THEN
Ni00(ikxs:ikxe,ikys:ikye,izs:ize) = moments_i(ips_e,1,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
ENDIF
- CALL manual_3D_bcast(Ni00(ikxs:ikxe,ikys:ikye,izs:ize))
+ ! CALL manual_3D_bcast(Ni00(ikxs:ikxe,ikys:ikye,izs:ize))
CALL write_field3d(Ni00(ikxs:ikxe,ikys:ikye,izs:ize), 'Ni00')
ENDIF
+ !! Fuid moments
+ IF (write_dens .OR. write_fvel .OR. write_temp) &
+ CALL compute_fluid_moments
+
IF (write_dens) THEN
- CALL compute_density
IF(KIN_E)&
CALL write_field3d(dens_e(ikxs:ikxe,ikys:ikye,izs:ize), 'dens_e')
CALL write_field3d(dens_i(ikxs:ikxe,ikys:ikye,izs:ize), 'dens_i')
ENDIF
+ IF (write_fvel) THEN
+ IF(KIN_E)&
+ CALL write_field3d(upar_e(ikxs:ikxe,ikys:ikye,izs:ize), 'upar_e')
+ CALL write_field3d(upar_i(ikxs:ikxe,ikys:ikye,izs:ize), 'upar_i')
+ IF(KIN_E)&
+ CALL write_field3d(uper_e(ikxs:ikxe,ikys:ikye,izs:ize), 'uper_e')
+ CALL write_field3d(uper_i(ikxs:ikxe,ikys:ikye,izs:ize), 'uper_i')
+ ENDIF
+
IF (write_temp) THEN
- CALL compute_temperature
+ IF(KIN_E)&
+ CALL write_field3d(Tpar_e(ikxs:ikxe,ikys:ikye,izs:ize), 'Tpar_e')
+ CALL write_field3d(Tpar_i(ikxs:ikxe,ikys:ikye,izs:ize), 'Tpar_i')
+ IF(KIN_E)&
+ CALL write_field3d(Tper_e(ikxs:ikxe,ikys:ikye,izs:ize), 'Tper_e')
+ CALL write_field3d(Tper_i(ikxs:ikxe,ikys:ikye,izs:ize), 'Tper_i')
IF(KIN_E)&
CALL write_field3d(temp_e(ikxs:ikxe,ikys:ikye,izs:ize), 'temp_e')
CALL write_field3d(temp_i(ikxs:ikxe,ikys:ikye,izs:ize), 'temp_i')
@@ -459,7 +510,7 @@ CONTAINS
IF (num_procs .EQ. 1) THEN ! no data distribution
CALL putarr(fidres, dset_name, field(ikxs:ikxe, ikys:ikye, izs:ize), ionode=0)
ELSE
- CALL putarrnd(fidres, dset_name, field(ikxs:ikxe, ikys:ikye, izs:ize), (/1, 1/))
+ CALL putarrnd(fidres, dset_name, field(ikxs:ikxe, ikys:ikye, izs:ize), (/1, 3/))
ENDIF
CALL attach(fidres, dset_name, "time", time)
@@ -516,7 +567,7 @@ SUBROUTINE diagnose_5d
IF (num_procs .EQ. 1) THEN
CALL putarr(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,ikys:ikye,izs:ize), ionode=0)
ELSE
- CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,ikys:ikye,izs:ize), (/1,3/))
+ CALL putarrnd(fidres, dset_name, field(ips_e:ipe_e,ijs_e:ije_e,ikxs:ikxe,ikys:ikye,izs:ize), (/1,3,5/))
ENDIF
CALL attach(fidres, dset_name, 'cstep', cstep)
CALL attach(fidres, dset_name, 'time', time)
diff --git a/src/diagnostics_par_mod.F90 b/src/diagnostics_par_mod.F90
index 36f41076..59a7578f 100644
--- a/src/diagnostics_par_mod.F90
+++ b/src/diagnostics_par_mod.F90
@@ -9,7 +9,7 @@ MODULE diagnostics_par
LOGICAL, PUBLIC, PROTECTED :: write_gamma, write_hf ! output particle transport and heat flux
LOGICAL, PUBLIC, PROTECTED :: write_phi, write_Na00
LOGICAL, PUBLIC, PROTECTED :: write_Napj, write_Sapj
- LOGICAL, PUBLIC, PROTECTED :: write_dens, write_temp
+ LOGICAL, PUBLIC, PROTECTED :: write_dens, write_fvel, write_temp
INTEGER, PUBLIC, PROTECTED :: nsave_0d, nsave_1d, nsave_2d, nsave_3d, nsave_5d
@@ -36,7 +36,7 @@ CONTAINS
NAMELIST /OUTPUT_PAR/ nsave_0d, nsave_1d, nsave_2d, nsave_3d, nsave_5d
NAMELIST /OUTPUT_PAR/ write_doubleprecision, write_gamma, write_hf, write_phi
NAMELIST /OUTPUT_PAR/ write_Na00, write_Napj, write_Sapj
- NAMELIST /OUTPUT_PAR/ write_dens, write_temp
+ NAMELIST /OUTPUT_PAR/ write_dens, write_fvel, write_temp
NAMELIST /OUTPUT_PAR/ job2load
READ(lu_in,output_par)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index d8bc22b7..460529a4 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -20,7 +20,8 @@ contains
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
REAL(dp) :: kx,ky
- COMPLEX(dp), DIMENSION(izs:ize) :: integrant
+ COMPLEX(dp), DIMENSION(izgs:izge) :: integrant
+ INTEGER :: fid
!
IF( (Ny .EQ. 1) .AND. (Nz .EQ. 1)) THEN !1D perp linear run
IF( my_id .eq. 0 ) WRITE(*,*) '1D perpendicular geometry'
@@ -34,11 +35,11 @@ contains
! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
! normalized to rhos_
DO eo = 0,1
- DO iz = izs,ize
- DO iky = ikys, ikye
- ky = kyarray(iky)
- DO ikx = ikxs, ikxe
- kx = kxarray(ikx)
+ DO iky = ikys, ikye
+ ky = kyarray(iky)
+ DO ikx = ikxs, ikxe
+ kx = kxarray(ikx)
+ DO iz = izgs,izge
kparray(ikx, iky, iz, eo) = &
SQRT( gxx(iz,eo)*kx**2 + 2._dp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/hatB(iz,eo)
! there is a factor 1/B from the normalization; important to match GENE
@@ -49,9 +50,10 @@ contains
two_third_kpmax = 2._dp/3._dp * MAXVAL(kparray)
!
! Compute the inverse z integrated Jacobian (useful for flux averaging)
- integrant = Jacobian(:,0) ! Convert into complex array
+ integrant = Jacobian(izgs:izge,0) ! Convert into complex array
CALL simpson_rule_z(integrant,iInt_Jacobian)
iInt_Jacobian = 1._dp/iInt_Jacobian ! reverse it
+
END SUBROUTINE eval_magnetic_geometry
!
!--------------------------------------------------------------------------------
@@ -60,19 +62,31 @@ contains
subroutine eval_salphaB_geometry
! evaluate s-alpha geometry model
implicit none
- REAL(dp) :: z, kx, ky
+ REAL(dp) :: z, kx, ky, alpha_MHD
+ alpha_MHD = 0._dp
parity: DO eo = 0,1
- zloop: DO iz = izs,ize
+ zloop: DO iz = izgs,izge
z = zarray(iz,eo)
! metric
gxx(iz,eo) = 1._dp
- gxy(iz,eo) = shear*z
- gyy(iz,eo) = 1._dp + (shear*z)**2
+ gxy(iz,eo) = shear*z - alpha_MHD*SIN(z)
+ gxz(iz,eo) = 0._dp
+ gyy(iz,eo) = 1._dp + (shear*z - alpha_MHD*SIN(z))**2
+ gyz(iz,eo) = 1._dp/(eps + EPSILON(eps)) !avoid 1/0 in Zpinch config
+ gzz(iz,eo) = 0._dp
+ dxdR(iz,eo)= COS(z)
+ dxdZ(iz,eo)= SIN(z)
! Relative strengh of radius
hatR(iz,eo) = 1._dp + eps*COS(z)
+ hatZ(iz,eo) = 1._dp + eps*SIN(z)
+
+ ! toroidal coordinates
+ Rc (iz,eo) = hatR(iz,eo)
+ phic(iz,eo) = z
+ Zc (iz,eo) = hatZ(iz,eo)
! Jacobian
Jacobian(iz,eo) = q0*hatR(iz,eo)
@@ -81,15 +95,16 @@ contains
hatB(iz,eo) = 1._dp / hatR(iz,eo)
! Derivative of the magnetic field strenght
- gradxB(iz,eo) = -COS(z)
- gradzB(iz,eo) = eps * SIN(z) / hatR(iz,eo)
+ gradxB(iz,eo) = -COS(z) ! Gene put a factor hatB^2 or 1/hatR^2 in this
+ gradyB(iz,eo) = 0._dp
+ gradzB(iz,eo) = eps * SIN(z) / hatR(iz,eo) ! Gene put a factor hatB or 1/hatR in this
! Curvature operator
DO iky = ikys, ikye
ky = kyarray(iky)
DO ikx= ikxs, ikxe
kx = kxarray(ikx)
- Ckxky(ikx, iky, iz,eo) = (-SIN(z)*kx - (COS(z) + shear* z* SIN(z))*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
+ Ckxky(ikx, iky, iz,eo) = (-SIN(z)*kx - (COS(z) + (shear*z - alpha_MHD*SIN(z))* SIN(z))*ky) * hatB(iz,eo) ! .. multiply by hatB to cancel the 1/ hatB factor in moments_eqs_rhs.f90 routine
ENDDO
ENDDO
! coefficient in the front of parallel derivative
diff --git a/src/ghosts_mod.F90 b/src/ghosts_mod.F90
index fa0757f9..a50e0fc3 100644
--- a/src/ghosts_mod.F90
+++ b/src/ghosts_mod.F90
@@ -1,6 +1,6 @@
module ghosts
USE basic
-USE fields, ONLY : moments_e, moments_i
+USE fields, ONLY : moments_e, moments_i, phi
USE grid
USE time_integration
USE model, ONLY : KIN_E
@@ -9,11 +9,11 @@ IMPLICIT NONE
INTEGER :: status(MPI_STATUS_SIZE), source, dest, count, ipg
-PUBLIC :: update_ghosts
+PUBLIC :: update_ghosts_p_moments, update_ghosts_z_phi, update_ghosts_z_moments
CONTAINS
-SUBROUTINE update_ghosts
+SUBROUTINE update_ghosts_p_moments
CALL cpu_time(t0_ghost)
IF (num_procs_p .GT. 1) THEN ! Do it only if we share the p
@@ -24,9 +24,10 @@ SUBROUTINE update_ghosts
CALL update_ghosts_p_i
ENDIF
+ CALL update_ghosts_z_moments
CALL cpu_time(t1_ghost)
tc_ghost = tc_ghost + (t1_ghost - t0_ghost)
-END SUBROUTINE update_ghosts
+END SUBROUTINE update_ghosts_p_moments
!Communicate p+1, p+2 moments to left neighboor and p-1, p-2 moments to right one
@@ -44,25 +45,25 @@ END SUBROUTINE update_ghosts
SUBROUTINE update_ghosts_p_e
IMPLICIT NONE
- count = (ijeg_e-ijsg_e+1)*(ikxe-ikxs+1)*(ikye-ikys+1)*(ize-izs+1)
+ count = (ijge_e-ijgs_e+1)*(ikxe-ikxs+1)*(ikye-ikys+1)*(izge-izgs+1)
!!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!!
! Send the last local moment to fill the -1 neighbour ghost
- CALL mpi_sendrecv(moments_e(ipe_e ,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 10, & ! Send to right
- moments_e(ips_e-1,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 10, & ! Recieve from left
+ CALL mpi_sendrecv(moments_e(ipe_e ,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 10, & ! Send to right
+ moments_e(ips_e-1,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 10, & ! Recieve from left
comm0, status, ierr)
IF (deltape .EQ. 1) & ! If we have odd Hermite degrees we need a 2nd order stencil
- CALL mpi_sendrecv(moments_e(ipe_e-1,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 11, & ! Send to right
- moments_e(ips_e-2,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 11, & ! Recieve from left
+ CALL mpi_sendrecv(moments_e(ipe_e-1,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 11, & ! Send to right
+ moments_e(ips_e-2,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 11, & ! Recieve from left
comm0, status, ierr)
!!!!!!!!!!! Send ghost to left neighbour !!!!!!!!!!!!!!!!!!!!!!
- CALL mpi_sendrecv(moments_e(ips_e ,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 12, & ! Send to left
- moments_e(ipe_e+1,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 12, & ! Recieve from right
+ CALL mpi_sendrecv(moments_e(ips_e ,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 12, & ! Send to left
+ moments_e(ipe_e+1,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 12, & ! Recieve from right
comm0, status, ierr)
IF (deltape .EQ. 1) & ! If we have odd Hermite degrees we need a 2nd order stencil
- CALL mpi_sendrecv(moments_e(ips_e+1,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 13, & ! Send to left
- moments_e(ipe_e+2,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 13, & ! Recieve from right
+ CALL mpi_sendrecv(moments_e(ips_e+1,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 13, & ! Send to left
+ moments_e(ipe_e+2,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 13, & ! Recieve from right
comm0, status, ierr)
END SUBROUTINE update_ghosts_p_e
@@ -72,27 +73,144 @@ SUBROUTINE update_ghosts_p_i
IMPLICIT NONE
- count = (ijeg_i-ijsg_i+1)*(ikxe-ikxs+1)*(ikye-ikys+1)*(ize-izs+1) ! Number of elements sent
+ count = (ijge_i-ijgs_i+1)*(ikxe-ikxs+1)*(ikye-ikys+1)*(izge-izgs+1) ! Number of elements sent
!!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!!
- CALL mpi_sendrecv(moments_i(ipe_i ,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 14, &
- moments_i(ips_i-1,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 14, &
+ CALL mpi_sendrecv(moments_i(ipe_i ,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 14, &
+ moments_i(ips_i-1,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 14, &
comm0, status, ierr)
IF (deltapi .EQ. 1) & ! If we have odd Hermite degrees we need a 2nd order stencil
- CALL mpi_sendrecv(moments_i(ipe_i-1,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 15, &
- moments_i(ips_i-2,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 15, &
+ CALL mpi_sendrecv(moments_i(ipe_i-1,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 15, &
+ moments_i(ips_i-2,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 15, &
comm0, status, ierr)
!!!!!!!!!!! Send ghost to left neighbour !!!!!!!!!!!!!!!!!!!!!!
CALL mpi_cart_shift(comm0, 0, -1, source , dest , ierr)
- CALL mpi_sendrecv(moments_i(ips_i ,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 16, &
- moments_i(ipe_i+1,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 16, &
+ CALL mpi_sendrecv(moments_i(ips_i ,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 16, &
+ moments_i(ipe_i+1,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 16, &
comm0, status, ierr)
IF (deltapi .EQ. 1) & ! If we have odd Hermite degrees we need a 2nd order stencil
- CALL mpi_sendrecv(moments_i(ips_i+1,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 17, &
- moments_i(ipe_i+2,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 17, &
+ CALL mpi_sendrecv(moments_i(ips_i+1,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_L, 17, &
+ moments_i(ipe_i+2,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izgs:izge,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_R, 17, &
comm0, status, ierr)
END SUBROUTINE update_ghosts_p_i
+SUBROUTINE update_ghosts_z_moments
+ IMPLICIT NONE
+
+ IF(KIN_E) &
+ CALL update_ghosts_z_e
+ CALL update_ghosts_z_i
+
+END SUBROUTINE update_ghosts_z_moments
+
+!Communicate z+1, z+2 moments to left neighboor and z-1, z-2 moments to right one
+! [a b|C D|e f] : proc n has moments a to f where a,b,e,f are ghosts
+!
+!proc 0: [0 1 2 3 4|5 6]
+! V V ^ ^
+!proc 1: [3 4|5 6 7 8|9 10]
+! V V ^ ^
+!proc 2: [7 8|9 10 11 12|13 14]
+! V V ^ ^
+!proc 3: [11 12|13 14 15 16|17 18]
+! ^ ^
+!Periodic: 0 1
+SUBROUTINE update_ghosts_z_e
+
+ IMPLICIT NONE
+
+ CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
+ IF (num_procs_z .GT. 1) THEN
+ count = (ipge_e-ipgs_e+1)*(ijge_e-ijgs_e+1)*(ikxe-ikxs+1)*(ikye-ikys+1)
+
+ !!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!!
+ ! Send the last local moment to fill the -1 neighbour ghost
+ CALL mpi_sendrecv(moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,ize ,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 20, & ! Send to right
+ moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,ize-1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 20, & ! Recieve from left
+ comm0, status, ierr)
+ CALL mpi_sendrecv(moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,ize-1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 21, & ! Send to right
+ moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,ize-2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 21, & ! Recieve from left
+ comm0, status, ierr)
+
+ !!!!!!!!!!! Send ghost to left neighbour !!!!!!!!!!!!!!!!!!!!!!
+ CALL mpi_sendrecv(moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izs ,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 22, & ! Send to left
+ moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izs+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 22, & ! Recieve from right
+ comm0, status, ierr)
+ CALL mpi_sendrecv(moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izs+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 23, & ! Send to left
+ moments_e(ipgs_e:ipge_e,ijgs_e:ijge_e,ikxs:ikxe,ikys:ikye,izs+2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 23, & ! Recieve from right
+ comm0, status, ierr)
+ ELSE ! still need to perform periodic boundary conditions
+ moments_e(:,:,:,:, 0,updatetlevel) = moments_e(:,:,:,:, Nz,updatetlevel)
+ moments_e(:,:,:,:, -1,updatetlevel) = moments_e(:,:,:,:,Nz-1,updatetlevel)
+ moments_e(:,:,:,:,Nz+1,updatetlevel) = moments_e(:,:,:,:, 1,updatetlevel)
+ moments_e(:,:,:,:,Nz+2,updatetlevel) = moments_e(:,:,:,:, 2,updatetlevel)
+ ENDIF
+
+END SUBROUTINE update_ghosts_z_e
+
+SUBROUTINE update_ghosts_z_i
+
+ IMPLICIT NONE
+ CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
+ IF (num_procs_z .GT. 1) THEN
+ count = (ipge_i-ipgs_i+1)*(ijge_i-ijgs_i+1)*(ikxe-ikxs+1)*(ikye-ikys+1)
+
+ !!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!!
+ ! Send the last local moment to fill the -1 neighbour ghost
+ CALL mpi_sendrecv(moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,ize ,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 30, & ! Send to right
+ moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,ize-1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 30, & ! Recieve from left
+ comm0, status, ierr)
+ CALL mpi_sendrecv(moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,ize-1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 31, & ! Send to right
+ moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,ize-2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 31, & ! Recieve from left
+ comm0, status, ierr)
+
+ !!!!!!!!!!! Send ghost to left neighbour !!!!!!!!!!!!!!!!!!!!!!
+ CALL mpi_sendrecv(moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izs ,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 32, & ! Send to left
+ moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izs+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 32, & ! Recieve from right
+ comm0, status, ierr)
+ CALL mpi_sendrecv(moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izs+1,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_D, 33, & ! Send to left
+ moments_i(ipgs_i:ipge_i,ijgs_i:ijge_i,ikxs:ikxe,ikys:ikye,izs+2,updatetlevel), count, MPI_DOUBLE_COMPLEX, nbr_U, 33, & ! Recieve from right
+ comm0, status, ierr)
+ ELSE ! still need to perform periodic boundary conditions
+ moments_i(:,:,:,:, 0,updatetlevel) = moments_i(:,:,:,:, Nz,updatetlevel)
+ moments_i(:,:,:,:, -1,updatetlevel) = moments_i(:,:,:,:,Nz-1,updatetlevel)
+ moments_i(:,:,:,:,Nz+1,updatetlevel) = moments_i(:,:,:,:, 1,updatetlevel)
+ moments_i(:,:,:,:,Nz+2,updatetlevel) = moments_i(:,:,:,:, 2,updatetlevel)
+ ENDIF
+END SUBROUTINE update_ghosts_z_i
+
+SUBROUTINE update_ghosts_z_phi
+
+ IMPLICIT NONE
+ CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
+ IF (num_procs_z .GT. 1) THEN
+ count = (ikxe-ikxs+1)*(ikye-ikys+1)
+
+ !!!!!!!!!!! Send ghost to right neighbour !!!!!!!!!!!!!!!!!!!!!!
+ ! Send the last local moment to fill the -1 neighbour ghost
+ CALL mpi_sendrecv(phi(ikxs:ikxe,ikys:ikye, ize), count, MPI_DOUBLE_COMPLEX, nbr_U, 40, & ! Send to right
+ phi(ikxs:ikxe,ikys:ikye,ize-1), count, MPI_DOUBLE_COMPLEX, nbr_D, 40, & ! Recieve from left
+ comm0, status, ierr)
+ CALL mpi_sendrecv(phi(ikxs:ikxe,ikys:ikye,ize-1), count, MPI_DOUBLE_COMPLEX, nbr_U, 41, & ! Send to right
+ phi(ikxs:ikxe,ikys:ikye,ize-2), count, MPI_DOUBLE_COMPLEX, nbr_D, 41, & ! Recieve from left
+ comm0, status, ierr)
+ !!!!!!!!!!! Send ghost to left neighbour !!!!!!!!!!!!!!!!!!!!!!
+
+ CALL mpi_sendrecv(phi(ikxs:ikxe,ikys:ikye, izs), count, MPI_DOUBLE_COMPLEX, nbr_U, 42, & ! Send to right
+ phi(ikxs:ikxe,ikys:ikye,izs+1), count, MPI_DOUBLE_COMPLEX, nbr_D, 42, & ! Recieve from left
+ comm0, status, ierr)
+ CALL mpi_sendrecv(phi(ikxs:ikxe,ikys:ikye,izs+1), count, MPI_DOUBLE_COMPLEX, nbr_U, 43, & ! Send to right
+ phi(ikxs:ikxe,ikys:ikye,izs+2), count, MPI_DOUBLE_COMPLEX, nbr_D, 43, & ! Recieve from left
+ comm0, status, ierr)
+
+ ELSE ! still need to perform periodic boundary conditions
+ phi(:,:, 0) = phi(:,:, Nz)
+ phi(:,:, -1) = phi(:,:, Nz-1)
+ phi(:,:,Nz+1) = phi(:,:, 1)
+ phi(:,:,Nz+2) = phi(:,:, 2)
+ ENDIF
+END SUBROUTINE update_ghosts_z_phi
+
END MODULE ghosts
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index 6bb65b00..4f67a86e 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -19,6 +19,8 @@ MODULE grid
INTEGER, PUBLIC, PROTECTED :: Ny = 16 ! Number of total internal grid points in y
REAL(dp), PUBLIC, PROTECTED :: Ly = 1._dp ! vertical length of the spatial box
INTEGER, PUBLIC, PROTECTED :: Nz = 1 ! Number of total perpendicular planes
+ REAL(dp), PUBLIC, PROTECTED :: Npol = 1._dp ! number of poloidal turns
+ INTEGER, PUBLIC, PROTECTED :: Odz = 4 ! order of z interp and derivative schemes
REAL(dp), PUBLIC, PROTECTED :: q0 = 1._dp ! safety factor
REAL(dp), PUBLIC, PROTECTED :: shear = 0._dp ! magnetic field shear
REAL(dp), PUBLIC, PROTECTED :: eps = 0._dp ! inverse aspect ratio
@@ -42,21 +44,29 @@ MODULE grid
! Local and global z grids, 2D since it has to store odd and even grids
REAL(dp), DIMENSION(:,:), ALLOCATABLE, PUBLIC :: zarray
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: zarray_full
- INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: izarray
+ ! local z weights for computing simpson rule
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: zweights_SR
REAL(dp), PUBLIC, PROTECTED :: deltax, deltay, deltaz, inv_deltaz
INTEGER, PUBLIC, PROTECTED :: ixs, ixe, iys, iye, izs, ize
INTEGER, PUBLIC, PROTECTED :: izgs, izge ! ghosts
LOGICAL, PUBLIC, PROTECTED :: SG = .true.! shifted grid flag
INTEGER, PUBLIC :: ir,iz ! counters
+ ! Data about parallel distribution for kx
integer(C_INTPTR_T), PUBLIC :: local_nkx, local_nky
integer(C_INTPTR_T), PUBLIC :: local_nkx_offset, local_nky_offset
INTEGER, PUBLIC :: local_nkp
+ ! "" for p
INTEGER, PUBLIC :: local_np_e, local_np_i
INTEGER, PUBLIC :: total_np_e, total_np_i
integer(C_INTPTR_T), PUBLIC :: local_np_e_offset, local_np_i_offset
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_np_e, counts_np_i
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_np_e, displs_np_i
-
+ ! "" for z
+ INTEGER, PUBLIC :: local_nz
+ INTEGER, PUBLIC :: total_nz
+ integer(C_INTPTR_T), PUBLIC :: local_nz_offset
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: counts_nz
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: displs_nz
! Grids containing position in fourier space
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kxarray, kxarray_full
REAL(dp), DIMENSION(:), ALLOCATABLE, PUBLIC :: kyarray, kyarray_full
@@ -78,13 +88,15 @@ MODULE grid
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC :: jarray_i, jarray_i_full
INTEGER, PUBLIC, PROTECTED :: ips_e,ipe_e, ijs_e,ije_e ! Start and end indices for pol. deg.
INTEGER, PUBLIC, PROTECTED :: ips_i,ipe_i, ijs_i,ije_i
- INTEGER, PUBLIC, PROTECTED :: ipsg_e,ipeg_e, ijsg_e,ijeg_e ! Ghosts start and end indices
- INTEGER, PUBLIC, PROTECTED :: ipsg_i,ipeg_i, ijsg_i,ijeg_i
+ INTEGER, PUBLIC, PROTECTED :: ipgs_e,ipge_e, ijgs_e,ijge_e ! Ghosts start and end indices
+ INTEGER, PUBLIC, PROTECTED :: ipgs_i,ipge_i, ijgs_i,ijge_i
INTEGER, PUBLIC, PROTECTED :: deltape, ip0_e, ip1_e, ip2_e ! Pgrid spacing and moment 0,1,2 index
INTEGER, PUBLIC, PROTECTED :: deltapi, ip0_i, ip1_i, ip2_i
-
+ LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip0_e, CONTAINS_ip0_i
+ LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip1_e, CONTAINS_ip1_i
+ LOGICAL, PUBLIC, PROTECTED :: CONTAINS_ip2_e, CONTAINS_ip2_i
! Usefull inverse numbers
- REAL(dp), PUBLIC, PROTECTED :: inv_Nx, inv_Ny
+ REAL(dp), PUBLIC, PROTECTED :: inv_Nx, inv_Ny, inv_Nz
! Public Functions
PUBLIC :: init_1Dgrid_distr
@@ -106,7 +118,7 @@ CONTAINS
INTEGER :: lu_in = 90 ! File duplicated from STDIN
NAMELIST /GRID/ pmaxe, jmaxe, pmaxi, jmaxi, &
- Nx, Lx, Ny, Ly, Nz, q0, shear, eps, SG
+ Nx, Lx, Ny, Ly, Nz, Npol, q0, shear, eps, SG
READ(lu_in,grid)
!! Compute the maximal degree of full GF moments set
@@ -162,8 +174,8 @@ CONTAINS
local_np_e = ipe_e - ips_e + 1
local_np_i = ipe_i - ips_i + 1
! Ghosts boundaries
- ipsg_e = ips_e - 2/deltape; ipeg_e = ipe_e + 2/deltape;
- ipsg_i = ips_i - 2/deltapi; ipeg_i = ipe_i + 2/deltapi;
+ ipgs_e = ips_e - 2/deltape; ipge_e = ipe_e + 2/deltape;
+ ipgs_i = ips_i - 2/deltapi; ipge_i = ipe_i + 2/deltapi;
! List of shift and local numbers between the different processes (used in scatterv and gatherv)
ALLOCATE(counts_np_e (1:num_procs_p))
ALLOCATE(counts_np_i (1:num_procs_p))
@@ -179,20 +191,45 @@ CONTAINS
ENDDO
! local grid computation
- ALLOCATE(parray_e(ipsg_e:ipeg_e))
- ALLOCATE(parray_i(ipsg_i:ipeg_i))
- DO ip = ipsg_e,ipeg_e
+ CONTAINS_ip0_e = .FALSE.
+ CONTAINS_ip1_e = .FALSE.
+ CONTAINS_ip2_e = .FALSE.
+ CONTAINS_ip0_i = .FALSE.
+ CONTAINS_ip1_i = .FALSE.
+ CONTAINS_ip2_i = .FALSE.
+ ALLOCATE(parray_e(ipgs_e:ipge_e))
+ ALLOCATE(parray_i(ipgs_i:ipge_i))
+ DO ip = ipgs_e,ipge_e
parray_e(ip) = (ip-1)*deltape
- ! Storing indices of particular degrees for DG and fluid moments computations
- IF(parray_e(ip) .EQ. 0) ip0_e = ip
- IF(parray_e(ip) .EQ. 1) ip1_e = ip
- IF(parray_e(ip) .EQ. 2) ip2_e = ip
+ ! Storing indices of particular degrees for fluid moments computations
+ IF(parray_e(ip) .EQ. 0) THEN
+ ip0_e = ip
+ CONTAINS_ip0_e = .TRUE.
+ ENDIF
+ IF(parray_e(ip) .EQ. 1) THEN
+ ip1_e = ip
+ CONTAINS_ip1_e = .TRUE.
+ ENDIF
+ IF(parray_e(ip) .EQ. 2) THEN
+ ip2_e = ip
+ CONTAINS_ip2_e = .TRUE.
+ ENDIF
END DO
- DO ip = ipsg_i,ipeg_i
+ DO ip = ipgs_i,ipge_i
parray_i(ip) = (ip-1)*deltapi
- IF(parray_i(ip) .EQ. 0) ip0_i = ip
- IF(parray_i(ip) .EQ. 1) ip1_i = ip
- IF(parray_i(ip) .EQ. 2) ip2_i = ip
+ ! Storing indices of particular degrees for fluid moments computations
+ IF(parray_i(ip) .EQ. 0) THEN
+ ip0_i = ip
+ CONTAINS_ip0_i = .TRUE.
+ ENDIF
+ IF(parray_i(ip) .EQ. 1) THEN
+ ip1_i = ip
+ CONTAINS_ip1_i = .TRUE.
+ ENDIF
+ IF(parray_i(ip) .EQ. 2) THEN
+ ip2_i = ip
+ CONTAINS_ip2_i = .TRUE.
+ ENDIF
END DO
!DGGK operator uses moments at index p=2 (ip=3) for the p=0 term so the
! process that contains ip=1 MUST contain ip=3 as well for both e and i.
@@ -218,12 +255,12 @@ CONTAINS
ijs_e = 1; ije_e = jmaxe + 1
ijs_i = 1; ije_i = jmaxi + 1
! Ghosts boundaries
- ijsg_e = ijs_e - 1; ijeg_e = ije_e + 1;
- ijsg_i = ijs_i - 1; ijeg_i = ije_i + 1;
- ALLOCATE(jarray_e(ijsg_e:ijeg_e))
- ALLOCATE(jarray_i(ijsg_i:ijeg_i))
- DO ij = ijsg_e,ijeg_e; jarray_e(ij) = ij-1; END DO
- DO ij = ijsg_i,ijeg_i; jarray_i(ij) = ij-1; END DO
+ ijgs_e = ijs_e - 1; ijge_e = ije_e + 1;
+ ijgs_i = ijs_i - 1; ijge_i = ije_i + 1;
+ ALLOCATE(jarray_e(ijgs_e:ijge_e))
+ ALLOCATE(jarray_i(ijgs_i:ijge_i))
+ DO ij = ijgs_e,ijge_e; jarray_e(ij) = ij-1; END DO
+ DO ij = ijgs_i,ijge_i; jarray_i(ij) = ij-1; END DO
! Precomputations
maxj = MAX(jmaxi, jmaxe)
jmaxe_dp = real(jmaxe,dp)
@@ -236,7 +273,6 @@ CONTAINS
USE model, ONLY: LINEARITY
IMPLICIT NONE
INTEGER :: i_
-
Nkx = Nx/2+1 ! Defined only on positive kx since fields are real
! Grid spacings
IF (Nx .EQ. 1) THEN
@@ -248,21 +284,16 @@ CONTAINS
kx_max = Nkx*deltakx
kx_min = deltakx
ENDIF
-
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(kxarray_full(1:Nkx))
DO ikx = 1,Nkx
kxarray_full(ikx) = REAL(ikx-1,dp) * deltakx
END DO
-
!! Parallel distribution
! Start and END indices of grid
- ! ikxs = 1
- ! ikxe = Nkx
ikxs = local_nkx_offset + 1
ikxe = ikxs + local_nkx - 1
ALLOCATE(kxarray(ikxs:ikxe))
-
local_kxmax = 0._dp
! Creating a grid ordered as dk*(0 1 2 3)
DO ikx = ikxs,ikxe
@@ -282,7 +313,6 @@ CONTAINS
contains_kxmax = .true.
ENDIF
END DO
-
! Orszag 2/3 filter
two_third_kxmax = 2._dp/3._dp*deltakx*(Nkx-1)
ALLOCATE(AA_x(ikxs:ikxe))
@@ -363,60 +393,76 @@ CONTAINS
SUBROUTINE set_zgrid
USE prec_const
IMPLICIT NONE
- INTEGER :: i_, ngz
- REAL :: grids_shift
- ! Start and END indices of grid
- izs = 1
- ize = Nz
- ALLOCATE(zarray(izs:ize,0:1))
- IF (Nz .EQ. 1) THEN ! full perp case
- deltaz = 1._dp
- zarray(1,0) = 0._dp
- zarray(1,1) = 0._dp
- SG = .false. ! unique perp plane at z=0
- izgs = izs
- izge = ize
- ELSE ! Parallel dimension exists
- deltaz = 2._dp*PI/REAL(Nz,dp)
- inv_deltaz = 1._dp/deltaz
- IF(SG) THEN ! Shifted grids option
- grids_shift = deltaz/2._dp ! we shift both z grid
- ELSE
- grids_shift = 0._dp
- ENDIF
- DO iz = izs,ize
- ! Even z grid (Napj with p even and phi)
- zarray(iz,0) = REAL((iz-1),dp)*deltaz - PI
- ! Odd z grid (Napj with p odd)
- zarray(iz,1) = REAL((iz-1),dp)*deltaz - PI + grids_shift
- ENDDO
- ! 2 ghosts cell for four point stencil
- izgs = izs-2
- izge = ize+2
+ INTEGER :: i_, fid
+ REAL :: grid_shift, Lz
+ INTEGER :: ip, istart, iend, in
+ total_nz = Nz
+ ! Length of the flux tube (in ballooning angle)
+ Lz = 2_dp*pi*Npol
+ ! Z stepping (#interval = #points since periodic)
+ deltaz = Lz/REAL(Nz,dp)
+ inv_deltaz = 1._dp/deltaz
+ IF (SG) THEN
+ grid_shift = deltaz/2._dp
+ ELSE
+ grid_shift = 0._dp
ENDIF
- if(my_id.EQ.0) write(*,*) '#parallel planes = ', Nz
! Build the full grids on process 0 to diagnose it without comm
ALLOCATE(zarray_full(1:Nz))
- ! z from -pi to pi
- IF (Nz .GT. 1) THEN
- DO iz = 1,Nz
- zarray_full(iz) = deltaz*(iz-1) - PI
- END DO
+ DO iz = 1,total_nz
+ zarray_full(iz) = REAL(iz-1,dp)*deltaz - PI*REAL(Npol,dp)
+ END DO
+ !! Parallel data distribution
+ ! Local data distribution
+ CALL decomp1D(total_nz, num_procs_z, rank_z, izs, ize)
+ local_nz = ize - izs + 1
+ ! Ghosts boundaries (depend on the order of z operators)
+ IF(Nz .EQ. 1) THEN
+ izgs = izs; izge = ize;
+ ELSEIF(Nz .GE. 4) THEN
+ izgs = izs - 2; izge = ize + 2;
ELSE
- zarray_full(1) = 0
+ ERROR STOP 'Error stop: Nz is not appropriate!!'
ENDIF
-
- ! Boundary conditions for FDF ddz derivative
- ! 4 stencil deritative -> 2 ghosts each sides
- ngz = 2
- ALLOCATE(izarray((1-ngz):(Nz+ngz)))
- DO iz = 1,Nz
- izarray(iz) = iz !points to usuall indices
- END DO
- ! Periodic BC for parallel centered finite differences
- izarray(-1) = Nz-1; izarray(0) = Nz;
- izarray(Nz+1) = 1; izarray(Nz+2) = 2;
-
+ ! List of shift and local numbers between the different processes (used in scatterv and gatherv)
+ ALLOCATE(counts_nz (1:num_procs_z))
+ ALLOCATE(displs_nz (1:num_procs_z))
+ DO in = 0,num_procs_z-1
+ CALL decomp1D(total_nz, num_procs_z, in, istart, iend)
+ counts_nz(in+1) = iend-istart+1
+ displs_nz(in+1) = istart-1
+ ENDDO
+ ! Local z array
+ ALLOCATE(zarray(izgs:izge,0:1))
+ DO iz = izgs,izge
+ IF(iz .EQ. 0) THEN
+ zarray(iz,0) = zarray_full(total_nz)
+ zarray(iz,1) = zarray_full(total_nz) + grid_shift
+ ELSEIF(iz .EQ. -1) THEN
+ zarray(iz,0) = zarray_full(total_nz-1)
+ zarray(iz,1) = zarray_full(total_nz-1) + grid_shift
+ ELSEIF(iz .EQ. total_nz + 1) THEN
+ zarray(iz,0) = zarray_full(1)
+ zarray(iz,1) = zarray_full(1) + grid_shift
+ ELSEIF(iz .EQ. total_nz + 2) THEN
+ zarray(iz,0) = zarray_full(2)
+ zarray(iz,1) = zarray_full(2) + grid_shift
+ ELSE
+ zarray(iz,0) = zarray_full(iz)
+ zarray(iz,1) = zarray_full(iz) + grid_shift
+ ENDIF
+ ENDDO
+ ! Weitghs for Simpson rule
+ ALLOCATE(zweights_SR(izs:ize))
+ DO iz = izs,ize
+ IF((iz .EQ. 1) .OR. (iz .EQ. Nz)) THEN
+ zweights_SR(iz) = 1._dp
+ ELSEIF(MODULO(iz-1,2)) THEN
+ zweights_SR(iz) = 4._dp
+ ELSE
+ zweights_SR(iz) = 2._dp
+ ENDIF
+ ENDDO
END SUBROUTINE set_zgrid
SUBROUTINE grid_outputinputs(fidres, str)
diff --git a/src/inital.F90 b/src/inital.F90
index 43cd7a49..d527d74f 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -13,9 +13,10 @@ SUBROUTINE inital
USE closure
USE ghosts
USE restarts
- USE numerics, ONLY: play_with_modes, save_EM_ZF_modes
- USE processing, ONLY: compute_nadiab_moments
- USE model, ONLY: KIN_E, LINEARITY
+ USE numerics, ONLY: play_with_modes, save_EM_ZF_modes
+ USE processing, ONLY: compute_nadiab_moments, compute_fluid_moments
+ USE model, ONLY: KIN_E, LINEARITY
+ USE nonlinear, ONLY: compute_Sapj, nonlinear_init
IMPLICIT NONE
CALL set_updatetlevel(1)
@@ -25,7 +26,7 @@ SUBROUTINE inital
IF ( job2load .GE. 0 ) THEN
IF (my_id .EQ. 0) WRITE(*,*) 'Load moments'
CALL load_moments ! get N_0
-
+ CALL update_ghosts_z_moments
CALL poisson ! compute phi_0=phi(N_0)
! through initialization
ELSE
@@ -34,51 +35,60 @@ SUBROUTINE inital
CASE ('phi')
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy phi'
CALL init_phi
+ CALL update_ghosts_z_phi
! set moments_00 (GC density) with noise and compute phi afterwards
CASE('mom00')
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy gyrocenter density'
CALL init_gyrodens ! init only gyrocenter density
- ! CALL init_moments ! init all moments randomly (unadvised)
- CALL poisson ! get phi_0 = phi(N_0)
+ CALL update_ghosts_z_moments
+ CALL poisson
+ ! init all moments randomly (unadvised)
CASE('allmom')
IF (my_id .EQ. 0) WRITE(*,*) 'Init noisy moments'
CALL init_moments ! init all moments
- CALL poisson ! get phi_0 = phi(N_0)
+ CALL update_ghosts_z_moments
+ CALL poisson
+ ! init a gaussian blob in gyrodens
CASE('blob')
IF (my_id .EQ. 0) WRITE(*,*) '--init a blob'
CALL initialize_blob
- CALL poisson ! get phi_0 = phi(N_0)
+ CALL update_ghosts_z_moments
+ CALL poisson
+ ! init moments 00 with a power law similarly to GENE
CASE('ppj')
IF (my_id .EQ. 0) WRITE(*,*) 'ppj init ~ GENE'
call init_ppj
- CALL poisson ! get phi_0 = phi(N_0)
+ CALL update_ghosts_z_moments
+ CALL poisson
END SELECT
ENDIF
+ ! closure of j>J, p>P and j<0, p<0 moments
+ IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure'
+ CALL apply_closure_model
+ ! ghosts for p parallelization
+ IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication'
+ CALL update_ghosts_p_moments
+ CALL update_ghosts_z_moments
+ !! End of phi and moments initialization
- ! Save zonal and entropy modes
+ ! Save (kx,0) and (0,ky) modes for num exp
CALL save_EM_ZF_modes
-
! Freeze/Wipe some selected modes (entropy,zonal,turbulent)
CALL play_with_modes
- IF (my_id .EQ. 0) WRITE(*,*) 'Apply closure'
- CALL apply_closure_model
-
- IF (my_id .EQ. 0) WRITE(*,*) 'Ghosts communication'
- CALL update_ghosts
+ ! Load the COSOlver collision operator coefficients
+ IF(cosolver_coll) CALL load_COSOlver_mat
- IF (my_id .EQ. 0) WRITE(*,*) 'Computing non adiab moments'
- CALL compute_nadiab_moments
+ !! Preparing auxiliary arrays at initial state
+ ! particle density, fluid velocity and temperature (used in diagnose)
+ IF (my_id .EQ. 0) WRITE(*,*) 'Computing fluid moments'
+ CALL compute_fluid_moments
- !!!!!! Set Sepj, Sipj and dnjs coeff table !!!!!!
- IF (my_id .EQ. 0) WRITE(*,*) 'Init Sapj'
+ ! init auxval for nonlinear
+ CALL nonlinear_init
+ ! compute nonlinear for t=0 diagnostic
CALL compute_Sapj ! compute S_0 = S(phi_0,N_0)
- !!!!!! Load the COSOlver collision operator coefficients !!!!!!
- IF(cosolver_coll) CALL load_COSOlver_mat
- ! Compute collision
- CALL compute_TColl ! compute C_0 = C(N_0)
-
END SUBROUTINE inital
!******************************************************************************!
@@ -379,16 +389,16 @@ SUBROUTINE init_ppj
USE utility, ONLY: checkfield
USE initial_par
USE model, ONLY: KIN_E, LINEARITY
+ USE geometry, ONLY: Jacobian, iInt_Jacobian
IMPLICIT NONE
REAL(dp) :: noise
- REAL(dp) :: kx, ky, sigma, gain, ky_shift
+ REAL(dp) :: kx, ky, sigma_z, amplitude, ky_shift, z
INTEGER, DIMENSION(12) :: iseedarr
- ! Seed random number generator
- iseedarr(:)=iseed
- CALL RANDOM_SEED(PUT=iseedarr+my_id)
+ sigma_z = pi/4.0
+ amplitude = 0.1
!**** Broad noise initialization *******************************************
! Electrons
@@ -401,19 +411,24 @@ SUBROUTINE init_ppj
DO iky=ikys,ikye
ky = kyarray(iky)
DO iz=izs,ize
- IF (kx .NE. 0) THEN
- IF(ky .NE. 0) THEN
+ z = zarray(iz,0)
+ IF (kx .EQ. 0) THEN
+ IF(ky .EQ. 0) THEN
moments_e(ip,ij,ikx,iky,iz,:) = 0._dp
ELSE
moments_e(ip,ij,ikx,iky,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
ENDIF
ELSE
- IF(ky .NE. 0) THEN
+ IF(ky .GT. 0) THEN
moments_e(ip,ij,ikx,iky,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
ELSE
moments_e(ip,ij,ikx,iky,iz,:) = 0.5_dp*(kx_min/(ABS(kx)+kx_min))
ENDIF
ENDIF
+ ! z-dep
+ moments_e(ip,ij,ikx,iky,iz,:) = moments_e(ip,ij,ikx,iky,iz,:) * &
+ ! (1 + exp(-(z/sigma_z)**2/2.0)*sqrt(2.0*sqrt(pi)/sigma_z))
+ (Jacobian(iz,0)*iInt_Jacobian)**2
END DO
END DO
END DO
@@ -439,19 +454,24 @@ SUBROUTINE init_ppj
DO iky=ikys,ikye
ky = kyarray(iky)
DO iz=izs,ize
- IF (kx .NE. 0) THEN
- IF(ky .NE. 0) THEN
+ z = zarray(iz,0)
+ IF (kx .EQ. 0) THEN
+ IF(ky .EQ. 0) THEN
moments_i(ip,ij,ikx,iky,iz,:) = 0._dp
ELSE
moments_i(ip,ij,ikx,iky,iz,:) = 0.5_dp * ky_min/(ABS(ky)+ky_min)
ENDIF
ELSE
- IF(ky .NE. 0) THEN
+ IF(ky .GT. 0) THEN
moments_i(ip,ij,ikx,iky,iz,:) = (kx_min/(ABS(kx)+kx_min))*(ky_min/(ABS(ky)+ky_min))
ELSE
moments_i(ip,ij,ikx,iky,iz,:) = 0.5_dp*(kx_min/(ABS(kx)+kx_min))
ENDIF
ENDIF
+ ! z-dep
+ moments_i(ip,ij,ikx,iky,iz,:) = moments_i(ip,ij,ikx,iky,iz,:) * &
+ ! (1 + exp(-(z/sigma_z)**2/2.0)*sqrt(2.0*sqrt(pi)/sigma_z))
+ (Jacobian(iz,0)*iInt_Jacobian)**2
END DO
END DO
END DO
diff --git a/src/memory.F90 b/src/memory.F90
index b1508b08..24be2ae9 100644
--- a/src/memory.F90
+++ b/src/memory.F90
@@ -13,22 +13,26 @@ SUBROUTINE memory
IMPLICIT NONE
! Electrostatic potential
- CALL allocate_array( phi, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( phi, ikxs,ikxe, ikys,ikye, izgs,izge)
CALL allocate_array( phi_ZF, ikxs,ikxe, izs,ize)
CALL allocate_array( phi_EM, ikys,ikye, izs,ize)
- CALL allocate_array(inv_poisson_op, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array(inv_poisson_op, ikxs,ikxe, ikys,ikye, izgs,izge)
!Electrons arrays
IF(KIN_E) THEN
CALL allocate_array( Ne00, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( dens_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( upar_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( uper_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Tpar_e, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Tper_e, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( temp_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Kernel_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize, 0,1)
- CALL allocate_array( moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
- CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
- CALL allocate_array( nadiab_moments_e, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( moments_e_ZF, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikxs,ikxe, izs,ize)
- CALL allocate_array( moments_e_EM, ipsg_e,ipeg_e, ijsg_e,ijeg_e, ikys,ikye, izs,ize)
+ CALL allocate_array( Kernel_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge, 0,1)
+ CALL allocate_array( moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_e, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
+ CALL allocate_array( nadiab_moments_e, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, ikys,ikye, izgs,izge)
+ CALL allocate_array( moments_e_ZF, ipgs_e,ipge_e, ijgs_e,ijge_e, ikxs,ikxe, izs,ize)
+ CALL allocate_array( moments_e_EM, ipgs_e,ipge_e, ijgs_e,ijge_e, ikys,ikye, izs,ize)
CALL allocate_array( TColl_e, ips_e,ipe_e, ijs_e,ije_e , ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( Sepj, ips_e,ipe_e, ijs_e,ije_e, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( xnepj, ips_e,ipe_e, ijs_e,ije_e)
@@ -48,15 +52,19 @@ SUBROUTINE memory
ENDIF
!Ions arrays
- CALL allocate_array(Ni00, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array(dens_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array(temp_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( Kernel_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize, 0,1)
- CALL allocate_array( moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
- CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
- CALL allocate_array( nadiab_moments_i, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, ikys,ikye, izs,ize)
- CALL allocate_array( moments_i_ZF, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikxs,ikxe, izs,ize)
- CALL allocate_array( moments_i_EM, ipsg_i,ipeg_i, ijsg_i,ijeg_i, ikys,ikye, izs,ize)
+ CALL allocate_array( Ni00, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( dens_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( upar_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( uper_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Tpar_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Tper_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( temp_i, ikxs,ikxe, ikys,ikye, izs,ize)
+ CALL allocate_array( Kernel_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge, 0,1)
+ CALL allocate_array( moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge, 1,ntimelevel )
+ CALL allocate_array( moments_rhs_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize, 1,ntimelevel )
+ CALL allocate_array( nadiab_moments_i, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, ikys,ikye, izgs,izge)
+ CALL allocate_array( moments_i_ZF, ipgs_i,ipge_i, ijgs_i,ijge_i, ikxs,ikxe, izs,ize)
+ CALL allocate_array( moments_i_EM, ipgs_i,ipge_i, ijgs_i,ijge_i, ikys,ikye, izs,ize)
CALL allocate_array( TColl_i, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( Sipj, ips_i,ipe_i, ijs_i,ije_i, ikxs,ikxe, ikys,ikye, izs,ize)
CALL allocate_array( xnipj, ips_i,ipe_i, ijs_i,ije_i)
@@ -83,19 +91,27 @@ SUBROUTINE memory
CALL allocate_array( xphijm1, ips_i,ipe_i, ijs_i,ije_i)
! Curvature and geometry
- CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye, izs,ize, 0,1)
- CALL allocate_array( kparray, ikxs,ikxe, ikys,ikye, izs,ize, 0,1)
- CALL allocate_array( Jacobian, izs,ize, 0,1)
- CALL allocate_array( gxx, izs,ize, 0,1)
- CALL allocate_array( gxy, izs,ize, 0,1)
- CALL allocate_array( gyy, izs,ize, 0,1)
- CALL allocate_array( gyz, izs,ize, 0,1)
- CALL allocate_array( gxz, izs,ize, 0,1)
- CALL allocate_array( gradzB, izs,ize, 0,1)
- CALL allocate_array( gradxB, izs,ize, 0,1)
- CALL allocate_array( hatB, izs,ize, 0,1)
- CALL allocate_array( hatR, izs,ize, 0,1)
- call allocate_array(gradz_coeff, izs,ize, 0,1)
+ CALL allocate_array( Ckxky, ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
+ CALL allocate_array( kparray, ikxs,ikxe, ikys,ikye,izgs,izge,0,1)
+ CALL allocate_array( Jacobian,izgs,izge, 0,1)
+ CALL allocate_array( gxx,izgs,izge, 0,1)
+ CALL allocate_array( gxy,izgs,izge, 0,1)
+ CALL allocate_array( gxz,izgs,izge, 0,1)
+ CALL allocate_array( gyy,izgs,izge, 0,1)
+ CALL allocate_array( gyz,izgs,izge, 0,1)
+ CALL allocate_array( gzz,izgs,izge, 0,1)
+ CALL allocate_array( gradxB,izgs,izge, 0,1)
+ CALL allocate_array( gradyB,izgs,izge, 0,1)
+ CALL allocate_array( gradzB,izgs,izge, 0,1)
+ CALL allocate_array( hatB,izgs,izge, 0,1)
+ CALL allocate_array( hatR,izgs,izge, 0,1)
+ CALL allocate_array( hatZ,izgs,izge, 0,1)
+ CALL allocate_array( Rc,izgs,izge, 0,1)
+ CALL allocate_array( phic,izgs,izge, 0,1)
+ CALL allocate_array( Zc,izgs,izge, 0,1)
+ CALL allocate_array( dxdR,izgs,izge, 0,1)
+ CALL allocate_array( dxdZ,izgs,izge, 0,1)
+ call allocate_array(gradz_coeff,izgs,izge, 0,1)
!___________________ 2x5D ARRAYS __________________________
!! Collision matrices
diff --git a/src/moments_eq_rhs.F90 b/src/moments_eq_rhs_mod.F90
similarity index 61%
rename from src/moments_eq_rhs.F90
rename to src/moments_eq_rhs_mod.F90
index a7d1a71b..10c765b3 100644
--- a/src/moments_eq_rhs.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -1,3 +1,7 @@
+MODULE moments_eq_rhs
+ IMPLICIT NONE
+ PUBLIC :: moments_eq_rhs_e, moments_eq_rhs_i
+CONTAINS
!_____________________________________________________________________________!
!_____________________________________________________________________________!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -15,19 +19,19 @@ SUBROUTINE moments_eq_rhs_e
USE prec_const
USE collision
use geometry
- USE calculus, ONLY : interp_z, grad_z
+ USE calculus, ONLY : interp_z, grad_z, grad_z4
IMPLICIT NONE
INTEGER :: p_int, j_int ! loops indices and polynom. degrees
REAL(dp) :: kx, ky, kperp2, dzlnB_o_J
- COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1, Tpare, Tphi ! Terms from b x gradB and drives
- COMPLEX(dp) :: Tmir, Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments
- COMPLEX(dp) :: UNepm1j, UNepm1jp1, UNepm1jm1 ! Terms from mirror force with adiab moments
+ COMPLEX(dp) :: Tnepj, Tnepp2j, Tnepm2j, Tnepjp1, Tnepjm1 ! Terms from b x gradB and drives
+ COMPLEX(dp) :: Tnepp1j, Tnepm1j, Tnepp1jm1, Tnepm1jm1 ! Terms from mirror force with non adiab moments
+ COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi
+ COMPLEX(dp) :: Unepm1j, Unepm1jp1, Unepm1jm1 ! Terms from mirror force with adiab moments
COMPLEX(dp) :: i_ky
- INTEGER :: izm2, izm1, izp1, izp2 ! indices for centered FDF ddz
! To store derivatives and odd-even z grid interpolations
COMPLEX(dp), DIMENSION(izs:ize) :: ddznepp1j, ddznepm1j, &
- nepp1j, nepp1jm1, nepm1j, nepm1jm1, nepm1jp1
+ nepp1j, nepp1jm1, nepm1j, nepm1jm1, nepm1jp1, ddz4Nepj
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -46,22 +50,22 @@ SUBROUTINE moments_eq_rhs_e
i_ky = imagu * ky ! toroidal derivative
IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
! Compute z derivatives and odd-even z interpolations
- CALL grad_z(eo,nadiab_moments_e(ip+1,ij ,ikx,iky,:),ddznepp1j(:))
- CALL grad_z(eo,nadiab_moments_e(ip-1,ij ,ikx,iky,:),ddznepm1j(:))
- CALL interp_z(eo,nadiab_moments_e(ip+1,ij ,ikx,iky,:), nepp1j (:))
- CALL interp_z(eo,nadiab_moments_e(ip+1,ij-1,ikx,iky,:), nepp1jm1(:))
- CALL interp_z(eo,nadiab_moments_e(ip-1,ij ,ikx,iky,:), nepm1j (:))
- CALL interp_z(eo,nadiab_moments_e(ip-1,ij+1,ikx,iky,:), nepm1jp1(:))
- CALL interp_z(eo,nadiab_moments_e(ip-1,ij-1,ikx,iky,:), nepm1jm1(:))
+ CALL grad_z(eo,nadiab_moments_e(ip+1,ij ,ikx,iky,izgs:izge), ddznepp1j(izs:ize))
+ CALL grad_z(eo,nadiab_moments_e(ip-1,ij ,ikx,iky,izgs:izge), ddznepm1j(izs:ize))
+ CALL interp_z(eo,nadiab_moments_e(ip+1,ij ,ikx,iky,izgs:izge), nepp1j (izs:ize))
+ CALL interp_z(eo,nadiab_moments_e(ip+1,ij-1,ikx,iky,izgs:izge), nepp1jm1(izs:ize))
+ CALL interp_z(eo,nadiab_moments_e(ip-1,ij ,ikx,iky,izgs:izge), nepm1j (izs:ize))
+ CALL interp_z(eo,nadiab_moments_e(ip-1,ij+1,ikx,iky,izgs:izge), nepm1jp1(izs:ize))
+ CALL interp_z(eo,nadiab_moments_e(ip-1,ij-1,ikx,iky,izgs:izge), nepm1jm1(izs:ize))
+ ! Parallel hyperdiffusion
+ CALL grad_z4(moments_e(ip,ij,ikx,iky,:,updatetlevel), ddz4Nepj(:))
+
zloope : DO iz = izs,ize
- ! Obtain the index with an array that accounts for boundary conditions
- ! e.g. : 4 stencil with periodic BC, izarray(Nz+2) = 2, izarray(-1) = Nz-1
- izp1 = izarray(iz+1); izp2 = izarray(iz+2);
- izm1 = izarray(iz-1); izm2 = izarray(iz-2);
- !
+ ! kperp
kperp2= kparray(ikx,iky,iz,eo)**2
!! Compute moments mixing terms
+ Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
! Perpendicular dynamic
! term propto n_e^{p,j}
Tnepj = xnepj(ip,ij)* nadiab_moments_e(ip,ij,ikx,iky,iz)
@@ -74,21 +78,20 @@ SUBROUTINE moments_eq_rhs_e
! term propto n_e^{p,j-1}
Tnepjm1 = xnepjm1(ij) * nadiab_moments_e(ip,ij-1,ikx,iky,iz)
! Parallel dynamic
- Tpare = 0._dp; Tmir = 0._dp
IF(Nz .GT. 1) THEN
- ! ddz derivative for Landau damping term
- Tpare = xnepp1j(ip) * ddznepp1j(iz) + xnepm1j(ip) * ddznepm1j(iz)
- ! Mirror terms
- Tnepp1j = ynepp1j (ip,ij) * nepp1j (iz)
- Tnepp1jm1 = ynepp1jm1(ip,ij) * nepp1jm1(iz)
- Tnepm1j = ynepm1j (ip,ij) * nepm1j (iz)
- Tnepm1jm1 = ynepm1jm1(ip,ij) * nepm1jm1(iz)
- ! Trapping terms
- UNepm1j = znepm1j (ip,ij) * nepm1j (iz)
- UNepm1jp1 = znepm1jp1(ip,ij) * nepm1jp1(iz)
- UNepm1jm1 = znepm1jm1(ip,ij) * nepm1jm1(iz)
+ ! ddz derivative for Landau damping term
+ Tpar = xnepp1j(ip) * ddznepp1j(iz) + xnepm1j(ip) * ddznepm1j(iz)
+ ! Mirror terms (trapping)
+ Tnepp1j = ynepp1j (ip,ij) * nepp1j (iz)
+ Tnepp1jm1 = ynepp1jm1(ip,ij) * nepp1jm1(iz)
+ Tnepm1j = ynepm1j (ip,ij) * nepm1j (iz)
+ Tnepm1jm1 = ynepm1jm1(ip,ij) * nepm1jm1(iz)
+ ! Trapping terms
+ Unepm1j = znepm1j (ip,ij) * nepm1j (iz)
+ Unepm1jp1 = znepm1jp1(ip,ij) * nepm1jp1(iz)
+ Unepm1jm1 = znepm1jm1(ip,ij) * nepm1jm1(iz)
- Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + UNepm1j + UNepm1jp1 + UNepm1jm1
+ Tmir = Tnepp1j + Tnepp1jm1 + Tnepm1j + Tnepm1jm1 + Unepm1j + Unepm1jp1 + Unepm1jm1
ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
@@ -99,29 +102,22 @@ SUBROUTINE moments_eq_rhs_e
Tphi = 0._dp
ENDIF
- !! Collision
- ! IF (CO .EQ. 0) THEN ! Lenard Bernstein
- ! CALL LenardBernstein_e(ip,ij,ikx,iky,iz,TColl)
- ! ELSEIF (CO .EQ. 1) THEN ! GK Dougherty
- ! CALL DoughertyGK_e(ip,ij,ikx,iky,iz,TColl)
- ! ELSE ! COSOLver matrix
- ! TColl = TColl_e(ip,ij,ikx,iky,iz)
- ! ENDIF
-
!! Sum of all linear terms (the sign is inverted to match RHS)
moments_rhs_e(ip,ij,ikx,iky,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo)* (Tnepj + Tnepp2j + Tnepm2j + Tnepjp1 + Tnepjm1)&
! Parallel coupling (Landau Damping)
- - Tpare*inv_deltaz*gradz_coeff(iz,eo) &
+ - Tpar*gradz_coeff(iz,eo) &
! Mirror term (parallel magnetic gradient)
- gradzB(iz,eo)* Tmir *gradz_coeff(iz,eo) &
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Electrostatic background gradients
- i_ky * K_E * moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
- ! Numerical hyperdiffusion (totally artificial, for stability purpose)
+ ! Numerical perpendicular hyperdiffusion (totally artificial, for stability purpose)
- (mu_x*kx**4 + mu_y*ky**4)*moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
+ ! Numerical parallel hyperdiffusion "- (mu_z*kz**4)"
+ - mu_z * ddz4Nepj(iz) &
! Collision term
+ TColl_e(ip,ij,ikx,iky,iz)
@@ -156,19 +152,19 @@ SUBROUTINE moments_eq_rhs_i
USE model
USE prec_const
USE collision
- USE calculus, ONLY : interp_z, grad_z
+ USE calculus, ONLY : interp_z, grad_z, grad_z4
IMPLICIT NONE
INTEGER :: p_int, j_int ! loops indices and polynom. degrees
REAL(dp) :: kx, ky, kperp2
- COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1, Tpari, Tphi
- COMPLEX(dp) :: Tmir, Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments
- COMPLEX(dp) :: UNipm1j, UNipm1jp1, UNipm1jm1 ! Terms from mirror force with adiab moments
+ COMPLEX(dp) :: Tnipj, Tnipp2j, Tnipm2j, Tnipjp1, Tnipjm1
+ COMPLEX(dp) :: Tnipp1j, Tnipm1j, Tnipp1jm1, Tnipm1jm1 ! Terms from mirror force with non adiab moments
+ COMPLEX(dp) :: Unipm1j, Unipm1jp1, Unipm1jm1 ! Terms from mirror force with adiab moments
+ COMPLEX(dp) :: Tperp, Tpar, Tmir, Tphi
COMPLEX(dp) :: i_ky
- INTEGER :: izm2, izm1, izp1, izp2 ! indices for centered FDF ddz
! To store derivatives and odd-even z grid interpolations
COMPLEX(dp), DIMENSION(izs:ize) :: ddznipp1j, ddznipm1j, &
- nipp1j, nipp1jm1, nipm1j, nipm1jm1, nipm1jp1
+ nipp1j, nipp1jm1, nipm1j, nipm1jm1, nipm1jp1, ddz4Nipj
! Measuring execution time
CALL cpu_time(t0_rhs)
@@ -186,18 +182,21 @@ SUBROUTINE moments_eq_rhs_i
i_ky = imagu * ky ! toroidal derivative
IF (Nky .EQ. 1) i_ky = imagu * kxarray(ikx) ! If 1D simulation we put kx as ky
! Compute z derivatives and odd-even z interpolations
- CALL grad_z(eo,nadiab_moments_i(ip+1,ij ,ikx,iky,:),ddznipp1j(:))
- CALL grad_z(eo,nadiab_moments_i(ip-1,ij ,ikx,iky,:),ddznipm1j(:))
- CALL interp_z(eo,nadiab_moments_i(ip+1,ij ,ikx,iky,:), nipp1j (:))
- CALL interp_z(eo,nadiab_moments_i(ip+1,ij-1,ikx,iky,:), nipp1jm1(:))
- CALL interp_z(eo,nadiab_moments_i(ip-1,ij ,ikx,iky,:), nipm1j (:))
- CALL interp_z(eo,nadiab_moments_i(ip-1,ij+1,ikx,iky,:), nipm1jp1(:))
- CALL interp_z(eo,nadiab_moments_i(ip-1,ij-1,ikx,iky,:), nipm1jm1(:))
+ CALL grad_z(eo,nadiab_moments_i(ip+1,ij ,ikx,iky,izgs:izge),ddznipp1j(izs:ize))
+ CALL grad_z(eo,nadiab_moments_i(ip-1,ij ,ikx,iky,izgs:izge),ddznipm1j(izs:ize))
+ CALL interp_z(eo,nadiab_moments_i(ip+1,ij ,ikx,iky,izgs:izge), nipp1j (izs:ize))
+ CALL interp_z(eo,nadiab_moments_i(ip+1,ij-1,ikx,iky,izgs:izge), nipp1jm1(izs:ize))
+ CALL interp_z(eo,nadiab_moments_i(ip-1,ij ,ikx,iky,izgs:izge), nipm1j (izs:ize))
+ CALL interp_z(eo,nadiab_moments_i(ip-1,ij+1,ikx,iky,izgs:izge), nipm1jp1(izs:ize))
+ CALL interp_z(eo,nadiab_moments_i(ip-1,ij-1,ikx,iky,izgs:izge), nipm1jm1(izs:ize))
+ ! for hyperdiffusion
+ CALL grad_z4(moments_i(ip,ij,ikx,iky,:,updatetlevel), ddz4Nipj(:))
+
zloopi : DO iz = izs,ize
- ! kperp2= gxx(iz)*kx**2 + 2._dp*gxy(iz)*kx*ky + gyy(iz)*ky**2
kperp2= kparray(ikx,iky,iz,eo)**2
!! Compute moments mixing terms
+ Tperp = 0._dp; Tpar = 0._dp; Tmir = 0._dp
! Perpendicular dynamic
! term propto n_i^{p,j}
Tnipj = xnipj(ip,ij) * nadiab_moments_i(ip ,ij ,ikx,iky,iz)
@@ -205,27 +204,27 @@ SUBROUTINE moments_eq_rhs_i
Tnipp2j = xnipp2j(ip) * nadiab_moments_i(ip+pp2,ij ,ikx,iky,iz)
! term propto n_i^{p-2,j}
Tnipm2j = xnipm2j(ip) * nadiab_moments_i(ip-pp2,ij ,ikx,iky,iz)
- ! term propto n_e^{p,j+1}
+ ! term propto n_i^{p,j+1}
Tnipjp1 = xnipjp1(ij) * nadiab_moments_i(ip ,ij+1,ikx,iky,iz)
- ! term propto n_e^{p,j-1}
+ ! term propto n_i^{p,j-1}
Tnipjm1 = xnipjm1(ij) * nadiab_moments_i(ip ,ij-1,ikx,iky,iz)
+ ! Tperp
+ Tperp = Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1
! Parallel dynamic
- Tpari = 0._dp; Tmir = 0._dp
IF(Nz .GT. 1) THEN
- ! term propto N_i^{p,j+1}, centered FDF
- Tpari = xnipp1j(ip) * ddznipp1j(iz) + xnipm1j(ip) * ddznipm1j(iz)
-
- ! Mirror terms
- Tnipp1j = ynipp1j (ip,ij) * nipp1j (iz)
- Tnipp1jm1 = ynipp1jm1(ip,ij) * nipp1jm1(iz)
- Tnipm1j = ynipm1j (ip,ij) * nipm1j (iz)
- Tnipm1jm1 = ynipm1jm1(ip,ij) * nipm1jm1(iz)
- ! Trapping terms
- UNipm1j = znipm1j (ip,ij) * nipm1j (iz)
- UNipm1jp1 = znipm1jp1(ip,ij) * nipm1jp1(iz)
- UNipm1jm1 = znipm1jm1(ip,ij) * nipm1jm1(iz)
+ ! ddz derivative for Landau damping term
+ Tpar = xnipp1j(ip) * ddznipp1j(iz) + xnipm1j(ip) * ddznipm1j(iz)
+ ! Mirror terms
+ Tnipp1j = ynipp1j (ip,ij) * nipp1j (iz)
+ Tnipp1jm1 = ynipp1jm1(ip,ij) * nipp1jm1(iz)
+ Tnipm1j = ynipm1j (ip,ij) * nipm1j (iz)
+ Tnipm1jm1 = ynipm1jm1(ip,ij) * nipm1jm1(iz)
+ ! Trapping terms
+ Unipm1j = znipm1j (ip,ij) * nipm1j (iz)
+ Unipm1jp1 = znipm1jp1(ip,ij) * nipm1jp1(iz)
+ Unipm1jm1 = znipm1jm1(ip,ij) * nipm1jm1(iz)
- Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + UNipm1j + UNipm1jp1 + UNipm1jm1
+ Tmir = Tnipp1j + Tnipp1jm1 + Tnipm1j + Tnipm1jm1 + Unipm1j + Unipm1jp1 + Unipm1jm1
ENDIF
!! Electrical potential term
IF ( p_int .LE. 2 ) THEN ! kronecker p0 p1 p2
@@ -236,29 +235,22 @@ SUBROUTINE moments_eq_rhs_i
Tphi = 0._dp
ENDIF
- !! Collision
- ! IF (CO .EQ. 0) THEN ! Lenard Bernstein
- ! CALL LenardBernstein_i(ip,ij,ikx,iky,iz,TColl)
- ! ELSEIF (CO .EQ. 1) THEN ! GK Dougherty
- ! CALL DoughertyGK_i(ip,ij,ikx,iky,iz,TColl)
- ! ELSE! COSOLver matrix (Sugama, Coulomb)
- ! TColl = TColl_i(ip,ij,ikx,iky,iz)
- ! ENDIF
-
!! Sum of all linear terms (the sign is inverted to match RHS)
moments_rhs_i(ip,ij,ikx,iky,iz,updatetlevel) = &
! Perpendicular magnetic gradient/curvature effects
- - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo)*(Tnipj + Tnipp2j + Tnipm2j + Tnipjp1 + Tnipjm1)&
+ - imagu*Ckxky(ikx,iky,iz,eo)*hatR(iz,eo) * Tperp &
! Parallel coupling (Landau Damping)
- - Tpari*inv_deltaz*gradz_coeff(iz,eo) &
+ - gradz_coeff(iz,eo) * Tpar &
! Mirror term (parallel magnetic gradient)
- - gradzB(iz,eo)*Tmir*gradz_coeff(iz,eo) &
+ - gradzB(iz,eo) * gradz_coeff(iz,eo) * Tmir &
! Drives (density + temperature gradients)
- i_ky * Tphi &
! Electrostatic background gradients
- i_ky * K_E * moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
! Numerical hyperdiffusion (totally artificial, for stability purpose)
- (mu_x*kx**4 + mu_y*ky**4)*moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
+ ! Numerical parallel hyperdiffusion "- (mu_z*kz**4)"
+ - mu_z * ddz4Nipj(iz) &
! Collision term
+ TColl_i(ip,ij,ikx,iky,iz)
@@ -276,3 +268,5 @@ SUBROUTINE moments_eq_rhs_i
tc_rhs = tc_rhs + (t1_rhs-t0_rhs)
END SUBROUTINE moments_eq_rhs_i
+
+END MODULE moments_eq_rhs
diff --git a/src/compute_Sapj.F90 b/src/nonlinear_mod.F90
similarity index 96%
rename from src/compute_Sapj.F90
rename to src/nonlinear_mod.F90
index ecf8df56..b327b44b 100644
--- a/src/compute_Sapj.F90
+++ b/src/nonlinear_mod.F90
@@ -1,8 +1,4 @@
-SUBROUTINE compute_Sapj
- ! This routine is meant to compute the non linear term for each specie and degree
- !! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
- !! Sapj = Sum_n (ikx Kn phi)#(iky Sum_s d_njs Naps) - (iky Kn phi)#(ikx Sum_s d_njs Naps)
- !! where # denotes the convolution.
+MODULE nonlinear
USE array, ONLY : dnjs, Sepj, Sipj, kernel_i, kernel_e,&
moments_e_ZF, moments_i_ZF, phi_ZF
USE initial_par, ONLY : ACT_ON_MODES
@@ -16,15 +12,34 @@ SUBROUTINE compute_Sapj
IMPLICIT NONE
INCLUDE 'fftw3-mpi.f03'
- COMPLEX(dp), DIMENSION(ikxs:ikxe,ikys:ikye) :: F_cmpx, G_cmpx
- COMPLEX(dp), DIMENSION(ikxs:ikxe,ikys:ikye) :: Fx_cmpx, Gy_cmpx
- COMPLEX(dp), DIMENSION(ikxs:ikxe,ikys:ikye) :: Fy_cmpx, Gx_cmpx, F_conv_G
- REAL(dp), DIMENSION(ixs:ixe,iys:iye) :: fr_real, gz_real
- REAL(dp), DIMENSION(ixs:ixe,iys:iye) :: fz_real, gr_real, f_times_g
+ COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: F_cmpx, G_cmpx
+ COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: Fx_cmpx, Gy_cmpx
+ COMPLEX(dp), DIMENSION(:,:), ALLOCATABLE :: Fy_cmpx, Gx_cmpx, F_conv_G
INTEGER :: in, is, p_int, j_int
INTEGER :: nmax, smax ! Upper bound of the sums
REAL(dp):: kx, ky, kerneln
+ PUBLIC :: compute_Sapj, nonlinear_init
+
+CONTAINS
+
+SUBROUTINE nonlinear_init
+ ALLOCATE( F_cmpx(ikxs:ikxe,ikys:ikye))
+ ALLOCATE( G_cmpx(ikxs:ikxe,ikys:ikye))
+
+ ALLOCATE(Fx_cmpx(ikxs:ikxe,ikys:ikye))
+ ALLOCATE(Gy_cmpx(ikxs:ikxe,ikys:ikye))
+ ALLOCATE(Fy_cmpx(ikxs:ikxe,ikys:ikye))
+ ALLOCATE(Gx_cmpx(ikxs:ikxe,ikys:ikye))
+
+ ALLOCATE(F_conv_G(ikxs:ikxe,ikys:ikye))
+END SUBROUTINE nonlinear_init
+
+SUBROUTINE compute_Sapj
+ ! This routine is meant to compute the non linear term for each specie and degree
+ !! In real space Sapj ~ b*(grad(phi) x grad(g)) which in moments in fourier becomes
+ !! Sapj = Sum_n (ikx Kn phi)#(iky Sum_s d_njs Naps) - (iky Kn phi)#(ikx Sum_s d_njs Naps)
+ !! where # denotes the convolution.
! Execution time start
CALL cpu_time(t0_Sapj)
@@ -47,7 +62,7 @@ SUBROUTINE compute_Sapj
CALL cpu_time(t1_Sapj)
tc_Sapj = tc_Sapj + (t1_Sapj - t0_Sapj)
-CONTAINS
+END SUBROUTINE compute_Sapj
SUBROUTINE compute_nonlinear
IMPLICIT NONE
@@ -402,4 +417,4 @@ zloopi: DO iz = izs,ize
ENDDO zloopi
END SUBROUTINE compute_semi_linear_NZ
-END SUBROUTINE compute_Sapj
+END MODULE nonlinear
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index a925f3fb..4d3f5af8 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -60,23 +60,27 @@ SUBROUTINE evaluate_kernels
DO eo = 0,1
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
-DO iz = izs,ize
+DO iz = izgs,izge
!!!!! Electron kernels !!!!!
IF(KIN_E) THEN
- DO ij = ijsg_e, ijeg_e
+ DO ij = ijgs_e, ijge_e
j_int = jarray_e(ij)
j_dp = REAL(j_int,dp)
y_ = sigmae2_taue_o2 * kparray(ikx,iky,iz,eo)**2
kernel_e(ij,ikx,iky,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
+ IF (ijs_e .EQ. 1) &
+ kernel_e(ijgs_e,ikx,iky,iz,eo) = 0._dp
ENDIF
!!!!! Ion kernels !!!!!
- DO ij = ijsg_i, ijeg_i
+ DO ij = ijgs_i, ijge_i
j_int = jarray_i(ij)
j_dp = REAL(j_int,dp)
y_ = sigmai2_taui_o2 * kparray(ikx,iky,iz,eo)**2
kernel_i(ij,ikx,iky,iz,eo) = y_**j_int*EXP(-y_)/GAMMA(j_dp+1._dp)!factj
ENDDO
+ IF (ijs_i .EQ. 1) &
+ kernel_i(ijgs_i,ikx,iky,iz,eo) = 0._dp
ENDDO
ENDDO
ENDDO
@@ -99,7 +103,7 @@ SUBROUTINE evaluate_poisson_op
kxloop: DO ikx = ikxs,ikxe
kyloop: DO iky = ikys,ikye
- zloop: DO iz = izs,ize
+ zloop: DO iz = izgs,izge
! This term is evalued on the even z grid since poisson uses only p=0 and phi
IF( (kxarray(ikx).EQ.0._dp) .AND. (kyarray(iky).EQ.0._dp) ) THEN
inv_poisson_op(ikx, iky, iz) = 0._dp
diff --git a/src/poisson.F90 b/src/poisson.F90
index 95b5342a..18eac749 100644
--- a/src/poisson.F90
+++ b/src/poisson.F90
@@ -6,6 +6,7 @@ SUBROUTINE poisson
USE array
USE fields
USE grid
+ USE ghosts, ONLY: update_ghosts_z_phi
USE calculus, ONLY : simpson_rule_z
USE utility, ONLY : manual_3D_bcast
use model, ONLY : qe2_taue, qi2_taui, q_e, q_i, lambdaD, KIN_E
@@ -18,8 +19,7 @@ SUBROUTINE poisson
COMPLEX(dp) :: fa_phi, intf_ ! current flux averaged phi
INTEGER :: count !! mpi integer to broadcast the electric potential at the end
COMPLEX(dp) :: buffer(ikxs:ikxe,ikys:ikye)
- COMPLEX(dp), DIMENSION(izs:ize) :: rho_i, rho_e ! charge density q_a n_a
- COMPLEX(dp), DIMENSION(izs:ize) :: integrant ! charge density q_a n_a
+ COMPLEX(dp), DIMENSION(izgs:izge) :: rho_i, rho_e ! charge density q_a n_a
!! Poisson can be solved only for process containing ip=1
IF ( (ips_e .EQ. ip0_e) .AND. (ips_i .EQ. ip0_i) ) THEN
@@ -31,38 +31,38 @@ SUBROUTINE poisson
!!!! Compute ion gyro charge density
rho_i = 0._dp
DO ini=1,jmaxi+1
- rho_i = rho_i &
- +q_i*kernel_i(ini,ikx,iky,:,0)*moments_i(ip0_i,ini,ikx,iky,:,updatetlevel)
+ rho_i(izs:ize) = rho_i(izs:ize) &
+ +q_i*kernel_i(ini,ikx,iky,izs:ize,0)*moments_i(ip0_i,ini,ikx,iky,izs:ize,updatetlevel)
END DO
!!!! Compute electron gyro charge density
rho_e = 0._dp
IF (KIN_E) THEN ! Kinetic electrons
DO ine=1,jmaxe+1
- rho_e = rho_e &
- +q_e*kernel_e(ine,ikx,iky,:,0)*moments_e(ip0_e,ine, ikx,iky,:,updatetlevel)
+ rho_e(izs:ize) = rho_e(izs:ize) &
+ +q_e*kernel_e(ine,ikx,iky,izs:ize,0)*moments_e(ip0_e,ine, ikx,iky,izs:ize,updatetlevel)
END DO
ELSE ! Adiabatic electrons
! Adiabatic charge density (linked to flux averaged phi)
fa_phi = 0._dp
- IF(kyarray(iky).EQ.0._dp) THEN
- DO ini=1,jmaxi+1
- rho_e(:) = Jacobian(:,0)*moments_i(ip0_i,ini,ikx,iky,:,updatetlevel)&
- *kernel_i(ini,ikx,iky,:,0)*(inv_poisson_op(ikx,iky,:)-1._dp)
- call simpson_rule_z(rho_e,intf_)
+ IF(kyarray(iky).EQ.0._dp) THEN ! take ky=0 mode (average)
+ DO ini=1,jmaxi+1 ! some over FLR contributions
+ rho_e(izgs:izge) = Jacobian(izgs:izge,0)*moments_i(ip0_i,ini,ikx,iky,izgs:izge,updatetlevel)&
+ *kernel_i(ini,ikx,iky,izgs:izge,0)*(inv_poisson_op(ikx,iky,izgs:izge)-1._dp)
+ call simpson_rule_z(rho_e,intf_) ! integrate over z
fa_phi = fa_phi + intf_
ENDDO
- rho_e = fa_phi*iInt_Jacobian !Normalize by 1/int(Jxyz)dz
+ rho_e(izs:ize) = fa_phi*iInt_Jacobian !Normalize by 1/int(Jxyz)dz
ENDIF
ENDIF
!!!!!!!!!!!!!!! Inverting the poisson equation !!!!!!!!!!!!!!!!!!!!!!!!!!
- phi(ikx, iky, :) = (rho_e + rho_i)*inv_poisson_op(ikx,iky,:)
+ phi(ikx, iky, izs:ize) = (rho_e(izs:ize) + rho_i(izs:ize))*inv_poisson_op(ikx,iky,izs:ize)
END DO kyloop
END DO kxloop
! Cancel origin singularity
- IF (contains_kx0 .AND. contains_ky0) phi(ikx_0,iky_0,:) = 0._dp
+ IF (contains_kx0 .AND. contains_ky0) phi(ikx_0,iky_0,izs:ize) = 0._dp
ENDIF
diff --git a/src/ppinit.F90 b/src/ppinit.F90
index c9418310..a51d666a 100644
--- a/src/ppinit.F90
+++ b/src/ppinit.F90
@@ -7,7 +7,7 @@ SUBROUTINE ppinit
INTEGER :: version_prov=-1
! Variables for cartesian domain decomposition
- INTEGER, PARAMETER :: ndims=2 ! p and kx
+ INTEGER, PARAMETER :: ndims=3 ! p, kx and z
INTEGER, DIMENSION(ndims) :: dims=0, coords=0, coords_L=0, coords_R=0
LOGICAL :: periods(ndims) = .FALSE., reorder=.FALSE.
CHARACTER(len=32) :: str
@@ -33,16 +33,20 @@ SUBROUTINE ppinit
! CALL MPI_DIMS_CREATE(num_procs, ndims, dims, ierr)
dims(1) = 1
dims(2) = num_procs
+ dims(3) = 1
END IF
- num_procs_p = dims(1) ! Number of processes along p
+ num_procs_p = dims(1) ! Number of processes along p
num_procs_kx = dims(2) ! Number of processes along kx
+ num_procs_z = dims(3) ! Number of processes along z
!
!periodicity in p
periods(1)=.FALSE.
!periodicity in kx
periods(2)=.FALSE.
+ !periodicity in z
+ periods(3)=.TRUE.
CALL MPI_CART_CREATE(MPI_COMM_WORLD, ndims, dims, periods, reorder, comm0, ierr)
@@ -53,13 +57,16 @@ SUBROUTINE ppinit
!
! Partitions 2-dim cartesian topology of comm0 into 1-dim cartesian subgrids
!
- CALL MPI_CART_SUB (comm0, (/.TRUE.,.FALSE./), comm_p, ierr)
- CALL MPI_CART_SUB (comm0, (/.FALSE.,.TRUE./), comm_kx, ierr)
+ CALL MPI_CART_SUB (comm0, (/.TRUE.,.FALSE.,.FALSE./), comm_p, ierr)
+ CALL MPI_CART_SUB (comm0, (/.FALSE.,.TRUE.,.FALSE./), comm_kx, ierr)
+ CALL MPI_CART_SUB (comm0, (/.FALSE.,.FALSE.,.TRUE./), comm_z, ierr)
! Find id inside the sub communicators
- CALL MPI_COMM_RANK(comm_p, rank_p, ierr)
+ CALL MPI_COMM_RANK(comm_p, rank_p, ierr)
CALL MPI_COMM_RANK(comm_kx, rank_kx, ierr)
+ CALL MPI_COMM_RANK(comm_z, rank_z, ierr)
! Find neighbours
CALL MPI_CART_SHIFT(comm0, 0, 1, nbr_L, nbr_R, ierr)
CALL MPI_CART_SHIFT(comm0, 1, 1, nbr_B, nbr_T, ierr)
+ CALL MPI_CART_SHIFT(comm0, 2, 1, nbr_D, nbr_U, ierr)
END SUBROUTINE ppinit
diff --git a/src/prec_const_mod.F90 b/src/prec_const_mod.F90
index da75b190..90388443 100644
--- a/src/prec_const_mod.F90
+++ b/src/prec_const_mod.F90
@@ -26,6 +26,7 @@ MODULE prec_const
REAL(dp), PARAMETER :: PI=3.141592653589793238462643383279502884197_dp
REAL(dp), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_dp
REAL(dp), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_dp
+ REAL(dp), PARAMETER :: ONEOPI=0.3183098861837906912164442019275156781077_dp
REAL(dp), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_dp
REAL(dp), PARAMETER :: INVSQRT2=0.7071067811865475244008443621048490392848359377_dp
REAL(dp), PARAMETER :: SQRT3=1.73205080756887729352744634150587236694281_dp
diff --git a/src/processing_mod.F90 b/src/processing_mod.F90
index 5419e3b7..06430661 100644
--- a/src/processing_mod.F90
+++ b/src/processing_mod.F90
@@ -3,111 +3,138 @@ MODULE processing
USE basic
USE prec_const
USE grid
+ USE geometry
USE utility
implicit none
REAL(dp), PUBLIC, PROTECTED :: pflux_ri, gflux_ri
REAL(dp), PUBLIC, PROTECTED :: hflux_x
- PUBLIC :: compute_density, compute_temperature, compute_nadiab_moments
+ PUBLIC :: compute_nadiab_moments
+ PUBLIC :: compute_density, compute_upar, compute_uperp
+ PUBLIC :: compute_Tpar, compute_Tperp, compute_fluid_moments
PUBLIC :: compute_radial_ion_transport, compute_radial_heatflux
CONTAINS
! 1D diagnostic to compute the average radial particle transport <n_i v_ExB>_r
SUBROUTINE compute_radial_ion_transport
USE fields, ONLY : moments_i, phi
- USE array, ONLY : kernel_i
+ USE array, ONLY : kernel_i, Jacobian
USE time_integration, ONLY : updatetlevel
IMPLICIT NONE
- COMPLEX(dp) :: pflux_local, gflux_local
+ COMPLEX(dp) :: pflux_local, gflux_local, integral
REAL(dp) :: ky_, buffer(1:2)
INTEGER :: i_, world_rank, world_size, root
+ COMPLEX(dp), DIMENSION(izgs:izge) :: integrant
pflux_local = 0._dp ! particle flux
gflux_local = 0._dp ! gyrocenter flux
+ integrant = 0._dp ! auxiliary variable for z integration
IF(ips_i .EQ. 1) THEN
- ! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Ni00 * phi
- DO iz = izs,ize
- DO iky = ikys,ikye
- ky_ = kyarray(iky)
- DO ikx = ikxs,ikxe
- gflux_local = gflux_local + &
- imagu * ky_ * moments_i(ip0_i,1,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
- DO ij = ijs_i, ije_i
- pflux_local = pflux_local + &
- imagu * ky_ * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
- ENDDO
- ENDDO
- ENDDO
+ !! Gyro center flux (drift kinetic)
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ DO iz = izgs,izge
+ integrant(iz) = Jacobian(iz,0)*imagu * ky_ * moments_i(ip0_i,1,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
+ ENDDO
+ ! Integrate over z
+ call simpson_rule_z(integrant,integral)
+ ! add contribution
+ gflux_local = gflux_local + integral*iInt_Jacobian
ENDDO
- gflux_local = gflux_local/Nz ! Average over parallel planes
- pflux_local = pflux_local/Nz
-
- buffer(1) = REAL(gflux_local)
- buffer(2) = REAL(pflux_local)
- root = 0
- !Gather manually among the rank_p=0 processes and perform the sum
- gflux_ri = 0
- pflux_ri = 0
- IF (num_procs_kx .GT. 1) THEN
- !! Everyone sends its local_sum to root = 0
- IF (rank_kx .NE. root) THEN
- CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_kx, ierr)
- ELSE
- ! Recieve from all the other processes
- DO i_ = 0,num_procs_kx-1
- IF (i_ .NE. rank_kx) &
- CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_kx, MPI_STATUS_IGNORE, ierr)
- gflux_ri = gflux_ri + buffer(1)
- pflux_ri = pflux_ri + buffer(2)
+ ENDDO
+ !! Particle flux
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ integrant = 0._dp ! auxiliary variable for z integration
+ DO ij = ijs_i, ije_i
+ DO iz = izgs,izge
+ integrant(iz) = integrant(iz) + &
+ Jacobian(iz,0)*imagu * ky_ * kernel_i(ij,ikx,iky,iz,0) &
+ *moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel) * CONJG(phi(ikx,iky,iz))
ENDDO
- ENDIF
+ ENDDO
+ ! Integrate over z
+ call simpson_rule_z(integrant,integral)
+ ! add contribution
+ pflux_local = pflux_local + integral*iInt_Jacobian
+ ENDDO
+ ENDDO
+ ENDIF
+ buffer(1) = REAL(gflux_local)
+ buffer(2) = REAL(pflux_local)
+ root = 0
+ !Gather manually among the rank_p=0 processes and perform the sum
+ gflux_ri = 0
+ pflux_ri = 0
+ IF (num_procs_kx .GT. 1) THEN
+ !! Everyone sends its local_sum to root = 0
+ IF (rank_kx .NE. root) THEN
+ CALL MPI_SEND(buffer, 2 , MPI_DOUBLE_PRECISION, root, 1234, comm_kx, ierr)
+ ELSE
+ ! Recieve from all the other processes
+ DO i_ = 0,num_procs_kx-1
+ IF (i_ .NE. rank_kx) &
+ CALL MPI_RECV(buffer, 2 , MPI_DOUBLE_PRECISION, i_, 1234, comm_kx, MPI_STATUS_IGNORE, ierr)
+ gflux_ri = gflux_ri + buffer(1)
+ pflux_ri = pflux_ri + buffer(2)
+ ENDDO
ENDIF
+ ELSE
+ gflux_ri = gflux_local
+ pflux_ri = pflux_local
ENDIF
+ ! if(my_id .eq. 0) write(*,*) 'pflux_ri = ',pflux_ri
END SUBROUTINE compute_radial_ion_transport
! 1D diagnostic to compute the average radial particle transport <n_i v_ExB>_r
SUBROUTINE compute_radial_heatflux
USE fields, ONLY : moments_i, moments_e, phi
- USE array, ONLY : dens_e, dens_i, kernel_e, kernel_i
+ USE array, ONLY : kernel_e, kernel_i, Jacobian
USE time_integration, ONLY : updatetlevel
- USE model, ONLY : q_e, q_i, tau_e, tau_i, KIN_E
+ USE model, ONLY : qe_taue, qi_taui, KIN_E
IMPLICIT NONE
- COMPLEX(dp) :: hflux_local
+ COMPLEX(dp) :: hflux_local, integral
REAL(dp) :: ky_, buffer(1:2), j_dp
INTEGER :: i_, world_rank, world_size, root
+ COMPLEX(dp), DIMENSION(izgs:izge) :: integrant ! charge density q_a n_a
hflux_local = 0._dp ! particle flux
IF(ips_i .EQ. 1 .AND. ips_e .EQ. 1) THEN
! Loop to compute gamma_kx = sum_ky sum_j -i k_z Kernel_j Ni00 * phi
- DO iz = izs,ize
- DO iky = ikys,ikye
- ky_ = kyarray(iky)
- DO ikx = ikxs,ikxe
- DO ij = ijs_i, ije_i
- j_dp = REAL(ij-1,dp)
- hflux_local = hflux_local - imagu*ky_*CONJG(phi(ikx,iky,iz))&
- *(2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky,iz,0) - (j_dp+1)*kernel_i(ij+1,ikx,iky,iz,0) - j_dp*kernel_i(ij-1,ikx,iky,iz,0))&
- * (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
- + SQRT2/3._dp * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip2_i,ij,ikx,iky,iz,updatetlevel))
+ DO iky = ikys,ikye
+ ky_ = kyarray(iky)
+ DO ikx = ikxs,ikxe
+ integrant = 0._dp
+ DO ij = ijs_i, ije_i
+ DO iz = izgs,izge
+ integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(ikx,iky,iz))&
+ *(twothird * ( 2._dp*j_dp * kernel_i(ij ,ikx,iky,iz,0) &
+ - (j_dp+1) * kernel_i(ij+1,ikx,iky,iz,0) &
+ - j_dp * kernel_i(ij-1,ikx,iky,iz,0))&
+ * (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+qi_taui*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
+ + SQRT2*onethird * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip2_i,ij,ikx,iky,iz,updatetlevel))
ENDDO
- IF(KIN_E) THEN
- DO ij = ijs_e, ije_e
- j_dp = REAL(ij-1,dp)
- hflux_local = hflux_local - imagu*ky_*CONJG(phi(ikx,iky,iz))&
- *(2._dp/3._dp *(2._dp*j_dp * kernel_e(ij ,ikx,iky,iz,0) &
- -(j_dp+1) * kernel_e(ij+1,ikx,iky,iz,0) &
- -j_dp * kernel_e(ij-1,ikx,iky,iz,0))&
- *(moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)&
- +q_e/tau_e * kernel_e(ij ,ikx,iky,iz,0) * phi(ikx,iky,iz)) &
- +SQRT2/3._dp * kernel_e(ij,ikx,iky,iz,0) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel))
+ ENDDO
+ IF(KIN_E) THEN
+ DO ij = ijs_e, ije_e
+ DO iz = izgs,izge
+ integrant(iz) = integrant(iz) + Jacobian(iz,0)*imagu*ky_*CONJG(phi(ikx,iky,iz))&
+ *(twothird * ( 2._dp*j_dp * kernel_e(ij ,ikx,iky,iz,0) &
+ - (j_dp+1) * kernel_e(ij+1,ikx,iky,iz,0) &
+ - j_dp * kernel_e(ij-1,ikx,iky,iz,0))&
+ * (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+qe_taue*kernel_e(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
+ + SQRT2*onethird * kernel_e(ij,ikx,iky,iz,0) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel))
ENDDO
- ENDIF
ENDDO
+ ENDIF
+ ! Integrate over z
+ call simpson_rule_z(integrant,integral)
+ hflux_local = hflux_local + integral*iInt_Jacobian
ENDDO
ENDDO
- hflux_local = hflux_local/Nz
-
buffer(2) = REAL(hflux_local)
root = 0
!Gather manually among the rank_p=0 processes and perform the sum
@@ -124,6 +151,8 @@ SUBROUTINE compute_radial_heatflux
hflux_x = hflux_x + buffer(2)
ENDDO
ENDIF
+ ELSE
+ hflux_x = hflux_local
ENDIF
ENDIF
END SUBROUTINE compute_radial_heatflux
@@ -139,115 +168,302 @@ SUBROUTINE compute_nadiab_moments
USE model, ONLY : qe_taue, qi_taui, KIN_E
implicit none
- ! Add non-adiabatique term
+ ! Electron non adiab moments
+ DO ikx = ikxs,ikxe
+ DO iky = ikys,ikye
+ DO iz = izgs,izge
IF(KIN_E) THEN
- DO ip=ipsg_e,ipeg_e
+ DO ip=ipgs_e,ipge_e
IF(parray_e(ip) .EQ. 0) THEN
- DO ij=ijsg_e,ijeg_e
- nadiab_moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize)&
- = moments_e(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) &
- + qe_taue*kernel_e(ij,ikxs:ikxe,ikys:ikye,izs:ize,0)*phi(ikxs:ikxe,ikys:ikye,izs:ize)
+ DO ij=ijgs_e,ijge_e
+ nadiab_moments_e(ip,ij,ikx,iky,iz) = moments_e(ip,ij,ikx,iky,iz,updatetlevel) &
+ + qe_taue*kernel_e(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)
ENDDO
ELSE
- nadiab_moments_e(ip,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize) &
- = moments_e(ip,ijsg_e:ijeg_e,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
+ DO ij=ijgs_e,ijge_e
+ nadiab_moments_e(ip,ij,ikx,iky,iz) = moments_e(ip,ij,ikx,iky,iz,updatetlevel)
+ ENDDO
ENDIF
ENDDO
ENDIF
- ! Add non-adiabatique term
- DO ip=ipsg_i,ipeg_i
+ ! Ions non adiab moments
+ DO ip=ipgs_i,ipge_i
IF(parray_i(ip) .EQ. 0) THEN
- DO ij=ijsg_i,ijeg_i
- nadiab_moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize)&
- = moments_i(ip,ij,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel) &
- + qi_taui*kernel_i(ij,ikxs:ikxe,ikys:ikye,izs:ize,0)*phi(ikxs:ikxe,ikys:ikye,izs:ize)
+ DO ij=ijgs_i,ijge_i
+ nadiab_moments_i(ip,ij,ikx,iky,iz) = moments_i(ip,ij,ikx,iky,iz,updatetlevel) &
+ + qi_taui*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)
ENDDO
ELSE
- nadiab_moments_i(ip,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize) &
- = moments_i(ip,ijsg_i:ijeg_i,ikxs:ikxe,ikys:ikye,izs:ize,updatetlevel)
+ DO ij=ijgs_i,ijge_i
+ nadiab_moments_i(ip,ij,ikx,iky,iz) = moments_i(ip,ij,ikx,iky,iz,updatetlevel)
+ ENDDO
ENDIF
ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
!
END SUBROUTINE compute_nadiab_moments
! Compute the 2D particle density for electron and ions (sum over Laguerre)
SUBROUTINE compute_density
- USE fields, ONLY : moments_i, moments_e, phi
- USE array, ONLY : dens_e, dens_i, kernel_e, kernel_i
- USE time_integration, ONLY : updatetlevel
- USE model, ONLY : q_e, q_i, tau_e, tau_i, KIN_E
+ USE array, ONLY : dens_e, dens_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i
+ USE model, ONLY : KIN_E
IMPLICIT NONE
+ COMPLEX(dp) :: dens
- IF ( (ips_e .EQ. ip0_e) .AND. (ips_i .EQ. ip0_i) ) THEN
+ IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i ) THEN
! Loop to compute dens_i = sum_j kernel_j Ni0j at each k
DO iky = ikys,ikye
DO ikx = ikxs,ikxe
DO iz = izs,ize
IF(KIN_E) THEN
! electron density
- dens_e(ikx,iky,iz) = 0._dp
+ dens = 0._dp
DO ij = ijs_e, ije_e
- dens_e(ikx,iky,iz) = dens_e(ikx,iky,iz) + kernel_e(ij,ikx,iky,iz,0) * &
- (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz,0)*phi(ikx,iky,iz))
+ dens = dens + kernel_e(ij,ikx,iky,iz,0) * nadiab_moments_e(ip0_e,ij,ikx,iky,iz)
ENDDO
+ dens_e(ikx,iky,iz) = dens
ENDIF
! ion density
- dens_i(ikx,iky,iz) = 0._dp
+ dens = 0._dp
DO ij = ijs_i, ije_i
- dens_i(ikx,iky,iz) = dens_i(ikx,iky,iz) + kernel_i(ij,ikx,iky,iz,0) * &
- (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz))
+ dens = dens + kernel_i(ij,ikx,iky,iz,0) * nadiab_moments_i(ip0_e,ij,ikx,iky,iz)
ENDDO
+ dens_i(ikx,iky,iz) = dens
ENDDO
ENDDO
ENDDO
ENDIF
- IF(KIN_E)&
- CALL manual_3D_bcast(dens_e(ikxs:ikxe,ikys:ikye,izs:ize))
- CALL manual_3D_bcast(dens_i(ikxs:ikxe,ikys:ikye,izs:ize))
+ ! IF(KIN_E)&
+ ! CALL manual_3D_bcast(dens_e(ikxs:ikxe,ikys:ikye,izs:ize))
+ ! CALL manual_3D_bcast(dens_i(ikxs:ikxe,ikys:ikye,izs:ize))
END SUBROUTINE compute_density
+! Compute the 2D particle fluid perp velocity for electron and ions (sum over Laguerre)
+SUBROUTINE compute_uperp
+ USE array, ONLY : uper_e, uper_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i
+ USE model, ONLY : KIN_E
+ IMPLICIT NONE
+ COMPLEX(dp) :: uperp
+
+ IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i ) THEN
+ DO iky = ikys,ikye
+ DO ikx = ikxs,ikxe
+ DO iz = izs,ize
+ IF(KIN_E) THEN
+ ! electron
+ uperp = 0._dp
+ DO ij = ijs_e, ije_e
+ uperp = uperp + kernel_e(ij,ikx,iky,iz,0) *&
+ 0.5_dp*(nadiab_moments_e(ip0_e,ij,ikx,iky,iz) - nadiab_moments_e(ip0_e,ij-1,ikx,iky,iz))
+ ENDDO
+ uper_e(ikx,iky,iz) = uperp
+ ENDIF
+ ! ion
+ uperp = 0._dp
+ DO ij = ijs_i, ije_i
+ uperp = uperp + kernel_i(ij,ikx,iky,iz,0) *&
+ 0.5_dp*(nadiab_moments_i(ip0_i,ij,ikx,iky,iz) - nadiab_moments_i(ip0_e,ij-1,ikx,iky,iz))
+ ENDDO
+ uper_i(ikx,iky,iz) = uperp
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! IF(KIN_E)&
+ ! CALL manual_3D_bcast(uper_e(ikxs:ikxe,ikys:ikye,izs:ize))
+ ! CALL manual_3D_bcast(uper_i(ikxs:ikxe,ikys:ikye,izs:ize))
+END SUBROUTINE compute_uperp
+
+! Compute the 2D particle fluid par velocity for electron and ions (sum over Laguerre)
+SUBROUTINE compute_upar
+ USE array, ONLY : upar_e, upar_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i
+ USE model, ONLY : KIN_E
+ IMPLICIT NONE
+ COMPLEX(dp) :: upar
+
+ IF ( CONTAINS_ip1_e .AND. CONTAINS_ip1_i ) THEN
+ DO iky = ikys,ikye
+ DO ikx = ikxs,ikxe
+ DO iz = izs,ize
+ IF(KIN_E) THEN
+ ! electron
+ upar = 0._dp
+ DO ij = ijs_e, ije_e
+ upar = upar + kernel_e(ij,ikx,iky,iz,1)*nadiab_moments_e(ip1_e,ij,ikx,iky,iz)
+ ENDDO
+ upar_e(ikx,iky,iz) = upar
+ ENDIF
+ ! ion
+ upar = 0._dp
+ DO ij = ijs_i, ije_i
+ upar = upar + kernel_i(ij,ikx,iky,iz,1)*nadiab_moments_i(ip1_i,ij,ikx,iky,iz)
+ ENDDO
+ upar_i(ikx,iky,iz) = upar
+ ENDDO
+ ENDDO
+ ENDDO
+ ELSE
+ IF(KIN_E)&
+ upar_e = 0
+ upar_i = 0
+ ENDIF
+ ! IF(KIN_E)&
+ ! CALL manual_3D_bcast(upar_e(ikxs:ikxe,ikys:ikye,izs:ize))
+ ! CALL manual_3D_bcast(upar_i(ikxs:ikxe,ikys:ikye,izs:ize))
+END SUBROUTINE compute_upar
+
! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
-SUBROUTINE compute_temperature
- USE fields, ONLY : moments_i, moments_e, phi
- USE array, ONLY : temp_e, temp_i, kernel_e, kernel_i
- USE time_integration, ONLY : updatetlevel
- USE model, ONLY : q_e, q_i, tau_e, tau_i, KIN_E
+SUBROUTINE compute_tperp
+ USE array, ONLY : Tper_e, Tper_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i
+ USE model, ONLY : KIN_E
+ IMPLICIT NONE
+ REAL(dp) :: j_dp
+ COMPLEX(dp) :: Tperp
+
+ IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i .AND. &
+ CONTAINS_ip2_e .AND. CONTAINS_ip2_i ) THEN
+ ! Loop to compute T = 1/3*(Tpar + 2Tperp)
+ DO iky = ikys,ikye
+ DO ikx = ikxs,ikxe
+ DO iz = izs,ize
+ ! electron temperature
+ IF(KIN_E) THEN
+ Tperp = 0._dp
+ DO ij = ijs_e, ije_e
+ j_dp = REAL(ij-1,dp)
+ Tperp = Tperp + kernel_e(ij,ikx,iky,iz,0)*&
+ ((2_dp*j_dp+1)*nadiab_moments_e(ip0_e,ij ,ikx,iky,iz)&
+ -j_dp *nadiab_moments_e(ip0_e,ij-1,ikx,iky,iz)&
+ -j_dp+1 *nadiab_moments_e(ip0_e,ij+1,ikx,iky,iz))
+ ENDDO
+ Tper_e(ikx,iky,iz) = Tperp
+ ENDIF
+ ! ion temperature
+ Tperp = 0._dp
+ DO ij = ijs_i, ije_i
+ j_dp = REAL(ij-1,dp)
+ Tperp = Tperp + kernel_i(ij,ikx,iky,iz,0)*&
+ ((2_dp*j_dp+1)*nadiab_moments_i(ip0_i,ij ,ikx,iky,iz)&
+ -j_dp *nadiab_moments_i(ip0_i,ij-1,ikx,iky,iz)&
+ -j_dp+1 *nadiab_moments_i(ip0_i,ij+1,ikx,iky,iz))
+ ENDDO
+ Tper_i(ikx,iky,iz) = Tperp
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+ ! IF(KIN_E)&
+ ! CALL manual_3D_bcast(Tper_e(ikxs:ikxe,ikys:ikye,izs:ize))
+ ! CALL manual_3D_bcast(Tper_i(ikxs:ikxe,ikys:ikye,izs:ize))
+END SUBROUTINE compute_Tperp
+
+! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
+SUBROUTINE compute_Tpar
+ USE array, ONLY : Tpar_e, Tpar_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i
+ USE model, ONLY : KIN_E
IMPLICIT NONE
REAL(dp) :: j_dp
- COMPLEX(dp) :: Tperp, Tpar
+ COMPLEX(dp) :: tpar
- IF( ((ips_i .EQ. 1) .AND. (ips_e .EQ. 1)) ) THEN
+ IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i .AND. &
+ CONTAINS_ip2_e .AND. CONTAINS_ip2_i ) THEN
! Loop to compute T = 1/3*(Tpar + 2Tperp)
DO iky = ikys,ikye
DO ikx = ikxs,ikxe
DO iz = izs,ize
! electron temperature
IF(KIN_E) THEN
- temp_e(ikx,iky,iz) = 0._dp
+ Tpar = 0._dp
DO ij = ijs_e, ije_e
j_dp = REAL(ij-1,dp)
- temp_e(ikx,iky,iz) = temp_e(ikx,iky,iz) + &
- 2._dp/3._dp * (2._dp*j_dp*kernel_e(ij,ikx,iky,iz,0) - (j_dp+1)*kernel_e(ij+1,ikx,iky,iz,0) - j_dp*kernel_e(ij-1,ikx,iky,iz,0))&
- * (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
- + SQRT2/3._dp * kernel_e(ij,ikx,iky,iz,0) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel)
+ Tpar = Tpar + kernel_e(ij,ikx,iky,iz,0)*&
+ (SQRT2 * nadiab_moments_e(ip2_e,ij,ikx,iky,iz) &
+ + nadiab_moments_e(ip0_e,ij,ikx,iky,iz))
ENDDO
+ Tpar_e(ikx,iky,iz) = Tpar
ENDIF
! ion temperature
- temp_i(ikx,iky,iz) = 0._dp
+ Tpar = 0._dp
DO ij = ijs_i, ije_i
j_dp = REAL(ij-1,dp)
- temp_i(ikx,iky,iz) = temp_i(ikx,iky,iz) + &
- 2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky,iz,0) - (j_dp+1)*kernel_i(ij+1,ikx,iky,iz,0) - j_dp*kernel_i(ij-1,ikx,iky,iz,0))&
- * (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
- + SQRT2/3._dp * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip2_i,ij,ikx,iky,iz,updatetlevel)
+ Tpar = Tpar + kernel_i(ij,ikx,iky,iz,0)*&
+ (SQRT2 * nadiab_moments_i(ip2_i,ij,ikx,iky,iz) &
+ + nadiab_moments_i(ip0_i,ij,ikx,iky,iz))
ENDDO
+ Tpar_i(ikx,iky,iz) = Tpar
ENDDO
ENDDO
ENDDO
ENDIF
- IF(KIN_E)&
- CALL manual_3D_bcast(temp_e(ikxs:ikxe,ikys:ikye,izs:ize))
- CALL manual_3D_bcast(temp_i(ikxs:ikxe,ikys:ikye,izs:ize))
-END SUBROUTINE compute_temperature
+ ! IF(KIN_E)&
+ ! CALL manual_3D_bcast(Tpar_e(ikxs:ikxe,ikys:ikye,izs:ize))
+ ! CALL manual_3D_bcast(Tpar_i(ikxs:ikxe,ikys:ikye,izs:ize))
+END SUBROUTINE compute_Tpar
+
+! Compute the 2D particle fluid moments for electron and ions (sum over Laguerre)
+SUBROUTINE compute_fluid_moments
+ USE array, ONLY : dens_e, Tpar_e, Tper_e, dens_i, Tpar_i, Tper_i
+ USE model, ONLY : KIN_E
+ IMPLICIT NONE
+ CALL compute_density
+ CALL compute_upar
+ CALL compute_uperp
+ CALL compute_Tpar
+ CALL compute_Tperp
+ ! Temperature
+ temp_e = (Tpar_e - 2._dp * Tper_e)/3._dp - dens_e
+ temp_i = (Tpar_i - 2._dp * Tper_i)/3._dp - dens_i
+END SUBROUTINE compute_fluid_moments
+
+! Compute the 2D particle temperature for electron and ions (sum over Laguerre)
+! SUBROUTINE compute_temperature
+! USE array, ONLY : temp_e, temp_i, kernel_e, kernel_i, nadiab_moments_e, nadiab_moments_i
+! USE model, ONLY : KIN_E
+! IMPLICIT NONE
+! REAL(dp) :: j_dp
+! COMPLEX(dp) :: Tperp, Tpar, dens
+!
+! IF ( CONTAINS_ip0_e .AND. CONTAINS_ip0_i .AND. &
+! CONTAINS_ip2_e .AND. CONTAINS_ip2_i ) THEN
+! ! Loop to compute T = 1/3*(Tpar + 2Tperp)
+! DO iky = ikys,ikye
+! DO ikx = ikxs,ikxe
+! DO iz = izs,ize
+! ! electron temperature
+! IF(KIN_E) THEN
+! dens = 0._dp
+! Tpar = 0._dp
+! Tperp = 0._dp
+! DO ij = ijs_e, ije_e
+! j_dp = REAL(ij-1,dp)
+! temp_e(ikx,iky,iz) = temp_e(ikx,iky,iz) + &
+! 2._dp/3._dp * (2._dp*j_dp*kernel_e(ij,ikx,iky,iz,0) - (j_dp+1)*kernel_e(ij+1,ikx,iky,iz,0) - j_dp*kernel_e(ij-1,ikx,iky,iz,0))&
+! * (moments_e(ip0_e,ij,ikx,iky,iz,updatetlevel)+q_e/tau_e*kernel_e(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
+! + SQRT2/3._dp * kernel_e(ij,ikx,iky,iz,0) * moments_e(ip2_e,ij,ikx,iky,iz,updatetlevel)
+! ENDDO
+! temp_e(ikx,iky,iz) = (Tpar - 2._dp * Tperp)/3._dp - dens
+! ENDIF
+! ! ion temperature
+! dens = 0._dp
+! Tpar = 0._dp
+! Tperp = 0._dp
+! DO ij = ijs_i, ije_i
+! j_dp = REAL(ij-1,dp)
+! temp_i(ikx,iky,iz) = temp_i(ikx,iky,iz) + &
+! 2._dp/3._dp * (2._dp*j_dp*kernel_i(ij,ikx,iky,iz,0) - (j_dp+1)*kernel_i(ij+1,ikx,iky,iz,0) - j_dp*kernel_i(ij-1,ikx,iky,iz,0))&
+! * (moments_i(ip0_i,ij,ikx,iky,iz,updatetlevel)+q_i/tau_i*kernel_i(ij,ikx,iky,iz,0)*phi(ikx,iky,iz)) &
+! + SQRT2/3._dp * kernel_i(ij,ikx,iky,iz,0) * moments_i(ip2_i,ij,ikx,iky,iz,updatetlevel)
+! ENDDO
+! temp_i(ikx,iky,iz) = (Tpar - 2._dp * Tperp)/3._dp - dens
+! ENDDO
+! ENDDO
+! ENDDO
+! ENDIF
+! IF(KIN_E) CALL manual_3D_bcast(temp_e(ikxs:ikxe,ikys:ikye,izs:ize))
+! CALL manual_3D_bcast(temp_i(ikxs:ikxe,ikys:ikye,izs:ize))
+! END SUBROUTINE compute_temperature
+
END MODULE processing
diff --git a/src/srcinfo.h b/src/srcinfo.h
index b1faef28..e4739288 100644
--- a/src/srcinfo.h
+++ b/src/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='65a7746-dirty')
+parameter (VERSION='aa401ea-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Mon Mar 14 16:05:57 CET 2022')
+parameter (EXECDATE='Tue Apr 5 17:06:15 CEST 2022')
parameter (HOST ='spcpc606')
diff --git a/src/srcinfo/srcinfo.h b/src/srcinfo/srcinfo.h
index b1faef28..e4739288 100644
--- a/src/srcinfo/srcinfo.h
+++ b/src/srcinfo/srcinfo.h
@@ -3,8 +3,8 @@ character(len=40) BRANCH
character(len=20) AUTHOR
character(len=40) EXECDATE
character(len=40) HOST
-parameter (VERSION='65a7746-dirty')
+parameter (VERSION='aa401ea-dirty')
parameter (BRANCH='master')
parameter (AUTHOR='ahoffman')
-parameter (EXECDATE='Mon Mar 14 16:05:57 CET 2022')
+parameter (EXECDATE='Tue Apr 5 17:06:15 CEST 2022')
parameter (HOST ='spcpc606')
diff --git a/src/stepon.F90 b/src/stepon.F90
index 79f5dfbb..5aef1d2e 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -4,7 +4,6 @@ SUBROUTINE stepon
USE array , ONLY: moments_rhs_e, moments_rhs_i, Sepj, Sipj
USE basic
USE closure
- USE collision, ONLY : compute_TColl
USE fields, ONLY: moments_e, moments_i, phi
USE initial_par, ONLY: ACT_ON_MODES
USE ghosts
@@ -13,9 +12,8 @@ SUBROUTINE stepon
use prec_const
USE time_integration
USE numerics, ONLY: play_with_modes
- USE processing, ONLY: compute_nadiab_moments
+ USE processing, ONLY: compute_nadiab_moments, compute_fluid_moments
USE utility, ONLY: checkfield
-
IMPLICIT NONE
INTEGER :: num_step
@@ -25,29 +23,30 @@ SUBROUTINE stepon
!----- BEFORE: All fields are updated for step = n
! Compute right hand side from current fields
! N_rhs(N_n, nadia_n, phi_n, S_n, Tcoll_n)
- IF(KIN_E) CALL moments_eq_rhs_e
- CALL moments_eq_rhs_i
+ CALL compute_RHS
! ---- step n -> n+1 transition
! Advance from updatetlevel to updatetlevel+1 (according to num. scheme)
CALL advance_time_level
+ ! ----
+
! Update moments with the hierarchy RHS (step by step)
! N_n+1 = N_n + N_rhs(n)
CALL advance_moments
- ! Closure enforcement of N_n+1
+ ! Closure enforcement of moments
CALL apply_closure_model
! Exchanges the ghosts values of N_n+1
- CALL update_ghosts
+ CALL update_ghosts_p_moments
+ CALL update_ghosts_z_moments
+
! Update electrostatic potential phi_n = phi(N_n+1)
CALL poisson
- ! Update non adiabatic moments n -> n+1
- CALL compute_nadiab_moments
- ! Update collision C_n+1 = C(N_n+1)
- CALL compute_TColl
- ! Update nonlinear term S_n -> S_n+1(phi_n+1,N_n+1)
- CALL compute_Sapj
+ CALL update_ghosts_z_phi
+
+ ! Numerical experiments
! Store or cancel/maintain zonal modes artificially
CALL play_with_modes
+
!- Check before next step
CALL checkfield_all()
IF( nlend ) EXIT ! exit do loop
@@ -57,6 +56,23 @@ SUBROUTINE stepon
END DO
CONTAINS
+ !!!! Basic structure to simplify stepon
+ SUBROUTINE compute_RHS
+ USE moments_eq_rhs, ONLY: moments_eq_rhs_e,moments_eq_rhs_i
+ USE collision, ONLY: compute_TColl
+ USE nonlinear, ONLY: compute_Sapj
+ IMPLICIT NONE
+ ! compute auxiliary non adiabatic moments arrays
+ CALL compute_nadiab_moments
+ ! compute nonlinear term ("if linear" is included inside)
+ CALL compute_Sapj
+ ! compute collision term ("if coll, if nu >0" is included inside)
+ CALL compute_TColl
+ ! compute the moments equation rhs for electrons and ions
+ IF (KIN_E) &
+ CALL moments_eq_rhs_e
+ CALL moments_eq_rhs_i
+ END SUBROUTINE compute_RHS
SUBROUTINE checkfield_all ! Check all the fields for inf or nan
! Execution time start
diff --git a/src/utility_mod.F90 b/src/utility_mod.F90
index a295575e..577a130f 100644
--- a/src/utility_mod.F90
+++ b/src/utility_mod.F90
@@ -114,8 +114,10 @@ SUBROUTINE manual_3D_bcast(field_)
IMPLICIT NONE
COMPLEX(dp), INTENT(INOUT) :: field_(ikxs:ikxe,ikys:ikye,izs:ize)
COMPLEX(dp) :: buffer(ikxs:ikxe,ikys:ikye,izs:ize)
- INTEGER :: i_, root, world_rank, world_size
+ INTEGER :: i_, root, world_rank, world_size, count
root = 0;
+ count = (ikxe-ikxs+1) * (ikye-ikys+1) * (ize-izs+1);
+
CALL MPI_COMM_RANK(comm_p,world_rank,ierr)
CALL MPI_COMM_SIZE(comm_p,world_size,ierr)
IF (world_size .GT. 1) THEN
@@ -132,11 +134,11 @@ SUBROUTINE manual_3D_bcast(field_)
! Send it to all the other processes
DO i_ = 0,num_procs_p-1
IF (i_ .NE. world_rank) &
- CALL MPI_SEND(buffer, local_nkx * Nky * Nz, MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
+ CALL MPI_SEND(buffer, count, MPI_DOUBLE_COMPLEX, i_, 0, comm_p, ierr)
ENDDO
ELSE
! Recieve buffer from root
- CALL MPI_RECV(buffer, local_nkx * Nky * Nz, MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
+ CALL MPI_RECV(buffer, count, MPI_DOUBLE_COMPLEX, root, 0, comm_p, MPI_STATUS_IGNORE, ierr)
! Write it in phi
DO ikx = ikxs,ikxe
DO iky = ikys,ikye
diff --git a/wk/analysis_3D.m b/wk/analysis_3D.m
index cb59d386..7c194520 100644
--- a/wk/analysis_3D.m
+++ b/wk/analysis_3D.m
@@ -1,16 +1,17 @@
-addpath(genpath('../matlab')) % ... add
-addpath(genpath('../matlab/plot')) % ... add
-addpath(genpath('../matlab/compute')) % ... add
-addpath(genpath('../matlab/load')) % ... add
+helazdir = '/home/ahoffman/HeLaZ/';
+addpath(genpath([helazdir,'matlab'])) % ... add
+addpath(genpath([helazdir,'matlab/plot'])) % ... add
+addpath(genpath([helazdir,'matlab/compute'])) % ... add
+addpath(genpath([helazdir,'matlab/load'])) % ... add
%% Load the results
-LOCALDIR = ['../results/',outfile,'/'];
+LOCALDIR = [helazdir,'results/',outfile,'/'];
MISCDIR = ['/misc/HeLaZ_outputs/results/',outfile,'/'];
system(['mkdir -p ',MISCDIR]);
CMD = ['rsync ', LOCALDIR,'outputs* ',MISCDIR]; disp(CMD);
system(CMD);
% Load outputs from jobnummin up to jobnummax
-JOBNUMMIN = 00; JOBNUMMAX = 20;
+JOBNUMMIN = 00; JOBNUMMAX = 10;
data = compile_results(MISCDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from first output found to JOBNUMMAX if existing
@@ -19,17 +20,18 @@ default_plots_options
disp('Plots')
FMT = '.fig';
-if 1
+if 0
%% Space time diagramm (fig 11 Ivanov 2020)
-TAVG_0 = 1000; TAVG_1 = 11000; % Averaging times duration
+TAVG_0 = 0; TAVG_1 = 400; % Averaging times duration
+compz = 'avg';
% chose your field to plot in spacetime diag (uzf,szf,Gx)
-fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'v_y',1);
+fig = plot_radial_transport_and_spacetime(data,TAVG_0,TAVG_1,'phi',1,compz);
save_figure(data,fig)
end
if 0
%% statistical transport averaging
-options.T = [1500 2500];
+options.T = [16000 17000];
fig = statistical_transport_averaging(data,options);
end
@@ -39,16 +41,18 @@ if 0
options.INTERP = 1;
options.POLARPLOT = 0;
options.NAME = '\phi';
-% options.NAME = 'v_x';
+% options.NAME = 'N_i^{00}';
+% options.NAME = 'v_y';
% options.NAME = 'n_i^{NZ}';
% options.NAME = '\Gamma_x';
% options.NAME = 'n_i';
-% options.PLAN = 'xy';
+options.PLAN = 'yz';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 1;
-% options.TIME = data.Ts5D;
-options.TIME = 900:10:2000;
+% options.TIME = dat.Ts5D;
+options.TIME = 0:0.5:400;
+data.EPS = 0.1;
data.a = data.EPS * 2000;
create_film(data,options,'.gif')
end
@@ -56,19 +60,20 @@ end
if 0
%% 2D snapshots
% Options
-options.INTERP = 1;
+options.INTERP = 0;
options.POLARPLOT = 0;
options.AXISEQUAL = 1;
-% options.NAME = '\phi';
-options.NAME = 'v_y';
-% options.NAME = 'n_e^{NZ}';
+options.NAME = '\phi';
+% options.NAME = 'v_y';
+% options.NAME = 'N_i^{00}';
+% options.NAME = 'T_i';
% options.NAME = '\Gamma_x';
% options.NAME = 'k^2n_e';
-options.PLAN = 'xy';
+options.PLAN = 'yz';
% options.NAME = 'f_e';
% options.PLAN = 'sx';
options.COMP = 'avg';
-options.TIME = [000:200:1000];
+options.TIME = [10 20 100];
data.a = data.EPS * 2000;
fig = photomaton(data,options);
save_figure(data,fig)
@@ -78,9 +83,10 @@ if 0
%% 3D plot on the geometry
options.INTERP = 1;
options.NAME = 'n_i';
-options.PLANES = 1;
-options.TIME = [0 500 1000];
-data.rho_o_R = 1e-3; % Sound larmor radius over Machine size ratio
+options.PLANES = 1:3:30;
+options.TIME = [100];
+options.PLT_MTOPO = 1;
+data.rho_o_R = 2e-3; % Sound larmor radius over Machine size ratio
fig = show_geometry(data,options);
save_figure(data,fig)
end
@@ -102,25 +108,24 @@ end
if 0
%% Hermite-Laguerre spectrum
% options.TIME = 'avg';
-options.TIME = 1000:4000;
options.P2J = 1;
-options.ST = 0;
+options.ST = 1;
options.NORMALIZED = 0;
fig = show_moments_spectrum(data,options);
save_figure(data,fig)
end
if 0
-%% 1D spectrum
-options.TIME = 5000:10:5050;
-options.NORM = 1;
+%% Time averaged spectrum
+options.TIME = 1000:5000;
+options.NORM = 0;
options.NAME = '\phi';
-% options.NAME = 'n_i';
-% options.NAME ='\Gamma_x';
+options.NAME = 'n_i';
+options.NAME ='\Gamma_x';
options.PLAN = 'kxky';
-options.COMPZ = 1;
+options.COMPZ = 'avg';
options.OK = 0;
-options.COMPXY = 'sum';
+options.COMPXY = 0;
options.COMPT = 'avg';
options.PLOT = 'semilogy';
fig = spectrum_1D(data,options);
@@ -145,9 +150,9 @@ end
if 0
%% Mode evolution
-options.NORMALIZED = 0;
+options.NORMALIZED = 1;
options.K2PLOT = 0.01:0.01:1.0;
-options.TIME = 4900:1:5100;
+options.TIME = 20:1:600;
options.NMA = 1;
options.NMODES = 20;
fig = mode_growth_meter(data,options);
@@ -155,9 +160,21 @@ save_figure(data,fig)
end
if 0
-%% ZF caracteristics (space time diagrams
+%% ZF caracteristics (space time diagrams)
TAVG_0 = 200; TAVG_1 = 3000; % Averaging times duration
% chose your field to plot in spacetime diag (uzf,szf,Gx)
fig = ZF_spacetime(data,TAVG_0,TAVG_1);
save_figure(data,fig)
end
+
+if 0
+%% Metric infos
+fig = plot_metric(data);
+end
+
+if 1
+%% linear growth rate for 3D fluxtube
+trange = [0 100];
+nplots = 1;
+lg = compute_fluxtube_growth_rate(data,trange,nplots);
+end
\ No newline at end of file
diff --git a/wk/analysis_header.m b/wk/analysis_header.m
index 4ac8664b..18b07606 100644
--- a/wk/analysis_header.m
+++ b/wk/analysis_header.m
@@ -1,166 +1,10 @@
%% Directory of the simulation
% if 1% Local results
outfile ='';
-outfile ='';
-outfile ='';
-outfile ='';
-% outfile ='debug/ppj_init';
-%% nu = 5e-1
-% Sugama
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_5e-01_SGGK';% also in 7x4
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_5e-01_SGGK';
-% outfile ='Hallenbert_nu_5e-01/200x32_7x4_L_120_kN_1.7_kT_0.425_nu_5e-01_SGGK';% also in 7x4
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_5e-01_SGGK';
-% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_5e-01_SGGK';%also in 7x4
-
-%% nu = 1e-1
-% Landau
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_LDGK';
-
-% Sugama
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_SGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_SGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_SGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_SGGK';
-
-% Dougherty
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_DGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_DGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_DGGK';
-% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_DGGK';
-
-%% nu = 5e-2
-% outfile ='Hallenbert_nu_5e-02/200x32_11x6_L_120_kN_1.8_kT_0.45_nu_5e-02_SGGK';%For GENE benchmark % to analyse (added HD)
-% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
-
-% testing various NL closures
-% outfile ='Hallenbert_nu_5e-02/200x32_7x4_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
-
-% outfile ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
-% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
-
-% outfile ='Hallenbert_nu_5e-02/256x64_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
-% outfile ='Hallenbert_nu_5e-02/256x64_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
-% outfile ='Hallenbert_nu_5e-02/200x32_21x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
-% outfile ='Hallenbert_nu_5e-02/200x32_17x9_L_120_kN_1.8_kT_0.45_nu_5e-02_SGDK';
-
-% outfile ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_DGGK';
-% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_DGGK';
-% outfile ='Hallenbert_nu_5e-02/128x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_FCGK';
-
-%% nu = 1e-2
-% Landau
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_LDGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_LDGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_LDGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_LDGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_LDGK';
-
-% Sugama
-% outfile ='kobayashi_2015_fig1/150x150_5x3_L_100_kN_1.4_nu_5e-03_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/300x64_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_SGGK'; % To analyse (added HD)
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
-% Dougherty
-% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_DGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_DGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_8x5_L_120_kN_1.6_kT_0.4_nu_0e+00_DGGK';
-% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_DGGK';
-
-%% nu = 5e-3
-% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
-% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
-% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
-% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
-
-%% nu = 0
-
-% outfile ='Hallenbert_fig2a/200x32_21x11_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
-% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
-% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
-% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2';
-% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2';
-% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuSGGK_0.1_muxy_1e-2';
-% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nuSGGK_0.1_muxy_1e-2';
-% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuLDGK_0.1_muxy_1e-2';
-% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuLRGK_0.1_muxy_1e-2';
-
-
-% outfile ='Hallenbert_fig2b/200x32_11x6_Lx_240_Ly_120_kN_2.5_eta_0.25_nu_0_muxy_1e-1';
-% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nu_0_muxy_1e-1';
-% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuSGGK_0.1_muxy_1e-1';
-% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuDGGK_0.1_muxy_1e-1';
-% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuLDGK_0.1_muxy_1e-1';
-% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuLRGK_0.1_muxy_1e-1';
-
-% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuSGGK_0.1_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nuSGGK_0.1_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuDGGK_0.1_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nuDGGK_0.1_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLDGK_0.1_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLRGK_0.1_muxy_5e-2';
-% outfile ='Hallenbert_fig2c/200x32_9x5_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLRGK_0.1_muxy_5e-2';
-
-%% Transport scan
-% outfile = 'nu_0.1_transport_scan/colless_kn_1.7_to_2.0';
-% outfile = 'nu_0.1_transport_scan/colless_kn_2.1_to_2.5';
-
-% outfile = 'nu_0.1_transport_scan/LB_kn_2.0';
-
-% outfile = 'nu_0.1_transport_scan/DG_kn_1.8_to_2.1';
-% outfile = 'nu_0.1_transport_scan/DG_kn_2.2_to_2.5';
-% outfile = 'nu_0.1_transport_scan/DG_conv_kN_1.9';
-
-% outfile = 'nu_0.1_transport_scan/SG_kn_1.7_to_2.0';
-% outfile = 'nu_0.1_transport_scan/SG_10x5_conv_test';
-% outfile = 'nu_0.1_transport_scan/SG_kn_2.2_to_2.5';
-
-% outfile = 'nu_0.1_transport_scan/LD_kn_2.0_to_2.5';
-% outfile = 'nu_0.1_transport_scan/LD_kn_1.7_to_2.5';
-
-% outfile = 'nu_0.1_transport_scan/LR_kn_1.7_to_2.0';
-% outfile = 'nu_0.1_transport_scan/LR_kn_2.1_to_2.5';
-
-% outfile = 'nu_0.1_transport_scan/colless_kn_2.2_Lx1.5';
-% outfile = 'nu_0.1_transport_scan/colless_kn_2.2_HD';
-
-% outfile = 'nu_0.1_transport_scan/colless_kn_1.6_HD';
-
-% outfile = 'nu_0.1_transport_scan/large_box_kN_2.1_nu_0.1';
-% outfile = 'nu_0.1_transport_scan/large_box_kN_2.0_nu_0.1';
-
-outfile = 'predator_prey_nu_scan/DG_Kn_1.7_nu_0.01';
-
-% outfile = 'ZF_damping_linear_nu_0_20x10_kn_1.6_GK/LR_4x2_nu_0.1';
-% outfile = 'ZF_damping_nu_0_20x10_kn_1.6_GK/HSG_4x2_nu_0.1';
-% outfile = 'ZF_damping_nu_0_5x3_kn_2.5_GK/LR_4x2_nu_0.1';
-% outfile = 'hacked_sugama/hacked_B_kn_1.6_200x32_L_120x60_nu_0.1';
+outfile ='linear_shearless_cyclone/dbg';
+% outfile ='debug/test_zpar/';
+% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_1.0_colless/';
+% outfile ='shearless_cyclone/128x128x16x4x2_L_120_CTC_0.8/';
-% outfile = 'shearless_cyclone/200x32x24_5x4_Lx_120_Ly_60_q0_1.4_e_0.18_kN_2.22_kT_6.9_nu_0';
-% else% Marconi results
-% outfile ='';
-% outfile ='';
-% outfile ='';
-% outfile ='';
-% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A_new/300x300_5x3_L_120_kN_1.6667_nu_1e-01_SGGK/out.txt';
-% % outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
-% % BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
-% MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
-% end
analysis_3D
\ No newline at end of file
diff --git a/wk/analysis_header_2DZP.m b/wk/analysis_header_2DZP.m
new file mode 100644
index 00000000..cd4f356a
--- /dev/null
+++ b/wk/analysis_header_2DZP.m
@@ -0,0 +1,171 @@
+%% Directory of the simulation
+% if 1% Local results
+outfile ='';
+outfile ='';
+outfile ='';
+outfile ='';
+% outfile ='debug/ppj_init';
+%% nu = 5e-1
+% Sugama
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_5e-01_SGGK';% also in 7x4
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_5e-01_SGGK';
+% outfile ='Hallenbert_nu_5e-01/200x32_7x4_L_120_kN_1.7_kT_0.425_nu_5e-01_SGGK';% also in 7x4
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_5e-01_SGGK';
+% outfile ='Hallenbert_nu_5e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_5e-01_SGGK';%also in 7x4
+
+%% nu = 1e-1
+% Landau
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/150x50_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_LDGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_LDGK';
+
+% Sugama
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_SGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-01_SGGK';
+
+% Dougherty
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-01_DGGK';
+% outfile ='Hallenbert_nu_1e-01/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-01_DGGK';
+
+%% nu = 5e-2
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_L_120_kN_1.8_kT_0.45_nu_5e-02_SGGK';%For GENE benchmark % to analyse (added HD)
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
+
+% testing various NL closures
+% outfile ='Hallenbert_nu_5e-02/200x32_7x4_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGGK';
+
+% outfile ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+
+% outfile ='Hallenbert_nu_5e-02/256x64_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+% outfile ='Hallenbert_nu_5e-02/256x64_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+% outfile ='Hallenbert_nu_5e-02/200x32_21x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+% outfile ='Hallenbert_nu_5e-02/200x32_17x9_L_120_kN_1.8_kT_0.45_nu_5e-02_SGDK';
+
+% outfile ='Hallenbert_nu_5e-02/200x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_DGGK';
+% outfile ='Hallenbert_nu_5e-02/200x32_11x6_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_DGGK';
+% outfile ='Hallenbert_nu_5e-02/128x32_5x3_Lx_120_Ly_60_kN_1.8_kT_0.45_nu_5e-02_FCGK';
+
+%% nu = 1e-2
+% Landau
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_LDGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_LDGK';
+
+% Sugama
+% outfile ='kobayashi_2015_fig1/150x150_5x3_L_100_kN_1.4_nu_5e-03_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.5_kT_0.375_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.6_kT_0.4_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/300x64_5x3_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.7_kT_0.425_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_SGGK'; % To analyse (added HD)
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_11x6_L_120_kN_1.9_kT_0.475_nu_1e-02_SGGK';
+% Dougherty
+% outfile ='Hallenbert_nu_1e-02/200x32_7x4_L_120_kN_1.5_kT_0.375_nu_1e-02_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.6_kT_0.4_nu_1e-02_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_8x5_L_120_kN_1.6_kT_0.4_nu_0e+00_DGGK';
+% outfile ='Hallenbert_nu_1e-02/200x32_5x3_L_120_kN_1.8_kT_0.45_nu_1e-02_DGGK';
+
+%% nu = 5e-3
+% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
+% outfile ='Hallenbert_nu_5e-03/200x32_5x3_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
+% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_mux_5e-2_muy_6e-1';
+% outfile ='Hallenbert_nu_5e-03/200x32_11x6_Lx_120_Ly_60_kN_1.8_eta_0.25_nuSG_5e-03_muxy_5e-2';
+
+%% nu = 0
+
+% outfile ='Hallenbert_fig2a/200x32_21x11_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nu_0_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nuDGGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuSGGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_11x6_Lx_120_Ly_60_kN_1.6_eta_0.4_nuSGGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuLDGK_0.1_muxy_1e-2';
+% outfile ='Hallenbert_fig2a/200x32_5x3_Lx_120_Ly_60_kN_1.6_eta_0.4_nuLRGK_0.1_muxy_1e-2';
+
+
+% outfile ='Hallenbert_fig2b/200x32_11x6_Lx_240_Ly_120_kN_2.5_eta_0.25_nu_0_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nu_0_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuSGGK_0.1_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuDGGK_0.1_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuLDGK_0.1_muxy_1e-1';
+% outfile ='Hallenbert_fig2b/200x32_5x3_Lx_240_Ly_120_kN_2.5_eta_0.25_nuLRGK_0.1_muxy_1e-1';
+
+% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nu_0_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuSGGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nuSGGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuDGGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_11x6_Lx_120_Ly_60_kN_2.0_eta_0.25_nuDGGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLDGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_5x3_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLRGK_0.1_muxy_5e-2';
+% outfile ='Hallenbert_fig2c/200x32_9x5_Lx_120_Ly_60_kN_2.0_eta_0.25_nuLRGK_0.1_muxy_5e-2';
+
+%% Transport scan
+% outfile = 'nu_0.1_transport_scan/colless_kn_1.7_to_2.0';
+% outfile = 'nu_0.1_transport_scan/colless_kn_2.1_to_2.5';
+
+% outfile = 'nu_0.1_transport_scan/LB_kn_2.0';
+
+% outfile = 'nu_0.1_transport_scan/DG_kn_1.8_to_2.1';
+% outfile = 'nu_0.1_transport_scan/DG_kn_2.2_to_2.5';
+% outfile = 'nu_0.1_transport_scan/DG_conv_kN_1.9';
+
+% outfile = 'nu_0.1_transport_scan/SG_kn_1.7_to_2.0';
+% outfile = 'nu_0.1_transport_scan/SG_10x5_conv_test';
+% outfile = 'nu_0.1_transport_scan/SG_kn_2.2_to_2.5';
+
+% outfile = 'nu_0.1_transport_scan/LD_kn_2.0_to_2.5';
+% outfile = 'nu_0.1_transport_scan/LD_kn_1.7_to_2.5';
+
+% outfile = 'nu_0.1_transport_scan/LR_kn_1.7_to_2.0';
+% outfile = 'nu_0.1_transport_scan/LR_kn_2.1_to_2.5';
+
+% outfile = 'nu_0.1_transport_scan/colless_kn_2.2_Lx1.5';
+% outfile = 'nu_0.1_transport_scan/colless_kn_2.2_HD';
+
+% outfile = 'nu_0.1_transport_scan/colless_kn_1.6_HD';
+
+% outfile = 'nu_0.1_transport_scan/large_box_kN_2.1_nu_0.1';
+% outfile = 'nu_0.1_transport_scan/large_box_kN_2.0_nu_0.1';
+
+% outfile = 'predator_prey_nu_scan/DG_Kn_1.7_nu_0.01';
+
+% outfile = 'ZF_damping_linear_nu_0_20x10_kn_1.6_GK/LR_4x2_nu_0.1';
+% outfile = 'ZF_damping_nu_0_20x10_kn_1.6_GK/HSG_4x2_nu_0.1';
+% outfile = 'ZF_damping_nu_0_5x3_kn_2.5_GK/LR_4x2_nu_0.1';
+% outfile = 'hacked_sugama/hacked_B_kn_1.6_200x32_L_120x60_nu_0.1';
+
+% outfile = 'shearless_cyclone/200x32x24_5x4_Lx_120_Ly_60_q0_1.4_e_0.18_kN_2.22_kT_6.9_nuLR_0.01_adiab_e';
+% outfile = 'shearless_cyclone/no_sg_128x32x36_6x3_Lx_120_Ly_60_q0_1.4_e_0.18_kN_2.22_kT_6.9_adiab_e';
+% outfile = 'shearless_cyclone/sgrid_128x64x32_4x2_Lx_100_Ly_120_q0_1.4_e_0.18_kN_2.22_kT_6.96_adiab_e';
+% outfile = 'shearless_cyclone/sgrid_128x64x32_4x2_Lx_100_Ly_120_q0_1.4_e_0.18_kN_1.78_kT_5.52_adiab_e';
+% outfile = 'linear_shearless_cyclone/4_2_cyclone_1.0';
+% outfile = 'linear_shearless_cyclone/test_fmom';
+% else% Marconi results
+% outfile ='';
+% outfile ='';
+% outfile ='';fd
+% outfile ='';
+% outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A_new/300x300_5x3_L_120_kN_1.6667_nu_1e-01_SGGK/out.txt';
+% % outfile ='/marconi_scratch/userexternal/ahoffman/HeLaZ/results/simulation_A/300x150_L_120_P_8_J_4_eta_0.6_nu_1e-01_SGGK_mu_0e+00/out.txt';
+% % BASIC.RESDIR = ['../',outfile(46:end-8),'/'];
+% MISCDIR = ['/misc/HeLaZ_outputs/',outfile(46:end-8),'/'];
+% end
+
+analysis_3D
\ No newline at end of file
diff --git a/wk/benchmark_HeLaZ_gene_transport_zpinch.m b/wk/benchmark_HeLaZ_gene_transport_zpinch.m
index f0e5ddc6..b10312db 100644
--- a/wk/benchmark_HeLaZ_gene_transport_zpinch.m
+++ b/wk/benchmark_HeLaZ_gene_transport_zpinch.m
@@ -119,9 +119,9 @@ xlim([ 1.55 2.55]);
%% Burst study Kn = 1.7
Kn = 1.7;
NU = [0.0e+0 1.0e-2 2.0e-2 3.0e-2 4.0e-2 5.0e-2 6.0e-2 7.0e-2 8.0e-2 9.0e-2 1.0e-1];
-SGGK_200x32x05x03 = [1.0e-2 2.4e-2 2.8e-2 3.5e-2 4.5e-2 5.5e-2 6.5e-2 7.5e-2 8.3e-2 8.8e-2 1.0e-1];
-DGGK_200x32x05x03 = [1.0e-2 7.0e-3 7.3e-3 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
-LDGK_200x32x05x03 = [1.0e-2 8.5e-2 1.4e-1 2.0e-1 0.0e-2 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0 0.0e+0];
+SGGK_200x32x05x03 = [1.0e-2 1.1e-2 1.7e-2 2.3e-2 3.1e-2 4.1e-2 5.1e-2 6.1e-2 8.3e-2 8.8e-2 1.0e-1];
+DGGK_200x32x05x03 = [1.0e-2 1.0e-3 1.1e-3 2.3e-3 3.0e-3 3.0e-3 3.7e-3 4.0e-3 5.7e-3 5.8e-3 9.2e-3];
+LDGK_200x32x05x03 = [1.0e-2 9.0e-2 1.5e-1 2.2e-1 2.8e-1 3.5e-1 4.0e-1 4.5e-1 5.0e-1 5.6e-1 6.2e-1];
% NOCO_200x32x05x03 = [0.0e+0 0.0e+0 0.0e+0 0.0e+0];
figure
plot(NU,SGGK_200x32x05x03/Kn,'s','Color',clr_(1,:)); hold on
diff --git a/wk/gene_analysis_3D.m b/wk/gene_analysis_3D.m
new file mode 100644
index 00000000..2289b64f
--- /dev/null
+++ b/wk/gene_analysis_3D.m
@@ -0,0 +1,84 @@
+% folder = '/misc/gene_results/shearless_cyclone/miller_output_1.0/';
+% folder = '/misc/gene_results/shearless_cyclone/miller_output_0.8/';
+% folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_1.0/';
+folder = '/misc/gene_results/shearless_cyclone/s_alpha_output_0.8/';
+% folder = '/misc/gene_results/HP_fig_2b_mu_5e-2/';
+% folder = '/misc/gene_results/HP_fig_2c_mu_5e-2/';
+gene_data = load_gene_data(folder);
+gene_data = rotate_c_plane_nxnky_to_nkxny(gene_data);
+if 1
+%% Space time diagramm (fig 11 Ivanov 2020)
+TAVG_0 = 500; TAVG_1 = 700; % Averaging times duration
+% chose your field to plot in spacetime diag (uzf,szf,Gx)
+field = 'phi';
+compz = 'avg';
+nmvm = 1;
+fig = plot_radial_transport_and_spacetime(gene_data,TAVG_0,TAVG_1,field,nmvm,compz);
+% save_figure(data,fig)
+end
+
+if 0
+%% 2D snapshots
+% Options
+options.INTERP = 0;
+options.POLARPLOT = 0;
+options.AXISEQUAL = 1;
+options.NAME = '\phi';
+% options.NAME = 'v_y';
+% options.NAME = 'T_i';
+% options.NAME = '\Gamma_x';
+% options.NAME = 'k^2n_e';
+options.PLAN = 'xy';
+% options.NAME ='f_e';
+% options.PLAN = 'sx';
+options.COMP = 9;
+options.TIME = [100 200 400];
+data.a = data.EPS * 2000;
+fig = photomaton(gene_data,options);
+save_figure(gene_data,fig)
+end
+
+if 0
+%% MOVIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Options
+options.INTERP = 1;
+options.POLARPLOT = 0;
+options.NAME = '\phi';
+% options.NAME = 'v_y';
+% options.NAME = 'n_i^{NZ}';
+% options.NAME = '\Gamma_x';
+% options.NAME = 'n_i';
+options.PLAN = 'xy';
+% options.NAME = 'f_e';
+% options.PLAN = 'sx';
+options.COMP = 9;
+options.TIME = 0:700;
+gene_data.a = data.EPS * 2000;
+create_film(gene_data,options,'.gif')
+end
+
+if 0
+%% Geometry
+names = {'$g^{xx}$','$g^{xy}$','$g^{xz}$','$g^{yy}$','$g^{yz}$','$g^{zz}$',...
+ '$B_0$','$\partial_x B_0$','$\partial_y B_0$','$\partial_z B_0$',...
+ '$J$','$R$','$\phi$','$Z$','$\partial_R x$','$\partial_Z x$'};
+figure;
+subplot(311)
+ for i = 1:6
+ plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(gene_data.z),max(gene_data.z)]); legend('show'); title('GENE geometry');
+
+subplot(312)
+ for i = 7:10
+ plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
+
+subplot(313)
+ for i = 11:16
+ plot(gene_data.z, gene_data.geo_arrays(:,i),'DisplayName',names{i}); hold on;
+ end
+ xlim([min(gene_data.z),max(gene_data.z)]); legend('show');
+
+end
diff --git a/wk/tmp.txt b/wk/tmp.txt
deleted file mode 100644
index e69de29b..00000000
--
GitLab