diff --git a/README.md b/README.md
index c042ae5e5b2f603d6405b31216844560d1dbe79c..0a8bfaed7acd400663b2ce30422ab9f77b745054 100644
--- a/README.md
+++ b/README.md
@@ -16,7 +16,7 @@ You should have received a copy of the GNU General Public License along with thi
Author: Antoine C.D. Hoffmann
-Contact: antoine.hoffmann@epfl.ch (do not hesitate!)
+Contact: antoine.hoffmann@epfl.ch
##### Citing GYACOMO
If you use GYACOMO in your work, please cite (at least) one of the following paper:
@@ -60,7 +60,7 @@ To compile and run GYACOMO, follow these steps:
1. Make sure the correct library paths are set in local/dirs.inc. Refer to INSTALATION.txt for instructions on installing the required libraries.
2. Go to the /gyacomo directory and run make to compile the code. The resulting binary will be located in /gyacomo/bin. You can also compile a debug version by running make dbg.
-3. Once the compilation is done, you can test the executable in the directory /testcases/cyclone_example (see README there).
+3. Once the compilation is done, you can test the executable in the directory /testcases/cyclone_example (see README there) or /testcases/zpinch_example (see wiki).
4. To stop the simulation without corrupting the output file, create a blank file called "mystop" using $touch mystop in the directory where the simulation is running. The code looks for this file every 100 time steps and the file will be removed once the simulation is stopped.
5. It is possible to chain simulations by using the parameter "Job2load" in the fort_XX.90 file. For example, to restart a simulation from the latest 5D state saved in outputs_00.h5, create a new fort.90 file called fort_01.90 and set "Job2load" to 0. Then run GYACOMO with the command ./gyacomo 0 or mpirun -np N ./gyacomo Np Ny Nz 0. This will create a new output file called output_01.h5.
6. To generate plots and gifs using the simulation results, use the script gyacomo/wk/gyacomo_analysis.m and specify the directory where the results are located. Note that this script is not currently a function.
diff --git a/matlab/compute/compute_coll.m b/matlab/compute/compute_coll.m
new file mode 100644
index 0000000000000000000000000000000000000000..c0f182ed033f1504772240b48f27f5a0fb6ea39b
--- /dev/null
+++ b/matlab/compute/compute_coll.m
@@ -0,0 +1,9 @@
+function [nuGENE, nuGYAC] = compute_coll(Lref,n_e,T_e,T_i)
+
+ tau = T_i/T_e;
+
+ nuGENE = 2.3031E-5*Lref*(n_e)/(T_e)^2*(24.-log(sqrt(n_e*1.0E13)/T_e*0.001));
+
+ nuGYAC = 3/8*sqrt(pi/2)*tau.^(3/2)*nuGENE;
+
+end
\ No newline at end of file
diff --git a/matlab/compute/compute_growth.m b/matlab/compute/compute_growth.m
new file mode 100644
index 0000000000000000000000000000000000000000..9ef8818b18b9c9f2dc4a2c4394459992f4458f52
--- /dev/null
+++ b/matlab/compute/compute_growth.m
@@ -0,0 +1,15 @@
+function [gamma,err] = compute_growth(t,y)
+gamma = zeros(size(y));
+ for it = 2:numel(t)
+ y_n = y(it);
+ y_nm1 = y(it-1);
+ dt = t(it)-t(it-1);
+ ZS = sum(y_nm1,1);
+ wl = log(y_n./y_nm1)/dt;
+ gamma(it) = squeeze(sum(wl.*y_nm1,1)./ZS);
+ end
+ % error estimation
+ n = 5;
+ [gamma, ~, err] = sliceAverage(gamma, n);
+ err = mean(err);
+end
\ No newline at end of file
diff --git a/matlab/compute/mode_growth_meter.m b/matlab/compute/mode_growth_meter.m
index b97ab56be1247f12890edd7c71f7e0c6068f9fe2..ab5736b4f47592164e2158576a322d4c088c1d71 100644
--- a/matlab/compute/mode_growth_meter.m
+++ b/matlab/compute/mode_growth_meter.m
@@ -1,4 +1,4 @@
-function [FIGURE, kykx, wkykx, ekykx] = mode_growth_meter(DATA,OPTIONS)
+function [FIGURE, wkykx, ekykx] = mode_growth_meter(DATA,OPTIONS)
NORMALIZED = OPTIONS.NORMALIZED;
Nma = OPTIONS.NMA; %Number moving average
@@ -41,7 +41,6 @@ TW = [OPTIONS.KY_TW; OPTIONS.KX_TW];
[~,it1] = min(abs(t-TW(1,1)));
[~,it2] = min(abs(t-TW(1,2)));
[wkykx, ekykx] = compute_growth_rates(FIELD(:,:,:,it1:it2),DATA.Ts3D(it1:it2));
-kykx = meshgrid(DATA.grids.ky,DATA.grids.kx)';
FIGURE = struct();
if OPTIONS.SHOWFIG
diff --git a/matlab/load/compile_results_low_mem.m b/matlab/load/compile_results_low_mem.m
index a90409a871c3274503ef7662d614733c9076dafe..33fd76b5eb6086530d02faf9553cbc2b7dfddc2f 100644
--- a/matlab/load/compile_results_low_mem.m
+++ b/matlab/load/compile_results_low_mem.m
@@ -23,6 +23,7 @@ DATA.outfilenames = {};
ii = 1;
while(CONTINUE)
filename = sprintf([DIRECTORY,'outputs_%.2d.h5'],JOBNUM);
+ fortname = sprintf([DIRECTORY,'fort_%.2d.90'],JOBNUM);
% Check presence and jobnummax
if (exist(filename, 'file') == 2 && JOBNUM <= JOBNUMMAX)
DATA.outfilenames{ii} = filename;
@@ -35,6 +36,9 @@ while(CONTINUE)
if openable
%% load results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Loading ID %.2d (%s)\n',JOBNUM,filename);
+ % loading input parameters
+ DATA.(sprintf('fort_%.2d',JOBNUM)) = struct();
+ DATA.(sprintf('fort_%.2d',JOBNUM)) = read_namelist(fortname);
% Loading from output file
CPUTIME = h5readatt(filename,'/data/input','cpu_time');
DT_SIM = h5readatt(filename,'/data/input/basic','dt');
diff --git a/matlab/plot/plot_ballooning.m b/matlab/plot/plot_ballooning.m
index 5d8f47a4c4d39ce199529ce938d989fb5a81bc0e..f43866e82dd9cc6cf5cb02dad988bbd859b9f083 100644
--- a/matlab/plot/plot_ballooning.m
+++ b/matlab/plot/plot_ballooning.m
@@ -1,5 +1,6 @@
function [FIG, b_angle, phib, psib, ky_] = plot_ballooning(DATA,OPTIONS)
FIG.fig = figure; iplot = 1;
+ psib = 0;
[~,it0] = min(abs(DATA.Ts3D - OPTIONS.time_2_plot(1)));
[~,it1] = min(abs(DATA.Ts3D - OPTIONS.time_2_plot(end)));
[~,ikyarray] = min(abs(DATA.grids.ky - OPTIONS.kymodes));
diff --git a/matlab/plot/plot_metric.m b/matlab/plot/plot_metric.m
index 39d3855107422f1215465ced486610a4c5fc94a5..1267999b25d0121511e955858f9b6b53c6509e4f 100644
--- a/matlab/plot/plot_metric.m
+++ b/matlab/plot/plot_metric.m
@@ -1,4 +1,4 @@
-function [ fig, arrays, R, Z ] = plot_metric( data, options )
+function [ fig, arrays, R, Z] = plot_metric( data, options )
fig = 0;
% names = {'Jacobian','gradxB','gradyB','gradzB','gradz_coeff',...
% 'gxx','gxy','gxz','gyy','gyz','gzz','hatB','hatR','hatZ'};
@@ -11,12 +11,25 @@ for i_ = 1:numel(names)
namae = names{i_};
geo_arrays(:,i_) = h5read(data.outfilenames{end},['/data/metric/',namae])';
end
+gxx = geo_arrays(:,1);
+gxy = geo_arrays(:,2);
+gxz = geo_arrays(:,3);
+gyy = geo_arrays(:,4);
+gyz = geo_arrays(:,5);
+% gzz = geo_arrays(:,6);
+B = geo_arrays(:,7);
+dxB = geo_arrays(:,8);
+dyB = geo_arrays(:,9);
+dzB = geo_arrays(:,10);
+Jac = geo_arrays(:,11);
+
try
- NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS;
+ NPLOT = options.SHOW_FLUXSURF + options.SHOW_METRICS + options.SHOW_CURVOP;
catch
NPLOT = 2;
options.SHOW_FLUXSURF = 1;
- options.SHOW_METRICS = 1;
+ options.SHOW_METRICS = 0;
+ options.SHOW_CURVOP = 1;
end
if NPLOT > 0
fig = figure;
@@ -41,12 +54,46 @@ if NPLOT > 0
legend('Interpreter','none')
end
if options.SHOW_FLUXSURF
- subplot(1,2,1)
+ subplot(1,NPLOT,1)
R = [geo_arrays(:,12);geo_arrays(1,12)];
Z = [geo_arrays(:,13);geo_arrays(1,13)];
- plot(R,Z);
+ plot(R./data.inputs.EPS,Z./data.inputs.EPS,'-k');
+ xlabel('R [m]'); ylabel('Z [m]');
+ axis tight
axis equal
end
+ if options.SHOW_CURVOP
+ G1 = gxx.*gyy - gxy.*gxy;
+ G2 = gxx.*gyz - gxy.*gxz;
+ G3 = gxy.*gyz - gyy.*gxz;
+ Ckx0 = -(dyB./B + G2./G1.*dzB./B);
+ C0ky = (dxB./B - G3./G1.*dzB./B);
+ % subplot(1,NPLOT,2);
+ % plot(data.grids.z,Ckx0,'DisplayName','$\mathcal C(k_x,k_y=0)/ik_x$'); hold on
+ % plot(data.grids.z,C0ky,'DisplayName','$\mathcal C(k_x=0,k_y)/ik_y$');
+ % xlabel('$z$'); legend("show");
+ subplot(1,NPLOT,2);
+ % z = [data.grids.z; data.grids.z(1)]+pi;
+ z = data.grids.z;
+ % rho = 1+0.5*(Ckx0)./max(abs(Ckx0));
+ rho = Ckx0;
+ % rho = [rho; rho(1)];
+ minr = min(rho); maxr = max(rho);
+ plot(z,rho,'DisplayName','$\mathcal C(k_x,k_y=0)/ik_x$'); hold on
+ % rho = 1+0.5*(C0ky)./max(abs(C0ky));
+ rho = C0ky;
+ % rho = [rho; rho(1)];
+ minr = min(minr,min(rho)); maxr = max(maxr,max(rho));
+ plot(z,rho,'DisplayName','$\mathcal C(k_x=0,k_y)/ik_y$');
+ % rho = 1+0.5*(B.^(-1)./Jac)./max(abs(B.^(-1)./Jac));
+ rho = B.^(-1)./Jac;
+ % rho = [rho; rho(1)];
+ minr = min(minr,min(rho)); maxr = max(maxr,max(rho));
+ plot(z,rho,'DisplayName','$\hat B^{-1}/J_{xyz}$');
+ % polarplot(z,zeros(size(rho)),'-k','DisplayName','$0$')
+ legend("show",'Location','south','FontSize',14);
+ % rlim(1.1*[minr maxr]);
+ end
end
%outputs
arrays = squeeze(geo_arrays(:,:));
diff --git a/matlab/process_field.m b/matlab/process_field.m
index 06290b4c8af281f769afe8aad7ff0da7f2b91bd9..9390b1b2cd33e40186f458a292373df32dad64bd 100644
--- a/matlab/process_field.m
+++ b/matlab/process_field.m
@@ -230,31 +230,31 @@ switch OPTIONS.NAME
FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
for it = 1:numel(FRAMES)
- tmp = zeros(DATA.Ny,DATA.Nx,Nz);
- for iz = 1:DATA.Nz
+ tmp = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
+ for iz = 1:DATA.grids.Nz
vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
gx_ = vx_.*ni_;
-% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
- tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))));
+% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))),2));
+ tmp(:,:,iz) = abs((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))));
end
FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)));
end
case 'Q_x' % ion heat flux
NAME = 'Qx';
- FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
+ % FLD_ = 0.*DATA.PHI(:,:,:,FRAMES);
OPE_ = 1;
for it = 1:numel(FRAMES)
- tmp = zeros(DATA.Ny,DATA.Nx,Nz);
- for iz = 1:DATA.Nz
+ tmp = zeros(DATA.grids.Ny,DATA.grids.Nx,Nz);
+ for iz = 1:DATA.grids.Nz
vx_ = real((ifourier_GENE(-1i*KY.*(DATA.PHI (:,:,iz,FRAMES(it))))));
ni_ = real((ifourier_GENE( DATA.DENS_I(:,:,iz,FRAMES(it)))));
Ti_ = real((ifourier_GENE( DATA.TEMP_I(:,:,iz,FRAMES(it)))));
qx_ = vx_.*ni_.*Ti_;
-% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.Ny,DATA.Nx))),2));
- tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.Ny,DATA.Nx))));
+% tmp(:,:,iz) = abs(fftshift((squeeze(fft2(gx_,DATA.grids.Ny,DATA.grids.Nx))),2));
+ tmp(:,:,iz) = abs((squeeze(fft2(qx_,DATA.grids.Ny,DATA.grids.Nx))));
end
- FLD_(:,:,it)= squeeze(compr(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:)));
+ FLD_(:,:,:,it)= squeeze(tmp(1:DATA.grids.Nky,1:DATA.grids.Nkx,:));
end
case 'f_i'
SKIP_COMP = 1;
diff --git a/src/solve_EM_fields.F90 b/src/solve_EM_fields.F90
index 936ac1ee4df88ee0854a5fe370c44b89a06542c0..2c1e20b72cc0247b99338647ce0ddc8d70958ef3 100644
--- a/src/solve_EM_fields.F90
+++ b/src/solve_EM_fields.F90
@@ -51,7 +51,7 @@ CONTAINS
! Adiabatic charge density (linked to flux surface averaged phi)
! We compute the flux surface average solving a flux surface averaged
! Poisson equation, i.e.
- ! [qi^2(1-sum_j K_j^2)/tau_i] <phi>_psi = <q_i n_i >_psi
+ ! [qi^2/taui (1-sum_j K_j^2)/tau_i] <phi>_psi = <q_i n_i >_psi
! inv_pol_ion^-1 fsa_phi = simpson(Jacobian rho_i ) * iInt_Jacobian
IF (ADIAB_E) THEN
fsa_phi = 0._xp
diff --git a/wk/get_param_from_profiles.m b/wk/get_param_from_profiles.m
index 6376d03511eb590356d6d43c633f33b129535e5b..2a8f461ecf6d421a7d4e83bf4cd77f3d344c3765 100644
--- a/wk/get_param_from_profiles.m
+++ b/wk/get_param_from_profiles.m
@@ -35,4 +35,7 @@ end
params.BETA = 403.0E-5*n_e*T_e/(Bref*Bref);
cref= sqrt(T_e*1000*1.6e-19/mref);
params.EXBRATE= K_w*wExB*1000*Lref/cref;
+ params.ne = n_e;
+ params.Te = T_e;
+ params.Ti = T_i;
end
\ No newline at end of file
diff --git a/wk/parameters/lin_DIIID_data.m b/wk/parameters/lin_DIIID_data.m
index 297f11e3c8565ea031fff4d18e0bcd3850fc0f0e..043de849c780e5a852bc4cff278a69f36461b7fd 100644
--- a/wk/parameters/lin_DIIID_data.m
+++ b/wk/parameters/lin_DIIID_data.m
@@ -46,6 +46,8 @@ Bref = abs(geom.Bref);
mref = m_i;
% Get values at a given location :
[params,profiles] = get_param_from_profiles(prof_folder,rho,Lref,mref,Bref,READPROF);
+params.NUeff_e = params.nuGENE*16/3/sqrt(pi)*geom.q0/geom.trpeps^(3/2);
+params.NUeff_i = params.nuGENE*8/3*sqrt(2)/sqrt(pi)*geom.q0/geom.trpeps^(3/2);
K_Ne = params.K_Ne;
K_Ni = params.K_Ne;
K_Te = params.K_Te;