diff --git a/.vscode/launch.json b/.vscode/launch.json index afe555a..ace68c0 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -1,59 +1,59 @@ { // Use IntelliSense to learn about possible attributes. // Hover to view descriptions of existing attributes. // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 "version": "0.0.1", "configurations": [ { "name": "Python: Current File", "type": "python", "request": "launch", "program": "${file}", "args": [ "/misc/espic2dsimus/openmp/trajectories.h5" ], "console": "integratedTerminal" }, { "name": "Fortran Launch (GDB)", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/espic2d", "args": ["${workspaceFolder}/wk/espic2d.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", }, { "name": "Fortran test.in (GDB)", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/espic2d", - "args": ["${workspaceFolder}/wk/test_short.in"], + "args": ["${workspaceFolder}/wk/job.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, - "preLaunchTask": "make debug", + "preLaunchTask": "", "miDebuggerPath": "gdb", }, { "name": "Fortran Launch stable.in", "type": "cppdbg", "request": "launch", "program": "${workspaceFolder}/src/espic2d", "args": ["${workspaceFolder}/wk/stable.in"], "stopAtEntry": false, "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", }, { "name": "Intel Debug Attach", "type": "cppdbg", "request": "attach", "processId": "${command:pickProcess}", "program": "${workspaceFolder}/espic2d" } ] } diff --git a/espic2d.code-workspace b/espic2d.code-workspace index 74da27a..66f6e41 100644 --- a/espic2d.code-workspace +++ b/espic2d.code-workspace @@ -1,66 +1,67 @@ { "folders": [ { "path": "." } ], "settings": { "python.pythonPath": "dataanalysis/bin/python3", "workbench.colorCustomizations": { "activityBar.background": "#88d7ea", "activityBar.activeBorder": "#de41bf", "activityBar.foreground": "#15202b", "activityBar.inactiveForeground": "#15202b99", "activityBarBadge.background": "#de41bf", "activityBarBadge.foreground": "#e7e7e7", "titleBar.activeBackground": "#5dc9e2", "titleBar.inactiveBackground": "#5dc9e299", "titleBar.activeForeground": "#15202b", "titleBar.inactiveForeground": "#15202b99", "statusBar.background": "#5dc9e2", "statusBarItem.hoverBackground": "#32bbda", "statusBar.foreground": "#15202b", "panel.border": "#88d7ea", "sideBar.border": "#88d7ea", "editorGroup.border": "#88d7ea", - "tab.activeBorder": "#88d7ea" + "tab.activeBorder": "#88d7ea", + "activityBar.activeBackground": "#88d7ea" }, "peacock.color": "#5dc9e2", "cSpell.enabledLanguageIds": [ "asciidoc", "c", "cpp", "csharp", "css", "fortran", "fortran-modern", "git-commit", "go", "handlebars", "haskell", "html", "jade", "java", "javascript", "javascriptreact", "json", "jsonc", "latex", "less", "markdown", "php", "plaintext", "pug", "python", "restructuredtext", "rust", "scala", "scss", "text", "typescript", "typescriptreact", "yaml", "yml" ] } } \ No newline at end of file diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m index 230a3f1..ac9c157 100644 --- a/matlab/analyse_espicfields.m +++ b/matlab/analyse_espicfields.m @@ -1,343 +1,342 @@ % file='teststablegauss_13_traject2.h5'; % file='unigauss.h5'; % file='Dav_5e14fine_wr+.h5'; % file='stablegauss_8e19fine.h5' % % %file='teststable_Dav1.h5'; % % %close all hidden; % if (~exist('M','var')) % M.file=file; % end % M=load_espic2d(file,M,'all'); t2dlength=size(M.t2d,1) fieldstart=max(1,t2dlength-500);%floor(0.95*t2dlength); deltat=t2dlength-fieldstart; %[~, rgridend]=min(abs(M.rgrid-0.045)); rgridend=sum(M.nnr(1:2)); [~, name, ~] = fileparts(M.file); M.displayenergy %% Radial profile try M.displayrprofile(deltat) catch end %fieldstart=2300; fieldend=2500; dens=mean(M.N(:,:,fieldstart:end),3); pot=mean(M.pot(:,:,fieldstart:end),3); brillratio=2*dens*M.me./(M.eps_0*M.Bz'.^2); [R,Z]=meshgrid(M.rgrid,M.zgrid); Rinv=1./R; -Rinv(:,1)=0; +Rinv(isnan(Rinv))=0; VTHET=mean(M.fluidUTHET(:,:,fieldstart:end),3); omegare=(VTHET.*Rinv'); maxdens=max(max(dens)); % f=figure(); % ax3=gca; % surface(ax3,M.zgrid(1:end-1),M.rgrid(1:end-1),(squeeze(mean(M.Presstens(4,:,:,fieldstart:end),4)))*M.vlight) % xlabel(ax3,'z [m]') % ylabel(ax3,'r [m]') % xlim(ax3,[M.zgrid(1) M.zgrid(end)]) % ylim(ax3,[M.rgrid(1) M.rgrid(50)]) % colormap(ax3,'jet') % c = colorbar(ax3); % c.Label.String= 'thermal v_\theta [m/s]'; % %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; % %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) % title(ax3,'Azimuthal velocity') % %set(ax3,'colorscale','log') % grid(ax3, 'on') % view(ax3,2) % f.PaperOrientation='landscape'; % [~, name, ~] = fileparts(M.file); % f.PaperUnits='centimeters'; % papsize=[16 14]; % f.PaperSize=papsize; % print(f,sprintf('%sfluid_thermtheta',name),'-dpdf','-fillpage') %% Density f=figure('Name', sprintf('%sfluid_dens',name)); ax1=gca; h=surface(ax1,M.zgrid,M.rgrid,dens); set(h,'edgecolor','none'); ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)]) xlim(ax1,[M.zgrid(1) M.zgrid(end)]) xlabel(ax1,'z [m]') ylabel(ax1,'r [m]') title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens))) c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; view(ax1,2) grid on; hold on; f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_dens',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_dens',name)) %% Fieldswide f=figure('Name', sprintf('%s fields',name)); ax1=gca; h=contourf(ax1,M.zgrid,M.rgrid,dens,'Displayname','n_e [m^{-3}]'); hold on; %[c1,h]=contour(ax1,M.zgrid,M.rgrid,pot./1e3,5,'m--','ShowText','off','linewidth',1.5,'Displayname','Equipotentials [kV]'); %clabel(c1,h,'Color','white') Blines=zeros(size(M.Athet)); for i=1:length(M.zgrid) Blines(i,:)=M.Athet(i,:).*M.rgrid'; end zindex=floor(length(M.zgrid/2)); levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),8); contour(ax1,M.zgrid,M.rgrid,Blines',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines'); xlim(ax1,[M.zgrid(1) M.zgrid(end)]) ylim(ax1,[M.rgrid(1) M.rgrid(end)]) legend xlabel(ax1,'z [m]') ylabel(ax1,'r [m]') title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens))) c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; view(ax1,2) %set(h,'edgecolor','none'); grid on; %rectangle('position',[M.zgrid(5) M.rgrid(3) M.zgrid(end-5)-M.zgrid(5) (M.rgrid(rgridend)-M.rgrid(1))],'Edgecolor','white','linestyle','--','linewidth',2) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[18 10]; f.PaperSize=papsize; print(f,sprintf('%sFieldswide',name),'-dpdf','-fillpage') savefig(f,sprintf('%sFieldswide',name)) %% Fields f=figure('Name', sprintf('%s fields',name)); ax1=gca; h=contourf(ax1,M.zgrid,M.rgrid,dens,'Displayname','n_e [m^{-3}]'); hold on; [c1,h]=contour(ax1,M.zgrid,M.rgrid,pot./1e3,'m--','ShowText','on','linewidth',1.5,'Displayname','Equipotentials [kV]'); clabel(c1,h,'Color','white') Blines=zeros(size(M.Athet)); for i=1:length(M.zgrid) Blines(i,:)=M.Athet(i,:).*M.rgrid'; end zindex=floor(length(M.zgrid/2)); levels=logspace( log10(min(min(Blines(:,2:end)))), log10(max(max(Blines(:,:)))),10); contour(ax1,M.zgrid,M.rgrid,Blines',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines'); xlim(ax1,[M.zgrid(1) M.zgrid(end)]) ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)]) legend xlabel(ax1,'z [m]') ylabel(ax1,'r [m]') title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens))) c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; view(ax1,2) %set(h,'edgecolor','none'); grid on; f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[18 10]; f.PaperSize=papsize; print(f,sprintf('%sFields',name),'-dpdf','-fillpage') savefig(f,sprintf('%sFields',name)) %% Brillouin f=figure('Name', sprintf('%s Brillouin',name)); ax1=gca; h=surface(ax1,M.zgrid,M.rgrid,brillratio); set(h,'edgecolor','none'); xlim(ax1,[M.zgrid(1) M.zgrid(end)]) ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)]) xlabel(ax1,'z [m]') ylabel(ax1,'r [m]') title(ax1,'Position') c = colorbar(ax1); c.Label.String= 'Brillouin Ratio'; view(ax1,2) caxis([0 1.2]) title(ax1,sprintf('Brillouin Ratio t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens))) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_Brillouin',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_Brillouin',name)) %% Omegar f=figure('Name', sprintf('%s Omegar',name)); ax3=gca; surface(ax3,M.zgrid,M.rgrid,omegare,'edgecolor','None') xlabel(ax3,'z [m]') ylabel(ax3,'r [m]') xlim(ax3,[M.zgrid(1) M.zgrid(end)]) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) colormap(ax3,'jet') c = colorbar(ax3); c.Label.String= '|\omega_\theta| [1/s]'; %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) title(ax3,sprintf('Azimuthal frequency t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens))) %set(ax3,'colorscale','log') grid(ax3, 'on') view(ax3,2) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_omegar',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_omegar',name)) %% VTHET f=figure('Name', sprintf('%s Vthet',name)); ax3=gca; surface(ax3,M.zgrid,M.rgrid,VTHET,'edgecolor','None') xlabel(ax3,'z [m]') ylabel(ax3,'r [m]') xlim(ax3,[M.zgrid(1) M.zgrid(end)]) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) colormap(ax3,'jet') c = colorbar(ax3); c.Label.String= '|V_\theta| [m/s]'; %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) title(ax3,sprintf('Azimuthal velocity t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens))) %set(ax3,'colorscale','log') grid(ax3, 'on') view(ax3,2) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_vthet',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_vthet',name)) -%% VDRIFT +%% diff vdrift f=figure('Name', sprintf('%s V drift',name)); ax3=gca; vdrift=-mean(M.Er(:,:,fieldstart:end),3).*M.Bz'./(M.B.^2)'; -v=abs(vdrift-VTHET); +v=abs(vdrift-VTHET)./VTHET; surface(ax3,M.zgrid,M.rgrid,v,'edgecolor','None') xlabel(ax3,'z [m]') ylabel(ax3,'r [m]') xlim(ax3,[M.zgrid(1) M.zgrid(end)]) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) colormap(ax3,'jet') c = colorbar(ax3); set(ax3,'zscale','log') set(ax3,'colorscale','log') -c.Label.String= '|V_{drift}| [m/s]'; +c.Label.String= '(V_\theta-V_E)/V_\theta'; %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) -title(ax3,sprintf('Azimuthal drift velocity t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens))) +title(ax3,sprintf('(V_\\theta-V_E)/V_\\theta velocity t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(end),double(maxdens))) %set(ax3,'colorscale','log') grid(ax3, 'on') view(ax3,2) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%sfluid_EB',name),'-dpdf','-fillpage') savefig(f,sprintf('%sfluid_EB',name)) %% HP if(M.R.nt>=2) taveragestart=max(floor(0.8*size(M.tpart,1)),2); M.displayHP(taveragestart); end %% 0D diagnostics BrillouinRatioth=2*M.omepe^2/M.omece^2 BrillouinMax=max(max(brillratio)) - NpartsTrapped=mean(M.nbparts(end-10:end))*M.msim/M.me dblengthiter=fieldstart:length(M.t2d); Pr=squeeze(M.Presstens(1,:,:,dblengthiter)); Pz=squeeze(M.Presstens(6,:,:,dblengthiter)); enddens=M.N(:,:,dblengthiter); meanPr=mean(Pr,3); meanPz=mean(Pz,3); Tr=Pr./enddens; Tz=Pz./enddens; Debye_length=sqrt(abs(min(min(min(mean(Tr(Tr>0)./enddens(Tr>0),3))),min(min(mean(Tz(Tz>0)./enddens(Tz>0),3)))))*M.eps_0/M.qe^2) %% Debye length f=figure('Name', sprintf('%s Debye',name)); subplot(2,1,1) surface(M.zgrid,M.rgrid,sqrt(mean(abs(Tr./enddens*M.eps_0/M.qe^2),3))) colorbar f.PaperOrientation='landscape'; xlabel('z [m]') ylabel('r [m]') c = colorbar; c.Label.String= '\lambda_r [m]'; subplot(2,1,2) surface(M.zgrid,M.rgrid,sqrt(mean(abs(Tz./enddens*M.eps_0/M.qe^2),3))) colorbar f.PaperOrientation='landscape'; xlabel('z [m]') ylabel('r [m]') c = colorbar; c.Label.String= '\lambda_z [m]'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%s_dblength',name),'-dpdf','-fillpage') savefig(f,sprintf('%s_dblength',name)) %% Temperature f=figure('Name', sprintf('%s Temperature',name)); ax1=subplot(2,1,1); title('Radial temperature') surface(M.zgrid,M.rgrid,mean(Tr,3)/M.qe,'edgecolor','none') colormap('jet') c = colorbar; c.Label.String= 'T_{er} [eV]'; templimits1=caxis; ax2=subplot(2,1,2); title('Axial tempreature') surface(M.zgrid,M.rgrid,mean(Tz,3)/M.qe,'edgecolor','none') colormap('jet') c = colorbar; c.Label.String= 'T_{ez} [eV]'; maxtemp=max([templimits1,caxis]); caxis(ax1,[0 maxtemp]) caxis(ax2,[0 maxtemp]) f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; print(f,sprintf('%s_temp',name),'-dpdf','-fillpage') savefig(f,sprintf('%s_temp',name)) %% Wells M.displaypotentialwell(t2dlength); diff --git a/matlab/espic2dhdf5.m b/matlab/espic2dhdf5.m index d38f80d..409f37a 100644 --- a/matlab/espic2dhdf5.m +++ b/matlab/espic2dhdf5.m @@ -1,849 +1,858 @@ classdef espic2dhdf5 %espic2dhdf5 General class used to treat hdf5 result files of espic2d code % A result file is loaded with a call to M=espic2dhdf5(filename) where filename is the relative or absolute file path % after loading, several quantities and composite diagnostics such as moments of the distribution function or individual particles % quantities can be accessed. properties filename name folder fullpath timestamp info t0d t1d t2d tpart it0 it1 it2 %% Physical constants vlight=299792458; qe=1.60217662E-19; me=9.109383E-31; eps_0=8.85418781762E-12; kb=1.38064852E-23; %% Run parameters dt % simulation time step nrun % number of time steps simulated nlres nlsave nlclassical % Was the equation of motion solved in the classical framework nlPhis % Was the self-consistent electric field computed nz % number of intervals in the z direction for the grid nnr % number of intervals in the r direction for the grid for each of the 3 mesh regions lz % physical axial dimension of the simulation space nplasma % Number of initial macro particles potinn % Normalized electric potential at the coaxial insert potout % Normalized electric potential at the cylinder surface B0 % Normalization for the magnetic field Rcurv % Magnetic mirror ratio width % Magnetic mirror length n0 % Initial particle density in case of old particle loading temp % Initial particle temperature in case of old particle loading femorder % finite element method order in z and r direction ngauss % Order of the Gauss integration method for the FEM plasmadim % initial dimensions of the plasma for the old particle loading system radii % Radial limits of the three mesh regions coarse,fine,coarse H0 % Initial particle Energy for Davidsons distribution function P0 % Initial particle Angular momentum for Davidsons distribution function normalized % Are the parts quantities normalized in the h5 file nbspecies % Number of species simulated %% Frequencies omepe % Reference plasma frequency used for normalization omece % Reference cyclotronic frequency for normalization %% Normalizations tnorm % Time normalization rnorm % Dimension normalization bnorm % Magnetic field normalization enorm % Electric field normalization phinorm % Electric potential normalization vnorm % Velocity normalization %% Grid data rgrid % Radial grid position points zgrid % Axial grid position points dz % Axial grid step dr % Radial grid step for the three mesh regions CellVol % Volume of the cell used for density calculation %% Magnetic field Br % Radial magnetic field Bz % Axial magnetic field Athet % Azimuthal component of the Magnetic potential vector rAthet % r*Athet used for the representation of magnetic field lines B % Magnetic field amplitude %% Energies epot % Time evolution of the particles potential energy ekin % Time evolution of the particles kinetic energy etot % Time evolution of the particles total energy etot0 % Time evolution of the reference particle total energy eerr % Time evolution of the error on the energy conservation %% 2D time data evaluated on grid points N % main specie Density fluidUR % main specie radial fluid velocity fluidUZ % main specie axial fluid velocity fluidUTHET % main specie azimuthal fluid velocity pot % Electric potential Er % Radial electric field Ez % Axial electric field Presstens % Pressure tensor %% Splines knotsr % Spline radial knots knotsz % Spline axial knots %% Particle parameters weight % Macro particle numerical weight of the main specie qsim % Macro particle charge msim % Macro particle mass nbparts % Time evolution of the number of simulated particles partepot % Electric potential at the particles positions R % Particles radial position Z % Particles axial position Rindex % Particles radial grid index Zindex % Particles axial grid index partindex % Particles unique id for tracing trajectories VR % Particles radial velocity VZ % Particles axial velocity VTHET % Particles azimuthal velocity THET % Particles azimuthal position species % Array containing the other simulated species %% Celldiag celldiag % Array containing the cell diagnostic data nbcelldiag % Total number of cell diagnostics end methods function file=file(obj) % returns the h5 file name file=obj.filename; end function obj = espic2dhdf5(filename,readparts) % Reads the new result file filename and read the parts data if readparts==true % Try catch are there for compatibility with older simulation files filedata=dir(filename); if (isempty(filedata)) error("File: ""%s"" doesn't exist",filename) end obj.folder=filedata.folder; obj.filename=filename; [~, obj.name, ~] = fileparts(obj.filename); obj.fullpath=[obj.folder,'/',obj.filename]; obj.timestamp=filedata.date; if nargin==1 readparts=true; end %obj.info=h5info(filename); %% Read the run parameters obj.dt = h5readatt(obj.fullpath,'/data/input.00/','dt'); obj.nrun = h5readatt(obj.fullpath,'/data/input.00/','nrun'); obj.nlres = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlres'),'y'); obj.nlsave = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlsave'),'y'); obj.nlclassical =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlclassical'),'y'); obj.nlPhis =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlPhis'),'y'); obj.nz = h5readatt(obj.fullpath,'/data/input.00/','nz'); obj.nnr = h5read(obj.fullpath,'/data/input.00/nnr'); obj.lz = h5read(obj.fullpath,'/data/input.00/lz'); obj.qsim = h5readatt(obj.fullpath,'/data/input.00/','qsim'); obj.msim = h5readatt(obj.fullpath,'/data/input.00/','msim'); try obj.weight=h5readatt(obj.fullpath,'/data/part/','weight'); catch obj.weight=obj.msim/obj.me; end obj.nplasma = h5readatt(obj.fullpath,'/data/input.00/','nplasma'); obj.potinn = h5readatt(obj.fullpath,'/data/input.00/','potinn'); obj.potout = h5readatt(obj.fullpath,'/data/input.00/','potout'); obj.B0 = h5readatt(obj.fullpath,'/data/input.00/','B0'); obj.Rcurv = h5readatt(obj.fullpath,'/data/input.00/','Rcurv'); obj.width = h5readatt(obj.fullpath,'/data/input.00/','width'); obj.n0 = h5readatt(obj.fullpath,'/data/input.00/','n0'); obj.temp = h5readatt(obj.fullpath,'/data/input.00/','temp'); try obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0d'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2d'); obj.it2 = h5readatt(obj.fullpath,'/data/input.00/','itparts'); catch obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it1'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2'); end try obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies'); obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y'); catch obj.nbspecies=1; obj.normalized=false; end try obj.nbcelldiag=h5readatt(obj.fullpath,'/data/celldiag/','nbcelldiag'); catch obj.nbcelldiag=0; end obj.omepe=sqrt(abs(obj.n0)*obj.qe^2/(obj.me*obj.eps_0)); obj.omece=obj.qe*obj.B0/obj.me; obj.nbparts= h5read(obj.fullpath, '/data/var0d/nbparts'); try obj.H0 = h5read(obj.fullpath,'/data/input.00/H0'); obj.P0 = h5read(obj.fullpath,'/data/input.00/P0'); catch obj.H0=3.2e-14; obj.P0=8.66e-25; end % Normalizations obj.tnorm=1/obj.omepe; obj.rnorm=obj.vlight*obj.tnorm; obj.bnorm=obj.B0; obj.enorm=obj.vlight*obj.bnorm; obj.phinorm=obj.enorm*obj.rnorm; obj.vnorm=obj.vlight; % Grid data obj.rgrid= h5read(obj.fullpath, '/data/var1d/rgrid')*obj.rnorm; obj.zgrid= h5read(obj.fullpath, '/data/var1d/zgrid')*obj.rnorm; obj.dz=(obj.zgrid(end)-obj.zgrid(1))/double(obj.nz); rid=1; for i=1:length(obj.nnr) obj.dr(i)=(obj.rgrid(sum(obj.nnr(1:i))+1)-obj.rgrid(rid))/double(obj.nnr(i)); rid=rid+obj.nnr(i); end Br = h5read(obj.fullpath,'/data/fields/Br')*obj.bnorm; obj.Br= reshape(Br,length(obj.zgrid),length(obj.rgrid)); Bz = h5read(obj.fullpath,'/data/fields/Bz')*obj.bnorm; obj.Bz= reshape(Bz,length(obj.zgrid),length(obj.rgrid)); try Atheta = h5read(obj.fullpath,'/data/fields/Athet')*obj.bnorm; obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid)); [rmeshgrid,~]=meshgrid(obj.rgrid,obj.zgrid); obj.rAthet=(rmeshgrid.*obj.Athet)'; catch end obj.B=sqrt(obj.Bz.^2+obj.Br.^2); clear Br Bz try obj.t0d=h5read(obj.fullpath,'/data/var0d/time'); catch obj.t0d=obj.dt.*double(0:length(obj.epot)-1); end obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder'); obj.ngauss = h5read(obj.fullpath,'/data/input.00/ngauss'); obj.plasmadim = h5read(obj.fullpath,'/data/input.00/plasmadim'); obj.radii = h5read(obj.fullpath,'/data/input.00/radii'); obj.epot = h5read(obj.fullpath,'/data/var0d/epot'); obj.ekin = h5read(obj.fullpath,'/data/var0d/ekin'); obj.etot = h5read(obj.fullpath,'/data/var0d/etot'); try obj.etot0 = h5read(obj.fullpath,'/data/var0d/etot0'); obj.eerr = obj.etot-obj.etot0; catch obj.eerr = obj.etot-obj.etot(2); end if(obj.normalized) obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,1); obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,1); obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,1); else obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,obj.phinorm); obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,obj.enorm); obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,obj.enorm); end try obj.t2d = h5read(obj.fullpath,'/data/fields/time'); catch info=h5info(obj.fullpath,'/data/fields/partdensity'); obj.t2d=obj.dt*(0:info.objspace.Size(2)-1)*double(obj.it1); end try info=h5info(obj.fullpath,'/data/fields/moments'); obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder'); kr=obj.femorder(2)+1; obj.knotsr=augknt(obj.rgrid,kr); kz=obj.femorder(1)+1; obj.knotsz=augknt(obj.zgrid,kz); try obj.CellVol= reshape(h5read(obj.fullpath,'/data/fields/volume'),length(obj.knotsz)-kz,length(obj.knotsr)-kr); obj.CellVol=permute(obj.CellVol,[2,1,3])*obj.rnorm^3; catch zvol=fnder(spmak(obj.knotsz,ones(1,length(obj.knotsz)-kz)), -1 ); rvol=fnder(spmak(obj.knotsr,2*pi*[obj.rgrid' 2*obj.rgrid(end)-obj.rgrid(end-1)]), -1 ); ZVol=diff(fnval(zvol,obj.knotsz)); RVol=diff(fnval(rvol,obj.knotsr)); obj.CellVol=RVol(3:end-1)*ZVol(3:end-1)'; obj.CellVol=padarray(obj.CellVol,[1,1],'replicate','post'); end if(obj.normalized) obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, 1); else obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), 1); end obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 2); obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 3); obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 4); if(obj.normalized) obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me); else obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim); end catch obj.CellVol=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)'; obj.CellVol=obj.CellVol'; obj.N=griddensity(obj.fullpath, '/data/fields/partdensity', sum(obj.nnr)+1, obj.nz+1, obj.CellVol, abs(obj.qsim/obj.qe), true); obj.fluidUR=gridquantity(obj.fullpath, '/data/fields/fluidur', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); obj.fluidUTHET=gridquantity(obj.fullpath, '/data/fields/fluiduthet', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); obj.fluidUZ=gridquantity(obj.fullpath, '/data/fields/fluiduz', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); end if(readparts) if(obj.normalized) obj.R = h5partsquantity(obj.fullpath,'/data/part','R'); obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z'); else obj.R = h5partsquantity(obj.fullpath,'/data/part','R',obj.rnorm); obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z',obj.rnorm); end try obj.THET = h5partsquantity(obj.fullpath,'/data/part','THET'); catch clear obj.THET end try obj.Rindex=h5partsquantity(obj.fullpath,'/data/part','Rindex'); obj.Zindex=h5partsquantity(obj.fullpath,'/data/part','Zindex'); catch clear obj.Rindex obj.Zindex end obj.VR = h5partsquantity(obj.fullpath,'/data/part','UR',obj.vnorm); obj.VZ = h5partsquantity(obj.fullpath,'/data/part','UZ',obj.vnorm); obj.VTHET= h5partsquantity(obj.fullpath,'/data/part','UTHET',obj.vnorm); if(obj.normalized) obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe); else obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe*obj.phinorm); end try obj.partindex = h5partsquantity(obj.fullpath,'/data/part/','partindex'); catch end end try obj.tpart = h5read(obj.fullpath,'/data/part/time'); catch obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2); end if(obj.nbspecies >1) obj.species=h5parts.empty; for i=2:obj.nbspecies nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/part/%2d',i))); if (nbparts(1)>0) obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj); end end end if(obj.nbcelldiag > 0) obj.celldiag=h5parts.empty; for i=1:obj.nbcelldiag nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/celldiag/%02d',i))); if (nbparts(1)>0) obj.celldiag(i)=h5parts(obj.fullpath,sprintf('/data/celldiag/%02d',i),obj); end end end end function Atheta=Atheta(obj,R,Z) halflz=(obj.zgrid(end)+obj.zgrid(1))/2; Atheta=0.5*obj.B0*(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)... .*besseli(1,2*pi*R/obj.width).*cos(2*pi*(Z-halflz)/obj.width)); end function quantity=H(obj,indices) %% computes the total energy for the main specie particle indices{1} at time indices{2} if indices{1}==':' p=1:obj.VR.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end if size(indices,1)>2 track=indices{3}; else track=false; end quantity=0.5*obj.me*(obj.VR(p,t,track).^2+obj.VTHET(p,t,track).^2+obj.VZ(p,t,track).^2)+obj.partepot(p,t,track); end function quantity=P(obj,indices) %% computes the canonical angular momentum for the main specie particle indices{1} at time indices{2} if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end if size(indices,1)>2 track=indices{3}; else track=false; end quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.me+sign(obj.qsim)*obj.qe*obj.Atheta(obj.R(p,t,track),obj.Z(p,t,track))); end function quantity=Vpar(obj,varargin) %Vpar Computes the parallel velocity for the main specie particle indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end if size(indices,1)>2 track=indices{3}; else track=false; end Zind=obj.Zindex(p,t,track)+1; Rind=obj.Rindex(p,t,track)+1; posind=sub2ind(size(obj.B),Zind,Rind); costhet=obj.Bz(posind)./obj.B(posind); sinthet=obj.Br(posind)./obj.B(posind); quantity=obj.VR(p,t,track).*sinthet+obj.VZ(p,t,track).*costhet; end function quantity=Vperp(obj,varargin) %Vperp Computes the perpendicular velocity in the guidind center reference frame, % for the main specie particle indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if indices{1}==':' p=1:obj.R.nparts; else p=indices{1}; end if indices{2}==':' t=1:length(obj.tpart); else t=indices{2}; end if size(indices,1)>2 track=indices{3}; else track=false; end Zind=obj.Zindex(p,t,track)+1; Rind=obj.Rindex(p,t,track)+1; posind=sub2ind(size(obj.B),Zind,Rind); costhet=obj.Bz(posind)./obj.B(posind); sinthet=obj.Br(posind)./obj.B(posind); Vdrift=zeros(size(Zind)); for j=1:length(t) tfield=find(obj.t2d==obj.tpart(t(j))); timeEr=obj.Er(:,:,tfield); timeEz=obj.Ez(:,:,tfield); posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); Vdrift(:,j)=(timeEz(posindE).*obj.Br(posind(:,j))-timeEr(posindE).*obj.Bz(posind(:,j)))./obj.B(posind(:,j)).^2; end quantity=sqrt((obj.VTHET(p,t,track)-Vdrift).^2+(obj.VR(p,t,track).*costhet-obj.VZ(p,t,track).*sinthet).^2); end function displaypsi(obj,deltat,Davidson) %% plot the initial and final radial profile at position z=0 and show the normalized enveloppe function Psi % relevant for Davidson annular distribution function f=figure('Name', sprintf('%s Psi',obj.name)); f.Name= sprintf('%s Psi',obj.name); zpos=floor(length(obj.zgrid)/2); tinit=1; tend=length(obj.t2d); deltat=cell2mat(deltat); if(obj.R.nt<2) H0=obj.H0; P0=obj.P0; else H0=mean(H(obj,{1:obj.VR.nparts,obj.VR.nt,false})); P0=mean(P(obj,{1:obj.VR.nparts,obj.VR.nt,false})); end lw=1.5; Mirrorratio=(obj.Rcurv-1)/(obj.Rcurv+1); locpot=mean(obj.pot(:,zpos,tend-deltat:tend),3); psi=1+obj.qe*locpot(:)/H0-1/(2*obj.me*H0)*(P0./obj.rgrid+obj.qe*0.5*obj.B0.*(obj.rgrid-obj.width/pi*Mirrorratio*cos(2*pi*obj.zgrid(zpos)/obj.width)*besseli(1,2*pi*obj.rgrid/obj.width))).^2; locdens=mean(obj.N(:,zpos,tend-deltat:tend),3); [maxn,In]=max(locdens);%M.N(:,zpos,tinit)); plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) hold on plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end)*maxn*obj.rgrid(In),'DisplayName','N=c*1/r','linewidth',lw) xlabel('r [m]') ylabel('n_e [m^{-3}]') I=find(psi>0); if (length(I)>1) I=[I(1)-2; I(1)-1; I; I(end)+1; I(end)+2]; else I=obj.nnr(1):length(psi); end rq=linspace(obj.rgrid(max(I(1),1)),obj.rgrid(I(end)),500); psiinterp=interp1(obj.rgrid(I),psi(I),rq,'pchip'); zeroindices=find(diff(psiinterp>=0),2); maxpsiinterp=max(psiinterp); plot(rq,maxn*psiinterp/abs(maxpsiinterp),'Displayname','normalized \Psi [a.u.]','linewidth',lw) ylim([0 inf]) for i=1:length(zeroindices) border=plot([rq(zeroindices(i)) rq(zeroindices(i))],[0 obj.rgrid(In)/obj.rgrid(In-2)*maxn],'k--','linewidth',lw); set(get(get(border,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); end legend xlim([0 0.02]) grid on title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%sPsi',obj.name),[15 10]); end function f=displayrprofile(obj,deltat,zpos) %% plot the initial and final radial profile at the middle of the simulation space f=figure('Name', sprintf('%s Prof',obj.name)); if nargin < 3 zpos=floor(length(obj.zgrid)/2); end tinit=1; tend=length(obj.t2d); if(iscell(deltat)) deltat=cell2mat(deltat); end lw=1.5; locdens=mean(obj.N(:,zpos,tend-deltat:tend),3); + Rinv=1./obj.rgrid; + Rinv(isnan(Rinv))=0; + VTHET=mean(obj.fluidUTHET(:,zpos,tend-deltat:tend),3); + omegare=(VTHET.*Rinv); plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) hold on plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) xlabel('r [m]') ylabel('n_e [m^{-3}]') legend xlim([obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))]) grid on + + yyaxis right + plot(obj.rgrid,omegare,'DisplayName',sprintf('<\\omega_{re}> t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) + ylabel('\omega_{re} [1/s]') + title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%srProf',obj.name),[15 10]); end function changed=ischanged(obj) %ischanged Check if the file has been changed and some data must be reloaded try filedata=dir(obj.fullpath); checkedtimestamp=filedata.date; if (max(checkedtimestamp > obj.timestamp) ) changed=true; return end changed=false; return catch changed=true; return end end function displayenergy(obj) %% Plot the time evolution of the system energy and number of simulated macro particles tmin=2; tmax=length(obj.ekin); f=figure('Name', sprintf('%s Energy',obj.name)); subplot(2,1,1) plot(obj.t0d(tmin:tmax),obj.ekin(tmin:tmax),'o-',... obj.t0d(tmin:tmax),obj.epot(tmin:tmax),'d-',... obj.t0d(tmin:tmax),obj.etot(tmin:tmax),'h-',... obj.t0d(tmin:tmax),obj.eerr(tmin:tmax),'x--') legend('ekin', 'epot', 'etot','eerr') xlabel('Time [s]') ylabel('Energies [J]') grid on xlimits=xlim(); subplot(2,1,2) try semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'h-') catch semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)/obj.etot(2)),'h-') end xlabel('t [s]') ylabel('E_{err}/E_{tot}') xlim(xlimits) grid on try yyaxis right plot(obj.t0d(tmin:tmax),abs(obj.nbparts(tmin:tmax)./obj.nbparts(1)*100),'d--') ylabel('Nparts %') %ylim([0 110]) catch end obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name)); end - function trappedParticles(obj) - %% Plot the time evolution of the number of physical trapped particles in the main specie + function SimParticles(obj) + %% Plot the time evolution of the number of simulated physical particles in the main specie f=figure('Name', sprintf('%s Trapped particles',obj.name)); plot(obj.t0d,obj.nbparts*obj.weight) xlabel('t [s]') ylabel('N particles') obj.savegraph(f,sprintf('%s/%sntrapped',obj.folder,obj.name)); end function fighandle=savegraph(obj, fighandle, name, papsize) %% Saves the given figure as a pdf and a .fig fighandle.PaperUnits='centimeters'; if (nargin < 4) papsize=[14 16]; end fighandle.PaperSize=papsize; print(fighandle,name,'-dpdf','-fillpage') savefig(fighandle,name) end function displayHP(obj,tstart) % Plot the histogramm of the total energy and canonical angular momentum at time tstart and % end time of the simulation over the full simulation space for the main specie if(iscell(tstart)) tstart=cell2mat(tstart); end if(obj.R.nt>=2) tstart=obj.R.nt; f=figure('Name', sprintf('%s HP',obj.name)); legtext=sprintf("t=%2.1f - %2.1f [ns]",obj.tpart(tstart)*1e9,obj.tpart(end)*1e9); subplot(1,2,1) partsmax=min(obj.nbparts(end),obj.R.nparts); Hloc=H(obj,{1:obj.nbparts(1),1,false}); h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9)); hold on Hloc=H(obj,{1:partsmax,obj.R.nt,false}); %,'Binwidth',h1.BinWidth h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',legtext); ylabel('counts') xlabel('H [J]') legend subplot(1,2,2) Ploc=P(obj,{1:obj.nbparts(1),1,false}); h2=histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9)); hold on Ploc=P(obj,{1:partsmax,obj.R.nt,false}); histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',legtext); ylabel('counts') xlabel('P [kg\cdotm^2\cdots^{-1}]') %clear P %clear H legend %xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)]) obj.savegraph(f,sprintf('%s/%sParts_HP',obj.folder,obj.name)); end end function displayaveragetemp(obj) % Computes and show the particles average temperature as a function of time f=figure('Name',sprintf('%s potinn=%f part temperature',obj.name,obj.potinn*obj.phinorm)); vr2=obj.VR(:,:,false); vr2=mean(vr2.^2,1)-mean(vr2,1).^2; vz2=obj.VZ(:,:,false); vz2=mean(vz2.^2,1)-mean(vz2,1).^2; vthet2=obj.VTHET(:,:,false); vthet2=mean(vthet2.^2,1)-mean(vthet2,1).^2; plot(obj.tpart,0.5*obj.me*vr2/obj.qe,'displayname','T_r') hold on plot(obj.tpart,0.5*obj.me*vz2/obj.qe,'displayname','T_z') plot(obj.tpart,0.5*obj.me*vthet2/obj.qe,'displayname','T_{thet}') xlabel('time [s]') ylabel('T [eV]') title(sprintf('\\phi_a=%.1f kV \\phi_b=%.1f kV R=%.1f',obj.potinn*obj.phinorm/1e3,obj.potout*obj.phinorm/1e3,obj.Rcurv)) legend grid obj.savegraph(f,sprintf('%s/%s_partstemp',obj.folder,obj.name)); end function Gamma=Axialflux(obj,timestep,zpos) % Computes the axial particle flux n*Uz at timestep timestep and axial position zpos Gamma=obj.fluidUZ(:,zpos,timestep).*obj.N(:,zpos,timestep); end function I=OutCurrents(obj,timestep) % Computes the Outgoing currens at the simulation axial boundaries at timestep timestep % This is simply the surface integral of the axial flux I=zeros(2,length(timestep)); flux=obj.Axialflux(timestep,[1 obj.nz+1]); I=squeeze(trapz(obj.rgrid,flux.*obj.rgrid)*2*pi*obj.qsim/obj.weight); end function displayCurrentsevol(obj,timesteps) % Computes and display the time evolution of the outgoing currents at timesteps timesteps if nargin<2 timesteps=1:length(obj.t2d); end currents=obj.OutCurrents(timesteps); f=figure('Name',sprintf('%s Currents',obj.name)); plot(obj.t2d(timesteps),currents(1,:),'Displayname',sprintf('z=%.2f cm',obj.zgrid(1)*100)); hold on plot(obj.t2d(timesteps),-currents(2,:),'Displayname',sprintf('z=%.2f cm',obj.zgrid(end)*100)); legend('location','east') xlabel('time [s]') ylabel('I [A]') grid on title(sprintf('\\phi_b-\\phi_a=%.2g kV, R=%.1f',(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv)) obj.savegraph(f,sprintf('%s/%s_outCurrents',obj.folder,obj.name)); end function model=potentialwellmodel(obj,timestep) % Computes the potential well at the given timestep and return the model to be able to % interpolate either using grid coordinates or magnetic field line coordinates if iscell(timestep) timestep=cell2mat(timestep); end Phi=-obj.pot(:,:,timestep); contpoints=contourc(obj.zgrid,obj.rgrid,obj.rAthet,obj.rAthet(:,1)'); [x,y,zcont]=C2xyz(contpoints); k=1; zdiff=[diff(zcont),0]; potfinal=zeros(numel(cell2mat(x)),length(timestep)); pot=cell(1,size(x,2)); rmin=obj.rgrid(1); rmax=obj.rgrid(end); rathet=x; for i=1:length(timestep) locPhi=Phi(:,:,i); for j=1:size(x,2) %for i=1:size(pot,1) xloc=x{j}; yloc=y{j}; % xloc=xloc(ylocrmin); % yloc=yloc(ylocrmin); % x{j}=xloc; % y{j}=yloc; rathet{j}=zcont(j)*ones(1,length(xloc)); if(length(xloc)>=1) pot{j}=interp2(obj.zgrid,obj.rgrid,locPhi,xloc,yloc,'makima'); pot{j}=pot{j}-max(pot{j}); end k=k+(zdiff(j)~=0); end potfinal(:,i)=cell2mat(pot); end model.z=cell2mat(x); model.r=cell2mat(y); model.pot=potfinal; model.rathet=cell2mat(rathet); end function [pot] = PotentialWell(obj,timestep) % PotentialWell Computes the potential well at the given timestep on the grid points model=obj.potentialwellmodel(timestep); z=model.z; r=model.r; modpot=model.pot; rathet=model.rathet; [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); pot=zeros(length(obj.zgrid),length(obj.rgrid),length(fieldstep)); for i=1:length(fieldstep) pot(:,:,i)=griddata(z,r,modpot(i,:),Zmesh,Rmesh); end end function displaypotentialwell(obj,timestep,rcoord) % Display the potential well at timestep timestep % if rcoord is true, the potential is evaluated at grid points in r,z coordinates % if false, the potential is evaluated at grid points in magnetic field line coordinates if iscell(timestep) timestep=cell2mat(timestep); end if nargin <3 rcoord=true; end f=figure('Name',sprintf('%f Potential well',obj.name)); model=obj.potentialwellmodel(timestep); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; title(sprintf('Potential well t=%1.2f [ns]',obj.t2d(timestep)*1e9)) if rcoord [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); pot=griddata(z,r,pot,Zmesh,Rmesh); surface(obj.zgrid(1:end),obj.rgrid(1:end),pot(1:end,1:end),'edgecolor','none') xlabel('z [m]') ylabel('r [m]') xlim([obj.zgrid(1) obj.zgrid(end)]) ylim([obj.rgrid(1) obj.rgrid(end)]) else [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rAthet(:,1)); pot=griddata(z,rathet,pot,Zmesh,Rmesh); surface(obj.zgrid(1:end),obj.rAthet(:,1),pot(1:end,1:end),'edgecolor','none') ylabel('rA_\theta [Tm^2]') xlabel('z [m]') xlim([obj.zgrid(1) obj.zgrid(end)]) ylim([obj.rAthet(1,1) obj.rAthet(sum(obj.nnr(1:2)),1)]) end c=colorbar; colormap('jet'); c.Label.String= 'depth [eV]'; obj.savegraph(f,sprintf('%s/%s_well',obj.folder,obj.name)); end function Epar = Epar(obj,fieldstep) % Computes the electric field component parallel to the magnetic field line Epar=obj.Er(:,:,fieldstep).*(obj.Br./obj.B)' + (obj.Bz./obj.B)'.*obj.Ez(:,:,fieldstep); end function Eperp = Eperp(obj,fieldstep) % Computes the electric field component perpendicular to the magnetic field line Eperp=obj.Er(:,:,fieldstep).*(obj.Bz./obj.B)' - (obj.Br./obj.B)'.*obj.Ez(:,:,fieldstep); end end end diff --git a/src/.depend b/src/.depend index cce7a7c..6f06fa2 100644 --- a/src/.depend +++ b/src/.depend @@ -1,36 +1,37 @@ auxval.o : auxval.f90 beam_mod.o fields_mod.o basic_mod.o constants.o auxval__genmod.o : auxval__genmod.f90 basic_mod.o : basic_mod.f90 mpihelper_mod.o constants.o beam_mod.o : beam_mod.f90 distrib_mod.o basic_mod.o mpihelper_mod.o constants.o -chkrst.o : chkrst.f90 fields_mod.o beam_mod.o basic_mod.o +celldiag_mod.o : celldiag_mod.f90 beam_mod.o basic_mod.o mpihelper_mod.o constants.o +chkrst.o : chkrst.f90 constants.o fields_mod.o beam_mod.o basic_mod.o chkrst__genmod.o : chkrst__genmod.f90 constants.o : constants.f90 cp2bk__genmod.o : cp2bk__genmod.f90 -diagnose.o : diagnose.f90 fields_mod.o beam_mod.o basic_mod.o +diagnose.o : diagnose.f90 celldiag_mod.o fields_mod.o xg_mod.o beam_mod.o basic_mod.o diagnose__genmod.o : diagnose__genmod.f90 distrib_mod.o : distrib_mod.f90 constants.o endrun.o : endrun.f90 fields_mod.o beam_mod.o basic_mod.o endrun__genmod.o : endrun__genmod.f90 energies.o : energies.f90 basic_mod.o fields_mod.o : fields_mod.f90 mpihelper_mod.o beam_mod.o basic_mod.o constants.o inital.o : inital.f90 maxwellsource_mod.o mpihelper_mod.o fields_mod.o beam_mod.o basic_mod.o inital__genmod.o : inital__genmod.f90 main.o : main.f90 basic_mod.o maxwellsource_mod.o : maxwellsource_mod.f90 distrib_mod.o beam_mod.o basic_mod.o mpihelper_mod.o constants.o mpihelper_mod.o : mpihelper_mod.f90 constants.o mv2bk.o : mv2bk.f90 mv2bk__genmod.o : mv2bk__genmod.f90 newrun.o : newrun.f90 newrun__genmod.o : newrun__genmod.f90 restart.o : restart.f90 restart__genmod.o : restart__genmod.f90 resume.o : resume.f90 maxwellsource_mod.o sort_mod.o fields_mod.o basic_mod.o beam_mod.o resume__genmod.o : resume__genmod.f90 sort_mod.o : sort_mod.f90 constants.o beam_mod.o start.o : start.f90 basic_mod.o start__genmod.o : start__genmod.f90 stepon.o : stepon.f90 maxwellsource_mod.o beam_mod.o fields_mod.o constants.o basic_mod.o stepon__genmod.o : stepon__genmod.f90 tesend.o : tesend.f90 basic_mod.o tesend__genmod.o : tesend__genmod.f90 xg_mod.o : xg_mod.f90 fields_mod.o beam_mod.o basic_mod.o constants.o