From c36a00672b4e29efef32d2228ac693145abb3d0f Mon Sep 17 00:00:00 2001
From: Antoine <antoine.hoffmann@epfl.ch>
Date: Fri, 16 Feb 2024 16:13:23 +0100
Subject: [PATCH] update matlab
---
matlab/profiler.m | 62 ++++++++++++++---------------
wk/HEL_lin_scan.m | 16 ++++++--
wk/fast_run.m | 68 ++++++++++++++++++++++++++++++++
wk/get_miller_GENE_py.m | 2 +-
wk/p3_geom_scan_analysis.m | 15 +++----
wk/parameters/lin_Ivanov.m | 25 +++++++-----
wk/plot_toroidal_fluxtube_geom.m | 12 +++---
7 files changed, 142 insertions(+), 58 deletions(-)
create mode 100644 wk/fast_run.m
diff --git a/matlab/profiler.m b/matlab/profiler.m
index 817c939a..1f895d1c 100644
--- a/matlab/profiler.m
+++ b/matlab/profiler.m
@@ -1,23 +1,19 @@
-function [] = profiler(data)
+function [] = profiler(DATADIR,JID)
%% load profiling
% 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=[]; EXBTC=[]; STETC=[]; TS0TC=[];
-for i = 1:numel(data.outfilenames)
- outfilename = data.outfilenames{i};
+% for i = 1:numel(data.outfilenames)
+ outfilename = [DATADIR,'/',sprintf('outputs_%02d.h5',JID)];
CPUTI = [ CPUTI; double(h5readatt(outfilename,'/data/input','cpu_time'))];
DTSIM = [ DTSIM; h5readatt(outfilename,'/data/input/basic','dt')];
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')];
@@ -28,18 +24,18 @@ for i = 1:numel(data.outfilenames)
DIATC = [ DIATC; h5read(outfilename,'/profiler/Tc_diag')];
STETC = [ STETC; h5read(outfilename,'/profiler/Tc_step')];
TS0TC = [ TS0TC; h5read(outfilename,'/profiler/time')];
-end
-CPUTI = CPUTI(end);
+% end
+CPUTI = sum(CPUTI);
DTSIM = mean(DTSIM);
N_T = 13;
-missing_Tc = STETC - RHSTC - ADVTC - GHOTC - CLOTC - EXBTC ...
+MISTC = 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(EXBTC); ...
- diff(CHKTC); diff(DIATC); diff(GRATC); diff(NADTC); diff(missing_Tc)];
+ diff(CHKTC); diff(DIATC); diff(GRATC); diff(NADTC); diff(MISTC)];
TIME_PER_FCT = reshape(TIME_PER_FCT,[numel(TIME_PER_FCT)/N_T,N_T]);
TIME_PER_STEP = sum(TIME_PER_FCT,2);
@@ -57,7 +53,7 @@ checkfield_Ta = mean(diff(CHKTC));
grad_Ta = mean(diff(GRATC));
nadiab_Ta = mean(diff(NADTC));
diag_Ta = mean(diff(DIATC));
-miss_Ta = mean(diff(missing_Tc));
+miss_Ta = mean(diff(MISTC));
total_Ta = mean(diff(total_Tc));
names = {...
'Mrhs';
@@ -79,9 +75,10 @@ Ts_A = [rhs_Ta adv_field_Ta ghost_Ta clos_Ta coll_Ta poisson_Ta...
NSTEP_PER_SAMP= mean(diff(TS0TC))/DTSIM;
%% Plots
+fig = figure;
if 1
+ % subplot(121)
%% Area plot
-fig = figure;
% colors = rand(N_T,3);
% colors = lines(N_T);
colors = distinguishable_colors(N_T);
@@ -100,21 +97,21 @@ for i = 2:numel(x_)
yy_(2*i-1,:) = y_(i,:)/(dx/DTSIM);
yy_(2*i ,:) = y_(i,:)/(dx/DTSIM);
end
-p1 = area(xx_,yy_,'LineStyle','none');
+p1 = area(xx_/DTSIM,yy_,'LineStyle','none');
for i = 1:N_T; p1(i).FaceColor = colors(i,:);
% LEGEND{i} = sprintf('%s t=%1.1e[s] (%0.1f %s)',names{i},Ts_A(i),Ts_A(i)/total_Tc(end)*100,'\%');
LEGEND{i} = [names{i},' $\hat t=$',sprintf('%1.1e[s] (%0.1f %s)',Ts_A(i)/NSTEP_PER_SAMP,Ts_A(i)/total_Ta*100,'\%')];
-end;
+end
legend(LEGEND,'Location','bestoutside');
% legend('Compute RHS','Adv. fields','ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Process', 'Missing')
-xlabel('Sim. Time [$\rho_s/c_s$]'); ylabel('Step Comp. Time [s]')
-xlim([TS0TC(2),TS0TC(end)]);
-ylim([0, 1.1*CPUTI/(TS0TC(end)/DTSIM)])
+xlabel('Simulation step'); ylabel('Comp. time per step [s]')
+xlim([TS0TC(2),TS0TC(end)]./DTSIM);
+ylim([0, 1.5*total_Ta/(dx/DTSIM)])
h_ = floor(CPUTI/3600);
m_ = floor(floor(CPUTI/60)-60*h_);
s_ = CPUTI - 3600*h_ - 60*m_;
-title(sprintf('Gyacomo 2 (%.0f [h] ~%.0f [min] ~%.0f [s])',...
- h_,m_,s_))
+title([DATADIR,sprintf(' (%.0f [h] ~%.0f [min] ~%.0f [s])',...
+ h_,m_,s_)])
hold on
FIGNAME = 'profiler';
% save_figure
@@ -145,20 +142,21 @@ FIGNAME = 'profiler';
end
if 0
+ subplot(122)
%% Histograms
-fig = figure;
-histogram(diff(rhs_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(adv_field_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(ghost_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(clos_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(coll_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(poisson_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(Sapj_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(process_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(checkfield_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
-histogram(diff(diag_Tc)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+% fig = figure;
+histogram(diff(RHSTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(ADVTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(GHOTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(CLOTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(COLTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(POITC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(SAPTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(CHKTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(DIATC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
+histogram(diff(MISTC)/NSTEP_PER_SAMP,'Normalization','probability');hold on
grid on;
-legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Process','Check+sym', 'Diagnos.', 'Missing')
+legend('Compute RHS','Adv. fields','Ghosts comm', 'closure', 'collision','Poisson','Nonlin','Check+sym', 'Diagnos.', 'Missing')
xlabel('Step Comp. Time [s]'); ylabel('')
set(gca,'Xscale','log')
FIGNAME = 'profiler';
diff --git a/wk/HEL_lin_scan.m b/wk/HEL_lin_scan.m
index 9c725524..719fbb66 100644
--- a/wk/HEL_lin_scan.m
+++ b/wk/HEL_lin_scan.m
@@ -1,14 +1,24 @@
-Nsim = 32;
+gyacomodir = pwd; gyacomodir = gyacomodir(1:end-2); % get code directory
+addpath(genpath([gyacomodir,'matlab'])) % ... add
+addpath(genpath([gyacomodir,'matlab/plot'])) % ... add
+addpath(genpath([gyacomodir,'matlab/compute'])) % ... add
+addpath(genpath([gyacomodir,'matlab/load'])) % ... add
+addpath(genpath([gyacomodir,'wk/parameters'])) % ... add
+
+Nsim = 48;
gr_ = zeros(61,Nsim);
i = 1;
-tau_a = logspace(0,-3,Nsim);
+tau_a = logspace(0,-3.5,Nsim);
for i = 1:Nsim
run lin_Ivanov;
TAU = tau_a(i);
NU = 0.1*3/2/TAU/4; % Collision frequency
K_Ti = 0.35*2/TAU; % ion Temperature
- DT = 1e-2;
+ DT = 5e-3;
+ PMAX = 2; JMAX = 1;
+ % PMAX = 2; JMAX = 0;
+ % PMAX = 0; JMAX = 1;
run fast_run;
% growth rates
[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
diff --git a/wk/fast_run.m b/wk/fast_run.m
new file mode 100644
index 00000000..498e00fb
--- /dev/null
+++ b/wk/fast_run.m
@@ -0,0 +1,68 @@
+%% QUICK RUN SCRIPT
+% This script runs a simulation directly from the Matlab framework.
+% It is meant to run only small problems in linear for benchmarking and
+% debugging purposes since it makes Matlab "busy".
+
+%% Set up the paths for the necessary Matlab modules
+gyacomodir = pwd;
+gyacomodir = gyacomodir(1:end-2);
+% mpirun = 'mpirun';
+mpirun = '/opt/homebrew/bin/mpirun'; % for macos
+addpath(genpath([gyacomodir,'matlab'])) % Add matlab folder
+addpath(genpath([gyacomodir,'matlab/plot'])) % Add plot folder
+addpath(genpath([gyacomodir,'matlab/compute'])) % Add compute folder
+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
+default_plots_options
+% EXECNAME = 'gyacomo23_sp_save'; % single precision
+% EXECNAME = 'gyacomo23_sp'; % single precision
+EXECNAME = 'gyacomo23_dp'; % double precision
+% EXECNAME = 'gyacomo23_debug'; % double precision
+%% RUN
+setup
+% system(['rm fort*.90']);
+% Run linear simulation
+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 6 ',gyacomodir,'bin/',EXECNAME,' 3 2 1 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]);
+end
+
+%% quick load of basic data
+filename = [SIMID,'/',PARAMS,'/']; % Create the filename based on SIMID and PARAMS
+LOCALDIR = [gyacomodir,'results/',filename,'/']; % Create the local directory path based on gyacomodir, results directory, and filename
+FIGDIR = LOCALDIR; % Set FIGDIR to the same path as LOCALDIR
+% Load outputs from jobnummin up to jobnummax
+J0 = 0; J1 = 0;
+data = {}; % Initialize data as an empty cell array
+% load grids, inputs, and time traces
+data = compile_results_low_mem(data,LOCALDIR,J0,J1);
+
+if 0 % Activate or not
+%% plot mode evolution and growth rates
+% Load phi
+[data.PHI, data.Ts3D] = compile_results_3D(LOCALDIR,J0,J1,'phi');
+options.NORMALIZED = 0;
+options.TIME = data.Ts3D;
+ % Time window to measure the growth of kx/ky modes
+options.KY_TW = [0.5 1.0]*data.Ts3D(end);
+options.KX_TW = [0.5 1.0]*data.Ts3D(end);
+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
+options.FIELD = 'phi';
+options.SHOWFIG = 1;
+[fig, 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
diff --git a/wk/get_miller_GENE_py.m b/wk/get_miller_GENE_py.m
index 056eae65..4c2c3931 100644
--- a/wk/get_miller_GENE_py.m
+++ b/wk/get_miller_GENE_py.m
@@ -2,7 +2,7 @@ function [variables] = get_miller_GENE_py(prof_folder,rho)
filePath = [prof_folder,'/equilibrium.txt'];
-command = ['python extract_miller_from_eqdsk.py ',...
+command = ['python3 extract_miller_from_eqdsk.py ',...
filePath,' ',num2str(rho),' > tmp.txt'];
system(command);
diff --git a/wk/p3_geom_scan_analysis.m b/wk/p3_geom_scan_analysis.m
index f364ab56..0ae25549 100644
--- a/wk/p3_geom_scan_analysis.m
+++ b/wk/p3_geom_scan_analysis.m
@@ -14,7 +14,7 @@ GETFLUXSURFACE = 0;
% partition= '../results/paper_3/';
% Get the scan directory
-switch 5
+switch 6
case 1 % delta kappa scan
casename = 'DTT rho85';
partition= '/misc/gyacomo23_outputs/paper_3/DTT_rho85_geom_scan/';
@@ -56,20 +56,21 @@ switch 5
t1 = 50; t2 = 150;
case 5 % delta K_T tau=1
casename = 'DIIID rho95 $\tau=1$';
- partition= '/misc/gyacomo23_outputs/paper_3/DIIID_tau_1_rho95_geom_scan/';
- scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
+ partition= '../results/paper_3/DIIID_tau_1_rho95_geom_scan/';
+ % scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
+ scandir = '2_1_delta_RT_scan'; scanname= '(2,1)';
% scandir = '5x3x192x48x32_delta_RT_scan'; scanname= '(4,2)';
% scandir = 'delta_RT_scan_PJ_21'; scanname= '(2,1)';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1.0;
- t1 = 200; t2 = 800;
+ t1 = 200; t2 = 500;
case 6 % delta K_T cold ions
casename = 'DIIID rho95 $\tau=10^{-3}$';
- partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
+ partition= '../results/paper_3/DIIID_cold_ions_rho95_geom_scan/';
scandir = '3x2x192x48x32_delta_RT_scan'; scanname= '(2,1)';
nml1 = 'GEOMETRY'; pnam1 = '$\delta$'; attr1 = 'delta'; pref1 = 0; scale1 =1.0;
- nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1e3;
- t1 = 200; t2 = 480;
+ nml2 = 'SPECIES'; pnam2 = '$R_0/L_T\times T_i/T_e$'; attr2 = 'K_T_'; pref2 = 0.8; scale2 =1e3/2;
+ t1 = 150; t2 = 280;
case 7 % delta s_delta
casename = 'DIIID rho95 $\tau=10^{-3}$';
partition= '/misc/gyacomo23_outputs/paper_3/DIIID_cold_ions_rho95_geom_scan/';
diff --git a/wk/parameters/lin_Ivanov.m b/wk/parameters/lin_Ivanov.m
index a638d755..80dccb8e 100644
--- a/wk/parameters/lin_Ivanov.m
+++ b/wk/parameters/lin_Ivanov.m
@@ -3,12 +3,12 @@ SIMID = 'lin_Ivanov'; % Name of the simulation
%% Set up physical parameters
CLUSTER.TIME = '99:00:00'; % Allocation time hh:mm:ss
-TAU = 1e-3; % e/i temperature ratio
-NU = 0.1/TAU; % Collision frequency
+TAU = 0.001; % e/i temperature ratio
+NU = 1.0*3/2/TAU/4; % Collision frequency
K_Ne = 0*2.22; % ele Density
K_Te = 0*6.96; % ele Temperature
K_Ni = 0*2.22; % ion Density gradient drive
-K_Ti = 0.36*2/TAU; % ion Temperature
+K_Ti = 1.0*2/TAU+0*5*TAU; % ion Temperature
SIGMA_E = 0.0233380; % mass ratio sqrt(m_a/m_i) (correct = 0.0233380)
NA = 1; % number of kinetic species
ADIAB_E = (NA==1); % adiabatic electron model
@@ -20,9 +20,9 @@ J = 1;%P/2;
PMAX = P; % Hermite basis size
JMAX = J; % Laguerre basis size
NX = 2; % real space x-gridpoints
-NY = 40; % real space y-gridpoints
+NY = 120; % real space y-gridpoints
LX = 2*pi/0.1; % Size of the squared frequency domain in x direction
-LY = 2*pi/0.15; % Size of the squared frequency domain in y direction
+LY = 2*pi/0.05; % Size of the squared frequency domain in y direction
NZ = 1; % number of perpendicular planes (parallel grid)
SG = 0; % Staggered z grids option
NEXC = 0; % To extend Lx if needed (Lx = Nexc/(kymin*shear))
@@ -34,15 +34,18 @@ EPS = 0.0; % inverse aspect ratio
Q0 = 0.0; % safety factor
SHEAR = 0.0; % magnetic shear
KAPPA = 1.0; % elongation
+S_KAPPA = 0.0;
DELTA = 0.0; % triangularity
+S_DELTA = 0.0;
ZETA = 0.0; % squareness
+S_ZETA = 0.0;
PARALLEL_BC = 'dirichlet'; % Boundary condition for parallel direction ('dirichlet','periodic','shearless','disconnected')
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
+DT = 1e-3; % Time step
DTSAVE0D = 1; % Sampling time for 0D arrays
DTSAVE2D = -1; % Sampling time for 2D arrays
DTSAVE3D = 1; % Sampling time for 3D arrays
@@ -51,12 +54,12 @@ JOB2LOAD = -1; % Start a new simulation serie
%% OPTIONS
LINEARITY = 'linear'; % activate non-linearity (is cancelled if KXEQ0 = 1)
-CO = 'SG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
-GKCO = 0; % Gyrokinetic operator
+CO = 'DG'; % Collision operator (LB:L.Bernstein, DG:Dougherty, SG:Sugama, LR: Lorentz, LD: Landau)
+GKCO = 1; % Gyrokinetic operator
ABCO = 1; % INTERSPECIES collisions
INIT_ZF = 0; % Initialize zero-field quantities
% HRCY_CLOS = 'max_degree'; % Closure model for higher order moments
-HRCY_CLOS = 'truncation'; % Closure model for higher order moments
+HRCY_CLOS = 'max_degree'; % Closure model for higher order moments
DMAX = -1;
NLIN_CLOS = 'truncation'; % Nonlinear closure model for higher order moments
NMAX = 0;
@@ -93,3 +96,7 @@ BCKGD0 = 0.0; % Initial background
k_gB = 1.0; % Magnetic gradient strength
k_cB = 1.0; % Magnetic curvature strength
COLL_KCUT = 1.0; % Cutoff for collision operator
+PB_PHASE = 0.0;
+ADIAB_I = 0;
+MHD_PD = 0;
+EXBRATE = 0;
\ No newline at end of file
diff --git a/wk/plot_toroidal_fluxtube_geom.m b/wk/plot_toroidal_fluxtube_geom.m
index a43c212c..b6aa8e4e 100644
--- a/wk/plot_toroidal_fluxtube_geom.m
+++ b/wk/plot_toroidal_fluxtube_geom.m
@@ -8,13 +8,13 @@ default_plots_options
Nx = 128;
Ny = 128;
-Nz = 64;
Nturns = 1;
+Nz = 128;
x = linspace(-60,60,Nx)*0.001;
y = linspace(-60,60,Ny)*0.001;
FIELD = ones(Nx,Ny,Nz);
z = linspace(-Nturns*pi,Nturns*pi,Nz);
-N_field_lines = 120;
+N_field_lines = 2;
N_magn_flux_surf = 1;
openangle = pi/3;
phi0 = openangle/2; phi1= 2*pi-openangle/2;
@@ -28,9 +28,9 @@ select = 3;
switch select
case 1
% CBC
- rho = 0.50; drho = 0.1;% Torus minor radius
+ rho = 0.3; drho = 0.1;% Torus minor radius
eps = 0.3;
- q0 = 0.000001; % Flux tube safety factor
+ q0 = 1.4; % Flux tube safety factor
shear = 0.8;
kappa = 1.0;
s_kappa= 0.0;
@@ -171,7 +171,7 @@ figure; set(gcf, 'Position', DIMENSIONS)
end
%plot field lines
if N_field_lines > 0
- theta = linspace(-Nturns*pi, Nturns*pi, 256) ; % Poloidal angle
+ theta = linspace(-Nturns*pi, Nturns*pi, Nturns*256) ; % Poloidal angle
vecfl = Rvec(theta);
plot3(vecfl(1,:),vecfl(2,:),vecfl(3,:),'-b'); hold on;
% Multiple field lines
@@ -183,7 +183,7 @@ figure; set(gcf, 'Position', DIMENSIONS)
y_ = [y_ sin(t)*vecfl(1,:) + cos(t)*vecfl(2,:)];
z_ = [z_ vecfl(3,:)];
end
- plot3(x_,y_,z_,'.','color',[1.0 0.6 0.6]*0.8); hold on;
+ plot3(x_,y_,z_,'-','color',[1.0 0.6 0.6]*0.8); hold on;
end
%plot fluxe tube
if PLT_FTUBE
--
GitLab