diff --git a/matlab/FBspline.m b/matlab/FBspline.m index 36308f5..5b91d5d 100644 --- a/matlab/FBspline.m +++ b/matlab/FBspline.m @@ -1,262 +1,350 @@ -function [Forces, Continuity]=FBspline(M,it,Forces, Continuity, Density,norm,log,contourscale, zlims) +function [Forces, Continuity]=FBspline(M,it,Forces, Continuity, Density,norm,log,parper,contourscale, zlims) %ForceBalance Show the radial force balance % Plot the three radial forces for the given time-step it or averaged % over the range of time steps defined in it if it==':' it=floor(0.95*size(M.t2d)):size(M.t2d); end switch nargin case 5 + parper=false; contourscale=0.6; norm=false; log=false; zlims=[-inf inf]; case 6 + parper=false; log=false; contourscale=0.6; zlims=[-inf inf]; case 7 + parper=false; contourscale=0.6; zlims=[-inf inf]; case 8 + contourscale=0.6; zlims=[-inf inf]; case 9 + zlims=[-inf inf]; + case 10 otherwise error("Invalid number of arguments") end contourcolor=[0 0 0];%[255 20 147]/255; N=Density.N; maxdens=max(N(:)); densitycontour=contour(M.zgrid,M.rgrid,mean(N,3),[contourscale contourscale]*maxdens); -Fmaxr=max([max(Forces.Eforcer(:)),max(Forces.Bforcer(:)),max(Forces.inertforce(:)),max(Forces.Pressforcer(:))]); -Fmaxz=max([max(Forces.Eforcez(:)),max(Forces.Bforcez(:)),max(Forces.Pressforcez(:))]); +Fmaxr=max([max(Forces.Eforcer(:)),max(Forces.Bforcer(:)),max(Forces.inertforcer(:)),max(Forces.Pressforcer(:))]); +Fmaxz=max([max(Forces.Eforcez(:)),max(Forces.Bforcez(:)),max(Forces.Pressforcez(:)),max(Forces.inertforcez(:))]); +FEr=mean(Forces.Eforcer,3); +FBr=mean(Forces.Bforcer,3); +FPr=mean(Forces.Pressforcer,3); +FIr=mean(Forces.inertforcer,3); +FEz=mean(Forces.Eforcez,3); +FBz=mean(Forces.Bforcez,3); +FPz=mean(Forces.Pressforcez,3); +FIz=mean(Forces.inertforcez,3); + if log Er=Forces.Eforcer(Forces.Eforcer~=0); Br=Forces.Bforcer(Forces.Bforcer~=0); - inert=Forces.inertforce(Forces.inertforce~=0); + inertr=Forces.inertforcer(Forces.inertforcer~=0); + inertz=Forces.inertforcez(Forces.inertforcez~=0); pressr=Forces.Pressforcer(Forces.Pressforcer~=0); - Fminr=min([min(abs(Er(:))),min(abs(Br(:))),min(abs(inert(:))),min(abs(pressr(:)))]); + Fminr=min([min(abs(Er(:))),min(abs(Br(:))),min(abs(inertr(:))),min(abs(pressr(:)))]); else - Fminr=min([min(Forces.Eforcer(:)),min(Forces.Bforcer(:)),min(Forces.inertforce(:)),min(Forces.Pressforcer(:))]); + Fminr=min([min(Forces.Eforcer(:)),min(Forces.Bforcer(:)),min(Forces.inertforcer(:)),min(Forces.Pressforcer(:))]); end rgridmax=sum(M.nnr(1:2)); +if parper + costhet=(M.Br./M.B)'; + sinthet=(M.Bz./M.B)'; + FErd=FEr.*sinthet-FEz.*costhet; + FBrd=FBr.*sinthet-FBz.*costhet; + FPrd=FPr.*sinthet-FPz.*costhet; + FIrd=FIr.*sinthet-FIz.*costhet; + FEzd=FEr.*costhet+FEz.*sinthet; + FBzd=FBr.*costhet+FBz.*sinthet; + FPzd=FPr.*costhet+FPz.*sinthet; + FIzd=FIr.*costhet+FIz.*sinthet; + +else + FErd=FEr; + FBrd=FBr; + FPrd=FPr; + FIrd=FIr; + FEzd=FEz; + FBzd=FBz; + FPzd=FPz; + FIzd=FIz; +end + %% Radial forces f=figure(); ax1=subplot(4,1,1); -plotvalue(ax1,mean(Forces.Eforcer,3),'Electric force density', 'F_{Er} [Nm^{-3}]',true); +plotvalue(ax1,FErd,'Electric force density', 'F_{E,r} [Nm^{-3}]',true); %zlim=[Fmin Fmax]; zlim=[-inf inf]; ax2=subplot(4,1,2); -plotvalue(ax2,mean(Forces.Bforcer,3),'Magnetic force density','F_{Br} [Nm^{-3}]',true,zlim) +plotvalue(ax2,FBrd,'Magnetic force density','F_{B,r} [Nm^{-3}]',true,zlim) ax3=subplot(4,1,3); -plotvalue(ax3,mean(Forces.inertforce,3),'Centrifugal force density','F_{i} [Nm^{-3}]',true,zlim) +plotvalue(ax3,FIrd,'Inertial force density','F_{i,r} [Nm^{-3}]',true,zlim) ax4=subplot(4,1,4); -plotvalue(ax4,mean(Forces.Pressforcer,3),'Pressure force density','F_{pr} [Nm^{-3}]',true,zlim) +plotvalue(ax4,FPrd,'Pressure force density','F_{p,r} [Nm^{-3}]',true,zlim) linkaxes([ax1 ax2 ax3 ax4],'xy') +if parper + sgtitle(sprintf('Perp forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) +else sgtitle(sprintf('Radial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) - +end f.PaperOrientation='portrait'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[20 29]; f.PaperSize=papsize; print(f,sprintf('%sforce_balancer',name),'-dpdf','-fillpage') savefig(f,sprintf('%sforce_balancer',name)) %% Radial relative forces f=figure(); ax1=subplot(4,1,1); color_limits=[-1 1]; -reference=Forces.Bforcer; -plotvalue(ax1,mean(Forces.Eforcer./reference,3),'Relative Electric force density', 'F_{Er}/F_{Br}',true,color_limits); +referencer=FBrd; +plotvalue(ax1,FErd./referencer,'Relative Electric force density', 'F_{Er}/F_{Br}',true,color_limits); %zlim=[Fmin Fmax]; zlim=[-inf inf]; ax2=subplot(4,1,2); -plotvalue(ax2,mean(Forces.Bforcer./reference,3),'Relative Magnetic force density','F_{Br}/F_{Br}',true,color_limits) +plotvalue(ax2,FBrd./referencer,'Relative Magnetic force density','F_{Br}/F_{Br}',true,color_limits) ax3=subplot(4,1,3); -plotvalue(ax3,mean(Forces.inertforce./reference,3),'Relative Centrifugal force density','F_i/F_{Br}',true,color_limits) +plotvalue(ax3,FIrd./referencer,'Relative Inertial force density','F_i/F_{Br}',true,color_limits) ax4=subplot(4,1,4); -plotvalue(ax4,mean(Forces.Pressforcer./reference,3),'Relative Pressure force density','F_{pr}/F_{Br}',true,color_limits) +plotvalue(ax4,FPrd./referencer,'Relative Pressure force density','F_{pr}/F_{Br}',true,color_limits) linkaxes([ax1 ax2 ax3 ax4],'xy') +if parper + sgtitle(sprintf('Relative perp forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) +else sgtitle(sprintf('Relative radial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) - +end f.PaperOrientation='portrait'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[20 29]; f.PaperSize=papsize; -print(f,sprintf('%s_relatforce_balance',name),'-dpdf','-fillpage') -savefig(f,sprintf('%s_relatforce_balance',name)) +print(f,sprintf('%s_relatforce_balancer',name),'-dpdf','-fillpage') +savefig(f,sprintf('%s_relatforce_balancer',name)) %% Axial forces f=figure(); -ax1=subplot(3,1,1); -plotvalue(ax1,mean(Forces.Eforcez,3),'Electric force density', 'F_{Ez} [Nm^{-3}]',true); +ax1=subplot(4,1,1); +plotvalue(ax1,FEzd,'Electric force density', 'F_{E,z} [Nm^{-3}]',true); +%zlim=[Fmin Fmax]; zlim=[-inf inf]; -ax2=subplot(3,1,2); -plotvalue(ax2,mean(Forces.Bforcez,3),'Magnetic force density','F_{Bz} [Nm^{-3}]',true,zlim) +ax2=subplot(4,1,2); +plotvalue(ax2,FBzd,'Magnetic force density','F_{B,z} [Nm^{-3}]',true,zlim) + +ax3=subplot(4,1,3); +plotvalue(ax3,FIzd,'Inertial force density','F_{i,z} [Nm^{-3}]',true,zlim) -ax3=subplot(3,1,3); -plotvalue(ax3,mean(Forces.Pressforcez,3),'Pressure force density','F_{pz} [Nm^{-3}]',true,zlim) +ax4=subplot(4,1,4); +plotvalue(ax4,FPzd,'Pressure force density','F_{p,z} [Nm^{-3}]',true,zlim) -linkaxes([ax1 ax2 ax3],'xy') +linkaxes([ax1 ax2 ax3 ax4],'xy') sgtitle(sprintf('Axial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) f.PaperOrientation='portrait'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[20 29]; f.PaperSize=papsize; print(f,sprintf('%sforce_balancez',name),'-dpdf','-fillpage') savefig(f,sprintf('%sforce_balancez',name)) +%% Axial relative forces +f=figure(); +ax1=subplot(4,1,1); +color_limits=[-1 1]; +referencez=FEzd; +plotvalue(ax1,FEzd./referencez,'Relative Electric force density', 'F_{Ez}/F_{Ez}',true,color_limits); + +%zlim=[Fmin Fmax]; +zlim=[-inf inf]; + +ax2=subplot(4,1,2); +plotvalue(ax2,FBzd./referencez,'Relative Magnetic force density','F_{Bz}/F_{Ez}',true,color_limits) + +ax3=subplot(4,1,3); +plotvalue(ax3,FIzd./referencez,'Relative Inertial force density','F_{Iz}/F_{Ez}',true,color_limits) + +ax4=subplot(4,1,4); +plotvalue(ax4,FPzd./referencez,'Relative Pressure force density','F_{Pz}/F_{Ez}',true,color_limits) + +linkaxes([ax1 ax2 ax3 ax4],'xy') +if parper + sgtitle(sprintf('Relative Parallel forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) +else +sgtitle(sprintf('Relative axial forces t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) +end +f.PaperOrientation='portrait'; +[~, name, ~] = fileparts(M.file); +f.PaperUnits='centimeters'; +papsize=[20 29]; +f.PaperSize=papsize; +print(f,sprintf('%s_relatforce_balancez',name),'-dpdf','-fillpage') +savefig(f,sprintf('%s_relatforce_balancez',name)) + %% Forces sums f2=figure(); ax4=subplot(3,1,1); -totforcer=mean(Forces.Eforcer+Forces.Bforcer+Forces.inertforce+Forces.Pressforcer,3); -Forcemaxr=cat(3,mean(Forces.Eforcer,3),mean(Forces.Bforcer,3),mean(Forces.inertforce,3),mean(Forces.Pressforcer,3)); +totforcer=FErd+FBrd+FIrd+FPrd; +Forcemaxr=cat(3,FErd,FBrd,FIrd,FPrd); Forcemaxr=max(Forcemaxr,[],3); -plotvalue(ax4,totforcer,'Total radial force density','F [Nm^{-3}]',true,zlims) - +if parper +plotvalue(ax4,totforcer,'Total Perp force density','F [Nm^{-3}]',true,zlims) +else + plotvalue(ax4,totforcer,'Total radial force density','F [Nm^{-3}]',true,zlims) +end ax5=subplot(3,1,2); -totforcez=mean(Forces.Eforcez+Forces.Bforcez+Forces.Pressforcez,3); -Forcemaxz=cat(3,mean(Forces.Eforcez,3),mean(Forces.Bforcez,3),mean(Forces.Pressforcez,3)); +totforcez=FEzd+FBzd+FIzd+FPzd; +Forcemaxz=cat(3,FEzd,FBzd,FIzd,FPzd); Forcemaxz=max(Forcemaxz,[],3); -plotvalue(ax5,totforcez,'Total axial force density','F [Nm^{-3}]',true,zlims) - +if parper + plotvalue(ax5,totforcez,'Total par force density','F [Nm^{-3}]',true,zlims) +else + plotvalue(ax5,totforcez,'Total axial force density','F [Nm^{-3}]',true,zlims) +end ax6=subplot(3,1,3); surface(ax6,M.zgrid,M.rgrid,mean(N,3),'edgecolor','none'); hold on; zlevel=interp2(M.zgrid, M.rgrid, mean(N,3), densitycontour(1,2:end), densitycontour(2,2:end)); plot3(ax6,densitycontour(1,2:end),densitycontour(2,2:end),zlevel,'-.','linewidth',2,'color',contourcolor) xlim(ax6,[M.zgrid(1) M.zgrid(end)]) if M.conformgeom - ylim(axis6,[M.rgrid(1) M.rgrid(rgridmax)]) + ylim(ax6,[M.rgrid(1) M.rgrid(rgridmax)]) else - ylim(axis6,[M.rgrid(1) M.rgrid(end)]) + ylim(ax6,[M.rgrid(1) M.rgrid(end)]) end colormap(ax6,'parula') xlabel(ax6,'z [m]') ylabel(ax6,'r [m]') title(ax6,'Density') c = colorbar(ax6); c.Label.String= 'n[m^{-3}]'; view(ax6,2) sgtitle(sprintf('Forces sum t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) f2.PaperOrientation='portrait'; [~, name, ~] = fileparts(M.file); f2.PaperUnits='centimeters'; papsize=[20 29]; f2.PaperSize=papsize; print(f2,sprintf('%sforce_dens',name),'-dpdf','-fillpage') savefig(f2,sprintf('%sforce_dens',name)) %% Continuity if( ~isempty(Continuity)) f3=figure(); t=M.t2d(it); ax7=subplot(4,1,1); plotvalue(ax7,mean(Continuity.dndt,3),'\partialn/\partialt','\partialn/\partialt [m^{-3}s^{-1}]',true) zlim=[max([min(min(mean(Continuity.ndivu(2:end,:,:),3))),min(min(mean(Continuity.ugradn(2:end,:,:),3)))]) min([max(max(mean(Continuity.ndivu(2:end,:,:),3))),max(max(mean(Continuity.ugradn(2:end,:,:),3)))])]; ax8=subplot(4,1,2); plotvalue(ax8,mean(Continuity.ndivu,3),'n\nabla(u)','n\nabla(u)[m^{-3}s^{-1}]',true,zlim) ax9=subplot(4,1,3); plotvalue(ax9,mean(Continuity.ugradn,3),'u\nabla(n)','u\nabla(n)[m^{-3}s^{-1}]',true,zlim) ax10=subplot(4,1,4); plotvalue(ax10,mean(Continuity.ndivu+Continuity.ugradn+Continuity.dndt,3),'Total','total[m^{-3}s^{-1}]',true,zlim) sgtitle(sprintf('Continuity t=[%1.2g-%1.2g]s n_e=%1.2g m^{-3}',M.t2d(min(it)),M.t2d(max(it)),double(maxdens))) f3.PaperOrientation='portrait'; [~, name, ~] = fileparts(M.file); f3.PaperUnits='centimeters'; papsize=[20 19]; f3.PaperSize=papsize; print(f3,sprintf('%scontinuity',name),'-dpdf','-fillpage') savefig(f,sprintf('%scontinuity',name)) end function plotvalue(axis,matrix,titl, clab, rm1strow, zlims) if nargin < 5 rm1strow=false; end if nargin>5 zmin=zlims(1); zmax=zlims(2); else zmin=-inf; zmax=inf; end if rm1strow if log h=surface(axis,M.zgrid,M.rgrid(12:end),abs(matrix(12:end,:))); else h=surface(axis,M.zgrid,M.rgrid(12:end),matrix(12:end,:)); end else if log h=surface(axis,M.zgrid,M.rgrid,abs(matrix)); else h=surface(axis,M.zgrid,M.rgrid,matrix); end end set(h,'edgecolor','none'); hold on; if log zpos=interp2(M.zgrid, M.rgrid, abs(matrix), densitycontour(1,2:end), densitycontour(2,2:end)); else zpos=interp2(M.zgrid, M.rgrid, matrix, densitycontour(1,2:end), densitycontour(2,2:end)); end plot3(axis,densitycontour(1,2:end),densitycontour(2,2:end),zpos,'-.','linewidth',2,'color',contourcolor) xlabel(axis,'z [m]') ylabel(axis,'r [m]') xlim(axis,[M.zgrid(1) M.zgrid(end)]) if M.conformgeom ylim(axis,[M.rgrid(1) M.rgrid(rgridmax)]) else ylim(axis,[M.rgrid(1) M.rgrid(end)]) end colormap(axis,'jet') c = colorbar(axis); c.Label.String=clab; caxis(axis,[zmin zmax]) if norm caxis(axis,[Fminr Fmaxr]); end if log set(axis,'colorscale','log') caxis(axis,'auto') end title(axis,titl) grid(axis, 'on') view(axis,2) end end diff --git a/matlab/ForceBalancespline.m b/matlab/ForceBalancespline.m index 1b528ef..e506f56 100644 --- a/matlab/ForceBalancespline.m +++ b/matlab/ForceBalancespline.m @@ -1,82 +1,72 @@ function [Forces, Density]=ForceBalancespline(M,it,norm,log,contourscale, zlims) %ForceBalance Show the radial force balance % Plot the three radial forces for the given time-step it or averaged % over the range of time steps defined in it if it==':' it=floor(0.95*size(M.t2d)):size(M.t2d); end switch nargin case 2 contourscale=0.6; norm=false; log=false; zlims=[-inf inf]; case 3 log=false; contourscale=0.6; zlims=[-inf inf]; case 4 contourscale=0.6; zlims=[-inf inf]; case 5 zlims=[-inf inf]; case 6 otherwise error("Invalid number of arguments") end contourcolor=[0 0 0];%[255 20 147]/255; fluiduThet=M.fluidUTHET(:,:,it); Density.fluiduThet=fluiduThet; N=M.N(:,:,it); [R,~]=meshgrid(M.rgrid,M.zgrid); Rinv=1./R'; Rinv(isinf(Rinv))=0; Density.N=N; maxdens=max(N(:)); densitycontour=contour(M.zgrid,M.rgrid,mean(N,3),[contourscale contourscale]*maxdens); Forces.Eforcer=-M.qe*N.*M.Er(:,:,it); Forces.Eforcez=-M.qe*N.*M.Ez(:,:,it); Bforcer=zeros(size(N,1),size(N,2),size(N,3)); Bforcez=zeros(size(N,1),size(N,2),size(N,3)); -inertforce=zeros(size(N,1),size(N,2),size(N,3)); +inertforcer=zeros(size(N,1),size(N,2),size(N,3)); +inertforcez=zeros(size(N,1),size(N,2),size(N,3)); Pressforcer=zeros(size(N,1),size(N,2),size(N,3)); Pressforcez=zeros(size(N,1),size(N,2),size(N,3)); parfor j=1:size(N,3) Bforcer(:,:,j)=-M.qe*N(:,:,j).*fluiduThet(:,:,j).*M.Bz'; Bforcez(:,:,j)=M.qe*N(:,:,j).*fluiduThet(:,:,j).*M.Br'; - inertforce(:,:,j)=M.me*N(:,:,j).*fluiduThet(:,:,j).^2.*Rinv; + inertforcer(:,:,j)=M.me*N(:,:,j).*(fluiduThet(:,:,j).^2.*Rinv+M.fluidUR(:,:,it(j)).*M.fluidUR.der(:,:,it(j),[1 0])+M.fluidUZ(:,:,it(j)).*M.fluidUR.der(:,:,it(j),[0 1])); + inertforcez(:,:,j)=M.me*N(:,:,j).*(M.fluidUR(:,:,it(j)).*M.fluidUZ.der(:,:,it(j),[1 0])+M.fluidUZ(:,:,it(j)).*M.fluidUZ.der(:,:,it(j),[0 1])); Pressforcer(:,:,j)=-( squeeze(M.Presstens.der(1,:,:,it(j),[1 0]))... + squeeze(M.Presstens(1,:,:,it(j)) - M.Presstens(4,:,:,it(j))).*Rinv... + squeeze(M.Presstens.der(3,:,:,it(j),[0 1])) ); Pressforcez(:,:,j)=-( squeeze(M.Presstens.der(3,:,:,it(j),[1 0]))... + squeeze(M.Presstens(3,:,:,it(j))).*Rinv... + squeeze(M.Presstens.der(6,:,:,it(j),[0 1])) ); end Forces.Bforcer=Bforcer; Forces.Bforcez=Bforcez; -Forces.inertforce=inertforce; +Forces.inertforcer=inertforcer; +Forces.inertforcez=inertforcez; Forces.Pressforcer=Pressforcer; Forces.Pressforcez=Pressforcez; - -Fmaxr=max([max(Forces.Eforcer(:)),max(Forces.Bforcer(:)),max(Forces.inertforce(:)),max(Forces.Pressforcer(:))]); -Fmaxz=max([max(Forces.Eforcez(:)),max(Forces.Bforcez(:)),max(Forces.Pressforcez(:))]); -if log - Er=Forces.Eforcer(Eforcer~=0); - Br=Forces.Bforcer(Bforcer~=0); - inert=Forces.inertforce(inertforce~=0); - pressr=Forces.Pressforcer(Pressforcer~=0); - Fminr=min([min(abs(Er(:))),min(abs(Br(:))),min(abs(inert(:))),min(abs(pressr(:)))]); -else - Fminr=min([min(Forces.Eforcer(:)),min(Forces.Bforcer(:)),min(Forces.inertforce(:)),min(Forces.Pressforcer(:))]); -end -rgridmax=sum(M.nnr(1:2)); end diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m index 9e6046d..cd39904 100644 --- a/matlab/analyse_espicfields.m +++ b/matlab/analyse_espicfields.m @@ -1,447 +1,446 @@ % 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-800);%floor(0.95*t2dlength); +fieldstart=max(1,t2dlength-400);%floor(0.95*t2dlength); fieldend=t2dlength;%23; -deltat=t2dlength-max(1,t2dlength-800); - +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:fieldend),3); pot=mean(M.pot(:,:,fieldstart:fieldend),3); brillratio=2*dens*M.me./(M.eps_0*M.Bz'.^2); [R,Z]=meshgrid(M.rgrid,M.zgrid); Rinv=1./R; Rinv(isnan(Rinv))=0; VTHET=mean(M.fluidUTHET(:,:,fieldstart:fieldend),3); omegare=(VTHET.*Rinv'); maxdens=max(max(dens)); geomw=M.geomweight(:,:,1); geomw(geomw<0)=0; geomw(geomw>0)=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 M.displaydensity(fieldstart,fieldend); %% 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'); [~, hContour]=contourf(ax1,M.zgrid,M.rgrid,geomw,2); drawnow; xlim(ax1,[M.zgrid(1) M.zgrid(end)]) ylim(ax1,[M.rgrid(1) M.rgrid(end)]) legend({'n_e [m^{-3}]','Magnetic field lines'},'location','southwest') 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(fieldend),double(maxdens))) c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; view(ax1,2) %set(h,'edgecolor','none'); grid on; hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' hFills(1).ColorData = uint8([255;255;255;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end %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*1000,M.rgrid*1000,dens,'Displayname','n_e [m^{-3}]'); hold on; pot(M.geomweight(:,:,1)<0)=0; 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); [~,h1]=contour(ax1,M.zgrid*1000,M.rgrid*1000,Blines',real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines'); [c1,h2]=contour(ax1,M.zgrid*1000,M.rgrid*1000,pot./1e3,'m--','ShowText','on','linewidth',1.2,'Displayname','Equipotentials [kV]'); clabel(c1,h2,'Color','white') %contour(ax1,M.zgrid*1000,M.rgrid*1000,M.geomweight(:,:,1),[0 0],'w-','linewidth',1.5); [c1,hContour]=contourf(ax1,M.zgrid*1000,M.rgrid*1000,geomw); drawnow; xlim(ax1,[M.zgrid(1)*1000 M.zgrid(end)*1000]) if(M.conformgeom) ylim(ax1,[M.rgrid(1)*1000 M.rgrid(rgridend)*1000]) else ylim(ax1,[M.rgrid(1)*1000 M.rgrid(end)*1000]) end legend([h1,h2],{'Magnetic field lines','Equipotentials [kV]'},'location','southwest') xlabel(ax1,'z [mm]') ylabel(ax1,'r [mm]') %title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; view(ax1,2) %set(h,'edgecolor','none'); grid on; hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' hFills(1).ColorData = uint8([150;150;150;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); f.PaperUnits='centimeters'; papsize=[20 16]; f.PaperSize=papsize; rectangle('Position',[M.zgrid(1) M.rgrid(1) M.zgrid(end)-M.zgrid(1) 0.064-M.rgrid(1)]*1000,'FaceColor',[150 150 150]/255) 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)]) if(M.conformgeom) ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)]) else ylim(ax1,[M.rgrid(1) M.rgrid(end)]) end 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]) +caxis([0 max(1.2,max(brillratio(:)))]) title(ax1,sprintf('Brillouin Ratio t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),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)]) if(M.conformgeom) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) else ylim(ax3,[M.rgrid(1) M.rgrid(end)]) end 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(fieldend),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)]) if(M.conformgeom) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) else ylim(ax3,[M.rgrid(1) M.rgrid(end)]) end 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(fieldend),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 f=figure('Name', sprintf('%s V drift',name)); ax3=gca; vdrift=(mean(M.Ez(:,:,fieldstart:end),3).*M.Br'-mean(M.Er(:,:,fieldstart:end),3).*M.Bz')./(M.B.^2)'; v=vdrift; 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)]) if(M.conformgeom) ylim(ax3,[M.rgrid(1) M.rgrid(rgridend)]) else ylim(ax3,[M.rgrid(1) M.rgrid(end)]) end colormap(ax3,'jet') c = colorbar(ax3); set(ax3,'zscale','log') set(ax3,'colorscale','log') c.Label.String= 'V_{ExB}'; %c.Limits=[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]; %caxis(ax3,[min(M.fluidUTHET(:)) max(M.fluidUTHET(:))]) title(ax3,sprintf('V_{ExB} velocity t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),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 %% VR, VTHET, VZ if(M.R.nt>=2) timesteppart=length(M.tpart); fig=figure; ax1=gca; nbp=min(M.nbparts(1),M.R.nparts); Vr=M.VR(1:nbp,1,false)/M.vlight; Vz=M.VZ(1:nbp,1,false)/M.vlight; Vthet=M.VTHET(1:nbp,1,false)/M.vlight; Vperp=sqrt(Vthet.^2+Vr.^2); nbp=min(M.nbparts(end),M.R.nparts); Vrend=M.VR(1:nbp,timesteppart,false)/M.vlight; Vzend=M.VZ(1:nbp,timesteppart,false)/M.vlight; Vthetend=M.VTHET(1:nbp,timesteppart,false)/M.vlight; Vperpend=sqrt(Vthetend.^2+Vrend.^2); binwidth=abs(max(Vrend)-min(Vrend))/sqrt(length(Vrend)); f=figure('Name',sprintf("%s v distrib",M.file)); tstudied=0; legtext=sprintf("t=%2.1f [ns]",M.tpart(end)*1e9); ax1=subplot(1,3,1); if(length(Vr)>1) h1=histogram(Vr,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); %h1=histfit(ax1,Vr); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end hold on h1=histogram(Vrend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); %h1=histfit(ax1,Vrend); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_r') grid on ax2=subplot(1,3,2); if(length(Vthet)>1) h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end %h1=histfit(ax2,Vthet(:,1)); hold on h1=histogram(Vthetend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); %h1=histfit(ax2,Vthetend); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_\theta') legend(ax2,'location','northoutside','orientation','vertical') grid on ax3=subplot(1,3,3); %h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); if(length(Vz)>1) h1=histogram(ax3,Vz(:,1)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end hold on %h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); h1=histogram(ax3,Vzend); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_z') grid on f=gcf; sgtitle(sprintf('dt=%1.2e[ns]',M.dt*1e9)) f.PaperOrientation='landscape'; f.PaperUnits='centimeters'; papsize=[16 10]; f.PaperSize=papsize; [~, name, ~] = fileparts(M.file); print(f,sprintf('%sParts_V_RZ',name),'-dpdf','-fillpage') savefig(f,sprintf('%sParts_V_RZ',name)) 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 +NpartsTrapped=mean(M.npart(end-10:end))*M.msim/M.me dblengthiter=fieldstart:fieldend; 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-1); M.displaypotentialwell(t2dlength-1,false); diff --git a/matlab/dispespicParts.m b/matlab/dispespicParts.m index c02edf7..1497112 100644 --- a/matlab/dispespicParts.m +++ b/matlab/dispespicParts.m @@ -1,171 +1,199 @@ -function dispespicParts(M) +function dispespicParts(M,fieldstep,parper) %dispespicParts Show on an interactive plot the phase space of the given specie in M % Detailed explanation goes here -fieldstep=1; - +if nargin<2 + fieldstep=1; +end +if nargin<3 + parper=false; +end f2=uifigure('Name','Particles data'); mf2=uipanel(f2,'Position',[5 50 f2.Position(3)-10 f2.Position(4)-55]); mf2.AutoResizeChildren='off'; m2=uipanel(f2,'Position',[5 5 f2.Position(3)-10 40]); sgtitle(mf2,sprintf('step=%d t=%0.5e s',fieldstep*M.it2,M.tpart(fieldstep))) sld2 = uislider(m2,'Position',[10 30 0.6*m2.Position(3) 3]); sld2.Value=fieldstep; sld2.Limits=[1 length(M.tpart)]; edt2 = uieditfield(m2,'numeric','Limits',[1 length(M.tpart)],'Value',1); edt2.Position=[sld2.Position(1)+sld2.Position(3)+25 5 40 20]; edt2.RoundFractionalValues='on'; Printbt2=uibutton(m2,'Position',[edt2.Position(1)+edt2.Position(3)+10 5 40 20],'Text', 'Save'); sld2.ValueChangingFcn={@updatefigpartdata,edt2,mf2}; edt2.ValueChangedFcn={@updatefigpartdata,sld2,mf2}; Printbt2.ButtonPushedFcn={@plotPartButtonPushed}; set(f2,'KeyPressFcn',{ @onKeyDown,sld2,edt2,mf2}) PlotEspic2dpartdata(mf2,M,fieldstep); function PlotEspic2dpartdata(fig,M,fieldstep) %PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it2,M.tpart(fieldstep))) + Z=M.Z(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); + R=M.R(1:min(M.nbparts(fieldstep),M.R.nparts),fieldstep); + + if parper + vpar=M.Vpar(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); + vper=M.Vperp(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); + else + vpar=M.VZ(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); + vper=M.VR(1:min(M.nbparts(fieldstep),M.Z.nparts),fieldstep); + end + ax1=subplot(3,2,1,'Parent',fig); if isa(M,'h5parts') contour(ax1,M.zgrid,M.rgrid,M.parent.geomweight(:,:,1),[0 0],'r--') else contour(ax1,M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'r--') end hold(ax1, 'on') - plot(ax1,M.Z(1:M.nbparts(fieldstep),fieldstep),M.R(1:M.nbparts(fieldstep),fieldstep),'.'); + plot(ax1,Z,R,'.'); xlim(ax1,[M.zgrid(1) M.zgrid(end)]) ylim(ax1,[M.rgrid(1) M.rgrid(end)]) xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') title(ax1,'Position') grid(ax1, 'on') - ax2=subplot(3,2,2,'Parent',fig); - plot(ax2,M.Vpar(1:M.nbparts(fieldstep),fieldstep),M.Vperp(1:M.nbparts(fieldstep),fieldstep),'.'); - xlabel(ax2,'V_{par} [m/s]') - ylabel(ax2,'V_{perp} [m/s]') + plot(ax2,vpar,vper,'.'); + xlabel(ax2,'V_Z,{par} [m/s]') + ylabel(ax2,'V_R,{perp} [m/s]') axis(ax2, 'equal') grid(ax2, 'on') ax3=subplot(3,2,3,'Parent',fig); - plot(ax3,M.Z(1:M.nbparts(fieldstep),fieldstep),M.VZ(1:M.nbparts(fieldstep),fieldstep),'.'); + plot(ax3,Z,vpar,'.'); xlim(ax3,[M.zgrid(1) M.zgrid(end)]) xlabel(ax3,'Z [m]') ylabel(ax3,'VZ [m/s]') grid(ax3, 'on') ax4=subplot(3,2,4,'Parent',fig); - plot(ax4,M.R(1:M.nbparts(fieldstep),fieldstep),M.VZ(1:M.nbparts(fieldstep),fieldstep),'.') + plot(ax4,R,vpar,'.') ylabel(ax4,'VZ [m/s]') xlabel(ax4,'R [m]') xlim(ax4,[M.rgrid(1) M.rgrid(end)]) grid(ax4, 'on') ax5=subplot(3,2,5,'Parent',fig); - plot(ax5,M.Z(1:M.nbparts(fieldstep),fieldstep),M.VR(1:M.nbparts(fieldstep),fieldstep),'.'); + plot(ax5,Z,vper,'.'); xlim(ax5,[M.zgrid(1) M.zgrid(end)]) xlabel(ax5,'Z [m]') ylabel(ax5,'VR [m/s]') grid(ax5, 'on') ax6=subplot(3,2,6,'Parent',fig); - plot(ax6,M.R(1:M.nbparts(fieldstep),fieldstep),M.VR(1:M.nbparts(fieldstep),fieldstep),'.'); + plot(ax6,R,vper,'.'); xlim(ax6,[M.rgrid(1) M.rgrid(end)]) ylabel(ax6,'VR [m/s]') xlabel(ax6,'R [m]') grid(ax6, 'on') end function onKeyDown(src,event,slider,editfield, fig) direction=0; if strcmp(event.Key,'leftarrow') direction=-1; elseif strcmp(event.Key,'rightarrow') direction=+1; elseif strcmp(event.Key,'uparrow') direction=+10; elseif strcmp(event.Key,'downarrow') direction=-10; end - + if(direction~=0) currval=slider.Value; slider.Value=max(slider.Limits(1),min(currval+direction,slider.Limits(2))); updatefigpartdata(slider, event, editfield ,fig) end end function updatefigpartdata(control, event, Othercontrol, fig) if strcmp(event.EventName,'ValueChanged') fieldstep=floor(control.Value); control.Value=fieldstep; elseif strcmp(event.EventName,'KeyPress') fieldstep=floor(control.Value); control.Value=fieldstep; else fieldstep=floor(event.Value); end Othercontrol.Value=fieldstep; sgtitle(fig,sprintf('t=%0.5e s',M.tpart(fieldstep))) - [~,t0dstep]=min(abs(M.tpart(fieldstep)-M.t0d)); + % if isa(M,'espic2dhdf5') + % [~,t0dstep]=min(abs(M.tpart(fieldstep)-M.t0d)); + % else + t0dstep=fieldstep; + % end + Z=M.Z(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); + R=M.R(1:min(M.nbparts(t0dstep),M.R.nparts),fieldstep); + + if parper + vpar=M.Vpar(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); + vper=M.Vperp(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); + else + vpar=M.VZ(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); + vper=M.VR(1:min(M.nbparts(t0dstep),M.Z.nparts),fieldstep); + end %% update Position plot ax1=fig.Children(end); - ax1.Children(1).XData=M.Z(1:M.nbparts(t0dstep),fieldstep); - ax1.Children(1).YData=M.R(1:M.nbparts(t0dstep),fieldstep); + ax1.Children(1).XData=Z; + ax1.Children(1).YData=R; view(ax1,2) - %% update VZ VPERP - fig.Children(end-1).Children(1).XData=M.Vpar(1:M.nbparts(t0dstep),fieldstep); - fig.Children(end-1).Children(1).YData=M.Vperp(1:M.nbparts(t0dstep),fieldstep); + %% update VPAR VPERP + fig.Children(end-1).Children(1).XData=vpar; + fig.Children(end-1).Children(1).YData=vper; axis(fig.Children(end-1),'equal') %% update Z VZ - fig.Children(end-2).Children(1).XData=M.Z(1:M.nbparts(t0dstep),fieldstep); - fig.Children(end-2).Children(1).YData=M.VZ(1:M.nbparts(t0dstep),fieldstep); + fig.Children(end-2).Children(1).XData=Z; + fig.Children(end-2).Children(1).YData=vpar; %% update VZ VR - fig.Children(end-3).Children(1).XData=M.R(1:M.nbparts(t0dstep),fieldstep); - fig.Children(end-3).Children(1).YData=M.VZ(1:M.nbparts(t0dstep),fieldstep); + fig.Children(end-3).Children(1).XData=R; + fig.Children(end-3).Children(1).YData=vpar; %% update R VR - fig.Children(end-4).Children(1).XData=M.Z(1:M.nbparts(t0dstep),fieldstep); - fig.Children(end-4).Children(1).YData=M.VR(1:M.nbparts(t0dstep),fieldstep); + fig.Children(end-4).Children(1).XData=Z; + fig.Children(end-4).Children(1).YData=vper; %% update Z VR - fig.Children(end-5).Children(1).XData=M.R(1:M.nbparts(t0dstep),fieldstep); - fig.Children(end-5).Children(1).YData=M.VR(1:M.nbparts(t0dstep),fieldstep); + fig.Children(end-5).Children(1).XData=R; + fig.Children(end-5).Children(1).YData=vper; %drawnow limitrate end function plotPartButtonPushed(btn,ax) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here f=figure(); PlotEspic2dpartdata(f,M,sld.Value); f.PaperOrientation='portrait'; [~, name, ~] = fileparts(M.file); print(f,sprintf('%sParts%d',name,sld.Value),'-dpdf','-fillpage') end end diff --git a/matlab/distribution_function.m b/matlab/distribution_function.m index ff56a28..b786d71 100644 --- a/matlab/distribution_function.m +++ b/matlab/distribution_function.m @@ -1,368 +1,368 @@ %% Show the particles velocity histogram at the position set using ginput % The histogram is compared between timesteppart and end timesteppart=length(M.tpart); fig=figure; ax1=gca; timestepN=find(M.tpart(end)==M.t2d,1); Ndistrib=M.N(:,:,timestepN); h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib); %set(h, 'EdgeColor', 'none'); hold on [r,z]=find(Ndistrib~=0); xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))]) ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))]) xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') c = colorbar(ax1); c.Label.String= 'n [m^{-3}]'; view(ax1,2) %set(ax1,'colorscale','log') [x,y]=ginput(1); -%Zindex=find(x>M.zgrid,1,'last'); -%Rindex=find(y>M.rgrid,1,'last'); - Rindex=155; - Zindex=344; +Zindex=find(x>M.zgrid,1,'last'); +Rindex=find(y>M.rgrid,1,'last'); +% Rindex=155; +% Zindex=344; Z=(M.zgrid(Zindex)); R=(M.rgrid(Rindex)); plot(Z,R,'rx','Markersize',12); % Rindex=16; Zindex=64; % Rindex=28; Zindex=floor(size(M.zgrid,1)/2)+1; nbp=min(M.R.nparts,M.nbparts(1)); Rp=M.R(1:nbp,1,false); Zp=M.Z(1:nbp,1,false); Indices=Rp>=M.rgrid(Rindex) & Rp=M.zgrid(Zindex) & Zp=M.rgrid(Rindex) & Rend=M.zgrid(Zindex) & Zend1) h1=histogram(Vr,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); %h1=histfit(ax1,Vr); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end hold on h1=histogram(Vrend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); %h1=histfit(ax1,Vrend); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_r') grid on ax2=subplot(1,3,2); binwidth=abs(max(Vthetend)-min(Vthetend))/sqrt(length(Vthetend)); if(length(Vthet)>1) h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end %h1=histfit(ax2,Vthet(:,1)); hold on h1=histogram(Vthetend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); %h1=histfit(ax2,Vthetend); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_\theta') legend(ax2,'location','northoutside','orientation','vertical') grid on ax3=subplot(1,3,3); %h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); binwidth=abs(max(Vzend)-min(Vzend))/sqrt(length(Vzend)); if(length(Vz)>1) h1=histogram(ax3,Vz(:,1),'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); end hold on %h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); h1=histogram(ax3,Vzend,'Binwidth',binwidth,'Normalization','count','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); set(h1,'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); ylabel('counts') xlabel('\beta_z') grid on f=gcf; sgtitle(sprintf('R=%1.2e[m] Z=%1.2e[m] dt=%1.2e[ns]',M.rgrid(Rindex+1),M.zgrid(Zindex+1),M.dt*1e9)) f.PaperOrientation='landscape'; f.PaperUnits='centimeters'; papsize=[16 10]; f.PaperSize=papsize; [~, name, ~] = fileparts(M.file); print(f,sprintf('%sParts_V_RZ',name),'-dpdf','-fillpage') savefig(f,sprintf('%sParts_V_RZ',name)) %% Shows the particles phase space at position (rindex, zindex) % and times tinit and timesteppart timeinit=1; timestepNinit=find(M.tpart(timeinit)==M.t2d); timesteppart=length(M.tpart); timestepNend=find(M.tpart(timesteppart)==M.t2d); Rp=M.R(:,timeinit,false); Zp=M.Z(:,timeinit,false); Indices=Rp>=M.rgrid(Rindex) & Rp=M.zgrid(Zindex) & Zp=M.rgrid(Rindex) & Rend=M.zgrid(Zindex) & Zend=M.rAthet(:,zleftlim),1,'last'); rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2)); rposleft=find(M.rgrid>=rleft,1,'first'); n2=mean(mean(N(rposleft,[1 end],:),3)); vpar2=(dN/V/n2)^2; if n2==0 vpar2=0; end b=M.rgrid(end); a=M.rgrid(1); vd1=((M.potout-M.potinn)*M.phinorm/M.rgrid(Rindex))/M.Bz(Zindex,Rindex)/log(b/a); vd1=-M.Er(Rindex,Zindex,1)/M.Bz(Zindex,Rindex); %vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex); vinit=sqrt(M.kb*10000/M.me); %vinit=sqrt(M.qe*80/M.me) %vinit=sqrt(120*M.qe/M.me) %vd2=((M.potout-M.potinn)*M.phinorm/M.rgrid(25))/M.Bz(1,25)/log(b/a); [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1)); Rcurv=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); deltaphicomp=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv)) ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv)) Ndistrib=M.N(:,:,timestepNend); model=M.potentialwellmodel(timestepNend); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; Zeval=[M.zgrid(zleftlim) M.zgrid(Zindex)]; Psieval=M.rAthet(Rindex,Zindex)*[1 1]; phis=griddata(Zmesh,M.rAthet,M.pot(:,:,end),Zeval,Psieval,'natural'); deltaphi=-diff(phis) %deltaphi=-215; R=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); vper=sqrt(2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight; vpar=vpar/M.vlight; vpar=vpar(real(vper)~=0); vper=vper(real(vper)~=0); xlimits=xlim; ylimits=ylim; if(length(vper)>0) p(3)=plot(vpar,vper,'b-','displayname','Loss parabolla simul'); hold on plot(-vpar,vper,'b-') end vpar=linspace(1,M.vlight,1000); vper=sqrt(-2*M.qe/M.me*deltaphicomp/(R-1)+vpar.^2/(R-1))/M.vlight; vpar=vpar/M.vlight; vpar=vpar(real(vper)~=0); vper=vper(real(vper)~=0); if(length(vper)>0) p(4)=plot(vpar,vper,'k--','displayname','Loss parabolla prediction'); hold on plot(-vpar,vper,'k--') end xlim(xlimits) ylim(ylimits) legend(p) axis equal % subplot(2,2,2) % scatter(Zp,Vpar,'.','displayname','Init') % hold on % scatter(Zend,Vparend,'.','displayname','End') % xlabel('z [m]') % ylabel('\beta_{par}') % % subplot(2,2,3) % scatter(Vperp,Rp,'.','displayname','Init') % hold on % scatter(Vperpend,Rend,'.','displayname','End') % xlabel('\beta_\perp') % ylabel('R [m]') ax1=subplot(1,2,2); h=contourf(ax1,M.zgrid,M.rgrid,Ndistrib); hold on [r,z]=find(Ndistrib~=0); xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))]) ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))]) xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') c = colorbar(ax1); c.Label.String= 'n [m^{-3}]'; Zx=(M.zgrid(Zindex)); Rx=(M.rgrid(Rindex)); plot(ax1,Zx,Rx,'rx','Markersize',12); sgtitle(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV] R=%1.3g',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn)*M.phinorm,M.Rcurv)) M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%dpos',M.folder,M.name,Rindex,Zindex),[16,12]) %% f=figure; p(1)=scatter(Vpar,Vperp,'.','displayname','Initial'); hold on p(2)=scatter(Vparend,Vperpend,'.','displayname','Final'); xlabel('v_{par}/c') ylabel('v_\perp/c') zleftlim=1; vpar=linspace(1,M.vlight,1000); dN=1e13*M.weight; dN=0; rp=9.85e-3; rm=7e-3; V=2*pi*(rp^2-rm^2); N=M.N(:,:,timestepNend); levels=M.rAthet(Rindex,Zindex)*[1 1]; rAthetposend=find(levels(1)>=M.rAthet(:,zleftlim),1,'last'); rleft=sqrt(levels(1)/(M.rAthet(rAthetposend,1)/M.rgrid(rAthetposend)^2)); rposleft=find(M.rgrid>=rleft,1,'first'); n2=mean(mean(N(rposleft,[1 end],:),3)); vpar2=(dN/V/n2)^2; if n2==0 vpar2=0; end b=M.rgrid(end); a=M.rgrid(1); vd1=((M.potout-M.potinn)*M.phinorm/M.rgrid(Rindex))/M.Bz(Zindex,Rindex)/log(b/a); vd1=-M.Er(Rindex,Zindex,1)/M.Bz(Zindex,Rindex); %vd1=-M.Er(Rindex,Zindex,timestepNend)/M.Bz(Zindex,Rindex); vinit=sqrt(M.kb*10000/M.me); %vinit=sqrt(M.qe*80/M.me) %vinit=sqrt(120*M.qe/M.me) %vd2=((M.potout-M.potinn)*M.phinorm/M.rgrid(25))/M.Bz(1,25)/log(b/a); [Zmesh,Rmesh]=meshgrid(M.zgrid,M.rAthet(:,floor(M.nz/2)+1)); Rcurv=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); deltaphicomp=-0.5*M.me/M.qe*(vpar2-(vd1+vinit)^2*(1-M.Rcurv)) ratioparper=vpar2/((vd1+vinit)^2*(1-M.Rcurv)) Ndistrib=M.N(:,:,timestepNend); model=M.potentialwellmodel(timestepNend); z=model.z; r=model.r; pot=model.pot; rathet=model.rathet; Zeval=[M.zgrid(zleftlim) M.zgrid(Zindex)]; Psieval=M.rAthet(Rindex,Zindex)*[1 1]; phis=griddata(Zmesh,M.rAthet,M.pot(:,:,end),Zeval,Psieval,'natural'); deltaphi=-diff(phis) %deltaphi=-215; R=griddata(Zmesh,M.rAthet,M.B',M.zgrid(zleftlim),M.rAthet(Rindex,Zindex),'natural')/M.B(Zindex,Rindex); vper=sqrt(2*M.qe/M.me*deltaphi/(R-1)+vpar.^2/(R-1))/M.vlight; vpar=vpar/M.vlight; vpar=vpar(real(vper)~=0); vper=vper(real(vper)~=0); xlimits=xlim; ylimits=ylim; if(length(vper)>0) p(3)=plot(vpar,vper,'b-','displayname','Simulation'); hold on plot(-vpar,vper,'b-') end vpar=linspace(1,M.vlight,1000); vper=sqrt(-2*M.qe/M.me*deltaphicomp/(R-1)+vpar.^2/(R-1))/M.vlight; vpar=vpar/M.vlight; vpar=vpar(real(vper)~=0); vper=vper(real(vper)~=0); if(length(vper)>0) p(4)=plot(vpar,vper,'k--','displayname','Prediction'); hold on plot(-vpar,vper,'k--') end xlim(xlimits) ylim(ylimits) legend(p) axis equal legend('location','northeast') title(sprintf('r=%1.3g [m] z=%1.3g [m] \\Delta\\phi=%1.3g[kV] R=%1.3g',M.rgrid(Rindex),M.zgrid(Zindex),(M.potout-M.potinn)*M.phinorm,M.Rcurv)) M.savegraph(f,sprintf('%s/%s_phasespaceR%dZ%d',M.folder,M.name,Rindex,Zindex),[12,8]) \ No newline at end of file diff --git a/matlab/espic2dhdf5.m b/matlab/espic2dhdf5.m index ea15419..33e6fe8 100644 --- a/matlab/espic2dhdf5.m +++ b/matlab/espic2dhdf5.m @@ -1,1158 +1,1171 @@ 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 + npart % Time evolution of the number of simulated particles %% 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 + 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 %% Curvilinear geometry conformgeom % stores if we use the conforming or nonconforming boundary conditions r_a r_b z_r z_0 r_0 r_r L_r L_z Interior above1 above2 interior walltype geomweight 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, ext] = fileparts(obj.filename); obj.filename=[obj.name,ext]; 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.r_a=h5readatt(obj.fullpath,'/data/input.00/geometry','r_a'); obj.r_b=h5readatt(obj.fullpath,'/data/input.00/geometry','r_b'); obj.z_r=h5readatt(obj.fullpath,'/data/input.00/geometry','z_r'); obj.r_r=h5readatt(obj.fullpath,'/data/input.00/geometry','r_r'); obj.r_0=h5readatt(obj.fullpath,'/data/input.00/geometry','r_0'); obj.z_0=h5readatt(obj.fullpath,'/data/input.00/geometry','z_0'); obj.above1=h5readatt(obj.fullpath,'/data/input.00/geometry','above1'); obj.above2=h5readatt(obj.fullpath,'/data/input.00/geometry','above2'); obj.interior=h5readatt(obj.fullpath,'/data/input.00/geometry','interior'); obj.walltype=h5readatt(obj.fullpath,'/data/input.00/geometry','walltype'); obj.conformgeom=false; catch obj.conformgeom=true; end 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'); + obj.npart= 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=min(abs(1/obj.omepe),abs(1/obj.omece)); 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 try obj.geomweight = h5read(obj.fullpath,'/data/input.00/geometry/geomweight'); obj.geomweight= reshape(obj.geomweight,length(obj.zgrid),length(obj.rgrid),[]); obj.geomweight = permute(obj.geomweight,[2,1,3]); catch obj.geomweight=ones(length(obj.rgrid),length(obj.zgrid),3); end geomweight=ones(length(obj.rgrid),length(obj.zgrid)); if(obj.normalized) obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, geomweight, 1); else obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), geomweight, 1); end obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 2); obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 3); obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 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, geomweight); else obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim, geomweight); 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 if obj.nlclassical vscale=obj.vnorm; else vscale=1; end obj.VR = h5partsquantity(obj.fullpath,'/data/part','UR',vscale); obj.VZ = h5partsquantity(obj.fullpath,'/data/part','UZ',vscale); obj.VTHET= h5partsquantity(obj.fullpath,'/data/part','UTHET',vscale); 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'); + obj.nbparts = h5read(obj.fullpath,'/data/part/Nparts'); catch + obj.nbparts=obj.npart; obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2); end if(obj.nbspecies >1) obj.species=h5parts.empty(obj.nbspecies,0); for i=2:obj.nbspecies obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj); 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) + if (sum(nbparts)>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 strcmp(indices{1},':') p=1:obj.VR.nparts; else p=indices{1}; end if strcmp(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 strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(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 strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(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 strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end - if size(indices,1)>2 + if size(indices,2)>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)); + if size(indices,2)>3 + gcs=indices{4}; + else + gcs=false; + end + Zp=obj.Z(p,t,track); + Rp=obj.R(p,t,track); + Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp); + Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp); + Bp=interp2(obj.zgrid,obj.rgrid,obj.B',Zp,Rp); + costhet=Bzp./Bp; + sinthet=Brp./Bp; + Vdrift=zeros(size(Zp)); + if gcs for j=1:length(t) - tfield=find(obj.t2d==obj.tpart(t(j))); + [~, tfield]=min(abs(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; + %posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); + timeErp=interp2(obj.zgrid,obj.rgrid,timeEr,Zp(:,j),Rp(:,j)); + timeEzp=interp2(obj.zgrid,obj.rgrid,timeEz,Zp(:,j),Rp(:,j)); + Vdrift(:,j)=(timeEzp.*Brp(:,j)-timeErp.*Bzp(:,j))./Bp(:,j).^2; + end 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) %% 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 + legend('location','Northwest') if obj.conformgeom xlim([obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))]) else xlim([obj.rgrid(1) obj.rgrid(end)]) end grid on ylimits=ylim(); if obj.conformgeom plot(obj.rgrid(1)*[1 1],ylimits,'k--') plot(obj.rgrid(end)*[1 1],ylimits,'k--') else plot(obj.r_a*[1 1],ylimits,'k--') if obj.walltype==0 plot(obj.r_b*[1 1],ylimits,'k--') elseif obj.walltype==1 rmax=obj.r_0-obj.r_r*sqrt(1-(obj.zgrid(zpos)-obj.z_0)^2/obj.z_r^2); plot(rmax*[1 1],ylimits,'k--') end end 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; + tmin=1; 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--') + plot(obj.t0d(tmin:tmax),abs(obj.npart(tmin:tmax)./obj.npart(1)*100),'d--') ylabel('Nparts %') %ylim([0 110]) catch end obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name)); end 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) + plot(obj.t0d,obj.npart*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); lvls=sort(unique([obj.rAthet(:,end)' obj.rAthet(1,:)])); contpoints=contourc(obj.zgrid,obj.rgrid,obj.rAthet,lvls(:)); [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('%s 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)) N0=obj.N(:,:,1); Nend=obj.N(:,:,timestep); geomw=obj.geomweight(:,:,1); if rcoord [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); pot=griddata(z,r,pot,Zmesh,Rmesh); pot(obj.geomweight(:,:,1)<=0)=0; surface(obj.zgrid(1:end),obj.rgrid(1:end),pot(1:end,1:end),'edgecolor','none','Displayname','Well') xlabel('z [m]') ylabel('r [m]') xlim([obj.zgrid(1) obj.zgrid(end)]) ylim([obj.rgrid(1) obj.rgrid(end)]) hold(gca, 'on') rdisp=obj.rgrid; else rdisp=sort(obj.rAthet(:,end)); [Zmesh,Rmesh]=meshgrid(obj.zgrid,rdisp); pot=griddata(z,rathet,pot,Zmesh,Rmesh); [Zinit,~]=meshgrid(obj.zgrid,obj.rAthet(:,1)); N0=griddata(Zinit,obj.rAthet,N0,Zmesh,Rmesh); Nend=griddata(Zinit,obj.rAthet,Nend,Zmesh,Rmesh); geomw=griddata(Zinit,obj.rAthet,geomw,Zmesh,Rmesh); pot(geomw<=0)=0; surface(obj.zgrid(1:end),rdisp,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(end,1)]) hold(gca, 'on') end maxdensend=max(Nend(:)); contourscale=0.1; Nend=(Nend-contourscale*maxdensend)/maxdensend; maxdens0=max(N0(:)); contourscale=0.1; N0=(N0-contourscale*maxdens0)/maxdens0; contour(obj.zgrid,rdisp,Nend,linspace(0,1-contourscale,5),'k--','linewidth',1.5,'Displayname','Cloud Boundaries'); contour(obj.zgrid,rdisp,N0,[0 0],'k-.','linewidth',1.5,'Displayname','Source boundaries'); contour(obj.zgrid,rdisp,geomw,[0 0],'w-.','linewidth',1.5,'Displayname','Vessel Boundaries'); c=colorbar; colormap('jet'); c.Label.String= 'depth [eV]'; if rcoord obj.savegraph(f,sprintf('%s/%s_wellr',obj.folder,obj.name)); else obj.savegraph(f,sprintf('%s/%s_wellpsi',obj.folder,obj.name)); end 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 function charge=totcharge(obj,fieldstep) % Integrates the density profile over the full volume to obtain % the total number of electrons in the volume N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder,ones(size(obj.CellVol)), 1, 1); charge=sum(sum(N(:,:,end))); end function [p, maxnb, c]=displayPhaseSpace(obj,type,partsstep, Rindex, Zindex,legendtext, figtitle, f, maxnb, c, gcs) if nargin<8 f=figure; f=gca; end if nargin<7 figtitle=sprintf('r=%1.2f [mm] z=%1.2f [mm] \\Delta\\phi=%1.1f[kV] R=%1.1f',obj.rgrid(Rindex)*1e3,obj.zgrid(Zindex)*1e3,(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv); end if nargin <6 legendtext=sprintf('t=%1.3g [s]',obj.tpart(partsstep)); end fieldstep=find(obj.tpart(partsstep(end))==obj.t2d,1); if nargin>=10 ctemp=c; n=zeros(length(c{1}),length(c{2})); else nbins=15; n=zeros(nbins); end if nargin <11 gcs=true; end for i=1:length(partsstep) odstep=find(obj.tpart(partsstep(i))==obj.t0d); - nbp=min(obj.R.nparts,obj.nbparts(odstep)); + nbp=min(obj.R.nparts,obj.nbparts(partsstep(i))); Rp=obj.R(1:nbp,partsstep(i),false); Zp=obj.Z(1:nbp,partsstep(i),false); deltar=obj.dr(2)/2; deltarm=obj.rgrid(Rindex)-sqrt(obj.rgrid(Rindex)^2-deltar^2-2*obj.rgrid(Rindex)*deltar); deltaz=obj.dz/2; Indices=Rp>=obj.rgrid(Rindex)-deltarm & Rp=obj.zgrid(Zindex)-deltaz & Zp0 + if sum(obj.nbparts)>0 obj.R = h5partsquantity(obj.fullpath,hdf5group,'R'); obj.Z = h5partsquantity(obj.fullpath,hdf5group,'Z'); obj.THET = h5partsquantity(obj.fullpath,hdf5group,'THET'); obj.VR = h5partsquantity(obj.fullpath,hdf5group,'UR',obj.vnorm); obj.VZ = h5partsquantity(obj.fullpath,hdf5group,'UZ',obj.vnorm); obj.VTHET= h5partsquantity(obj.fullpath,hdf5group,'UTHET',obj.vnorm); obj.q= h5readatt(obj.fullpath,hdf5group,'q'); obj.m= h5readatt(obj.fullpath,hdf5group,'m'); obj.weight= h5readatt(obj.fullpath,hdf5group,'weight'); obj.partepot = h5partsquantity(filename,hdf5group,'pot',obj.q); try obj.Rindex=h5partsquantity(obj.fullpath,hdf5group,'Rindex'); obj.Zindex=h5partsquantity(obj.fullpath,hdf5group,'Zindex'); obj.partindex=h5partsquantity(obj.fullpath,hdf5group,'partindex'); catch end end % try % obj.partindex = h5read(obj.fullpath,sprintf('%s/partindex',hdf5group)); % partindex=obj.partindex; % partindex(partindex==-1)=NaN; % %[~,Indices]=sort(partindex,'ascend'); % Indices=obj.partindex; % for i=1:size(Indices,2) % obj.H(Indices(:,i),i)=obj.H(:,i); % obj.P(Indices(:,i),i)=obj.P(:,i); % % obj.R(Indices(:,i),i)=obj.R(:,i); % obj.Z(Indices(:,i),i)=obj.Z(:,i); % obj.Rindex(Indices(:,i),i)=obj.Rindex(:,i); % obj.Zindex(Indices(:,i),i)=obj.Zindex(:,i); % obj.VR(Indices(:,i),i)=obj.VR(:,i); % obj.VZ(Indices(:,i),i)=obj.VZ(:,i); % obj.VTHET(Indices(:,i),i)=obj.VTHET(:,i); % obj.THET(Indices(:,i),i)=obj.THET(:,i); % end % catch % end % clear partindex; end function quantity=H(obj,varargin) if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.VR.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end if size(indices,2)>2 track=indices{3}; else track=false; end quantity=0.5*obj.m*(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,varargin) if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end if size(indices,2)>2 track=indices{3}; else track=false; end quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.m+obj.q*obj.parent.Atheta(obj.R(p,t,track),obj.Z(p,t,track))); end function quantity=Vpar(obj,varargin) %Vpar Computes the parallel velocity for the particle indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end if size(indices,2)>2 track=indices{3}; else track=false; end Zp=obj.Z(p,t,track); Rp=obj.R(p,t,track); Bzp=interp2(obj.zgrid,obj.rgrid,obj.parent.Bz',Zp,Rp); Brp=interp2(obj.zgrid,obj.rgrid,obj.parent.Br',Zp,Rp); Bp=interp2(obj.zgrid,obj.rgrid,obj.parent.B',Zp,Rp); costhet=Bzp./Bp; sinthet=Brp./Bp; 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 strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end if size(indices,2)>2 track=indices{3}; else track=false; end if size(indices,2)>3 gcs=indices{4}; else gcs=false; end Zp=obj.Z(p,t,track); Rp=obj.R(p,t,track); Bzp=interp2(obj.zgrid,obj.rgrid,obj.parent.Bz',Zp,Rp); Brp=interp2(obj.zgrid,obj.rgrid,obj.parent.Br',Zp,Rp); Bp=interp2(obj.zgrid,obj.rgrid,obj.parent.B',Zp,Rp); costhet=Bzp./Bp; sinthet=Brp./Bp; Vdrift=zeros(size(Zp)); if gcs for j=1:length(t) [~, tfield]=min(abs(obj.parent.t2d-obj.tpart(t(j)))); timeEr=obj.parent.Er(:,:,tfield); timeEz=obj.parent.Ez(:,:,tfield); %posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); timeErp=interp2(obj.zgrid,obj.rgrid,timeEr,Zp(:,j),Rp(:,j)); timeEzp=interp2(obj.zgrid,obj.rgrid,timeEz,Zp(:,j),Rp(:,j)); Vdrift(:,j)=(timeEzp.*Brp(:,j)-timeErp.*Bzp(:,j))./Bp(:,j).^2; end 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 sref = subsref(obj,s) % obj(i) is equivalent to obj.Data(i) switch s(1).type case '.' if(strcmp(s(1).subs,'H')) sref=H(obj,s(2).subs); elseif(strcmp(s(1).subs,'P')) sref=P(obj,s(2).subs); else sref=builtin('subsref',obj,s); end case '()' sref=builtin('subsref',obj,s); case '{}' error('MYDataClass:subsref',... 'Not a supported subscripted reference') end end end end \ No newline at end of file