diff --git a/Makefile b/Makefile
index ef0f1f81483d40d7962da7002b2016bdedd772f3..d8f9cd131c06c5fa675444dbd4904947c55b313c 100644
--- a/Makefile
+++ b/Makefile
@@ -118,6 +118,7 @@ FOBJ=$(OBJDIR)/advance_field_mod.o $(OBJDIR)/array_mod.o $(OBJDIR)/auxval.o \
$(OBJDIR)/basic_mod.o $(OBJDIR)/coeff_mod.o $(OBJDIR)/closure_mod.o $(OBJDIR)/circular_mod.o \
$(OBJDIR)/collision_mod.o $(OBJDIR)/nonlinear_mod.o $(OBJDIR)/control.o \
$(OBJDIR)/diagnose.o $(OBJDIR)/diagnostics_par_mod.o $(OBJDIR)/endrun.o \
+$(OBJDIR)/ExB_shear_flow_mod.o \
$(OBJDIR)/fields_mod.o $(OBJDIR)/fourier_mod.o $(OBJDIR)/geometry_mod.o \
$(OBJDIR)/ghosts_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/inital.o \
$(OBJDIR)/initial_par_mod.o $(OBJDIR)/lag_interp_mod.o $(OBJDIR)/main.o \
@@ -146,7 +147,8 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/auxval.o : src/auxval.F90 \
$(OBJDIR)/fourier_mod.o $(OBJDIR)/memory.o $(OBJDIR)/model_mod.o \
$(OBJDIR)/geometry_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/numerics_mod.o \
- $(OBJDIR)/parallel_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/CLA_mod.o
+ $(OBJDIR)/parallel_mod.o $(OBJDIR)/processing_mod.o $(OBJDIR)/CLA_mod.o \
+ $(OBJDIR)/ExB_shear_flow_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/auxval.F90 -o $@
$(OBJDIR)/basic_mod.o : src/basic_mod.F90 \
@@ -206,6 +208,11 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/prec_const_mod.o $(OBJDIR)/basic_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/endrun.F90 -o $@
+ $(OBJDIR)/ExB_shear_flow_mod.o : src/ExB_shear_flow_mod.F90 \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/grid_mod.o $(OBJDIR)/prec_const_mod.o\
+ $(OBJDIR)/geometry_mod.o $(OBJDIR)/numerics_mod.o
+ $(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ExB_shear_flow_mod.F90 -o $@
+
$(OBJDIR)/fields_mod.o : src/fields_mod.F90 \
$(OBJDIR)/prec_const_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/fields_mod.F90 -o $@
@@ -221,7 +228,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/geometry_mod.F90 -o $@
$(OBJDIR)/ghosts_mod.o : src/ghosts_mod.F90 \
- $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o \
+ $(OBJDIR)/basic_mod.o $(OBJDIR)/fields_mod.o $(OBJDIR)/grid_mod.o\
$(OBJDIR)/geometry_mod.o $(OBJDIR)/prec_const_mod.o $(OBJDIR)/time_integration_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/ghosts_mod.F90 -o $@
@@ -250,7 +257,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_mod.o
$(OBJDIR)/memory.o : src/memory.F90 $ \
$(OBJDIR)/array_mod.o $(OBJDIR)/basic_mod.o $(OBJDIR)/collision_mod.o\
- $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
+ $(OBJDIR)/fields_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
$(OBJDIR)/grid_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/memory.F90 -o $@
@@ -321,7 +328,7 @@ $(OBJDIR)/time_integration_mod.o $(OBJDIR)/utility_mod.o $(OBJDIR)/CLA_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)/solve_EM_fields.o\
$(OBJDIR)/utility_mod.o $(OBJDIR)/model_mod.o $(OBJDIR)/time_integration_mod.o \
- $(OBJDIR)/CLA_mod.o
+ $(OBJDIR)/CLA_mod.o $(OBJDIR)/ExB_shear_flow_mod.o
$(F90) -c $(F90FLAGS) $(FPPFLAGS) $(EXTMOD) $(EXTINC) src/stepon.F90 -o $@
$(OBJDIR)/tesend.o : src/tesend.F90 \
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index 9993b2653f967e3fab20fd595a2bc67ce0683e8e..b9aa6fe0b3d9c109569a0881eb77c47897d42823 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -1,4 +1,4 @@
-function [FIGURE] = mode_growth_meter(DATA,OPTIONS)
+function [FIGURE, kykx, wkykx, ekykx] = mode_growth_meter(DATA,OPTIONS)
NORMALIZED = OPTIONS.NORMALIZED;
Nma = OPTIONS.NMA; %Number moving average
@@ -30,7 +30,9 @@ TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
[~,ikzf] = max(mean(squeeze(abs(field(1,:,FRAMES))),2));
+% Amplitude ratio method to measure growth rates and freq.
[wkykx, ekykx] = compute_growth_rates(DATA.PHI(:,:,:,FRAMES),DATA.Ts3D(FRAMES));
+kykx = meshgrid(DATA.grids.ky,DATA.grids.kx)';
FIGURE.fig = figure; %set(gcf, 'Position', [100 100 1200 700])
FIGURE.FIGNAME = 'mode_growth_meter';
@@ -95,35 +97,13 @@ for i = 1:2
gamma = MODES;
amp = MODES;
- % w_ky = zeros(numel(MODES),numel(FRAMES)-1);
- % ce = zeros(numel(MODES),numel(FRAMES));
i_=1;
+ % Alternative method to measure growth
for ik = MODES
-
to_measure = log(plt(field,ik));
gr = polyfit(t(it1:it2),to_measure(it1:it2),1);
gamma(i_) = gr(1);
amp(i_) = gr(2);
-
- % % Second method for measuring growth rate (measures frequ. too)
- % if MODES_SELECTOR == 1
- % to_measure = reshape(DATA.PHI(ik,1,:,FRAMES),DATA.grids.Nz,numel(FRAMES));
- % else
- % to_measure = reshape(DATA.PHI(1,ik,:,FRAMES),DATA.grids.Nz,numel(FRAMES));
- % end
- % for it = 2:numel(FRAMES)
- % phi_n = to_measure(:,it);
- % phi_nm1 = to_measure(:,it-1);
- % dt = t(it)-t(it-1);
- % ZS = sum(phi_nm1,1); %sum over z
- %
- % wl = log(phi_n./phi_nm1)/dt;
- % w_ky(i_,it) = squeeze(sum(wl.*phi_nm1,1)./ZS);
- %
- % % for iky = 1:numel(w_ky(:,it))
- % % ce(iky,it) = abs(sum(abs(w_ky(iky,it)-wl(iky,:)).^2.*phi_nm1(iky,:),2)./ZS(iky,:));
- % % end
- % end
i_=i_+1;
end
mod2plot = 2:min(Nmax,OPTIONS.NMODES+1);
@@ -135,7 +115,7 @@ for i = 1:2
pclr = pcolor(XX,YY,abs(plt(fftshift(field,MODES_SELECTOR),1:numel(k))));set(pclr, 'edgecolor','none'); hold on;
set(gca,'YDir','normal')
% xlim([t(1) t(end)]); %ylim([1e-5 1])
- xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /\rho_s$');
+ xlabel(['$',kstr,'\rho_s$']); ylabel('$t c_s /R_0$');
title([MODESTR,', ',zstrcomp])
% subplot(2+d,3,2+3*(i-1))
@@ -151,13 +131,25 @@ for i = 1:2
plot(t(it1)*[1 1],[1e-10 1],'--k')
plot(t(it2)*[1 1],[1e-10 1],'--k')
xlim([t(1) t(end)]); %ylim([1e-5 1])
- xlabel('$t c_s /\rho_s$'); ylabel(['$|\phi_{',kstr,'}|$']);
+ xlabel('$t c_s /R_0$'); ylabel(['$|\phi_{',kstr,'}|$']);
title('Measure time window')
% subplot(2+d,3,3+3*(i-1))
FIGURE.axes(3+3*(i-1)) = subplot(2+d,3,3+3*(i-1),'parent',FIGURE.fig);
% yyaxis("left")
- errorbar(k(MODES),real(omega(ik)),real(err(ik)),'-k',...
+ if OPTIONS.GOK2
+ x_ = k(MODES);
+ y_ = real(omega(ik))./k(MODES).^2;
+ e_ = real(err(ik))./k(MODES).^2;
+ lb_= '\gamma/k^2';
+ else
+ x_ = k(MODES);
+ y_ = real(omega(ik));
+ e_ = real(err(ik));
+ lb_= '\gamma';
+ end
+ errorbar(x_,y_,e_,'-ok',...
+ 'LineWidth',1.5,...
'DisplayName',...
['$\gamma$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']);
hold on;
@@ -169,11 +161,12 @@ for i = 1:2
end
% ylabel('$\gamma$');
% yyaxis("right")
- errorbar(k(MODES),imag(omega(ik)),imag(err(ik)),'--k',...
+ errorbar(k(MODES),imag(omega(ik)),imag(err(ik)),'--ok',...
+ 'LineWidth',1.5,...
'DisplayName',...
['$\omega$, (',num2str(DATA.inputs.PMAX),',',num2str(DATA.inputs.JMAX),')']);
hold on;
- ylabel('$\gamma,\omega$');
+ ylabel(['$',lb_,',\omega$']);
xlabel(['$',kstr,'\rho_s$']);
title('Growth rates')
end
diff --git a/matlab/profiler.m b/matlab/profiler.m
index 43a3713a1c24cdb63d077c38b9d0e56dc72311fb..817c939a84212915d54cd9d3498fc386dcfbaa9c 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -3,8 +3,8 @@ function [] = profiler(data)
% filename = sprintf([BASIC.RESDIR,'outputs_%.2d.h5'],00);
% outfilename = ['/misc/HeLaZ_outputs',filename(3:end)];
CPUTI=[]; DTSIM=[]; RHSTC=[]; POITC=[]; SAPTC=[]; COLTC=[];
-GRATC=[]; NADTC=[]; ADVTC=[]; GHOTC=[]; CLOTC=[]; CHKTC=[];
-DIATC=[]; STETC=[]; TS0TC=[];
+GRATC=[]; NADTC=[]; ADVTC=[]; GHOTC=[]; CLOTC=[]; CHKTC=[];
+DIATC=[]; EXBTC=[]; STETC=[]; TS0TC=[];
for i = 1:numel(data.outfilenames)
outfilename = data.outfilenames{i};
CPUTI = [ CPUTI; double(h5readatt(outfilename,'/data/input','cpu_time'))];
@@ -13,6 +13,11 @@ for i = 1:numel(data.outfilenames)
RHSTC = [ RHSTC; h5read(outfilename,'/profiler/Tc_rhs')];
POITC = [ POITC; h5read(outfilename,'/profiler/Tc_poisson')];
SAPTC = [ SAPTC; h5read(outfilename,'/profiler/Tc_Sapj')];
+ try
+ EXBTC = [ EXBTC; h5read(outfilename,'/profiler/Tc_ExBshear')];
+ catch
+ EXBTC = 0.*SAPTC;
+ end
COLTC = [ COLTC; h5read(outfilename,'/profiler/Tc_coll')];
GRATC = [ GRATC; h5read(outfilename,'/profiler/Tc_grad')];
NADTC = [ NADTC; h5read(outfilename,'/profiler/Tc_nadiab')];
@@ -26,14 +31,14 @@ for i = 1:numel(data.outfilenames)
end
CPUTI = CPUTI(end);
DTSIM = mean(DTSIM);
-N_T = 12;
+N_T = 13;
-missing_Tc = STETC - RHSTC - ADVTC - GHOTC - CLOTC ...
+missing_Tc = STETC - RHSTC - ADVTC - GHOTC - CLOTC - EXBTC ...
-COLTC - POITC - SAPTC - CHKTC - DIATC - GRATC - NADTC;
total_Tc = STETC;
TIME_PER_FCT = [diff(RHSTC); diff(ADVTC); diff(GHOTC);...
- diff(CLOTC); diff(COLTC); diff(POITC); diff(SAPTC); ...
+ diff(CLOTC); diff(COLTC); diff(POITC); diff(SAPTC); diff(EXBTC); ...
diff(CHKTC); diff(DIATC); diff(GRATC); diff(NADTC); diff(missing_Tc)];
TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/N_T,N_T]);
@@ -47,6 +52,7 @@ clos_Ta = mean(diff(CLOTC));
coll_Ta = mean(diff(COLTC));
poisson_Ta = mean(diff(POITC));
Sapj_Ta = mean(diff(SAPTC));
+ExB_Ta = mean(diff(EXBTC));
checkfield_Ta = mean(diff(CHKTC));
grad_Ta = mean(diff(GRATC));
nadiab_Ta = mean(diff(NADTC));
@@ -61,14 +67,15 @@ names = {...
'Capj';
'Pois';
'Sapj';
+ 'ExBs';
'Chck';
'Diag';
'Grad';
'napj';
'Miss';
};
-Ts_A = [rhs_Ta adv_field_Ta ghost_Ta clos_Ta coll_Ta...
- poisson_Ta Sapj_Ta checkfield_Ta diag_Ta grad_Ta nadiab_Ta miss_Ta];
+Ts_A = [rhs_Ta adv_field_Ta ghost_Ta clos_Ta coll_Ta poisson_Ta...
+ Sapj_Ta ExB_Ta checkfield_Ta diag_Ta grad_Ta nadiab_Ta miss_Ta];
NSTEP_PER_SAMP= mean(diff(TS0TC))/DTSIM;
%% Plots
diff --git a/matlab/setup.m b/matlab/setup.m
index a5b6f420fb80ca48b15bc693a96a61445b9beb9f..89d826810aec2bfb4792c01c2937a23d9799909b 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -38,6 +38,7 @@ if ADIAB_E; MODEL.ADIAB_E = '.true.'; else; MODEL.ADIAB_E = '.false.';end;
if ADIAB_I; MODEL.ADIAB_I = '.true.'; else; MODEL.ADIAB_I = '.false.';end;
if MHD_PD; MODEL.MHD_PD = '.true.'; else; MODEL.MHD_PD = '.false.';end;
MODEL.beta = BETA;
+MODEL.ExBrate = EXBRATE;
MODEL.mu_x = MU_X;
MODEL.mu_y = MU_Y;
MODEL.N_HD = N_HD;
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index abba1afe5ec8545bd669767c9aef5661cfe3a3d1..0c0dc4e9822a405e99bbf15436ac40e850cb37f4 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -72,6 +72,7 @@ fprintf(fid,[' k_gB = ', num2str(MODEL.k_gB),'\n']);
fprintf(fid,[' k_cB = ', num2str(MODEL.k_cB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
fprintf(fid,[' beta = ', num2str(MODEL.beta),'\n']);
+fprintf(fid,[' ExBrate = ', num2str(MODEL.ExBrate),'\n']);
fprintf(fid,[' ADIAB_E = ', MODEL.ADIAB_E,'\n']);
fprintf(fid,[' ADIAB_I = ', MODEL.ADIAB_I,'\n']);
fprintf(fid,[' tau_i = ', num2str(MODEL.tau_i),'\n']);
diff --git a/src/ExB_shear_flow_mod.F90 b/src/ExB_shear_flow_mod.F90
new file mode 100644
index 0000000000000000000000000000000000000000..649f99664d72a5f4fee7e762ebd8f9b4eddef064
--- /dev/null
+++ b/src/ExB_shear_flow_mod.F90
@@ -0,0 +1,111 @@
+MODULE ExB_shear_flow
+ ! This module contains the necessary tools to implement ExB shearing flow effects.
+ ! The algorithm is taken from the presentation of Hammett et al. 2006 (APS) and
+ ! it the one used in GS2.
+ USE prec_const, ONLY: xp
+
+ IMPLICIT NONE
+ ! Variables
+ REAL(xp), PUBLIC, PROTECTED :: gamma_E = 0._xp ! ExB background shearing rate \gamma_E
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: sky_ExB ! shift of the kx modes, kx* = kx + s(ky)
+ INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: jump_ExB ! jump to do to shift the kx grids
+ LOGICAL, DIMENSION(:), ALLOCATABLE, PUBLIC, PROTECTED :: shiftnow_ExB ! Indicates if there is a line to shift
+
+ ! Routines
+ PUBLIC :: Setup_ExB_shear_flow, Apply_ExB_shear_flow, Update_ExB_shear_flow
+
+CONTAINS
+
+ ! Setup the variables for the ExB shear
+ SUBROUTINE Setup_ExB_shear_flow
+ USE grid, ONLY : local_nky
+ IMPLICIT NONE
+
+ ! Setup the ExB shift
+ ALLOCATE(sky_ExB(local_nky))
+ sky_ExB = 0._xp
+
+ ! Setup the jump array and shifting flag
+ ALLOCATE(jump_ExB(local_nky))
+ jump_ExB = 0
+ ALLOCATE(shiftnow_ExB(local_nky))
+ shiftnow_ExB = .FALSE.
+ END SUBROUTINE Setup_ExB_shear_flow
+
+ ! Update according to the current ExB shear value
+ ! -the grids
+ ! -the spatial operators
+ ! -the fields by imposing a shift on kx
+ SUBROUTINE Apply_ExB_shear_flow
+ USE basic, ONLY: chrono_ExBs, start_chrono, stop_chrono
+ USE grid, ONLY: local_nky, kxarray, update_grids, &
+ total_nkx, deltakx, kx_min, kx_max
+ USE prec_const, ONLY: PI
+ USE geometry, ONLY: gxx,gxy,inv_hatB2, evaluate_magn_curv
+ USE numerics, ONLY: evaluate_EM_op, evaluate_kernels
+ USE fields, ONLY: moments, phi, psi
+ IMPLICIT NONE
+ ! local var
+ INTEGER :: iky, ikx, ikx_s
+
+ CALL start_chrono(chrono_ExBs)
+
+ ! shift all fields and correct the shift value
+ DO iky = 1,local_Nky
+ IF(shiftnow_ExB(iky)) THEN
+ print*, "SHIFT ARRAYS"
+ ! shift all fields
+ DO ikx = 1,total_nkx
+ ikx_s = ikx + jump_ExB(iky)
+ ! We test if the shifted modes are still in contained in our resolution
+ ! IF( (kxarray(iky,ikx)-sky_ExB(iky) .GE. kx_min) .AND. (kxarray(iky,ikx)-sky_ExB(iky) .LE. kx_max) ) THEN
+ IF ( ((ikx_s .GT. 0 ) .AND. (ikx_s .LE. total_nkx )) .AND. &
+ (((ikx .LE. (total_nkx/2+1)) .AND. (ikx_s .LE. (total_nkx/2+1))) .OR. &
+ ((ikx .GT. (total_nkx/2+1)) .AND. (ikx_s .GT. (total_nkx/2+1)))) ) THEN
+ moments(:,:,:,iky,ikx,:,:) = moments(:,:,:,iky,ikx_s,:,:)
+ phi(iky,ikx,:) = phi(iky,ikx_s,:)
+ psi(iky,ikx,:) = psi(iky,ikx_s,:)
+ ELSE ! if it is not, it is lost (~dissipation for high modes)
+ moments(:,:,:,iky,ikx,:,:) = 0._xp
+ phi(iky,ikx,:) = 0._xp
+ psi(iky,ikx,:) = 0._xp
+ ENDIF
+ ENDDO
+ ! correct the shift value s(ky) for this row
+ sky_ExB(iky) = sky_ExB(iky) - jump_ExB(iky)*deltakx
+ ENDIF
+ ENDDO
+ ! After shifting and correction, we update the operators and grids
+ ! update the grids
+ CALL update_grids(sky_ExB,gxx,gxy,inv_hatB2)
+
+ ! update the EM op., the kernels and the curvature op.
+ ! CALL evaluate_kernels
+ ! CALL evaluate_EM_op
+ ! CALL evaluate_magn_curv
+
+ CALL stop_chrono(chrono_ExBs)
+ END SUBROUTINE Apply_ExB_shear_flow
+
+ ! update the ExB shear value for the next time step
+ SUBROUTINE Update_ExB_shear_flow
+ USE basic, ONLY: dt, chrono_ExBs, start_chrono, stop_chrono
+ USE grid, ONLY: local_nky, kyarray, inv_dkx
+ USE model, ONLY: ExBrate
+ IMPLICIT NONE
+ ! local var
+ INTEGER :: iky
+ CALL start_chrono(chrono_ExBs)
+ ! update the ExB shift, jumps and flags
+ shiftnow_ExB = .FALSE.
+ DO iky = 1,local_Nky
+ sky_ExB(iky) = sky_ExB(iky) - kyarray(iky)*ExBrate*dt
+ jump_ExB(iky) = NINT(sky_ExB(iky)*inv_dkx)
+ ! If the jump is 1 or more for a given ky, we flag the index
+ ! in shiftnow_ExB and will use it in Shift_fields to avoid
+ ! zero-shiftings that may be majoritary.
+ shiftnow_ExB(iky) = (abs(jump_ExB(iky)) .GT. 0)
+ ENDDO
+ CALL stop_chrono(chrono_ExBs)
+ END SUBROUTINE Update_ExB_shear_flow
+END MODULE ExB_shear_flow
\ No newline at end of file
diff --git a/src/auxval.F90 b/src/auxval.F90
index dd4f4a16d6b9cb27c3de193db4c965d50478fc59..6828d2aa879b778088aee66136bff0320e95d639 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -5,13 +5,15 @@ subroutine auxval
USE grid
USE array
USE model
- USE fourier, ONLY: init_grid_distr_and_plans
+ USE fourier, ONLY: init_grid_distr_and_plans
use prec_const
USE numerics
USE geometry
- USE closure, ONLY: set_closure_model, hierarchy_closure
- USE parallel, ONLY: init_parallel_var, my_id, num_procs, num_procs_p, num_procs_z, num_procs_ky, rank_p, rank_ky, rank_z
- USE processing, ONLY: init_process
+ USE closure, ONLY: set_closure_model, hierarchy_closure
+ USE parallel, ONLY: init_parallel_var, my_id, num_procs, &
+ num_procs_p, num_procs_z, num_procs_ky, rank_p, rank_ky, rank_z
+ USE processing, ONLY: init_process
+ USE ExB_shear_flow, ONLY: Setup_ExB_shear_flow
#ifdef TEST_SVD
USE CLA, ONLY: init_CLA
#endif
@@ -34,7 +36,7 @@ subroutine auxval
! precompute the kernels
CALL evaluate_kernels
! compute inverse of poisson and ampere operators
- CALL evaluate_EM_op
+ CALL evaluate_EM_op
! precompute the Laguerre nonlin product coeffs
IF ( LINEARITY .NE. 'linear' ) &
CALL build_dnjs_table
@@ -42,6 +44,8 @@ subroutine auxval
CALL build_dv4Hp_table
! set the closure scheme in use
CALL set_closure_model
+ ! Setup ExB shear variables
+ CALL Setup_ExB_shear_flow
#ifdef TEST_SVD
! If we want to test SVD decomposition etc.
CALL init_CLA(local_nky,local_np*local_nj)
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 0c45dd452978d58a64048cfb65693e292b7efb19..fb34690a6aad560d5cf31e7e4692f8ab9643a21c 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -34,7 +34,8 @@ MODULE basic
end type chrono
! Define the chronos for each relevant routines
type(chrono), PUBLIC, PROTECTED :: chrono_runt, chrono_mrhs, chrono_advf, chrono_pois, chrono_sapj,&
- chrono_diag, chrono_chck, chrono_step, chrono_clos, chrono_ghst, chrono_coll, chrono_napj, chrono_grad
+ chrono_diag, chrono_chck, chrono_step, chrono_clos, chrono_ghst, chrono_coll, chrono_napj, &
+ chrono_grad, chrono_ExBs
#ifdef TEST_SVD
! A chrono for SVD tests
type(chrono), PUBLIC, PROTECTED :: chrono_CLA
@@ -86,6 +87,7 @@ CONTAINS
chrono_chck%ttot = 0
chrono_diag%ttot = 0
chrono_step%ttot = 0
+ chrono_ExBs%ttot = 0
#ifdef TEST_SVD
chrono_CLA%ttot = 0
#endif
diff --git a/src/circular_mod.F90 b/src/circular_mod.F90
index cff23018af3d1453cd43a2f503691a6b5f6a1675..8a01fb7cbce6a4b3177fb3dac977a69499f9a882 100644
--- a/src/circular_mod.F90
+++ b/src/circular_mod.F90
@@ -8,7 +8,6 @@ MODULE circular
USE grid, ONLY: local_Nky, local_Nkx, local_nz, ngz, nzgrid, kyarray, kxarray, zarray, total_nz, local_nz_offset, iodd, ieven
! use discretization
USE lagrange_interpolation
- USE model
implicit none
public:: get_circ
diff --git a/src/collision_mod.F90 b/src/collision_mod.F90
index 5d842e4db52b432dcd5ff4c8213f8f74cc2b1b05..66df93ec942fb430b242f32df273819ea1825177 100644
--- a/src/collision_mod.F90
+++ b/src/collision_mod.F90
@@ -116,14 +116,14 @@ CONTAINS
!******************************************************************************!
SUBROUTINE compute_lenard_bernstein
USE grid, ONLY: ias,iae, ips,ipe, ijs,ije, parray, jarray, &
- ikys,ikye, ikxs,ikxe, izs,ize, kparray, ngp, ngj, ngz
+ ikys,ikye, ikxs,ikxe, izs,ize, kp2array, ngp, ngj, ngz
USE species, ONLY: sigma2_tau_o2, nu_ab
USE time_integration, ONLY: updatetlevel
USE array, ONLY: Capj
USE fields, ONLY: moments
IMPLICIT NONE
COMPLEX(xp) :: TColl_
- REAL(xp) :: j_xp, p_xp, ba_2, kp
+ REAL(xp) :: j_xp, p_xp, ba_2, kp2
INTEGER :: iz,ikx,iky,ij,ip,ia,eo,ipi,iji,izi
DO iz = izs,ize; izi = iz + ngz/2
DO ikx = ikxs, ikxe;
@@ -135,8 +135,8 @@ CONTAINS
eo = MODULO(parray(ipi),2)
p_xp = REAL(parray(ipi),xp)
j_xp = REAL(jarray(iji),xp)
- kp = kparray(iky,ikx,izi,eo)
- ba_2 = kp**2 * sigma2_tau_o2(ia) ! this is (ba/2)^2
+ kp2 = kp2array(iky,ikx,izi,eo)
+ ba_2 = kp2 * sigma2_tau_o2(ia) ! this is (ba/2)^2
!** Assembling collison operator **
! -nuee (p + 2j) Nepj
@@ -160,7 +160,7 @@ CONTAINS
USE grid, ONLY: local_na, local_np, local_nj, parray, jarray, ngp, ngj, &
ip0, ip2, ij0, ij1, ieven, &
local_nky, local_nkx, local_nz, ngz,&
- kparray
+ kp2array
USE species, ONLY: nu_ab, tau, q_tau
USE time_integration, ONLY: updatetlevel
USE array, ONLY: Capj
@@ -169,12 +169,12 @@ CONTAINS
IMPLICIT NONE
COMPLEX(xp) :: Tmp
INTEGER :: iz,ikx,iky,ij,ip,ia, ipi,iji,izi
- REAL(xp) :: j_xp, p_xp, kp_xp
+ REAL(xp) :: j_xp, p_xp, kp2_xp
DO iz = 1,local_nz
izi = iz + ngz/2
DO ikx = 1,local_nkx
DO iky = 1,local_nky
- kp_xp = kparray(iky,ikx,izi,ieven)
+ kp2_xp = kp2array(iky,ikx,izi,ieven)
DO ij = 1,local_nj
iji = ij + ngj/2
DO ip = 1,local_np
@@ -185,17 +185,17 @@ CONTAINS
j_xp = REAL(jarray(iji),xp)
!** Assembling collison operator **
IF( (p_xp .EQ. 0._xp) .AND. (j_xp .EQ. 0._xp)) THEN !Ca00
- Tmp = tau(ia)**2 * kp_xp**4*(&
+ Tmp = tau(ia)**2 * kp2_xp**2*(&
67_xp/120_xp *moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
+67_xp*SQRT2/240_xp*moments(ia,ip2,ij0,iky,ikx,izi,updatetlevel)&
-67_xp/240_xp *moments(ia,ip0,ij1,iky,ikx,izi,updatetlevel)&
-3_xp/10_xp *q_tau(ia)*phi(iky,ikx,izi))
ELSEIF( (p_xp .EQ. 2._xp) .AND. (j_xp .EQ. 0._xp)) THEN ! Ca20
- Tmp = tau(ia) * kp_xp**2*(&
+ Tmp = tau(ia) * kp2_xp*(&
-SQRT2*twothird*moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
-2._xp*twothird*moments(ia,ip2,ij0,iky,ikx,izi,updatetlevel))
ELSEIF( (p_xp .EQ. 0._xp) .AND. (j_xp .EQ. 1._xp)) THEN ! Ca01
- Tmp = tau(ia) * kp_xp**2*(&
+ Tmp = tau(ia) * kp2_xp*(&
2._xp*twothird*moments(ia,ip0,ij0,iky,ikx,izi,updatetlevel)&
-2._xp*twothird*moments(ia,ip0,ij1,iky,ikx,izi,updatetlevel))
ELSE
diff --git a/src/diagnose.F90 b/src/diagnose.F90
index 327ad696636784faad221a7be612e7fb08e760df..8807f076fa17bd74169421b4f7907081bc8b3e26 100644
--- a/src/diagnose.F90
+++ b/src/diagnose.F90
@@ -134,11 +134,12 @@ SUBROUTINE diagnose_full(kstep)
CALL creatd(fidres, 0, dims, "/profiler/Tc_coll", "cumulative collision computation time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_grad", "cumulative grad computation time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_nadiab", "cumulative nadiab moments computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_ExBshear", "cumulative ExB shear time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_adv_field", "cumulative adv. fields computation time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_ghost", "cumulative communication time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_clos", "cumulative closure computation time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_checkfield", "cumulative checkfield computation time")
- CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative sym computation time")
+ CALL creatd(fidres, 0, dims, "/profiler/Tc_diag", "cumulative diagnostic time")
CALL creatd(fidres, 0, dims, "/profiler/Tc_step", "cumulative total step computation time")
CALL creatd(fidres, 0, dims, "/profiler/time", "current simulation time")
#ifdef TEST_SVD
@@ -351,6 +352,7 @@ SUBROUTINE diagnose_0d
CALL append(fidres, "/profiler/Tc_diag", REAL(chrono_diag%ttot,dp),ionode=0)
CALL append(fidres, "/profiler/Tc_grad", REAL(chrono_grad%ttot,dp),ionode=0)
CALL append(fidres, "/profiler/Tc_nadiab", REAL(chrono_napj%ttot,dp),ionode=0)
+ CALL append(fidres, "/profiler/Tc_ExBshear", REAL(chrono_ExBs%ttot,dp),ionode=0)
CALL append(fidres, "/profiler/Tc_step", REAL(chrono_step%ttot,dp),ionode=0)
#ifdef TEST_SVD
CALL append(fidres, "/profiler/Tc_CLA", REAL(chrono_CLA%ttot,dp),ionode=0)
diff --git a/src/geometry_mod.F90 b/src/geometry_mod.F90
index 4d18d1c9b4ec00a6daad8f83db11a66922ed3495..91b83db374904c8290a598ac1009bfff3592811b 100644
--- a/src/geometry_mod.F90
+++ b/src/geometry_mod.F90
@@ -51,7 +51,7 @@ implicit none
! derivatives of magnetic field strength
REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: dBdx, dBdy, dBdz, dlnBdz
! Relative magnetic field strength
- REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB
+ REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatB, inv_hatB2 ! normalized bckg magnetic gradient amplitude and its inv squared
! Relative strength of major radius
REAL(xp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hatR, hatZ
! Some geometrical coefficients
@@ -71,7 +71,7 @@ implicit none
! Functions
PUBLIC :: geometry_readinputs, geometry_outputinputs,&
- eval_magnetic_geometry, set_ikx_zBC_map
+ eval_magnetic_geometry, set_ikx_zBC_map, evaluate_magn_curv
CONTAINS
@@ -99,7 +99,7 @@ CONTAINS
subroutine eval_magnetic_geometry
USE grid, ONLY: total_nky, total_nz, local_nkx, local_nky, local_nz, ngz,&
- kxarray, kyarray, set_kparray, nzgrid, zweights_SR, ieven, set_kxgrid
+ set_kparray, nzgrid, zweights_SR, ieven, set_kxgrid
USE basic, ONLY: speak
USE miller, ONLY: set_miller_parameters, get_miller
USE circular, ONLY: get_circ
@@ -108,10 +108,8 @@ CONTAINS
USE species, ONLY: Ptot
! evalute metrix, elementwo_third_kpmaxts, jacobian and gradient
implicit none
- REAL(xp) :: kx,ky
COMPLEX(xp), DIMENSION(local_nz) :: integrant
- real(xp) :: Cx, Cy, g_xz, g_yz, g_zz
- INTEGER :: eo,iz,iky,ikx
+
! Allocate arrays
CALL geometry_allocate_mem(local_nky,local_nkx,local_nz,ngz,nzgrid)
@@ -164,15 +162,39 @@ CONTAINS
ERROR STOP '>> ERROR << geometry not recognized!!'
END SELECT
ENDIF
+ ! inv squared of the magnetic gradient bckg amplitude (for fast kperp update)
+ inv_hatB2 = 1/hatB/hatB
! Reset kx grid (to account for Cyq0_x0 factor)
CALL set_kxgrid(shear,Npol,Cyq0_x0,LINEARITY,N_HD)
!
! Evaluate perpendicular wavenumber
! k_\perp^2 = g^{xx} k_x^2 + 2 g^{xy}k_x k_y + k_y^2 g^{yy}
! normalized to rhos_
- CALL set_kparray(gxx,gxy,gyy,hatB)
+ CALL set_kparray(gxx,gxy,gyy,inv_hatB2)
+ ! Evaluate magnetic curvature operator Ckxky
+ CALL evaluate_magn_curv
+ ! coefficient in the front of parallel derivative
+ gradz_coeff = 1._xp /(jacobian*hatB)
+ ! d/dz(ln B) to correspond to the formulation in Hoffmann et al. 2023
+ dlnBdz = dBdz/hatB
+ !
+ ! set the mapping for parallel boundary conditions
+ CALL set_ikx_zBC_map
+ !
+ ! Compute the inverse z integrated Jacobian (useful for flux averaging)
+ integrant = Jacobian((1+ngz/2):(local_nz+ngz/2),ieven) ! Convert into complex array
+ CALL simpson_rule_z(local_nz,zweights_SR,integrant,iInt_Jacobian)
+ iInt_Jacobian = 1._xp/iInt_Jacobian ! reverse it
+ END SUBROUTINE eval_magnetic_geometry
+
+ SUBROUTINE evaluate_magn_curv
+ USE grid, ONLY: local_nkx, local_nky, local_nz, ngz,&
+ kxarray, kyarray, nzgrid
+ IMPLICIT NONE
+ REAL(xp) :: kx,ky
+ real(xp) :: Cx, Cy, g_xz, g_yz, g_zz
+ INTEGER :: eo,iz,iky,ikx
DO eo = 1,nzgrid
- ! Curvature operator (Frei et al. 2022 eq 2.15)
DO iz = 1,local_nz+ngz
! !covariant metric
g_xz = jacobian(iz,eo)**2*(gxy(iz,eo)*gyz(iz,eo)-gyy(iz,eo)*gxz(iz,eo))
@@ -186,25 +208,13 @@ CONTAINS
DO iky = 1,local_nky
ky = kyarray(iky)
DO ikx= 1,local_nkx
- kx = kxarray(ikx)
+ kx = kxarray(iky,ikx)
Ckxky(iky, ikx, iz,eo) = (Cx*kx + Cy*ky)/C_xy
ENDDO
ENDDO
- ! coefficient in the front of parallel derivative
- gradz_coeff(iz,eo) = 1._xp /(jacobian(iz,eo)*hatB(iz,eo))
- ! d/dz(ln B) to correspond to the formulation in paper 2023
- dlnBdz(iz,eo) = dBdz(iz,eo)/hatB(iz,eo)
ENDDO
ENDDO
- !
- ! set the mapping for parallel boundary conditions
- CALL set_ikx_zBC_map
- !
- ! Compute the inverse z integrated Jacobian (useful for flux averaging)
- integrant = Jacobian((1+ngz/2):(local_nz+ngz/2),ieven) ! Convert into complex array
- CALL simpson_rule_z(local_nz,zweights_SR,integrant,iInt_Jacobian)
- iInt_Jacobian = 1._xp/iInt_Jacobian ! reverse it
- END SUBROUTINE eval_magnetic_geometry
+ END SUBROUTINE
!
!--------------------------------------------------------------------------------
!
@@ -382,10 +392,10 @@ CONTAINS
! Check if it has to be connected to the otherside of the kx grid
if(ikx_zBC_L(iky,ikx) .LE. 0) ikx_zBC_L(iky,ikx) = ikx_zBC_L(iky,ikx) + total_nkx
! Check if it points out of the kx domain
- ! print*, kxarray(ikx), shift, kx_min
- IF( (kxarray(ikx) - shift) .LT. kx_min ) THEN ! outside of the frequ domain
- ! print*,( ((kxarray(ikx) - shift)-kx_min) .LT. EPSILON(kx_min) )
- ! print*, kxarray(ikx)-kx_min, EPSILON(kx_min)
+ ! print*, kxarray(iky,ikx), shift, kx_min
+ IF( (kxarray(iky,ikx) - shift) .LT. kx_min ) THEN ! outside of the frequ domain
+ ! print*,( ((kxarray(iky,ikx) - shift)-kx_min) .LT. EPSILON(kx_min) )
+ ! print*, kxarray(iky,ikx)-kx_min, EPSILON(kx_min)
! IF( (ikx-mn_y*Nexc) .LT. 1 ) THEN
SELECT CASE(parallel_bc)
CASE ('dirichlet','disconnected') ! connected to 0
@@ -420,7 +430,7 @@ CONTAINS
! Check if it has to be connected to the otherside of the kx grid
if(ikx_zBC_R(iky,ikx) .GT. total_nkx) ikx_zBC_R(iky,ikx) = ikx_zBC_R(iky,ikx) - total_nkx
! Check if it points out of the kx domain
- IF( (kxarray(ikx) + shift) .GT. kx_max ) THEN ! outside of the frequ domain
+ IF( (kxarray(iky,ikx) + shift) .GT. kx_max ) THEN ! outside of the frequ domain
! IF( (ikx+mn_y*Nexc) .GT. total_nkx ) THEN ! outside of the frequ domain
SELECT CASE(parallel_bc)
CASE ('dirichlet','disconnected') ! connected to 0
@@ -428,7 +438,7 @@ CONTAINS
CASE ('periodic') ! connected to itself as for shearless
ikx_zBC_R(iky,ikx) = ikx
CASE ('cyclic')
- ! write(*,*) 'check',ikx,iky, kxarray(ikx) + shift, '>', kx_max
+ ! write(*,*) 'check',ikx,iky, kxarray(iky,ikx) + shift, '>', kx_max
ikx_zBC_R(iky,ikx) = MODULO(ikx_zBC_R(iky,ikx)-1,total_nkx)+1
CASE DEFAULT
ERROR STOP "The parallel BC is not recognized (dirichlet, periodic, cyclic or disconnected)"
@@ -517,6 +527,7 @@ END SUBROUTINE set_ikx_zBC_map
ALLOCATE( dBdy(local_nz+ngz,nzgrid))
ALLOCATE( dBdz(local_nz+ngz,nzgrid))
ALLOCATE( dlnBdz(local_nz+ngz,nzgrid))
+ ALLOCATE( inv_hatB2(local_nz+ngz,nzgrid))
ALLOCATE( hatB(local_nz+ngz,nzgrid))
ALLOCATE( hatR(local_nz+ngz,nzgrid))
ALLOCATE( hatZ(local_nz+ngz,nzgrid))
diff --git a/src/grid_mod.F90 b/src/grid_mod.F90
index dbeda1d1ba8c321c56bb4448322558320a6e7c6c..a3e8a7a5fc6e421a0efc23a66b5a8e226d3be422 100644
--- a/src/grid_mod.F90
+++ b/src/grid_mod.F90
@@ -25,11 +25,13 @@ MODULE grid
! Grid arrays
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: parray, parray_full
INTEGER, DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: jarray, jarray_full
- REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray, kxarray_full
+ REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray ! ExB shear makes it ky dependant
+ REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kxarray_full
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kyarray, kyarray_full
REAL(xp), DIMENSION(:,:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: zarray_full
REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kparray !kperp
+ REAL(xp), DIMENSION(:,:,:,:), ALLOCATABLE, PUBLIC,PROTECTED :: kp2array !kperp^2
! Kronecker delta for p=0, p=1, p=2, j=0, j=1
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kroneck_p0, kroneck_p1, kroneck_p2, kroneck_p3
REAL(xp), DIMENSION(:), ALLOCATABLE, PUBLIC,PROTECTED :: kroneck_j0, kroneck_j1
@@ -71,7 +73,7 @@ MODULE grid
integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nkx_ptr, local_nky_ptr
integer(C_INTPTR_T), PUBLIC,PROTECTED :: local_nkx_ptr_offset, local_nky_ptr_offset
! Grid spacing and limits
- REAL(xp), PUBLIC, PROTECTED :: deltap, deltaz, inv_deltaz
+ REAL(xp), PUBLIC, PROTECTED :: deltap, deltaz, inv_deltaz, inv_dkx
REAL(xp), PUBLIC, PROTECTED :: deltakx, deltaky, kx_max, ky_max, kx_min, ky_min!, kp_max
INTEGER , PUBLIC, PROTECTED :: local_pmin, local_pmax
INTEGER , PUBLIC, PROTECTED :: local_jmin, local_jmax
@@ -112,6 +114,7 @@ MODULE grid
PUBLIC :: set_grids, set_kxgrid, set_kparray
PUBLIC :: grid_readinputs, grid_outputinputs
PUBLIC :: bar
+ PUBLIC :: update_grids ! Update the kx and kperp grids for ExB shear flow
! Precomputations
real(xp), PUBLIC, PROTECTED :: pmax_xp, jmax_xp
@@ -214,7 +217,7 @@ CONTAINS
local_nkx_ptr = ikxe - ikxs + 1
local_nkx = ikxe - ikxs + 1
local_nkx_offset = ikxs - 1
- ALLOCATE(kxarray(local_nkx))
+ ALLOCATE(kxarray(local_nky,local_Nkx))
ALLOCATE(kxarray_full(total_nkx))
ALLOCATE(AA_x(local_nkx))
!!---------------- PARALLEL Z GRID (parallelized)
@@ -263,6 +266,7 @@ CONTAINS
ENDDO
!!---------------- Kperp grid
CALL allocate_array( kparray, 1,local_nky, 1,local_nkx, 1,local_nz+ngz, 1,nzgrid)
+ CALL allocate_array(kp2array, 1,local_nky, 1,local_nkx, 1,local_nz+ngz, 1,nzgrid)
END SUBROUTINE
!!!! GRID REPRESENTATION
! We define the grids that contain ghosts (p,j,z) with indexing from 1 to Nlocal + nghost
@@ -474,7 +478,7 @@ CONTAINS
REAL(xp), INTENT(IN) :: shear, Npol, Cyq0_x0
CHARACTER(len=*), INTENT(IN) ::LINEARITY
INTEGER, INTENT(IN) :: N_HD
- INTEGER :: ikx
+ INTEGER :: ikx, iky
REAL(xp):: Lx_adapted
IF(shear .GT. 0) THEN
IF(my_id.EQ.0) write(*,*) 'Magnetic shear detected: set up sheared kx grid..'
@@ -488,7 +492,8 @@ CONTAINS
! x length is adapted
Lx = Lx_adapted*Nexc
ENDIF
- deltakx = 2._xp*PI/Lx
+ deltakx = 2._xp*PI/Lx
+ inv_dkx = 1._xp/deltakx
IF(MODULO(total_nkx,2) .EQ. 0) THEN ! Even number of kx (-2 -1 0 1 2 3)
! Creating a grid ordered as dk*(0 1 2 3 -2 -1)
DO ikx = 1,total_nkx
@@ -499,18 +504,20 @@ CONTAINS
kx_min = MINVAL(kxarray_full)!-kx_max+deltakx
! Set local grid (not parallelized so same as full one)
local_kxmax = 0._xp
- DO ikx = 1,local_nkx
- kxarray(ikx) = kxarray_full(ikx+local_nkx_offset)
- ! Finding kx=0
- IF (kxarray(ikx) .EQ. 0) THEN
- ikx0 = ikx
- contains_kx0 = .true.
- ENDIF
- ! Finding local kxmax
- IF (ABS(kxarray(ikx)) .GT. local_kxmax) THEN
- local_kxmax = ABS(kxarray(ikx))
- ikx_max = ikx
- ENDIF
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ kxarray(iky,ikx) = kxarray_full(ikx+local_nkx_offset)
+ ! Finding kx=0
+ IF (kxarray(1,ikx) .EQ. 0) THEN
+ ikx0 = ikx
+ contains_kx0 = .true.
+ ENDIF
+ ! Finding local kxmax
+ IF (ABS(kxarray(iky,ikx)) .GT. local_kxmax) THEN
+ local_kxmax = ABS(kxarray(iky,ikx))
+ ikx_max = ikx
+ ENDIF
+ END DO
END DO
ELSE ! Odd number of kx (-2 -1 0 1 2)
ERROR STOP "Gyacomo is safer with an even Kx number"
@@ -518,13 +525,15 @@ CONTAINS
! Orszag 2/3 filter
two_third_kxmax = 2._xp/3._xp*kx_max;
! Antialiasing filter
- DO ikx = 1,local_nkx
- IF ( ((kxarray(ikx) .GT. -two_third_kxmax) .AND. &
- (kxarray(ikx) .LT. two_third_kxmax)) .OR. (LINEARITY .EQ. 'linear')) THEN
- AA_x(ikx) = 1._xp;
- ELSE
- AA_x(ikx) = 0._xp;
- ENDIF
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ IF ( ((kxarray(iky,ikx) .GT. -two_third_kxmax) .AND. &
+ (kxarray(iky,ikx) .LT. two_third_kxmax)) .OR. (LINEARITY .EQ. 'linear')) THEN
+ AA_x(ikx) = 1._xp;
+ ELSE
+ AA_x(ikx) = 0._xp;
+ ENDIF
+ END DO
END DO
! For hyperdiffusion
IF(LINEARITY.EQ.'linear') THEN
@@ -572,19 +581,6 @@ CONTAINS
! Find local extrema
local_zmax(eo) = zarray(local_nz+ngz/2,eo)
local_zmin(eo) = zarray(1+ngz/2,eo)
- ! Periodic z \in (-pi pi-dz)
- ! DO ig = 1,ngz/2 ! first ghost cells
- ! iglob = ig+local_nz_offset-ngz/2
- ! IF (iglob .LE. 0) &
- ! iglob = iglob + total_nz
- ! zarray(ig,eo) = zarray_full(iglob)
- ! ENDDO
- ! DO ig = local_nz+ngz/2,local_nz+ngz ! last ghost cells
- ! iglob = ig+local_nz_offset-ngz/2
- ! IF (iglob .GT. total_nz) &
- ! iglob = iglob - total_nz
- ! zarray(ig,eo) = zarray_full(iglob)
- ! ENDDO
! continue in z<pi and z>pi
DO ig = 1,ngz/2
zarray(local_nz+ngz/2+ig,eo) = zarray(local_nz+ngz/2,eo) + ig*deltaz
@@ -611,9 +607,9 @@ CONTAINS
ENDIF
END SUBROUTINE set_zgrid
- SUBROUTINE set_kparray(gxx, gxy, gyy,hatB)
+ SUBROUTINE set_kparray(gxx, gxy, gyy, inv_hatB2)
IMPLICIT NONE
- REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,hatB
+ REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,gyy,inv_hatB2
INTEGER :: eo,iz,iky,ikx
REAL(xp) :: kx, ky
DO eo = 1,nzgrid
@@ -621,11 +617,12 @@ CONTAINS
DO iky = 1,local_nky
ky = kyarray(iky)
DO ikx = 1,local_nkx
- kx = kxarray(ikx)
+ kx = kxarray(iky,ikx)
! there is a factor 1/B from the normalization; important to match GENE
! this factor comes from $b_a$ argument in the Bessel. Kperp is not used otherwise.
- kparray(iky, ikx, iz, eo) = &
- SQRT( gxx(iz,eo)*kx**2 + 2._xp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)/ hatB(iz,eo)
+ kp2array(iky, ikx, iz, eo) = &
+ (gxx(iz,eo)*kx**2 + 2._xp*gxy(iz,eo)*kx*ky + gyy(iz,eo)*ky**2)*inv_hatB2(iz,eo)
+ kparray(iky, ikx, iz, eo) = SQRT(kp2array(iky, ikx, iz, eo))
ENDDO
ENDDO
ENDDO
@@ -633,6 +630,34 @@ CONTAINS
two_third_kpmax = 2._xp/3._xp * MAXVAL(kparray)
END SUBROUTINE
+ SUBROUTINE update_grids (sky,gxx,gxy,inv_hatB2)
+ IMPLICIT NONE
+ REAL(xp), DIMENSION(local_nky),INTENT(IN) :: sky ! ExB grid shift
+ REAL(xp), DIMENSION(local_nz+ngz,nzgrid), INTENT(IN) :: gxx,gxy,inv_hatB2
+ INTEGER :: eo,iz,iky,ikx
+ REAL(xp) :: kx, ky, skp2
+ ! Update the kx grid
+ DO iky = 1,local_nky
+ DO ikx = 1,local_nkx
+ kxarray(iky,ikx) = kxarray(iky,ikx) + sky(iky)
+ ENDDO
+ ENDDO
+ ! Update the kperp grid
+ DO eo = 1,nzgrid
+ DO iz = 1,local_nz+ngz
+ DO iky = 1,local_nky
+ ky = kyarray(iky)
+ DO ikx = 1,local_nkx
+ kx = kxarray(iky,ikx)
+ ! shift in kp obtained from kp2* = kp2 + skp2
+ skp2 = sky(iky)*inv_hatB2(iz,eo)*(gxx(iz,eo)*(2._xp*kx+sky(iky)) +gxy(iz,eo)*ky)
+ kp2array(iky, ikx, iz, eo) = kp2array(iky, ikx, iz, eo) + skp2
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDDO
+ END SUBROUTINE
+
SUBROUTINE grid_outputinputs(fid)
! Write the input parameters to the results_xx.h5 file
USE futils, ONLY: attach, creatd
diff --git a/src/inital.F90 b/src/inital.F90
index cfc12ecb46382ff7a583f9bc8488bebea9904733..8871714afef4976ee6bf96f470170dc870f0303c 100644
--- a/src/inital.F90
+++ b/src/inital.F90
@@ -308,10 +308,10 @@ SUBROUTINE init_phi_ppj
INTEGER :: ikx, iky, iz
amp = 1.0_xp
!**** ppj initialization *******************************************
- DO ikx=1,total_nkx
- kx = kxarray(ikx)
- DO iky=1,local_nky
- ky = kyarray(iky)
+ DO iky=1,local_nky
+ ky = kyarray(iky)
+ DO ikx=1,total_nkx
+ kx = kxarray(iky,ikx)
DO iz=1,local_nz+ngz
z = zarray(iz,ieven)
IF (ky .NE. 0) THEN
@@ -381,7 +381,7 @@ SUBROUTINE initialize_blob
DO iz=1+ngz/2,local_nz+ngz/2
z = zarray(iz,ieven)
DO ikx=1,total_nkx
- kx = kxarray(ikx) + z*shear*ky
+ kx = kxarray(iky,ikx) + z*shear*ky
DO ip=1+ngp/2,local_np+ngp/2
p = parray(ip)
DO ij=1+ngj/2,local_nj+ngj/2
@@ -437,10 +437,10 @@ SUBROUTINE init_ppj
DO ip=1,local_np+ngp
DO ij=1,local_nj+ngj
IF ( (ip .EQ. 1) .AND. (ij .EQ. 1) ) THEN
- DO ikx=1,total_nkx
- kx = kxarray(ikx)
- DO iky=1,local_nky
- ky = kyarray(iky)
+ DO iky=1,local_nky
+ ky = kyarray(iky)
+ DO ikx=1,total_nkx
+ kx = kxarray(iky,ikx)
DO iz=1,local_nz+ngz
z = zarray(iz,ieven)
IF (kx .EQ. 0) THEN
diff --git a/src/miller_mod.F90 b/src/miller_mod.F90
index 463770f7b4bca8f28146f6b4a03cee4f82bd5af4..ad06e9448a8ef77ca711b06859bdb292d3f8a307 100644
--- a/src/miller_mod.F90
+++ b/src/miller_mod.F90
@@ -9,7 +9,6 @@ MODULE miller
USE grid, ONLY: local_Nky, local_Nkx, local_nz, ngz, nzgrid, kyarray, kxarray, zarray, total_nz, local_nz_offset, iodd, ieven
! use discretization
USE lagrange_interpolation
- USE model
implicit none
public:: get_miller, set_miller_parameters
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 4f4974ca1a20f352c96125481a5ad69572dc688d..f74659a8fddda82c533b699edfd916e03f4da199 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -17,14 +17,15 @@ MODULE model
CHARACTER(len=16), &
PUBLIC, PROTECTED :: HYP_V = 'hypcoll' ! hyperdiffusion model for velocity space ('none','hypcoll','dvpar4')
INTEGER, PUBLIC, PROTECTED :: Na = 1 ! number of evolved species
+ LOGICAL, PUBLIC :: ADIAB_E = .false. ! adiabatic electron model
+ LOGICAL, PUBLIC :: ADIAB_I = .false. ! adiabatic ion model
+ REAL(xp), PUBLIC, PROTECTED :: tau_i = 1.0 ! electron-ion temperature ratio for ion adiabatic model
REAL(xp), PUBLIC, PROTECTED :: nu = 0._xp ! collision frequency parameter
REAL(xp), PUBLIC, PROTECTED :: k_gB = 1._xp ! Magnetic gradient strength (L_ref/L_gB)
REAL(xp), PUBLIC, PROTECTED :: k_cB = 1._xp ! Magnetic curvature strength (L_ref/L_cB)
REAL(xp), PUBLIC, PROTECTED :: lambdaD = 0._xp ! Debye length
REAL(xp), PUBLIC, PROTECTED :: beta = 0._xp ! electron plasma Beta (8piNT_e/B0^2)
- LOGICAL, PUBLIC :: ADIAB_E = .false. ! adiabatic electron model
- LOGICAL, PUBLIC :: ADIAB_I = .false. ! adiabatic ion model
- REAL(xp), PUBLIC, PROTECTED :: tau_i = 1.0 ! electron-ion temperature ratio for ion adiabatic model
+ REAL(xp), PUBLIC, PROTECTED :: ExBrate = 0._xp ! ExB background shearing rate (radially constant shear flow)
! Auxiliary variable
LOGICAL, PUBLIC, PROTECTED :: EM = .true. ! Electromagnetic effects flag
LOGICAL, PUBLIC, PROTECTED :: MHD_PD = .true. ! MHD pressure drift
@@ -40,14 +41,15 @@ CONTAINS
SUBROUTINE model_readinputs
! Read the input parameters
- USE basic, ONLY: lu_in, speak
- USE parallel, ONLY: num_procs_p
- USE prec_const
+ USE basic, ONLY: lu_in, speak
+ USE parallel, ONLY: num_procs_p
+ USE prec_const, ONLY: xp
IMPLICIT NONE
NAMELIST /MODEL_PAR/ KERN, LINEARITY, RM_LD_T_EQ, FORCE_SYMMETRY, MHD_PD, &
- mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, Na,&
- nu, k_gB, k_cB, lambdaD, beta, ADIAB_E, ADIAB_I, tau_i
+ Na, ADIAB_E, ADIAB_I, tau_i, &
+ mu_x, mu_y, N_HD, HDz_h, mu_z, mu_p, mu_j, HYP_V, &
+ nu, k_gB, k_cB, lambdaD, beta, ExBrate
READ(lu_in,model_par)
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index ea1691cb4487c14bb17181b4cd4d28e52bb87c19..f7609cb1e7f8da4328e0ac08bf946d12db193e29 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -10,7 +10,7 @@ CONTAINS
USE grid, ONLY: local_na, local_np, local_nj, local_nkx, local_nky, local_nz,&
nzgrid,pp2,ngp,ngj,ngz,&
diff_dz_coeff,diff_kx_coeff,diff_ky_coeff,diff_p_coeff,diff_j_coeff,&
- parray,jarray,kxarray, kyarray, kparray
+ parray,jarray,kxarray, kyarray, kp2array
USE basic
USE closure, ONLY: evolve_mom
USE prec_const
@@ -35,10 +35,10 @@ CONTAINS
z:DO iz = 1,local_nz
izi = iz + ngz/2
x:DO ikx = 1,local_nkx
- kx = kxarray(ikx) ! radial wavevector
- i_kx = imagu * kx ! radial derivative
y:DO iky = 1,local_nky
+ kx = kxarray(iky,ikx) ! radial wavevector
ky = kyarray(iky) ! binormal wavevector
+ i_kx = imagu * kx ! radial derivative
i_ky = imagu * ky ! binormal derivative
phikykxz = phi(iky,ikx,izi)
psikykxz = psi(iky,ikx,izi)
@@ -50,7 +50,7 @@ CONTAINS
ipi = ip+ngp/2
p_int = parray(ipi) ! Hermite degree
eo = min(nzgrid,MODULO(p_int,2)+1) ! Indicates if we are on odd or even z grid
- kperp2= kparray(iky,ikx,izi,eo)**2
+ kperp2= kp2array(iky,ikx,izi,eo)
! Species loop
a:DO ia = 1,local_na
Napj = moments(ia,ipi,iji,iky,ikx,izi,updatetlevel)
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 0224a51566afd7ccc26b6b89572d832c21727377..6ec4d1de2ea60ef855242cf3f7923ad5c05a009d 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -74,7 +74,7 @@ SUBROUTINE evaluate_kernels
USE basic
USE prec_const, ONLY: xp
USE array, ONLY : kernel!, HF_phi_correction_operator
- USE grid, ONLY : local_na, local_nj,ngj, local_nkx, local_nky, local_nz, ngz, jarray, kparray,&
+ USE grid, ONLY : local_na, local_nj,ngj, local_nkx, local_nky, local_nz, ngz, jarray, kp2array,&
nzgrid
USE species, ONLY : sigma2_tau_o2
USE model, ONLY : ADIAB_I, tau_i
@@ -91,7 +91,7 @@ DO ia = 1,local_na
DO ij = 1,local_nj + ngj
j_int = jarray(ij)
j_xp = REAL(j_int,xp)
- y_ = sigma2_tau_o2(ia) * kparray(iky,ikx,iz,eo)**2
+ y_ = sigma2_tau_o2(ia) * kp2array(iky,ikx,iz,eo)
IF(j_int .LT. 0) THEN !ghosts values
kernel(ia,ij,iky,ikx,iz,eo) = 0._xp
ELSE
@@ -114,7 +114,7 @@ DO ia = 1,local_na
DO ij = 1,local_nj + ngj
j_int = jarray(ij)
j_xp = REAL(j_int,xp)
- y_ = sigma_i**2*tau_i/2._xp * kparray(iky,ikx,iz,eo)**2
+ y_ = sigma_i**2*tau_i/2._xp * kp2array(iky,ikx,iz,eo)
IF(j_int .LT. 0) THEN !ghosts values
kernel_i(ij,iky,ikx,iz,eo) = 0._xp
ELSE
@@ -174,7 +174,7 @@ SUBROUTINE evaluate_poisson_op
kxloop: DO ikx = 1,local_nkx
kyloop: DO iky = 1,local_nky
zloop: DO iz = 1,local_nz
- IF( (kxarray(ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
+ IF( (kxarray(iky,ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
inv_poisson_op(iky, ikx, iz) = 0._xp
inv_pol_ion (iky, ikx, iz) = 0._xp
ELSE
@@ -217,7 +217,7 @@ SUBROUTINE evaluate_ampere_op
USE prec_const, ONLY : xp
USE array, ONLY : kernel, inv_ampere_op
USE grid, ONLY : local_na, local_nkx, local_nky, local_nz, ngz, total_nj, ngj,&
- kparray, kxarray, kyarray, SOLVE_AMPERE, iodd
+ kp2array, kxarray, kyarray, SOLVE_AMPERE, iodd
USE model, ONLY : beta, ADIAB_I
USE species, ONLY : q, sigma
USE geometry, ONLY : hatB
@@ -232,8 +232,8 @@ SUBROUTINE evaluate_ampere_op
x:DO ikx = 1,local_nkx
y:DO iky = 1,local_nky
z:DO iz = 1,local_nz
- kperp2 = kparray(iky,ikx,iz+ngz/2,iodd)**2
- IF( (kxarray(ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
+ kperp2 = kp2array(iky,ikx,iz+ngz/2,iodd)
+ IF( (kxarray(iky,ikx).EQ.0._xp) .AND. (kyarray(iky).EQ.0._xp) ) THEN
inv_ampere_op(iky, ikx, iz) = 0._xp
ELSE
sum_jpol = 0._xp
diff --git a/src/prec_const_mod.F90 b/src/prec_const_mod.F90
index ce1efa47661cc1a0c5cb03aeaf34bf150b640350..f8fedacda57a3c3e1006e46c87f0a48140bad660 100644
--- a/src/prec_const_mod.F90
+++ b/src/prec_const_mod.F90
@@ -78,7 +78,7 @@ MODULE prec_const
! dp_r = range(b)
! dp_p = precision(b)
- REAL(xp) :: c = 1_xp
+ REAL(xp) :: c = 1._xp
xp_r = range(c)
xp_p = precision(c)
diff --git a/src/stepon.F90 b/src/stepon.F90
index 91a7298b1281a5072d2c7c0559e5e1d542b7e008..b511e179f0e616acca9f3b08e256b177de68f9e8 100644
--- a/src/stepon.F90
+++ b/src/stepon.F90
@@ -8,6 +8,7 @@ SUBROUTINE stepon
use mpi, ONLY: MPI_COMM_WORLD
USE time_integration, ONLY: ntimelevel
USE prec_const, ONLY: xp
+ USE ExB_shear_flow, ONLY: Apply_ExB_shear_flow, Update_ExB_shear_flow
#ifdef TEST_SVD
USE CLA, ONLY: test_svd,filter_sv_moments_ky_pj
#endif
@@ -16,6 +17,14 @@ SUBROUTINE stepon
INTEGER :: num_step, ierr
LOGICAL :: mlend
+ ! Update the ExB backg. shear flow before the step
+ ! This call includes :
+ ! - the kx grid update
+ ! - the kernel, poisson op. and ampere op update
+ ! - kx-shift of the fields at required ky values
+ ! - the ExB shear value (s(ky)) update for the next time step
+ CALL Apply_ExB_shear_flow
+
DO num_step=1,ntimelevel ! eg RK4 compute successively k1, k2, k3, k4
!----- BEFORE: All fields+ghosts are updated for step = n
! Compute right hand side from current fields
@@ -70,6 +79,9 @@ SUBROUTINE stepon
END DO
+ ! Update the ExB shear flow for the next step
+ CALL Update_ExB_shear_flow
+
CONTAINS
!!!! Basic structure to simplify stepon
SUBROUTINE assemble_RHS
diff --git a/wk/lin_run_script.m b/wk/lin_run_script.m
index 97644a1790f9d1067e49b2f388db569bb269b40b..b4832e7da8a2dc06bf8c925f311142352b0c82c2 100644
--- a/wk/lin_run_script.m
+++ b/wk/lin_run_script.m
@@ -15,34 +15,36 @@ addpath(genpath([gyacomodir,'matlab/load'])) % Add load folder
addpath(genpath([gyacomodir,'wk/parameters'])) % Add parameters folder
%% Setup run or load an executable
-RUN = 1; % To run or just to load
+RUN = 0; % To run or just to load
default_plots_options
+% EXECNAME = 'gyacomo23_sp_save'; % single precision
EXECNAME = 'gyacomo23_sp'; % single precision
% EXECNAME = 'gyacomo23_dp'; % double precision
+% EXECNAME = 'gyacomo23_debug'; % double precision
%% Setup parameters
% run lin_DTT_HM_rho85
% run lin_DTT_HM_rho98
% run lin_DTT_LM_rho90
% run lin_DTT_LM_rho95
-run lin_JET_rho97
+% run lin_JET_rho97
% run lin_Entropy
-% run lin_ITG
+run lin_ITG
% run lin_KBM
%% Change parameters
-NY = 2;
-NX = 6;
-NZ = 32;
-PMAX = 2;
+EXBRATE = 0.0; % Background ExB shear flow
+NY = 40;
+NX = 8;
+PMAX = 16;
JMAX = PMAX/2;
-ky = 0.5;
+ky = 0.05;
LY = 2*pi/ky;
-DT = 1e-3;
-TMAX = 10;
+DT = 5e-3;
+TMAX = 50;
% % SIGMA_E = 0.04;
% TMAX = 10;
% DTSAVE0D = 200*DT;
-DTSAVE3D = TMAX/50;
+% DTSAVE3D = TMAX/50;
%%-------------------------------------------------------------------------
%% RUN
setup
@@ -51,9 +53,9 @@ setup
if RUN
MVIN =['cd ../results/',SIMID,'/',PARAMS,'/;'];
% RUN =['time ',mpirun,' -np 2 ',gyacomodir,'bin/',EXECNAME,' 1 2 1 0;'];
- RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
- % RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
-% RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
+ % RUN =['time ',mpirun,' -np 4 ',gyacomodir,'bin/',EXECNAME,' 1 2 2 0;'];
+ RUN =['time ',mpirun,' -np 8 ',gyacomodir,'bin/',EXECNAME,' 2 2 2 0;'];
+ % RUN =['time ',mpirun,' -np 1 ',gyacomodir,'bin/',EXECNAME,' 1 1 1 0;'];
% RUN = ['./../../../bin/gyacomo23_sp 0;'];
MVOUT='cd ../../../wk;';
system([MVIN,RUN,MVOUT]);
@@ -84,8 +86,9 @@ options.NMA = 1; % Set NMA option to 1
options.NMODES = 999; % Set how much modes we study
options.iz = 'avg'; % Compressing z
options.ik = 1; %
+options.GOK2 = 0; % plot gamma/k^2
options.fftz.flag = 0; % Set fftz.flag option to 0
-fig = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
+[fig, kykx, wkykx, ekykx] = mode_growth_meter(data,options); % Call the function mode_growth_meter with data and options as input arguments, and store the result in fig
end
if (1 && NZ>4)
diff --git a/wk/lin_scan_script.m b/wk/lin_scan_script.m
index 8831740be18598ee3c075df149aef21272dbafd6..2dea07f623d52d7151946d1c7ed97270ef2d0c6c 100644
--- a/wk/lin_scan_script.m
+++ b/wk/lin_scan_script.m
@@ -30,6 +30,8 @@ run lin_JET_rho97
%% Change parameters
NY = 2;
+EXBRATE = 0;
+% SIGMA_E = 0.023;
%% Scan parameters
SIMID = [SIMID,'_scan'];
P_a = [2 4 6 8];
diff --git a/wk/load_and_replot.m b/wk/load_and_replot.m
new file mode 100644
index 0000000000000000000000000000000000000000..3d0e7ee51fd0f3156026196434d16f4f95057fb9
--- /dev/null
+++ b/wk/load_and_replot.m
@@ -0,0 +1,23 @@
+% Load the .fig file
+figFile = '/home/ahoffman/paper2_fig/31_CBC_linear_convergence.fig'; % Replace with the actual path to your .fig file
+% figFile = 'tmp.fig'; % Replace with the actual path to your .fig file
+figHandle = openfig(figFile);
+
+%% Extract data from the figure
+% Here, we assume that the plotted data is stored in a single line plot
+lineHandle = findobj(figHandle, 'Type', 'line');
+xData = get(lineHandle, 'XData');
+yData = get(lineHandle, 'YData');
+% Extract legend labels
+legendHandle = findobj(figHandle, 'Type', 'legend');
+legendEntries = flip(get(legendHandle, 'String'));
+% legendEntries = legendEntries{end:-1:1};
+% Close the loaded figure
+close(figHandle);
+%%
+ figure
+for i = numel(xData):-1:1
+ % plot(xData{i},yData{i},'DisplayName',legendEntries{i}); hold on
+ plot(xData{i},yData{i}./xData{i}.^2,'DisplayName',legendEntries{i}); hold on
+end
+legend show
\ No newline at end of file
diff --git a/wk/load_metadata_scan.m b/wk/load_metadata_scan.m
index f107d4afec12cbd35cd8844c282d73a5cc056d90..7ef57a678dfc60073e788d78e736c638f4d0a4be 100644
--- a/wk/load_metadata_scan.m
+++ b/wk/load_metadata_scan.m
@@ -15,7 +15,11 @@ addpath(genpath([gyacomodir,'matlab/load'])) % ... add% EXECNAME = 'gyacomo_1.0'
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_1_P_2_8_DGGK_0.05_be_0.0039.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/16x24_ky_0.1_0.9_P_2_8_kN_1.33_DGGK_0.56901_be_0.0039.mat';
% datafname = 'lin_DTT_AB_rho85_PT_scan/8x24_ky_0.05_1.5_P_2_8_kN_1.33_DGGK_0.1_be_0.0039.mat';
-datafname = 'lin_JET_rho97_scan/8x32_ky_0.031623_31.6228_P_1_1_kN_31.4286_DGGK_0.1_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/4x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.1_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/4x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.05_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/8x32_ky_0.01_10_P_2_8_kN_10_DGGK_0.05_be_0.0031.mat';
+% datafname = 'lin_JET_rho97_scan/8x32_ky_0.01_10_P_2_2_kN_10_DGGK_0.05_be_0.0031.mat';
+datafname = 'lin_JET_rho97_scan/8x32_ky_0.031623_31.6228_P_2_8_kN_10_DGGK_0.1_be_0.0031.mat';
%% Chose if we filter gamma>0.05
FILTERGAMMA = 1;
diff --git a/wk/parameters/lin_ITG.m b/wk/parameters/lin_ITG.m
index 7e2167a3f8646b08f56338a158002d9e295063ae..07d85ebc3aa01595ebd7e924280e8eaf25d875c9 100644
--- a/wk/parameters/lin_ITG.m
+++ b/wk/parameters/lin_ITG.m
@@ -3,16 +3,17 @@ SIMID = 'lin_ITG'; % Name of the simulation
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
-NU = 0.005; % Collision frequency
-TAU = 1.0; % e/i temperature ratio
-K_Ne = 0*2.22; % ele Density
-K_Te = 0*6.96; % ele Temperature
-K_Ni = 2.22; % ion Density gradient drive
-K_Ti = 6.96; % ion Temperature
+NU = 0.005; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 0*2.22; % ele Density
+K_Te = 0*6.96; % ele Temperature
+K_Ni = 2.22; % ion Density gradient drive
+K_Ti = 6.96; % ion Temperature
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
-NA = 1; % number of kinetic species
+NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model
-BETA = 0.0; % electron plasma beta
+BETA = 0.0; % electron plasma beta
+EXBRATE = 0.0; % Background ExB shear flow
%% Set up grid parameters
P = 4;
J = 2;%P/2;
@@ -40,8 +41,8 @@ SHIFT_Y = 0.0; % Shift in the periodic BC in z
NPOL = 1; % Number of poloidal turns
%% TIME PARAMETERS
-TMAX = 50; % Maximal time unit
-DT = 1e-2; % Time step
+TMAX = 20; % Maximal time unit
+DT = 2e-2; % Time step
DTSAVE0D = 1; % Sampling time for 0D arrays
DTSAVE2D = -1; % Sampling time for 2D arrays
DTSAVE3D = 1; % Sampling time for 3D arrays