From 6c5afda664ee7171deb0097596cb4f1b37dfc2c5 Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Fri, 29 Jul 2022 11:26:37 +0200
Subject: [PATCH] Ampere solver ready
---
matlab/compile_results.m | 28 ++++++-----
matlab/compute/compute_fluxtube_growth_rate.m | 15 ++++--
matlab/load/load_params.m | 15 +++---
matlab/plot/plot_ballooning.m | 47 ++++++++++++++-----
matlab/setup.m | 11 +++--
matlab/write_fort90.m | 7 +--
src/auxval.F90 | 2 +-
src/model_mod.F90 | 17 ++++---
src/moments_eq_rhs_mod.F90 | 26 +++++-----
src/numerics_mod.F90 | 43 ++++++++++-------
src/solve_EM_fields.F90 | 1 -
wk/quick_run.m | 28 +++++------
12 files changed, 142 insertions(+), 98 deletions(-)
diff --git a/matlab/compile_results.m b/matlab/compile_results.m
index da220097..08b7fc6c 100644
--- a/matlab/compile_results.m
+++ b/matlab/compile_results.m
@@ -7,7 +7,7 @@ DATA.NU_EVOL = []; % evolution of parameter nu between jobs
DATA.CO_EVOL = []; % evolution of CO
DATA.MUx_EVOL = []; % evolution of parameter mu between jobs
DATA.MUy_EVOL = []; % evolution of parameter mu between jobs
-DATA.K_N_EVOL = []; %
+DATA.K_Ni_EVOL = []; %
DATA.L_EVOL = []; %
DATA.DT_EVOL = []; %
% FIELDS
@@ -20,6 +20,7 @@ PGAMMAI_ = [];
GGAMMAE_ = [];
PGAMMAE_ = [];
PHI_ = [];
+PSI_ = [];
DENS_E_ = [];
DENS_I_ = [];
UPAR_E_ = [];
@@ -61,7 +62,7 @@ while(CONTINUE)
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');
-
+ BETA = h5readatt(filename,'/data/input','beta');
% Check polynomials degrees
Pe_new= numel(Pe); Je_new= numel(Je);
Pi_new= numel(Pi); Ji_new= numel(Ji);
@@ -129,6 +130,10 @@ while(CONTINUE)
if W_PHI
[ PHI, Ts3D, ~] = load_3D_data(filename, 'phi');
PHI_ = cat(4,PHI_,PHI); clear PHI
+ if BETA > 0
+ [ PSI, Ts3D, ~] = load_3D_data(filename, 'psi');
+ PSI_ = cat(4,PSI_,PSI); clear PSI
+ end
end
if W_NA00
if KIN_E
@@ -196,14 +201,14 @@ while(CONTINUE)
% Evolution of simulation parameters
load_params
- DATA.TJOB_SE = [DATA.TJOB_SE Ts0D(1) Ts0D(end)];
- DATA.NU_EVOL = [DATA.NU_EVOL DATA.NU DATA.NU];
- DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO];
- DATA.MUx_EVOL = [DATA.MUx_EVOL DATA.MUx DATA.MUx];
- DATA.MUy_EVOL = [DATA.MUy_EVOL DATA.MUy DATA.MUy];
- DATA.K_N_EVOL = [DATA.K_N_EVOL DATA.K_N DATA.K_N];
- DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L];
- DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM];
+ DATA.TJOB_SE = [DATA.TJOB_SE Ts0D(1) Ts0D(end)];
+ DATA.NU_EVOL = [DATA.NU_EVOL DATA.NU DATA.NU];
+ DATA.CO_EVOL = [DATA.CO_EVOL DATA.CO DATA.CO];
+ DATA.MUx_EVOL = [DATA.MUx_EVOL DATA.MUx DATA.MUx];
+ DATA.MUy_EVOL = [DATA.MUy_EVOL DATA.MUy DATA.MUy];
+ DATA.K_Ni_EVOL = [DATA.K_Ni_EVOL DATA.K_Ni DATA.K_Ni];
+ DATA.L_EVOL = [DATA.L_EVOL DATA.L DATA.L];
+ DATA.DT_EVOL = [DATA.DT_EVOL DATA.DT_SIM DATA.DT_SIM];
ii = ii + 1;
JOBFOUND = JOBFOUND + 1;
@@ -263,6 +268,7 @@ else
end
DATA.Ts5D = Ts5D_; DATA.Ts3D = Ts3D_; DATA.Ts0D = Ts0D_;
DATA.PHI = PHI_;
+ DATA.PSI = PSI_;
DATA.KIN_E=KIN_E;
% grids
DATA.Pe = Pe; DATA.Pi = Pi;
@@ -275,7 +281,7 @@ else
DATA.dir = DIRECTORY;
DATA.localdir = DIRECTORY;
DATA.param_title=['$\nu_{',DATA.CONAME,'}=$', num2str(DATA.NU), ...
- ', $\kappa_N=$',num2str(DATA.K_N),', $\kappa_T=$',num2str(DATA.K_T),...
+ ', $\kappa_{Ni}=$',num2str(DATA.K_Ni),', $\kappa_{Ti}=$',num2str(DATA.K_Ti),...
', $L=',num2str(DATA.L),'$, $N=',...
num2str(DATA.Nx),'$, $(P,J)=(',num2str(DATA.PMAXI),',',...
num2str(DATA.JMAXI),')$,',' $\mu_{hd}=$(',num2str(DATA.MUx),...
diff --git a/matlab/compute/compute_fluxtube_growth_rate.m b/matlab/compute/compute_fluxtube_growth_rate.m
index 7ee54f66..4cae40da 100644
--- a/matlab/compute/compute_fluxtube_growth_rate.m
+++ b/matlab/compute/compute_fluxtube_growth_rate.m
@@ -46,8 +46,8 @@ linear_gr.avg_w = mean(imag(w_ky(Is,:)),2);
if PLOT >0
figure
-if PLOT >1
- subplot(1,PLOT,1)
+if PLOT > 1
+ subplot(1,2,1)
end
plot(linear_gr.ky,linear_gr.g_ky(:,end),'-o','DisplayName','$Re(\omega_{k_y})$'); hold on;
@@ -67,14 +67,19 @@ end
title(DATA.param_title);
if PLOT > 1
- subplot(1,PLOT,2)
+ if PLOT == 2
+ subplot(1,2,2)
+ elseif PLOT == 3
+ subplot(2,2,2)
+ end
plot(DATA.Ts3D(its+1:ite),linear_gr.g_ky(Is,:)); hold on;
plot(DATA.Ts3D(its+1:ite),linear_gr.w_ky(Is,:));
- xlabel('t'); ylabel('$\gamma(t),\omega(t)$')
+ xlabel('t'); ylabel('$\gamma(t),\omega(t)$'); xlim([DATA.Ts3D(1) DATA.Ts3D(end)]);
end
if PLOT > 2
- subplot(1,PLOT,3)
+ xlabel([]); xticks([]);
+ subplot(2,2,4)
semilogy(DATA.Ts3D,squeeze(abs(DATA.PHI(2,1,DATA.Nz/2,:)))); hold on;
xlabel('t'); ylabel('$|\phi_{ky}|(t)$')
diff --git a/matlab/load/load_params.m b/matlab/load/load_params.m
index 613a62bd..89e09722 100644
--- a/matlab/load/load_params.m
+++ b/matlab/load/load_params.m
@@ -1,9 +1,8 @@
DATA.CO = h5readatt(filename,'/data/input','CO');
-% DATA.K_N = h5readatt(filename,'/data/input','eta_n');
-% DATA.K_T = h5readatt(filename,'/data/input','eta_T');
-DATA.K_N = h5readatt(filename,'/data/input','K_n');
-DATA.K_T = h5readatt(filename,'/data/input','K_T');
-DATA.K_E = h5readatt(filename,'/data/input','K_E');
+DATA.K_Ne = h5readatt(filename,'/data/input','K_Ne');
+DATA.K_Ni = h5readatt(filename,'/data/input','K_Ni');
+DATA.K_Te = h5readatt(filename,'/data/input','K_Te');
+DATA.K_Ti = h5readatt(filename,'/data/input','K_Ti');
DATA.Q0 = h5readatt(filename,'/data/input','q0');
DATA.SHEAR = h5readatt(filename,'/data/input','shear');
DATA.EPS = h5readatt(filename,'/data/input','eps');
@@ -22,7 +21,7 @@ DATA.DT_SIM = h5readatt(filename,'/data/input','dt');
DATA.MU = h5readatt(filename,'/data/input','mu');
DATA.MUx = h5readatt(filename,'/data/input','mu_x');
DATA.MUy = h5readatt(filename,'/data/input','mu_y');
-% MU = str2num(filename(end-18:end-14)); %bad...
+DATA.BETA = h5readatt(filename,'/data/input','beta');
DATA.W_GAMMA = h5readatt(filename,'/data/input','write_gamma') == 'y';
DATA.W_PHI = h5readatt(filename,'/data/input','write_phi') == 'y';
DATA.W_NA00 = h5readatt(filename,'/data/input','write_Na00') == 'y';
@@ -47,9 +46,9 @@ else
degngrad = ['Pe_',num2str(DATA.PMAXE),'_Je_',num2str(DATA.JMAXE),...
'_Pi_',num2str(DATA.PMAXI),'_Ji_',num2str(DATA.JMAXI)];
end
-degngrad = [degngrad,'_Kn_%1.1f_nu_%0.0e_',...
+degngrad = [degngrad,'_Kni_%1.1f_nu_%0.0e_',...
DATA.CONAME,'_CLOS_',num2str(DATA.CLOS),'_mu_%0.0e'];
-degngrad = sprintf(degngrad,[DATA.K_N,DATA.NU,DATA.MU]);
+degngrad = sprintf(degngrad,[DATA.K_Ni,DATA.NU,DATA.MU]);
% if ~DATA.LINEARITY; degngrad = ['lin_',degngrad]; end
resolution = [num2str(DATA.Nx),'x',num2str(DATA.Ny),'_'];
gridname = ['L_',num2str(DATA.L),'_'];
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index f4deff2a..acc3cc9c 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -7,6 +7,12 @@ function [FIG] = plot_ballooning(data,options)
% phi_imag=mean(imag(data.PHI(:,:,:,it0:it1)),4);
phi_real=real(data.PHI(:,:,:,it1));
phi_imag=imag(data.PHI(:,:,:,it1));
+ ncol = 1;
+ if data.BETA > 0
+ psi_real=real(data.PSI(:,:,:,it1));
+ psi_imag=imag(data.PSI(:,:,:,it1));
+ ncol = 2;
+ end
% Apply baollooning tranform
for iky=ikyarray
dims = size(phi_real);
@@ -29,19 +35,13 @@ function [FIG] = plot_ballooning(data,options)
phib_real = zeros( Nkx*dims(3) ,1);
phib_imag = phib_real;
b_angle = phib_real;
-
- kk_=[];
- ss_=[];
- ee_=[];
+
for i_ =1:Nkx
start_ = (i_ -1)*dims(3) +1;
end_ = i_*dims(3);
ikx = ordered_ikx(i_);
phib_real(start_:end_) = phi_real(iky,ikx,:);
phib_imag(start_:end_) = phi_imag(iky,ikx,:);
- kk_ = [kk_ data.kx(ikx)];
- ss_ = [ss_ start_];
- ee_ = [ee_ end_];
end
% Define ballooning angle
@@ -61,11 +61,10 @@ function [FIG] = plot_ballooning(data,options)
else
normalization = 1;
end
- phib_norm = phib / normalization ;
- phib_real_norm = real( phib_norm);%phib_real(:)/phib_real(idxLFS);
- phib_imag_norm = imag( phib_norm);%phib_imag(:)/ phib_imag(idxLFS);
-
- subplot(numel(ikyarray),1,iplot)
+ phib_norm = phib / normalization;
+ phib_real_norm = real( phib_norm);
+ phib_imag_norm = imag( phib_norm);
+ subplot(numel(ikyarray),ncol,ncol*(iplot-1)+1)
plot(b_angle/pi, phib_real_norm,'-b'); hold on;
plot(b_angle/pi, phib_imag_norm ,'-r');
plot(b_angle/pi, sqrt(phib_real_norm .^2 + phib_imag_norm.^2),'-k');
@@ -74,6 +73,30 @@ function [FIG] = plot_ballooning(data,options)
ylabel('$\phi_B (\chi)$');
title(['$k_y=',sprintf('%1.1f',data.ky(iky)),...
',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
+
+ if data.BETA > 0
+
+ psib_real = zeros( Nkx*dims(3) ,1);
+ psib_imag = psib_real;
+ for i_ =1:Nkx
+ start_ = (i_ -1)*dims(3) +1;
+ end_ = i_*dims(3);
+ ikx = ordered_ikx(i_);
+ psib_real(start_:end_) = psi_real(iky,ikx,:);
+ psib_imag(start_:end_) = psi_imag(iky,ikx,:);
+ end
+
+ subplot(numel(ikyarray),ncol,ncol*(iplot-1)+2)
+ plot(b_angle/pi, psib_real,'-b'); hold on;
+ plot(b_angle/pi, psib_imag ,'-r');
+ plot(b_angle/pi, sqrt(psib_real .^2 + psib_imag.^2),'-k');
+ legend('real','imag','norm')
+ xlabel('$\chi / \pi$')
+ ylabel('$\psi_B (\chi)$');
+ title(['$k_y=',sprintf('%1.1f',data.ky(iky)),...
+ ',t_{avg}\in [',sprintf('%1.1f',data.Ts3D(it0)),',',sprintf('%1.1f',data.Ts3D(it1)),']$']);
+ end
+
iplot = iplot + 1;
end
end
diff --git a/matlab/setup.m b/matlab/setup.m
index a11f4546..07f5a65a 100644
--- a/matlab/setup.m
+++ b/matlab/setup.m
@@ -42,9 +42,10 @@ MODEL.q_e =-1.0;
MODEL.q_i = 1.0;
if MODEL.q_e == 0; SIMID = [SIMID,'_i']; end;
% gradients L_perp/L_x
-MODEL.K_n = K_N; % source term kappa for HW
-MODEL.K_T = K_T; % Temperature
-MODEL.K_E = 0; % Electric
+MODEL.K_Ne = K_Ne; % source term kappa for HW
+MODEL.K_Ni = K_Ni; % source term kappa for HW
+MODEL.K_Te = K_Te; % Temperature
+MODEL.K_Ti = K_Te; % Temperature
MODEL.GradB = GRADB; % Magnetic gradient
MODEL.CurvB = CURVB; % Magnetic curvature
MODEL.lambdaD = LAMBDAD;
@@ -97,8 +98,8 @@ else
end
% temp. dens. drives
drives_ = [];
-if abs(K_N) > 0; drives_ = [drives_,'_kN_',num2str(K_N)]; end;
-if abs(K_T) > 0; drives_ = [drives_,'_kT_',num2str(K_T)]; end;
+if abs(K_Ni) > 0; drives_ = [drives_,'_kNi_',num2str(K_Ni)]; end;
+if abs(K_Ti) > 0; drives_ = [drives_,'_kTi_',num2str(K_Ti)]; end;
% collision
coll_ = ['_nu_%1.1e_',CONAME];
coll_ = sprintf(coll_,NU);
diff --git a/matlab/write_fort90.m b/matlab/write_fort90.m
index 90dee735..7dd03641 100644
--- a/matlab/write_fort90.m
+++ b/matlab/write_fort90.m
@@ -68,9 +68,10 @@ fprintf(fid,[' sigma_e = ', num2str(MODEL.sigma_e),'\n']);
fprintf(fid,[' sigma_i = ', num2str(MODEL.sigma_i),'\n']);
fprintf(fid,[' q_e = ', num2str(MODEL.q_e),'\n']);
fprintf(fid,[' q_i = ', num2str(MODEL.q_i),'\n']);
-fprintf(fid,[' K_n = ', num2str(MODEL.K_n),'\n']);
-fprintf(fid,[' K_T = ', num2str(MODEL.K_T),'\n']);
-fprintf(fid,[' K_E = ', num2str(MODEL.K_E),'\n']);
+fprintf(fid,[' K_Ne = ', num2str(MODEL.K_Ne),'\n']);
+fprintf(fid,[' K_Ni = ', num2str(MODEL.K_Ni),'\n']);
+fprintf(fid,[' K_Te = ', num2str(MODEL.K_Te),'\n']);
+fprintf(fid,[' K_Ti = ', num2str(MODEL.K_Ti),'\n']);
fprintf(fid,[' GradB = ', num2str(MODEL.GradB),'\n']);
fprintf(fid,[' CurvB = ', num2str(MODEL.CurvB),'\n']);
fprintf(fid,[' lambdaD = ', num2str(MODEL.lambdaD),'\n']);
diff --git a/src/auxval.F90 b/src/auxval.F90
index 5db2df05..586a9fa9 100644
--- a/src/auxval.F90
+++ b/src/auxval.F90
@@ -44,7 +44,7 @@ subroutine auxval
CALL evaluate_kernels ! precompute the kernels
- CALL evaluate_poisson_op ! precompute the kernels
+ CALL evaluate_EM_op ! compute inverse of poisson and ampere operators
IF ( LINEARITY .NE. 'linear' ) THEN;
CALL build_dnjs_table ! precompute the Laguerre nonlin product coeffs
diff --git a/src/model_mod.F90 b/src/model_mod.F90
index 12d16a58..d185f381 100644
--- a/src/model_mod.F90
+++ b/src/model_mod.F90
@@ -24,9 +24,10 @@ MODULE model
REAL(dp), PUBLIC, PROTECTED :: sigma_i = 1._dp !
REAL(dp), PUBLIC, PROTECTED :: q_e = -1._dp ! Charge
REAL(dp), PUBLIC, PROTECTED :: q_i = 1._dp !
- REAL(dp), PUBLIC, PROTECTED :: K_n = 1._dp ! Density drive
- REAL(dp), PUBLIC, PROTECTED :: K_T = 0._dp ! Temperature drive
- REAL(dp), PUBLIC, PROTECTED :: K_E = 0._dp ! Backg. electric field drive
+ REAL(dp), PUBLIC, PROTECTED :: K_Ne = 1._dp ! Density drive
+ REAL(dp), PUBLIC, PROTECTED :: K_ni = 1._dp ! Density drive
+ REAL(dp), PUBLIC, PROTECTED :: K_Te = 0._dp ! Temperature drive
+ REAL(dp), PUBLIC, PROTECTED :: K_Ti = 0._dp ! Backg. electric field drive
REAL(dp), PUBLIC, PROTECTED :: GradB = 1._dp ! Magnetic gradient
REAL(dp), PUBLIC, PROTECTED :: CurvB = 1._dp ! Magnetic curvature
REAL(dp), PUBLIC, PROTECTED :: lambdaD = 1._dp ! Debye length
@@ -62,7 +63,7 @@ CONTAINS
NAMELIST /MODEL_PAR/ CLOS, NL_CLOS, KERN, LINEARITY, KIN_E, &
mu_x, mu_y, N_HD, mu_z, mu_p, mu_j, nu,&
tau_e, tau_i, sigma_e, sigma_i, q_e, q_i,&
- K_n, K_T, K_E, GradB, CurvB, lambdaD, beta
+ K_Ne, K_Ni, K_Te, K_Ti, GradB, CurvB, lambdaD, beta
READ(lu_in,model_par)
@@ -132,10 +133,12 @@ CONTAINS
CALL attach(fidres, TRIM(str), "sigma_i", sigma_i)
CALL attach(fidres, TRIM(str), "q_e", q_e)
CALL attach(fidres, TRIM(str), "q_i", q_i)
- CALL attach(fidres, TRIM(str), "K_n", K_n)
- CALL attach(fidres, TRIM(str), "K_T", K_T)
- CALL attach(fidres, TRIM(str), "K_E", K_E)
+ CALL attach(fidres, TRIM(str), "K_Ne", K_Ne)
+ CALL attach(fidres, TRIM(str), "K_Ni", K_Ni)
+ CALL attach(fidres, TRIM(str), "K_Te", K_Te)
+ CALL attach(fidres, TRIM(str), "K_Ti", K_Ti)
CALL attach(fidres, TRIM(str), "lambdaD", lambdaD)
+ CALL attach(fidres, TRIM(str), "beta", beta)
END SUBROUTINE model_outputinputs
END MODULE model
diff --git a/src/moments_eq_rhs_mod.F90 b/src/moments_eq_rhs_mod.F90
index 2fe44a54..760f1505 100644
--- a/src/moments_eq_rhs_mod.F90
+++ b/src/moments_eq_rhs_mod.F90
@@ -300,7 +300,7 @@ SUBROUTINE add_Maxwellian_background_terms
! 40, 01,02, 21 with background gradient dependences.
USE prec_const
USE time_integration, ONLY : updatetlevel
- USE model, ONLY: taue_qe, taui_qi, K_N, K_T, KIN_E
+ USE model, ONLY: taue_qe, taui_qi, K_Ne, K_Ni, K_Te, K_Ti, KIN_E
USE array, ONLY: moments_rhs_e, moments_rhs_i
USE grid, ONLY: contains_kx0, contains_ky0, ikx_0, iky_0,&
ips_e,ipe_e,ijs_e,ije_e,ips_i,ipe_i,ijs_i,ije_i,&
@@ -320,28 +320,28 @@ SUBROUTINE add_Maxwellian_background_terms
SELECT CASE (ip-1)
CASE(0) ! Na00 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * (1.5_dp*K_N - 1.125_dp*K_T)
+ +taue_qe * sinz(izs:ize) * (1.5_dp*K_Ne - 1.125_dp*K_Te)
CASE(2) ! Na20 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*K_N - 2.75_dp*K_T)
+ +taue_qe * sinz(izs:ize) * (SQRT2*0.5_dp*K_Ne - 2.75_dp*K_Te)
CASE(4) ! Na40 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*K_T
+ +taue_qe * sinz(izs:ize) * SQRT6*0.75_dp*K_Te
END SELECT
CASE(1) ! j = 1
SELECT CASE (ip-1)
CASE(0) ! Na01 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taue_qe * sinz(izs:ize) * (K_N + 3.5_dp*K_T)
+ -taue_qe * sinz(izs:ize) * (K_Ne + 3.5_dp*K_Te)
CASE(2) ! Na21 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taue_qe * sinz(izs:ize) * SQRT2*K_T
+ -taue_qe * sinz(izs:ize) * SQRT2*K_Te
END SELECT
CASE(2) ! j = 2
SELECT CASE (ip-1)
CASE(0) ! Na02 term
moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_e(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taue_qe * sinz(izs:ize) * 2._dp*K_T
+ +taue_qe * sinz(izs:ize) * 2._dp*K_Te
END SELECT
END SELECT
ENDDO
@@ -355,28 +355,28 @@ SUBROUTINE add_Maxwellian_background_terms
SELECT CASE (ip-1)
CASE(0) ! Na00 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * (1.5_dp*K_N - 1.125_dp*K_T)
+ +taui_qi * sinz(izs:ize) * (1.5_dp*K_Ni - 1.125_dp*K_Ti)
CASE(2) ! Na20 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*K_N - 2.75_dp*K_T)
+ +taui_qi * sinz(izs:ize) * (SQRT2*0.5_dp*K_Ni - 2.75_dp*K_Ti)
CASE(4) ! Na40 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*K_T
+ +taui_qi * sinz(izs:ize) * SQRT6*0.75_dp*K_Ti
END SELECT
CASE(1) ! j = 1
SELECT CASE (ip-1)
CASE(0) ! Na01 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taui_qi * sinz(izs:ize) * (K_N + 3.5_dp*K_T)
+ -taui_qi * sinz(izs:ize) * (K_Ni + 3.5_dp*K_Ti)
CASE(2) ! Na21 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- -taui_qi * sinz(izs:ize) * SQRT2*K_T
+ -taui_qi * sinz(izs:ize) * SQRT2*K_Ti
END SELECT
CASE(2) ! j = 2
SELECT CASE (ip-1)
CASE(0) ! Na02 term
moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel) = moments_rhs_i(ip,ij,iky_0,ikx_0,izs:ize,updatetlevel)&
- +taui_qi * sinz(izs:ize) * 2._dp*K_T
+ +taui_qi * sinz(izs:ize) * 2._dp*K_Ti
END SELECT
END SELECT
ENDDO
diff --git a/src/numerics_mod.F90 b/src/numerics_mod.F90
index 1ede5ac6..61db0f4f 100644
--- a/src/numerics_mod.F90
+++ b/src/numerics_mod.F90
@@ -6,7 +6,7 @@ MODULE numerics
USE coeff
implicit none
- PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_poisson_op, evaluate_ampere_op
+ PUBLIC :: build_dnjs_table, evaluate_kernels, evaluate_EM_op
PUBLIC :: compute_lin_coeff, play_with_modes, save_EM_ZF_modes
CONTAINS
@@ -104,6 +104,13 @@ END SUBROUTINE evaluate_kernels
!******************************************************************************!
!******************************************************************************!
+SUBROUTINE evaluate_EM_op
+ IMPLICIT NONE
+
+ CALL evaluate_poisson_op
+ CALL evaluate_ampere_op
+
+END SUBROUTINE evaluate_EM_op
!!!!!!! Evaluate inverse polarisation operator for Poisson equation
!******************************************************************************!
SUBROUTINE evaluate_poisson_op
@@ -201,7 +208,7 @@ END SUBROUTINE evaluate_ampere_op
SUBROUTINE compute_lin_coeff
USE array
USE model, ONLY: taue_qe, taui_qi, sqrtTaue_qe, sqrtTaui_qi, &
- K_T, K_n, CurvB, GradB, KIN_E, tau_e, tau_i, sigma_e, sigma_i
+ K_Te, K_Ti, K_Ne, K_Ni, CurvB, GradB, KIN_E, tau_e, tau_i, sigma_e, sigma_i
USE prec_const
USE grid, ONLY: parray_e, parray_i, jarray_e, jarray_i, &
ip,ij, ips_e,ipe_e, ips_i,ipe_i, ijs_e,ije_e, ijs_i,ije_i
@@ -297,11 +304,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij_e(ip,ij) =+K_n + 2.*j_dp*K_T
- xphijp1_e(ip,ij) =-K_T*(j_dp+1._dp)
- xphijm1_e(ip,ij) =-K_T* j_dp
+ xphij_e(ip,ij) =+K_Ne + 2.*j_dp*K_Te
+ xphijp1_e(ip,ij) =-K_Te*(j_dp+1._dp)
+ xphijm1_e(ip,ij) =-K_Te* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
- xphij_e(ip,ij) =+K_T/SQRT2
+ xphij_e(ip,ij) =+K_Te/SQRT2
xphijp1_e(ip,ij) = 0._dp; xphijm1_e(ip,ij) = 0._dp;
ELSE
xphij_e(ip,ij) = 0._dp; xphijp1_e(ip,ij) = 0._dp
@@ -317,11 +324,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 0) THEN ! kronecker p0
- xphij_i(ip,ij) =+K_n + 2._dp*j_dp*K_T
- xphijp1_i(ip,ij) =-K_T*(j_dp+1._dp)
- xphijm1_i(ip,ij) =-K_T* j_dp
+ xphij_i(ip,ij) =+K_Ni + 2._dp*j_dp*K_Ti
+ xphijp1_i(ip,ij) =-K_Ti*(j_dp+1._dp)
+ xphijm1_i(ip,ij) =-K_Ti* j_dp
ELSE IF (p_int .EQ. 2) THEN ! kronecker p2
- xphij_i(ip,ij) =+K_T/SQRT2
+ xphij_i(ip,ij) =+K_Ti/SQRT2
xphijp1_i(ip,ij) = 0._dp; xphijm1_i(ip,ij) = 0._dp;
ELSE
xphij_i(ip,ij) = 0._dp; xphijp1_i(ip,ij) = 0._dp
@@ -339,11 +346,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 1) THEN ! kronecker p1
- xpsij_e(ip,ij) =+(K_n + (2._dp*j_dp+1._dp)*K_T) * SQRT(tau_e)/sigma_e
- xpsijp1_e(ip,ij) =- K_T*(j_dp+1._dp) * SQRT(tau_e)/sigma_e
- xpsijm1_e(ip,ij) =- K_T* j_dp * SQRT(tau_e)/sigma_e
+ xpsij_e(ip,ij) =+(K_Ne + (2._dp*j_dp+1._dp)*K_Te)* SQRT(tau_e)/sigma_e
+ xpsijp1_e(ip,ij) =- K_Te*(j_dp+1._dp) * SQRT(tau_e)/sigma_e
+ xpsijm1_e(ip,ij) =- K_Te* j_dp * SQRT(tau_e)/sigma_e
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
- xpsij_e(ip,ij) =+ K_T*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e
+ xpsij_e(ip,ij) =+ K_Te*SQRT3/SQRT2 * SQRT(tau_e)/sigma_e
xpsijp1_e(ip,ij) = 0._dp; xpsijm1_e(ip,ij) = 0._dp;
ELSE
xpsij_e(ip,ij) = 0._dp; xpsijp1_e(ip,ij) = 0._dp
@@ -359,11 +366,11 @@ SUBROUTINE compute_lin_coeff
j_dp = REAL(j_int,dp) ! REALof Laguerre degree
!! Electrostatic potential pj terms
IF (p_int .EQ. 1) THEN ! kronecker p1
- xpsij_i(ip,ij) =+(K_n + (2._dp*j_dp+1._dp)*K_T) * SQRT(tau_i)/sigma_i
- xpsijp1_i(ip,ij) =- K_T*(j_dp+1._dp) * SQRT(tau_i)/sigma_i
- xpsijm1_i(ip,ij) =- K_T* j_dp * SQRT(tau_i)/sigma_i
+ xpsij_i(ip,ij) =+(K_Ni + (2._dp*j_dp+1._dp)*K_Ti)* SQRT(tau_i)/sigma_i
+ xpsijp1_i(ip,ij) =- K_Ti*(j_dp+1._dp) * SQRT(tau_i)/sigma_i
+ xpsijm1_i(ip,ij) =- K_Ti* j_dp * SQRT(tau_i)/sigma_i
ELSE IF (p_int .EQ. 3) THEN ! kronecker p3
- xpsij_i(ip,ij) =+ K_T*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i
+ xpsij_i(ip,ij) =+ K_Ti*SQRT3/SQRT2 * SQRT(tau_i)/sigma_i
xpsijp1_i(ip,ij) = 0._dp; xpsijm1_i(ip,ij) = 0._dp;
ELSE
xpsij_i(ip,ij) = 0._dp; xpsijp1_i(ip,ij) = 0._dp
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index b102ac1a..d1de15b0 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -106,7 +106,6 @@ CONTAINS
CALL cpu_time(t0_poisson)
!! Ampere can be solved only with beta > 0 and for process containing p=1
IF ( SOLVE_AMPERE ) THEN
-
kxloop: DO ikx = ikxs,ikxe
kyloop: DO iky = ikys,ikye
psi(iky,ikx,izs:ize) = 0._dp
diff --git a/wk/quick_run.m b/wk/quick_run.m
index 0769f58e..2b061efb 100644
--- a/wk/quick_run.m
+++ b/wk/quick_run.m
@@ -4,7 +4,7 @@
% for benchmark and debugging purpose since it makes matlab "busy"
%
SIMID = 'Kinetic_ballooning_mode'; % Name of the simulation
-RUN = 0; % To run or just to load
+RUN = 1; % To run or just to load
addpath(genpath('../matlab')) % ... add
default_plots_options
HELAZDIR = '/home/ahoffman/HeLaZ/';
@@ -15,12 +15,12 @@ EXECNAME = 'helaz3';
CLUSTER.TIME = '99:00:00'; % allocation time hh:mm:ss
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PHYSICAL PARAMETERS
-NU = 0.05; % Collision frequency
-TAU = 1.0; % e/i temperature ratio
-K_N = 2.22;%2.0; % ion Density gradient drive
-K_Ne = K_N; % ele Density gradient drive
-K_T = 6.96;%0.25*K_N; % Temperature '''
-K_Te = K_T; % Temperature '''
+NU = 0.05; % Collision frequency
+TAU = 1.0; % e/i temperature ratio
+K_Ne = 3.0; % ele Density gradient drive
+K_Ni = 3.0;%2.0; % ion Density gradient drive
+K_Te = 4.5; % Temperature '''
+K_Ti = 8.0;%0.25*K_N; % Temperature '''
SIGMA_E = 0.05196152422706632; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
% SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
KIN_E = 1; % 1: kinetic electrons, 2: adiabatic electrons
@@ -32,11 +32,11 @@ PMAXE = P; % Hermite basis size of electrons
JMAXE = J; % Laguerre "
PMAXI = P; % " ions
JMAXI = J; % "
-NX = 12; % real space x-gridpoints
-NY = 10; % '' y-gridpoints
+NX = 6; % real space x-gridpoints
+NY = 2; % '' y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain
-LY = 2*pi/0.1; % Size of the squared frequency domain
-NZ = 16; % number of perpendicular planes (parallel grid)
+LY = 2*pi/0.25; % Size of the squared frequency domain
+NZ = 32; % number of perpendicular planes (parallel grid)
NPOL = 1;
SG = 0; % Staggered z grids option
%% GEOMETRY
@@ -47,7 +47,7 @@ Q0 = 1.4; % safety factor
SHEAR = 0.8; % magnetic shear (Not implemented yet)
EPS = 0.18; % inverse aspect ratio
%% TIME PARMETERS
-TMAX = 30; % Maximal time unit
+TMAX = 40; % Maximal time unit
DT = 5e-3; % Time step
SPS0D = 1; % Sampling per time unit for 2D arrays
SPS2D = 0; % Sampling per time unit for 2D arrays
@@ -115,7 +115,7 @@ data = compile_results(LOCALDIR,JOBNUMMIN,JOBNUMMAX); %Compile the results from
if 1
%% linear growth rate (adapted for 2D zpinch and fluxtube)
trange = [0.5 1]*data.Ts3D(end);
-nplots = 1;
+nplots = 3;
lg = compute_fluxtube_growth_rate(data,trange,nplots);
[gmax, kmax] = max(lg.g_ky(:,end));
[gmaxok, kmaxok] = max(lg.g_ky(:,end)./lg.ky);
@@ -128,7 +128,7 @@ if 0
options.time_2_plot = [120];
options.kymodes = [0.1 0.2 0.3];
options.normalized = 1;
-options.field = 'phi';
+% options.field = 'phi';
fig = plot_ballooning(data,options);
end
--
GitLab