diff --git a/matlab/CheckDiocStability.m b/matlab/CheckDiocStability.m index 78b5e69..61e76b2 100644 --- a/matlab/CheckDiocStability.m +++ b/matlab/CheckDiocStability.m @@ -1,129 +1,136 @@ t2dlength=size(M.t2d,1); fieldstart=max(1,t2dlength-200);%floor(0.95*t2dlength); deltat=t2dlength-fieldstart; maxl=12; %[~, rgridend]=min(abs(M.rgrid-0.045)); rgridend=sum(M.nnr(1:2)); [~, name, ~] = fileparts(M.file); dens=mean(M.N(:,:,fieldstart:end),3); [R,Z]=meshgrid(M.rgrid,M.zgrid); Rinv=1./R; Rinv(isnan(Rinv))=0; VTHET=mean(M.fluidUTHET(:,:,fieldstart:end),3); omegare=(VTHET.*Rinv'); Er=mean(M.Er(:,:,fieldstart:end),3); Ez=mean(M.Ez(:,:,fieldstart:end),3); vdrift=(M.Br'.*Ez-Er.*M.Bz')./M.B'.^2; omegadrift=(vdrift.*Rinv'); zindices=linspace(length(M.zgrid)/4,3/4*length(M.zgrid),11); for zid=4:8%1:length(zindices)%;%floor(length(M.zgrid)/2); zindex=floor(zindices(zid)); rindex=M.geomweight(:,zindex,1)>=0; n=dens(rindex,zindex); omegar=omegadrift(rindex,zindex); rgrid=M.rgrid(rindex); - Bz=M.Bz(zindex,M.nnr(1)+floor(M.nnr(2)/2)); + Bz=mean(M.Bz(zindex,rindex),2); zfolder=sprintf('diocz_%03d',zindex); if ~isfolder(zfolder) mkdir(zfolder) end - wce=M.qsim/M.msim*Bz; - wpe2=M.qsim^2/M.msim/M.eps_0*n; + wce=abs(M.qsim/M.msim*Bz); + wpe2=M.qsim^2/(M.msim*M.eps_0*M.weight)*n; wnum=complex(zeros(20,maxl)); rnumlength=max(length(rgrid)+2,4003); phinum=complex(zeros(20,rnumlength,maxl)); rnum=zeros(rnumlength,maxl); for l=1:maxl [w, phin, r] = diocotronstab(l, rgrid, n, omegar, wce, 20); % different profiles of dphi [~,ind]=sort(abs(imag(w)),'descend'); f=figure; subplot(3,1,2) hold on ax=gca; ax.LineStyleOrder={'-','--','-.',':'}; subplot(3,1,3) hold on ax=gca; ax.LineStyleOrder={'-','--','-.',':'}; for i=1:10 if abs(imag(w(ind(i)))/real(w(ind(i))))>1e-10 freq=w(ind(i)); disphinum=phin(:,ind(i)); disphinum=disphinum/max(abs(disphinum)); subplot(3,1,2) plot(r(2:end-1),real(disphinum),'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq))) subplot(3,1,3) plot(r(2:end-1),imag(disphinum),'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq))) end end - subplot(3,1,1) - if (~isreal(w)) - title(sprintf('Unstable case l=%d, \\omega_d=%1.3g',l,max(wpe2/wce))) - else - title(sprintf('Stable case l=%d, \\omega_d=%1.3g',l,max(wpe2/wce))) - end + wnum(:,l)=w(ind(1:20)); phinum(:,:,l)=padarray(phin(:,ind(1:20))',[0 rnumlength-size(phin,1)], 0, 'post'); rnum(:,l)=padarray(r,[0 rnumlength-length(r(:))],0, 'post')'; + + subplot(3,1,1) plots=gobjects(2,1); plots(1)=plot(rgrid,omegar,'displayname','\omega_r'); hold on ylabel('\omega_r [s^{-1}]') yyaxis right plots(2)=plot(rgrid,n,'displayname','n'); legend(plots,{'\omega_r','n'},'location','northwest') ylabel('n [m^{-3}]') ylabel('R [m]') + xlim([min(r) max(r)]) + if (abs(imag(w(ind(1)))/real(w(ind(1))))>1e-10) + title(sprintf('Unstable case l=%d, \\omega_d=%1.3g',l,max(wpe2/(2*wce)))) + else + title(sprintf('Stable case l=%d, \\omega_d=%1.3g',l,max(wpe2/(2*wce)))) + end + set(gca,'fontsize', 12) subplot(3,1,2) set(gca,'fontsize', 12) xlabel('R [m]') ylabel('Re(\delta\phi) [a.u.]') yyaxis right grid on + xlim([min(r) max(r)]) subplot(3,1,3) set(gca,'fontsize', 12) xlabel('R [m]') ylabel('Im(\delta\phi) [a.u.]') grid on + xlim([min(r) max(r)]) f.PaperUnits='centimeters'; f.PaperSize=[16,20]; legend('location','northwest') print(f,sprintf('%s/%s_dioc_l%d',zfolder,M.name,l),'-dpdf','-fillpage') savefig(f,sprintf('%s/%s_dioc_l%d',zfolder,M.name,l)) close(f) % values of w found numerically f=figure; hold on plot(w, 'rx') set(gca,'fontsize', 12) if (~isreal(w)) xlabel('\omega_r [rad s^{-1}]') ylabel('\omega_i [rad s^{-1}]') - title(sprintf('Unstable case l=%d, \\omega_d=%1.3g',l,max(wpe2/wce))) + title(sprintf('Unstable case l=%d, \\omega_d=%1.3g',l,max(wpe2/(2*wce)))) else xlabel('') ylabel('\omega_r [rad s^{-1}]') - title(sprintf('Stable case l=%d, \\omega_d=%1.3g',l,max(wpe2/wce))) + title(sprintf('Stable case l=%d, \\omega_d=%1.3g',l,max(wpe2/(2*wce)))) end grid on f.PaperUnits='centimeters'; f.PaperSize=[14,12]; print(f,sprintf('%s/%s_wdioc_l%d',zfolder,M.name,l),'-dpdf','-fillpage') savefig(f,sprintf('%s/%s_wdioc_l%d',zfolder,M.name,l)) drawnow close(f) end - savecase(zfolder,rnum,wnum,phinum,rgrid,wce,Er,Ez,n,omegare,vdrift,omegadrift,rgrid(1),rgrid(end),1:maxl); + lsteps=1:maxl; + savecase(zfolder,rnum,wnum,phinum,rgrid,wce,Er,Ez,n,omegare,vdrift,omegadrift,M.geomweight,rgrid(1),rgrid(end),lsteps); end -function savecase(folder,r,w,phinum,rgrid,wce,Er,Ez,n,omegare,vdrift,omegadrift,geomweight,a,b,l) - save(sprintf('%s/results',folder),'r','w','phinum','a','b','l','rgrid','wce','Er','Ez','n','omegare', 'vdrift','omegadrift','geomweight','phinum') +function savecase(folder,r,w,phinum,rgrid,wce,Er,Ez,n,omegare,vdrift,omegadrift,geomweight,a,b,lsteps) + save(sprintf('%s/results',folder),'r','w','phinum','a','b','lsteps','rgrid','wce','Er','Ez','n','omegare', 'vdrift','omegadrift','geomweight','phinum') end \ No newline at end of file diff --git a/matlab/CheckDiocStability_post.m b/matlab/CheckDiocStability_post.m index e96a06e..3a35219 100644 --- a/matlab/CheckDiocStability_post.m +++ b/matlab/CheckDiocStability_post.m @@ -1,98 +1,97 @@ zindices=linspace(length(M.zgrid)/4,3/4*length(M.zgrid),11); -for zid=5:7%1:length(zindices)%;%floor(length(M.zgrid)/2); +for zid=4:8%1:length(zindices)%;%floor(length(M.zgrid)/2); zindex=floor(zindices(zid)); zfolder=sprintf('diocz_%03d',zindex); if ~isfolder(zfolder) || ~isfile(sprintf('%s/results.mat',zfolder)) fprintf('Warning: result %s doesn''t exist\n',zfolder) continue end result=load(sprintf('%s/results.mat',zfolder)); if ~isfield(result,'geomweight') result.geomweight=M.geomweight(:,:,1); end rindex=M.geomweight(:,zindex,1)>=0; for l=result.l r=result.r(:,l); phin=result.phinum(:,:,l); w=result.w(:,l); n=result.n; omegar=result.omegadrift; rgrid=result.rgrid; omegar=omegar(rindex,zindex); a=min(rgrid); b=max(rgrid); idmin=find(r==a); idmax=find(r==b); f=figure; subplot(3,1,1) hold on set(gca,'fontsize', 18) xlabel('R [m]') ylabel('n [m^{-3}]') plots=gobjects(2,1); plots(1)=plot(rgrid,n,'displayname','n'); yyaxis right plots(2)=plot(rgrid,omegar,'displayname','\omega_r'); legend(plots,{'\omega_r','n'},'location','northwest') ylabel('\omega_r [Hz]') grid minor - ax=gca; ax.LineStyleOrder={'-','--','-.',':'}; subplot(3,1,[2,3]) hold on p=gobjects(); j=1; for i=1:size(phin,1) if imag(w(i))/abs(real(w(i)))>1e-10 freq=w(i); disphinum=phin(i,:); disphinum=disphinum/max(abs(disphinum)); yyaxis left p(j)=plot(r(idmin:idmax),[0 real(disphinum(idmin:idmax-2)) 0],'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq))); j=j+1; yyaxis right plot(r(idmin:idmax),[0 imag(disphinum(idmin:idmax-2)) 0],'displayname',sprintf('\\omega=%1.3g%+1.3gi',real(freq),imag(freq))) end end subplot(3,1,1) if (abs(imag(w(1))/real(w(1)))>1e-10) - title(sprintf('Unstable case l=%d',l)) + title(sprintf('Unstable case l=%d, \\omega_d=%1.3g',l,max(wpe2/wce))) else - title(sprintf('Stable case l=%d',l)) + title(sprintf('Stable case l=%d, \\omega_d=%1.3g',l,max(wpe2/wce))) end subplot(3,1,[2,3]) set(gca,'fontsize', 18) yyaxis left ylim([-1 1]) xlabel('R [m]') ylabel('real(\delta\phi) [a.u.]') grid on yyaxis right ylim([-1 1]) grid minor ylabel('imag(\delta\phi) [a.u.]') f.PaperUnits='centimeters'; f.PaperSize=[16,20]; legend('location','northwest') if length(p)>=1 legend(p) end subplot(3,1,1) ax1=gca; subplot(3,1,[2,3]) ax2=gca; drawnow posit1=get(ax1,'position'); posit2=get(ax2,'position'); posit2(3)=posit1(3); set(ax2,'position',posit2); print(f,sprintf('%s/%s_dioc_l%d',zfolder,M.name,l),'-dpdf','-fillpage') savefig(f,sprintf('%s/%s_dioc_l%d',zfolder,M.name,l)) close(f) end end \ No newline at end of file diff --git a/matlab/PotentialWell.m b/matlab/PotentialWell.m index de734ba..74da3ad 100644 --- a/matlab/PotentialWell.m +++ b/matlab/PotentialWell.m @@ -1,39 +1,39 @@ function [pot] = PotentialWell(Phi,Atheta,r,z,creategraph) -%PotentialWell Computes and display the potential well given the electrostatic potentila Phi +%PotentialWell Computes and display the potential well given the electrostatic potential Phi % and the magnetic vector potential Atheta if nargin <5 creategraph=true; end [rgrid,~]=meshgrid(r,z); rAthet=(rgrid.*Atheta)'; contpoints=contourc(z,r,rAthet,rAthet(:,1)'); [x,y,zcont]=C2xyz(contpoints); k=1; zdiff=[diff(zcont),0]; for j=1:size(x,2) %for i=1:size(pot,1) xloc=x{j}; yloc=y{j}; xloc=xloc(ylocr(1)); yloc=yloc(ylocr(1)); x{j}=xloc; y{j}=yloc; if(length(xloc)>1) pot{j}=interp2(z,r,Phi,xloc,yloc); pot{j}=pot{j}-min(pot{j}); end k=k+(zdiff(j)~=0); end x=cell2mat(x); y=cell2mat(y); pot=cell2mat(pot); [Z,R]=meshgrid(z,r); pot=griddata(x,y,pot,Z,R); if(creategraph) figure surface(z(2:end),r(2:end),-pot(2:end,2:end),'edgecolor','none') xlabel('r') ylabel('z') colorbar end end diff --git a/matlab/analyse_espicfields.m b/matlab/analyse_espicfields.m index 53e084a..5c71a18 100644 --- a/matlab/analyse_espicfields.m +++ b/matlab/analyse_espicfields.m @@ -1,387 +1,442 @@ % 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); fieldend=t2dlength;%23; deltat=t2dlength-max(1,t2dlength-800); %[~, 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 -f=figure('Name', sprintf('%sfluid_dens',name)); -ax1=gca; -h=surface(ax1,M.zgrid,M.rgrid,dens); -set(h,'edgecolor','none'); -if(M.conformgeom) - ylim(ax1,[M.rgrid(1) M.rgrid(rgridend)]) -else - ylim(ax1,[M.rgrid(1) M.rgrid(end)]) -end -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(fieldend),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)) +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'); -geomw=M.geomweight(:,:,1); -geomw(geomw<0)=0; -geomw(geomw>0)=max(max(dens)); [~, 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,M.rgrid,dens,'Displayname','n_e [m^{-3}]'); hold on; pot(M.geomweight(:,:,1)<0)=0; [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'); [~, hContour]=contourf(ax1,M.zgrid,M.rgrid,geomw,2); drawnow; 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 legend({'n_e [m^{-3}]','Equipotentials [kV]','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 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)]) 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]) 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 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/diocotronstab.m b/matlab/diocotronstab.m index ac222ae..66d786a 100644 --- a/matlab/diocotronstab.m +++ b/matlab/diocotronstab.m @@ -1,103 +1,103 @@ function [w, phin, r] = diocotronstab(l, R, n, wre, wce, nbeigs) if nargin<6 nbeigs=20; end % Necessary physical constants % electron charge and mass e = 1.602176634e-19; me = 9.1093837015e-31; % vacuum permittivity eps_0 = 8.854188e-12; % declaration of wpe2 (vector for numerical solution, renormalised on wce) wpe2 = n*(e^2)/(me*eps_0); % scale factor for the frequency - wscale=abs(wce); + wscale=max(abs(wpe2/(2*wce))); % construction of r, new vector of distances (chooses the thinner mesh and % interpolate all variables on that mesh delta = diff(R); delta=delta(delta>0); x = unique(delta); N = numel(x); count = zeros(N,1); for k = 1:N count(k) = sum(delta==x(k)); end x(count<2)=[]; rmin = min(x); nR = floor((max(R)-min(R))/rmin) + 1; - if nR<2000 - nR=2001; + if nR<3000 + nR=3001; end r = linspace(R(1),R(end),nR); b=max(R); rmin=rmin/b; r=r/b; R=R/b; % interpolation % of n, wre and wpe2 --> nin, wrin, wpe2in wrin = interp1(R(2:end),wre(2:end),r,'spline'); wpe2in = interp1(R,wpe2,r); wpe2in=wpe2in/wscale^2; wrin=wrin/wscale; wce=wce/wscale; % figure % hold on % plot(r*b,wpe2in) % plot(R*b,wpe2) % title('wpe2') % % figure % hold on % plot(r*b,wrin) % plot(R*b,wre/wce) % title('wr') % finite difference to compute d/dr(wpe^2) Dwpe2 = zeros([nR 1]); Dwpe2(2:(nR-1)) = (wpe2in(3:end)-wpe2in(1:(nR-2)))/(2*rmin); % figure; % plot(r,Dwpe2) % construction of matrices A & B such that A*phi = w*B*phi A = zeros(nR,nR); B = zeros(nR,nR); % finite differences in A & B + other terms of the eigenvalue % equation for i = 2:1:(nR-1) A(i,i) = l*( - 2*r(i)*wrin(i)/(rmin)^2 - l^2*wrin(i)/r(i) - Dwpe2(i)/wce); A(i,i-1) = l*wrin(i)*(r(i)/(rmin)^2-1/(2*rmin)); A(i,i+1) = l*wrin(i)*(r(i)/(rmin)^2 + 1/(2*rmin)); B(i,i) = - 2*r(i)/(rmin)^2 - (l^2/r(i)); B(i,i-1) = r(i)/(rmin)^2-1/(2*rmin); B(i,i+1) = 1/(2*rmin) + r(i)/(rmin)^2; end A=complex(A(2:(end-1),2:(end-1))); B=complex(B(2:(end-1),2:(end-1))); % Solution og the eigenvalue problem [phin,w] = eig(A, B); % A=sparse(complex(A(2:(end-1),2:(end-1)))); % B=sparse(complex(B(2:(end-1),2:(end-1)))); % [phin, w]=eigs(A,B,nbeigs,'bothendsimag','subspacedimension',3*nbeigs+1); % renormalisation of w and r w=diag(w); w = w*wscale; r=r*b; end diff --git a/matlab/dispespicFields.m b/matlab/dispespicFields.m index 300170e..fd2bd24 100644 --- a/matlab/dispespicFields.m +++ b/matlab/dispespicFields.m @@ -1,244 +1,272 @@ -function dispespicFields(M,logdensity,showgrid,fixed) +function dispespicFields(M,logdensity,showgrid,fixed,parper) %dispespicFields Allows to display the time evolution of the density, electric potential and electric fields % M is of class espic2dhdf5 and contains the simulation results fieldstep=1; if nargin <2 logdensity=false; showgrid=false; fixed=false; end if nargin <3 showgrid=false; fixed=false; end if nargin <4 fixed=false; end +if nargin <5 + parper=false; +end fixed=fi(fixed); f=uifigure('Name',sprintf('Grid data %s',M.name)); mf=uipanel(f,'Position',[5 50 f.Position(3)-10 f.Position(4)-55]); mf.AutoResizeChildren='off'; m=uipanel(f,'Position',[5 5 f.Position(3)-10 40]); sgtitle(mf,sprintf('step=%d t=%0.5e s',fieldstep*M.it1,M.t2d(fieldstep))) sld = uislider(m,'Position',[10 30 0.6*m.Position(3) 3]); sld.Value=fieldstep; sld.Limits=[1 size(M.t2d,1)]; edt = uieditfield(m,'numeric','Limits',[1 size(M.t2d,1)],'Value',1); edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20]; edt.RoundFractionalValues='on'; MaxN=0; Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save'); Play=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play'); Pause=uibutton(m,'Position',[Play.Position(1)+Play.Position(3)+10 5 40 20],'Text', 'Pause'); %Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause'); stop=false; sld.ValueChangingFcn={@updatefigdata,edt,mf}; edt.ValueChangedFcn={@updatefigdata,sld,mf}; Printbt.ButtonPushedFcn={@plotGridButtonPushed}; Play.ButtonPushedFcn={@plotPlayButtonPushed}; Pause.ButtonPushedFcn={@PauseButtonPushed}; set(f,'KeyPressFcn',{ @onKeyDown,sld,edt,mf}) PlotEspic2dgriddata(mf,M,fieldstep); function plotPlayButtonPushed(btn,ax) stop=false; i=sld.Value; while ~stop edt.Value=i; sld.Value=i; updatesubplotsdata(i,mf); pause(0.01) i=sld.Value; i=i+10; if(i>sld.Limits(2)) stop=true; end end end function PauseButtonPushed(btn,ax) stop = true; 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))); updatefigdata(slider, event, editfield ,fig) end end function PlotEspic2dgriddata(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.it1,M.t2d(fieldstep))) ax1=subplot(2,2,1,'Parent',fig); sf=surface(ax1,M.zgrid,M.rgrid,M.N(:,:,fieldstep),'edgecolor','none'); if(showgrid) set(sf,'edgecolor','black') end 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') c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; %c.Limits=[0 max(M.N(:))]; if(isboolean(fixed) && fixed) climits=caxis(ax1); MaxN=climits(2); elseif(~isboolean(fixed)) MaxN=fixed.data; caxis(ax1,[-Inf MaxN]); end view(ax1,2) hold(ax1, 'on') contour(ax1,M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'r-.','linewidth',1.5,'Displayname','Boundaries'); if logdensity set(ax1,'zscale','log') set(ax1,'colorscale','log') end ax2=subplot(2,2,2,'Parent',fig); pot=M.pot(:,:,fieldstep); pot(M.geomweight(:,:,1)<0)=0; contourf(ax2,M.zgrid,M.rgrid,pot); %plot(ax2,M.zgrid,data.pot(:,5)) xlabel(ax2,'z [m]') ylabel(ax2,'r [m]') colormap(ax2,'jet') c = colorbar(ax2); c.Label.String= '\Phi [V]'; title(ax2,'Es potential') grid(ax2, 'on') hold(ax2, 'on') + contour(ax2,M.zgrid,M.rgrid,M.geomweight(:,:,1),[0 0],'w-.','linewidth',1.5,'Displayname','Boundaries'); 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(ax2,M.zgrid,M.rgrid,Blines',real(levels),'m-.','linewidth',1.5,'Displayname','Magnetic field lines'); ax3=subplot(2,2,3,'Parent',fig); - Ez=M.Epar(fieldstep); + if parper + Ez=M.Epar(fieldstep); + else + Ez=M.Ez(:,:,fieldstep); + end Ez(M.geomweight(:,:,1)<0)=0; contourf(ax3,M.zgrid,M.rgrid,Ez) xlabel(ax3,'z [m]') ylabel(ax3,'r [m]') c = colorbar(ax3); c.Label.String= 'E_{z} [V/m]'; - title(ax3,'Parallel Electric field') + if parper + title(ax3,'Parallel Electric field') + else + title(ax3,'Axial Electric field') + end grid(ax3, 'on') ax4=subplot(2,2,4,'Parent',fig); - Er=M.Eperp(fieldstep); + if parper + Er=M.Eperp(fieldstep); + else + Er=M.Er(:,:,fieldstep); + end Er(M.geomweight(:,:,1)<0)=0; contourf(ax4,M.zgrid,M.rgrid,Er) xlabel(ax4,'z [m]') ylabel(ax4,'r [m]') c = colorbar(ax4); c.Label.String= 'E_{r} [V/m]'; - title(ax4,'Perpendicular Electric field') + if parper + title(ax4,'Perpendicular Electric field') + else + title(ax4,'Radial Electric field') + end grid(ax4, 'on') linkaxes([ax1,ax2,ax3,ax4]); end function plotGridButtonPushed(btn,ax) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here f=figure(); PlotEspic2dgriddata(f,M,sld.Value); f.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); print(f,sprintf('%sGrid%d',name,sld.Value),'-dpdf','-fillpage') end function updatefigdata(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; updatesubplotsdata(fieldstep, fig); end function updatesubplotsdata(fieldstep, fig) sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,M.t2d(floor(fieldstep)))) %% update Position histogram ax1=fig.Children(end); dens=M.N(:,:,fieldstep); ax1.Children(end).ZData=dens; ax1.Children(end).CData=dens; if(isboolean(fixed) && fixed) climits=caxis(ax1); MaxN=climits(2); nmax=max(dens(:)); if(nmax>MaxN) MaxN=nmax; end caxis(ax1,[0 MaxN]); elseif(~isboolean(fixed)) MaxN=fixed.data; caxis(ax1,[-Inf MaxN]); end % view(ax1,2) %% update ES Potential pot=M.pot(:,:,fieldstep); pot(M.geomweight(:,:,1)<0)=0; fig.Children(end-2).Children(end).ZData=pot; %% update Radial electric field - Ez=M.Epar(fieldstep); + if parper + Ez=M.Epar(fieldstep); + else + Ez=M.Ez(:,:,fieldstep); + end Ez(M.geomweight(:,:,1)<0)=0; fig.Children(end-4).Children(1).ZData=Ez; %% update Axial electric field - Er=M.Eperp(fieldstep); + if parper + Er=M.Eperp(fieldstep); + else + Er=M.Er(:,:,fieldstep); + end Er(M.geomweight(:,:,1)<0)=0; fig.Children(end-6).Children(1).ZData=Er; drawnow limitrate end end diff --git a/matlab/distribution_function.m b/matlab/distribution_function.m index ec0c14e..ff56a28 100644 --- a/matlab/distribution_function.m +++ b/matlab/distribution_function.m @@ -1,365 +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=29 -% Zindex=193; +%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; - -Rp=M.R(:,1,false); -Zp=M.Z(:,1,false); +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)); +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); +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 5687d8f..24987ec 100644 --- a/matlab/espic2dhdf5.m +++ b/matlab/espic2dhdf5.m @@ -1,1098 +1,1152 @@ 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 + % 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 + 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 %% 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 + % 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'); 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 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; + obj.species=h5parts.empty(obj.nbspecies,0); 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 + 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) 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}==':' + if strcmp(indices{1},':') p=1:obj.VR.nparts; else p=indices{1}; end - if indices{2}==':' + 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 indices{1}==':' + if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end - if indices{2}==':' + 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 indices{1}==':' + if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end - if indices{2}==':' + 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); + 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}==':' + if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end - if indices{2}==':' + 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); + 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) %% 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 + %% 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 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 + 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 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 + %% 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 + 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 + % 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}); + 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 - k=k+(zdiff(j)~=0); - end - potfinal(:,i)=cell2mat(pot); + 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 + 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); - zleftlim=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); - Rp=obj.R(1:obj.nbparts(odstep),partsstep(i),false); - Zp=obj.Z(1:obj.nbparts(odstep),partsstep(i),false); + nbp=min(obj.R.nparts,obj.nbparts(odstep)); + 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 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 indices{1}==':' + if strcmp(indices{1},':') p=1:obj.VR.nparts; else p=indices{1}; end - if indices{2}==':' + 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 indices{1}==':' + if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end - if indices{2}==':' + 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 indices{1}==':' + if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end - if indices{2}==':' + 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); - Rind=obj.Rindex(p,t,track); - posind=sub2ind(size(obj.parent.B),Zind,Rind); - costhet=obj.parent.Bz(posind)./obj.parent.B(posind); - sinthet=obj.parent.Br(posind)./obj.parent.B(posind); + 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 indices{1}==':' + if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end - if indices{2}==':' + 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); - Rind=obj.Rindex(p,t,track); - posind=sub2ind(size(obj.parent.B),Zind,Rind); - costhet=obj.parent.Bz(posind)./obj.parent.B(posind); - sinthet=obj.parent.Br(posind)./obj.parent.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.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)); - Vdrift(:,j)=(timeEz(posindE).*obj.parent.Br(posind(:,j))-timeEr(posindE).*obj.parent.Bz(posind(:,j)))./obj.parent.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 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 diff --git a/matlab/h5partsquantity.m b/matlab/h5partsquantity.m index accfdd3..e06c766 100644 --- a/matlab/h5partsquantity.m +++ b/matlab/h5partsquantity.m @@ -1,85 +1,85 @@ classdef h5partsquantity properties filename group dataset nparts nt scale end methods function obj=h5partsquantity(filename, group, dataset, scale) obj.filename=filename; obj.dataset=dataset; obj.group=group; obj.nparts=h5info(filename,sprintf('%s/%s',group,dataset)).Dataspace.Size(1); obj.nt=h5info(filename,sprintf('%s/%s',group,dataset)).Dataspace.Size(end); if nargin < 4 obj.scale=1; else obj.scale=scale; end end function ind=end(obj,k,n) switch k case 1 ind=obj.nparts; case 2 ind=obj.nt; case default error('Invalid number of dimensions'); end end function sref = subsref(obj,s) % obj(i) is equivalent to obj.Data(i) switch s(1).type case '.' if(strcmp(s(1).subs,'val')) sref = obj.(s(1).subs)(s(2).subs); else sref=builtin('subsref',obj,s); end case '()' sref = val(obj,s(1).subs); case '{}' error('MYDataClass:subsref',... 'Not a supported subscripted reference') end end function quantity=val(obj,indices) - if indices{1}==':' + if strcmp(indices{1},':') p=1:obj.nparts; else p=indices{1}; end - if indices{2}==':' + if strcmp(indices{2},':') t=1:obj.nt; else t=indices{2}; end if size(indices,2)>2 track=indices{3}; else track=false; end quantity=zeros(obj.nparts,length(t)); if track for i=1:length(t) indices = h5read(obj.filename, sprintf('%s/%s',obj.group,'partindex'),[1 t(i)],[Inf 1]); indices(indices<=0)=max(indices)+1; quantity(indices,i) = h5read(obj.filename, sprintf('%s/%s',obj.group,obj.dataset),[1 t(i)],[Inf 1]); end else for i=1:length(t) quantity(:,i) = h5read(obj.filename, sprintf('%s/%s',obj.group,obj.dataset),[1 t(i)],[Inf 1]); end end quantity=quantity(p,:)*obj.scale; end end end diff --git a/matlab/splinequantity.m b/matlab/splinequantity.m index adfe616..4aece42 100644 --- a/matlab/splinequantity.m +++ b/matlab/splinequantity.m @@ -1,123 +1,123 @@ classdef splinequantity < h5quantity properties knotsr knotsz femorder index postscale end methods function obj=splinequantity(filename, dataset, knotsr, knotsz, femorder, scale, postscale, index) if nargin<6 scale=1; end obj=obj@h5quantity(filename, dataset, length(knotsr)-2*femorder(2), length(knotsz)-2*femorder(1), scale); obj.knotsr=knotsr; obj.knotsz=knotsz; obj.femorder=double(femorder); if nargin < 7 obj.postscale=ones(length(knotsr)-2*femorder(2), length(knotsz)-2*femorder(1)); else obj.postscale=postscale end if nargin < 8 obj.index=-1; else obj.index=index; end end function quantity=coeffs(obj,indices) - if indices{1}==':' + if strcmp(indices{1},':') r=1:obj.nr+obj.femorder(2)-1; else r=indices{1}; end - if indices{2}==':' + if strcmp(indices{2},':')' z=1:obj.nz+obj.femorder(1)-1; else z=indices{2}; end - if indices{3}==':' + if strcmp(indices{3},':') t=1:obj.nt; else t=indices{3}; end temp=zeros((obj.nz+obj.femorder(1)-1)*(obj.nr+obj.femorder(2)-1),length(t)); if obj.index ~= -1 if(sum(diff(diff(t))) == 0 && length(t)>1) stride=t(2)-t(1); temp = h5read(obj.filename, obj.dataset,[obj.index 1 t(1)],[1 Inf length(t)],[1 1 stride]); else for i=1:length(t) temp(:,i) = h5read(obj.filename, obj.dataset,[obj.index 1 t(i)],[1 Inf 1]); end end else if(sum(diff(diff(t))) == 0 && length(t)>1) stride=t(2)-t(1); temp = h5read(obj.filename, obj.dataset,[obj.index 1 t(1)],[1 Inf length(t)],[1 1 stride]); else for i=1:length(t) temp(:,i) = h5read(obj.filename, obj.dataset,[1 t(i)],[Inf 1]) ; end end end temp=reshape(temp,obj.nz+obj.femorder(1)-1,obj.nr+obj.femorder(2)-1,[]); temp=permute(temp,[2,1,3]); quantity=temp(r,z,:)*obj.scale; end function quantity=val(obj,indices) - if indices{1}==':' + if strcmp(indices{1},':') r=1:obj.nr; else r=indices{1}; end - if indices{2}==':' + if strcmp(indices{2},':') z=1:obj.nz; else z=indices{2}; end - if indices{3}==':' + if strcmp(indices{3},':') t=1:obj.nt; else t=indices{3}; end count=length(t); temp=obj.coeffs({':',':',t}); quantity=zeros(length(r),length(z),count); for i=1:size(temp,3) quantity(:,:,i)=obj.postscale(r,z).*fnval(spmak({obj.knotsr,obj.knotsz},temp(:,:,i)),{obj.knotsr(r+obj.femorder(2)),obj.knotsz(z+obj.femorder(1))}); end end function quantity=der(obj,indices) - if indices{1}==':' + if strcmp(indices{1},':') r=1:obj.nr; else r=indices{1}; end - if indices{2}==':' + if strcmp(indices{2},':') z=1:obj.nz; else z=indices{2}; end - if indices{3}==':' + if strcmp(indices{3},':') t=1:obj.nt; else t=indices{3}; end order=indices{4}; count=length(t); temp=obj.coeffs({':',':',t}); quantity=zeros(length(r),length(z),count); for i=1:size(temp,3) quantity(:,:,i)=fnval(fnder(spmak({obj.knotsr,obj.knotsz},temp(:,:,i)), order),{obj.knotsr(r+obj.femorder(2)),obj.knotsz(z+obj.femorder(1))}); end end end end diff --git a/src/beam_mod.f90 b/src/beam_mod.f90 index 0a09af9..3019828 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,2095 +1,2088 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: beam ! !> @author !> Guillaume Le Bars EPFL/SPC !> Patryk Kaminski EPFL/SPC !> Trach Minh Tran EPFL/SPC ! ! DESCRIPTION: !> Module responsible for loading, advancing and computing the necessary diagnostics for the simulated particles. !------------------------------------------------------------------------------ MODULE beam ! USE constants USE mpi USE mpihelper USE basic, ONLY: mpirank, mpisize USE distrib IMPLICIT NONE !> Stores the particles properties for the run. TYPE particles INTEGER :: Nploc !< Local number of simulated particles INTEGER :: Nptot !< Total number of simulated particles REAL(kind=db) :: m !< Particle mass REAL(kind=db) :: q !< Particle charge REAL(kind=db) :: weight !< Number of particles represented by one macro-particle REAL(kind=db) :: qmRatio !< Charge over mass ratio REAL(kind=db) :: H0 REAL(kind=db) :: P0 REAL(kind=db) :: temperature LOGICAL :: Davidson=.false. LOGICAL :: is_test= .false. INTEGER, DIMENSION(4) :: nblost !< number of particles lost in (z_min,z_max,r_min,r_max) since last gather INTEGER :: nbadded !< number of particles added by source since last gather INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex !< Index in the electric potential grid for the R direction INTEGER, DIMENSION(:), ALLOCATABLE :: Zindex !< Index in the electric potential grid for the Z direction INTEGER, DIMENSION(:), ALLOCATABLE :: partindex !< Index of the particle to be able to follow it when it goes from one MPI host to the other REAL(kind=db), DIMENSION(:), ALLOCATABLE :: R !< radial coordinates of the particles REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Z !< longitudinal coordinates of the particles REAL(kind=db), DIMENSION(:), ALLOCATABLE :: THET !< azimuthal coordinates of the particles REAL(kind=db), DIMENSION(:), ALLOCATABLE :: BZ !< axial radial relative distances to the left grid line REAL(kind=db), DIMENSION(:), ALLOCATABLE :: BR !< radial relative distances to the bottom grid line REAL(kind=db), DIMENSION(:), ALLOCATABLE :: pot !< Electric potential REAL(kind=db), DIMENSION(:), ALLOCATABLE :: potxt !< External electric potential REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Er !< Radial Electric field REAL(kind=db), DIMENSION(:), ALLOCATABLE :: Ez !< Axial electric field REAL(kind=db), DIMENSION(:),POINTER:: UR !< normalized radial velocity at the current time step REAL(kind=db), DIMENSION(:),POINTER:: URold !< normalized radial velocity at the previous time step REAL(kind=db), DIMENSION(:),POINTER:: UTHET !< normalized azimuthal velocity at the current time step REAL(kind=db), DIMENSION(:),POINTER:: UTHETold !< normalized azimuthal velocity at the previous time step REAL(kind=db), DIMENSION(:),POINTER:: UZ !< normalized axial velocity at the current time step REAL(kind=db), DIMENSION(:),POINTER:: UZold !< normalized axial velocity at the previous time step REAL(kind=db), DIMENSION(:),POINTER:: Gamma !< Lorentz factor at the current time step REAL(kind=db), DIMENSION(:),POINTER:: Gammaold !< Lorentz factor at the previous time step Real(kind=db), Dimension(:,:),ALLOCATABLE:: geomweight !< geometric weight at the particle position LOGICAL:: collected !< Stores if the particles data have been collected to MPI root process during this timestep INTEGER, DIMENSION(:), ALLOCATABLE:: addedlist END TYPE particles ! !TYPE(particles) :: parts !< Storage for all the particles !SAVE :: parts TYPE(particles), DIMENSION(:), ALLOCATABLE, SAVE :: partslist ! Diagnostics (scalars) REAL(kind=db) :: ekin=0 !< Total kinetic energy (J) REAL(kind=db) :: epot=0 !< Total potential energy (J) REAL(kind=db) :: etot=0 !< Current total energy (J) REAL(kind=db) :: etot0=0 !< Initial total energy (J) REAL(kind=db) :: loc_etot0=0 !< theoretical local total energy (J) REAL(kind=db) :: Energies(4) !< (1) kinetic energy, (2) potential energy, (3) total energy and (4) gained/lossed energy due to gain or loss of particles (J) ! INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: Nplocs_all !< Array containing the local numbers of particles in each MPI process ! CONTAINS !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Allocate the memory for the particles variable storing the particles quantities. ! !> @param[inout] p the particles variable needing to be allocated. !> @param[in] nparts the maximum number of particles that will be stored in this variable !--------------------------------------------------------------------------- SUBROUTINE creat_parts(p, nparts) TYPE(particles) :: p INTEGER, INTENT(in) :: nparts IF (.NOT. ALLOCATED(p%Z) ) THEN p%Nploc = nparts p%Nptot = nparts ALLOCATE(p%Z(nparts)) ALLOCATE(p%R(nparts)) ALLOCATE(p%THET(nparts)) ALLOCATE(p%BZ(nparts)) ALLOCATE(p%BR(nparts)) ALLOCATE(p%UR(nparts)) ALLOCATE(p%UZ(nparts)) ALLOCATE(p%UTHET(nparts)) ALLOCATE(p%URold(nparts)) ALLOCATE(p%UZold(nparts)) ALLOCATE(p%UTHETold(nparts)) ALLOCATE(p%Gamma(nparts)) ALLOCATE(p%Rindex(nparts)) ALLOCATE(p%Zindex(nparts)) ALLOCATE(p%partindex(nparts)) ALLOCATE(p%pot(nparts)) ALLOCATE(p%potxt(nparts)) ALLOCATE(p%Er(nparts)) ALLOCATE(p%Ez(nparts)) ALLOCATE(p%GAMMAold(nparts)) Allocate(p%geomweight(nparts,0:2)) p%nblost=0 p%nbadded=0 p%partindex=-1 p%URold=0 p%UZold=0 p%UTHETold=0 p%rindex=0 p%zindex=0 p%BR=0 p%BZ=0 p%UR=0 p%UZ=0 p%UTHET=0 p%Z=0 p%R=0 p%THET=0 p%Gamma=1 p%Er=0 p%Ez=0 p%pot=0 p%potxt=0 p%gammaold=1 p%collected=.false. p%Davidson=.false. p%is_test=.false. p%m=me p%q=-elchar p%qmRatio=p%q/p%m p%weight=1.0_db p%H0=0 p%P0=0 p%temperature=0 p%geomweight=0 END IF END SUBROUTINE creat_parts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Loads the particles at the beginning of the simulation and create the parts variable if necessary !--------------------------------------------------------------------------- SUBROUTINE load_parts USE basic, ONLY: nplasma, mpirank, ierr, distribtype, mpisize, nlclassical, nbspecies, Zbounds USE mpi INTEGER:: i, j REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ALLOCATE(VZ(nplasma), VR(nplasma), VTHET(nplasma)) ! Select case to define the type of distribution SELECT CASE(distribtype) CASE(1) ! Gaussian distribution in V, uniform in Z and 1/R in R CALL loaduniformRZ(partslist(1), VR, VZ, VTHET) CASE(2) !Stable distribution from Davidson 4.95 p.119 CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir) CASE(3) !Stable distribution from Davidson 4.95 p.119 but with constant distribution in R CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodinvr) CASE(4) !Stable distribution from Davidson 4.95 p.119 but with gaussian distribution in R CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodgausr) CASE(5) !Stable distribution from Davidson 4.95 p.119 with gaussian in V computed from v_th given by temp CALL loadDavidson(partslist(1), VR, VZ, VTHET, lodunir) CASE(6) ! Uniform distribution in R and Z and Gaussian distribution in V with Vz @brief Checks for each particle if the z position is outside of the local/global simulation space. !> Depending on the boundary conditions, the leaving particles are sent to the correct neighbouring MPI process !> or deleted. ! !> @param[in] p particles structure ! !> @author Guillaume Le Bars EPFL/SPC !--------------------------------------------------------------------------- SUBROUTINE bound(p) USE basic, ONLY: zgrid, nz, Zbounds, mpirank, step, leftproc, rightproc USE IFPORT IMPLICIT NONE type(particles), INTENT(INOUT):: p INTEGER :: i, rsendnbparts, lsendnbparts, nblostparts INTEGER :: receivednbparts, partdiff, lostindex, sentindex INTEGER, DIMENSION(p%Nploc) :: sendhole INTEGER, DIMENSION(p%Nploc) :: losthole LOGICAL:: leftcomm, rightcomm, sent INTEGER, ALLOCATABLE:: partstoremove(:) rsendnbparts=0 lsendnbparts=0 nblostparts=0 losthole=0 sendhole=0 receivednbparts=0 ! We communicate with the left processus leftcomm = leftproc .ne. -1 ! We communicate with the right processus rightcomm = rightproc .ne. -1 IF (p%Nploc .gt. 0) THEN ! Boundary condition at z direction !$OMP PARALLEL DO DEFAULT(SHARED) DO i=1,p%Nploc ! If the particle is to the right of the local simulation space, it is sent to the right MPI process IF (p%Z(i) .ge. zgrid(Zbounds(mpirank+1))) THEN !$OMP CRITICAL (nbparts) IF(rightcomm) THEN rsendnbparts=rsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=i DO WHILE (p%Z(i) .GT. zgrid(nz)) p%Z(i) = p%Z(i) - zgrid(nz) + zgrid(0) END DO ELSE nblostparts=nblostparts+1 losthole(nblostparts)=i p%nblost(2)=p%nblost(2)+1 END IF !$OMP END CRITICAL (nbparts) ! If the particle is to the left of the local simulation space, it is sent to the left MPI process ELSE IF (p%Z(i) .lt. zgrid(Zbounds(mpirank))) THEN !$OMP CRITICAL (nbparts) IF(leftcomm) THEN ! We send the particle to the left process lsendnbparts=lsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=-i DO WHILE (p%Z(i) .LT. zgrid(0)) p%Z(i) = p%Z(i) + zgrid(nz) - zgrid(0) END DO ELSE ! we destroy the particle nblostparts=nblostparts+1 losthole(nblostparts)=i p%nblost(1)=p%nblost(1)+1 END IF !$OMP END CRITICAL (nbparts) END IF END DO !$OMP END PARALLEL DO END IF IF(mpisize .gt. 1) THEN ! We send the particles leaving the local simulation space to the closest neighbour CALL particlescommunication(p, lsendnbparts, rsendnbparts, sendhole, receivednbparts, (/leftproc,rightproc/)) END IF ! If the boundary conditions are not periodic, we delete the corresponding particles IF(nblostparts .gt. 0 .and. step .ne. 0) THEN DO i=nblostparts,1,-1 CALL delete_part(p, losthole(i), .false. ) END DO !WRITE(*,'(i8.2,a,i4.2)') nblostparts, " particles lost in z on process: ", mpirank END IF ! computes if we received less particles than we sent partdiff=max(lsendnbparts+rsendnbparts-receivednbparts,0) IF(nblostparts + partdiff .gt. 0) THEN ALLOCATE(partstoremove(nblostparts+partdiff)) partstoremove(1:partdiff)=abs(sendhole(receivednbparts+1:receivednbparts+partdiff)) partstoremove(partdiff+1:partdiff+nblostparts)=abs(losthole(1:nblostparts)) call qsort(partstoremove,size(partstoremove),sizeof(partstoremove(1)),compare_int) ! If we received less particles than we sent, or lost particles we fill the remaining holes with the particles from the end of the ! parts arrays DO i=1,nblostparts+partdiff CALL move_part(p, p%Nploc, partstoremove(i)) p%partindex(p%Nploc)=-1 p%Nploc = p%Nploc-1 END DO END IF CONTAINS INTEGER(2) function compare_int(a,b) INTEGER :: a, b, c c=b-a if(c.ne.0) c=1 compare_int=sign(c,b-a) end function END subroutine bound !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Compute the grid cell indices for each particle as well as the distance weight Wr, Wz. !> @param[in] p particles structure !--------------------------------------------------------------------------- SUBROUTINE localisation(p) USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, mpirank Use geometry, ONLY: r_a, geom_weight USE IFPORT type(particles), INTENT(INOUT):: p INTEGER :: i, nblostparts INTEGER, DIMENSION(p%Nploc) :: losthole INTEGER, DIMENSION(2):: nblost losthole=0 nblostparts=0 nblost=0 IF (p%Nploc .gt. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED), private(i) reduction(+:nblost) DO i=1,p%Nploc call geom_weight(p%Z(i), p%r(i), p%geomweight(i,:)) if(p%geomweight(i,0).le.0 .or. p%R(i) .gt. rgrid(nr) .or. p%R(i) .lt. rgrid(0)) then ! If the particle is outside of the simulation space in the r direction, it is deleted. !$OMP CRITICAL (lostparts) nblostparts=nblostparts+1 losthole(nblostparts)=i !$OMP END CRITICAL (lostparts) if(p%R(i) .le. r_a) then nblost(1)=nblost(1)+1 else nblost(2)=nblost(2)+1 end if else call p_calc_rzindex(p,i) end if END DO !$OMP END PARALLEL DO IF(nblostparts.gt.0) THEN p%nblost(3:4)=nblost+p%nblost(3:4) call qsort(losthole(1:nblostparts),nblostparts,sizeof(losthole(1)),compare_int) DO i=1,nblostparts !WRITE(*,'(a,i8.8,a,i4.2,a,3(1pe12.3))') "Particle: ", losthole(i) , " of process: ", mpirank, " is out of bound in r: ", p%R(losthole(i))*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm CALL delete_part(p,losthole(i)) END DO !WRITE(*,'(i6.2,a,i6.2)') nblostparts, " particles lost in r on process: ", mpirank END IF END IF CONTAINS INTEGER(2) function compare_int(a,b) INTEGER :: a, b, c c=b-a if(c.ne.0) c=1 compare_int=sign(c,b-a) end function END SUBROUTINE localisation subroutine p_calc_rzindex(p,i) use basic, only: rgrid,zgrid,invdz,invdr, nnr, nr integer::i type(particles)::p IF (p%R(i) .GT. rgrid(0) .AND. p%R(i) .LT. rgrid(nnr(1))) THEN p%rindex(i)=(p%R(i)-rgrid(0))*invdr(1) ELSE IF(p%R(i) .GE. rgrid(nnr(1)) .AND. p%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN p%rindex(i)=(p%R(i)-rgrid(nnr(1)))*invdr(2)+nnr(1) ELSE IF(p%R(i) .GE. rgrid(nnr(1)+nnr(2)) .AND. p%R(i) .LT. rgrid(nr)) THEN p%rindex(i)=(p%R(i)-rgrid(nnr(1)+nnr(2)))*invdr(3)+nnr(1)+nnr(2) End if p%zindex(i)=(p%Z(i)-zgrid(0))*invdz end subroutine p_calc_rzindex SUBROUTINE comp_mag_p(p) USE basic, ONLY: zgrid, rgrid, dz, BZ, BR, nz, invdz type(particles), INTENT(INOUT):: p INTEGER :: i, nblostparts Real(kind=db):: WZ,WR INTEGER:: j1,j2,j3,j4 !$OMP PARALLEL DO SIMD DEFAULT(SHARED) Private(J1,J2,J3,J4,WZ,WR) DO i=1,p%Nploc WZ=(p%Z(i)-zgrid(p%zindex(i)))*invdz; WR=(p%R(i)-rgrid(p%rindex(i)))/(rgrid(p%rindex(i)+1)-rgrid(p%rindex(i))); J1=(p%rindex(i))*(nz+1) + p%zindex(i)+1 J2=(p%rindex(i))*(nz+1) + p%zindex(i)+2 J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1 J4=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+2 ! Interpolation for magnetic field p%BZ(i)=(1-WZ)*(1-WR)*Bz(J4) & & +WZ*(1-WR)*Bz(J3) & & +(1-WZ)*WR*Bz(J2) & & +WZ*WR*Bz(J1) p%BR(i)=(1-WZ)*(1-WR)*Br(J4) & & +WZ*(1-WR)*Br(J3) & & +(1-WZ)*WR*Br(J2) & & +WZ*WR*Br(J1) END DO !$OMP END PARALLEL DO SIMD end subroutine comp_mag_p !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief General routine to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using a pointer to the routine computing gamma. This avoid the nlclassical flag check on each particle. ! !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE comp_velocity(p) ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : nlclassical type(particles), INTENT(INOUT):: p ! Store old Velocities CALL swappointer(p%UZold, p%UZ) CALL swappointer(p%URold, p%UR) CALL swappointer(p%UTHETold, p%UTHET) CALL swappointer(p%Gammaold, p%Gamma) IF (nlclassical) THEN CALL comp_velocity_fun(p, gamma_classical) ELSE CALL comp_velocity_fun(p, gamma_relativistic) END IF END SUBROUTINE comp_velocity !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine called by comp_velocity to compute the velocities at time t+1. !> This routine allows to treat the classical and relativistic case efficiently from a numerical standpoint, !> by using the routine computing gamma as an input. This avoid the nlclassical flag check on each particle. ! !> @param[in] gamma the function used to compute the value of the lorentz factor \f$\gamma\f$ !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE comp_velocity_fun(p, gamma) ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : bnorm, dt, BZ, BR, nz interface REAL(kind=db) FUNCTION gamma(UZ, UR, UTHET) USE constants REAL(kind=db), INTENT(IN):: UR,UZ,UTHET end FUNCTION end interface type(particles), INTENT(INOUT):: p REAL(kind=db) :: tau REAL(kind=db):: BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2 INTEGER:: J1, J2, J3, J4 INTEGER:: i ! Normalized time increment !tau=omegac/2/omegap*dt/tnorm tau=p%qmRatio*bnorm*0.5*dt IF (p%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2) DO i=1,p%Nploc !J1=(p%rindex(i))*(nz+1) + p%zindex(i)+1 !J2=(p%rindex(i))*(nz+1) + p%zindex(i)+2 !J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1 !J4=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+2 ! Interpolation for magnetic field !BRZ=(1-p%WZ(i))*(1-p%WR(i))*Bz(J4) & !& +p%WZ(i)*(1-p%WR(i))*Bz(J3) & !& +(1-p%WZ(i))*p%WR(i)*Bz(J2) & !& +p%WZ(i)*p%WR(i)*Bz(J1) !BRR=(1-p%WZ(i))*(1-p%WR(i))*Br(J4) & !& +p%WZ(i)*(1-p%WR(i))*Br(J3) & !& +(1-p%WZ(i))*p%WR(i)*Br(J2) & !& +p%WZ(i)*p%WR(i)*Br(J1) ! First half of electric pulse p%UZ(i)=p%UZold(i)+p%Ez(i)*tau p%UR(i)=p%URold(i)+p%ER(i)*tau p%Gamma(i)=gamma(p%UZ(i), p%UR(i), p%UTHETold(i)) ! Rotation along magnetic field !ZBZ=tau*BRZ/p%Gamma(i) !ZBR=tau*BRR/p%Gamma(i) ZBZ=tau*p%BZ(i)/p%Gamma(i) ZBR=tau*p%BR(i)/p%Gamma(i) ZPZ=p%UZ(i)-ZBR*p%UTHETold(i) !u'_{z} ZPR=p%UR(i)+ZBZ*p%UTHETold(i) !u'_{r} ZPTHET=p%UTHETold(i)+(ZBR*p%UZ(i)-ZBZ*p%UR(i)) !u'_{theta} SQR=1+ZBZ*ZBZ+ZBR*ZBR ZBZ2=2*ZBZ/SQR ZBR2=2*ZBR/SQR p%UZ(i)=p%UZ(i)-ZBR2*ZPTHET !u+_{z} p%UR(i)=p%UR(i)+ZBZ2*ZPTHET !u+_{r} p%UTHET(i)=p%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta} ! Second half of acceleration p%UZ(i)=p%UZ(i)+p%EZ(i)*tau p%UR(i)=p%UR(i)+p%ER(i)*tau ! Final computation of the Lorentz factor p%Gamma(i)=gamma(p%UZ(i), p%UR(i), p%UTHET(i)) END DO !$OMP END PARALLEL DO SIMD END IF p%collected=.false. END SUBROUTINE comp_velocity_fun !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the classical simulations. !> This routine systematically returns 1.0 to treat the system according to classical dynamic. ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- ELEMENTAL REAL(kind=db) FUNCTION gamma_classical(UZ, UR, UTHET) #if __INTEL_COMPILER > 1700 !$OMP declare simd(gamma_classical) #endif REAL(kind=db), INTENT(IN):: UR,UZ,UTHET gamma_classical=1.0 END FUNCTION gamma_classical !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Routine used to compute the lorentz factor \f$\gamma\f$ in the relativistic simulations. !> This routine computes the Lorentz factor \f$\gamma=\sqrt{1+\mathbf{\gamma\beta}^2}\f$ ! !> @param[out] gamma the lorentz factor \f$\gamma\f$ !> @param[in] UZ \f$\gamma\beta_z=\gamma v_z/c\f$ the normalized particle longitudinal velocity !> @param[in] UR \f$\gamma\beta_r=\gamma v_r/c\f$ the normalized particle radial velocity !> @param[in] UTHET \f$\gamma\beta_\theta=\gamma v_\theta/c\f$ the normalized particle azimuthal velocity !--------------------------------------------------------------------------- ELEMENTAL REAL(kind=db) FUNCTION gamma_relativistic(UZ, UR, UTHET) #if __INTEL_COMPILER > 1700 !$OMP declare simd(gamma_relativistic) #endif REAL(kind=db), INTENT(IN):: UR,UZ,UTHET gamma_relativistic=sqrt(1+UZ**2+UR**2+UTHET**2) END FUNCTION gamma_relativistic !--------------------------------------------------------------------------- !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes the particles position at time t+1 !> This routine computes the particles position at time t+1 according to the Bunemann algorithm. ! !> @param[in] p The particles structure being updated !--------------------------------------------------------------------------- SUBROUTINE push(p) Use basic, ONLY: dt, tnorm type(particles), INTENT(INOUT):: p REAL(kind=db):: XP, YP, COSA, SINA, U1, U2, ALPHA INTEGER :: i IF (p%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(XP, YP, COSA, SINA, U1, U2, ALPHA) DO i=1,p%Nploc ! Local Cartesian coordinates XP=p%R(i)+dt/tnorm*p%UR(i)/p%Gamma(i) YP=dt/tnorm*p%UTHET(i)/p%Gamma(i) ! Conversion to cylindrical coordiantes p%Z(i)=p%Z(i)+dt/tnorm*p%UZ(i)/p%Gamma(i) p%R(i)=sqrt(XP**2+YP**2) ! Computation of the rotation angle IF (p%R(i) .EQ. 0) THEN COSA=1 SINA=0 ALPHA=0 ELSE COSA=XP/p%R(i) SINA=YP/p%R(i) ALPHA=asin(SINA) END IF ! New azimuthal position p%THET(i)=MOD(p%THET(i)+ALPHA,2*pi) ! Velocity in rotated reference frame U1=COSA*p%UR(i)+SINA*p%UTHET(i) U2=-SINA*p%UR(i)+COSA*p%UTHET(i) p%UR(i)=U1 p%UTHET(i)=U2 END DO !$OMP END PARALLEL DO SIMD END IF p%collected=.false. END SUBROUTINE push !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Computes several diagnostic quantities !> This routine computes the total kinetic and electric potential energy. !> It keeps track of the reference energy and the number of particle per mpi node. ! !--------------------------------------------------------------------------- SUBROUTINE diagnostics ! ! Compute energies ! USE constants, ONLY: vlight USE basic, ONLY: phinorm, cstep, nlclassical, ierr, step, nlend,& & itparts INTEGER:: i ! Reset the quantities ekin=0 epot=0 etot=0 ! Computation of the kinetic and potential energy as well as fluid velocities and density !$OMP PARALLEL DO SIMD REDUCTION(+:epot, ekin) DEFAULT(SHARED) DO i=1,partslist(1)%Nploc ! Potential energy epot=epot+(partslist(1)%pot(i)+partslist(1)%potxt(i)) ! Kinetic energy IF(.not. nlclassical) THEN ekin=ekin+(partslist(1)%Gamma(i)-1) ELSE ekin=ekin+0.5*( partslist(1)%UR(i)**2 & & + partslist(1)%UZ(i)**2 & & + partslist(1)%UTHET(i)**2 ) END IF END DO !$OMP END PARALLEL DO SIMD epot=epot*phinorm*0.5*partslist(1)%q*partslist(1)%weight ekin=ekin*partslist(1)%m*partslist(1)%weight*vlight**2 ! Shift to Etot at cstep=1 (not valable yet at cstep=0!) IF(cstep.EQ. 0) THEN ! Compute the local total energy loc_etot0 = epot+ekin etot0=0 END IF !etot=loc_etot0 ! Compute the total energy etot=epot+ekin Energies=(/ekin,epot,etot,loc_etot0/) ! The computed energy is sent to the root process IF(mpisize .gt.1) THEN IF(mpirank .eq.0 ) THEN CALL MPI_REDUCE(MPI_IN_PLACE, Energies, 4, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ierr) etot0=etot0+Energies(4) ekin=Energies(1) epot=Energies(2) etot=Energies(3) ELSE CALL MPI_REDUCE(Energies, Energies, 4, db_type, db_sum_op, & & 0, MPI_COMM_WORLD, ierr) END IF ELSE etot0=etot0+loc_etot0 END IF loc_etot0=0 ! Send the local number of particles on each node to the root process IF(mpisize .gt. 1) THEN Nplocs_all(mpirank)=partslist(1)%Nploc IF(mpirank .eq.0 ) THEN CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) partslist(1)%Nptot=sum(Nplocs_all) ELSE CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) END IF ELSE partslist(1)%Nptot=partslist(1)%Nploc END IF !Calls the communications to send the particles data to the root process for diagnostic file every itparts time steps IF(mpisize .gt. 1 .and. (modulo(step,itparts) .eq. 0 .or. nlend)) THEN CALL collectparts(partslist(1)) END IF end subroutine diagnostics !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Collect the particles positions and velocities on the root process. !> If the collection has already been performed at this time step, the routine does nothing. ! !--------------------------------------------------------------------------- SUBROUTINE collectparts(p) USE basic, ONLY: mpirank, mpisize, ierr type(particles), INTENT(INOUT):: p INTEGER, DIMENSION(:), ALLOCATABLE :: displs, Nploc INTEGER:: i IF(p%collected) RETURN ! exit subroutine if particles have already been collected during this time step ALLOCATE(Nploc(0:mpisize-1)) ALLOCATE(displs(0:mpisize-1)) displs=0 Nploc(mpirank)=p%Nploc CALL MPI_Allgather(MPI_IN_PLACE, 1, MPI_INTEGER, Nploc, 1, MPI_INTEGER,& & MPI_COMM_WORLD, ierr) p%Nptot=sum(Nploc) IF(p%Nptot .eq. 0 ) THEN p%collected=.true. RETURN END IF IF(mpirank .ne. 0) THEN ! Send Particles informations to root process CALL MPI_Gatherv(p%Z, Nploc(mpirank), db_type, p%Z, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%R, Nploc(mpirank), db_type, p%R, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%THET, Nploc(mpirank), db_type, p%THET, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%UR, Nploc(mpirank), db_type, p%UR, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%UZ, Nploc(mpirank), db_type, p%UZ, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%UTHET, Nploc(mpirank), db_type, p%UTHET, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%pot, Nploc(mpirank), db_type, p%pot, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%Rindex, Nploc(mpirank), MPI_INTEGER, p%Rindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%Zindex, Nploc(mpirank), MPI_INTEGER, p%Zindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(p%partindex, Nploc(mpirank), MPI_INTEGER, p%partindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) ELSE Do i=1,mpisize-1 displs(i)=displs(i-1)+Nploc(i-1) END DO IF(p%Nptot .gt. size(p%R,1)) THEN CALL change_parts_allocation(p,p%Nptot-size(P%R,1)) END IF ! Receive particle information from all processes CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%Z, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%R, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%THET, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%UR, Nploc, displs,& & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%UZ, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%UTHET, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), db_type, p%pot, Nploc, displs, & & db_type, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), MPI_INTEGER, p%Rindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), MPI_INTEGER, p%Zindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) CALL MPI_Gatherv(MPI_IN_PLACE, Nploc(mpirank), MPI_INTEGER, p%partindex, Nploc, displs, & & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) p%partindex(sum(Nploc)+1:)=-1 END IF p%collected=.TRUE. END SUBROUTINE collectparts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Computes the velocities at time t-1/2 to keep the second order precision in time on the velocity. ! !--------------------------------------------------------------------------- SUBROUTINE adapt_vinit(p) !! Computes the velocity at time -dt/2 from velocities computed at time 0 ! USE basic, ONLY : bnorm, dt, BZ, BR, nlclassical, phinorm, nz, distribtype, vnorm type(particles), INTENT(INOUT):: p REAL(kind=db) :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & & SQR, Vperp, v2 INTEGER :: J1, J2, J3, J4, i REAL(kind=db), DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET ! In case Davidson distribution is used the longitudinal and radial velocities are adapted to take into account the ! electric potential. IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3 .OR. distribtype .EQ. 4 .or. p%Davidson) THEN ALLOCATE(VR(p%Nploc),VZ(p%Nploc),VTHET(p%Nploc)) CALL loduni(7,VZ) VZ=VZ*2*pi VTHET=p%UTHET/p%Gamma*vnorm DO i=1,p%Nploc Vperp=sqrt(MAX(2*p%H0/p%m-2*p%qmRatio*p%pot(i)*phinorm-VTHET(i)**2,0.0_db)) VR(i)=Vperp*sin(VZ(i)) VZ(i)=Vperp*cos(VZ(i)) IF(nlclassical) THEN p%Gamma(i)=1 ELSE v2=VR(i)**2+VZ(i)**2+VTHET(i)**2 p%Gamma(i)=sqrt(1/(1-v2/vnorm**2)) END IF p%UR(i)=p%Gamma(i)*VR(i)/vnorm p%UZ(i)=p%Gamma(i)*VZ(i)/vnorm p%UTHET(i)=p%Gamma(i)*VTHET(i)/vnorm END DO DEALLOCATE(VR,VZ,VTHET) END IF RETURN ! Normalised time increment !tau=-omegac/2/omegap*dt/tnorm tau=p%qmRatio*bnorm*0.5*dt ! Store old Velocities CALL swappointer(p%UZold, p%UZ) CALL swappointer(p%URold, p%UR) CALL swappointer(p%UTHETold, p%UTHET) CALL swappointer(p%Gammaold, p%Gamma) IF (p%Nploc .NE. 0) THEN !$OMP PARALLEL DO SIMD DEFAULT(SHARED) PRIVATE(J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR) DO i=1,p%Nploc ! Compute the particle linear indices for the magnetic field interpolation !J1=(p%rindex(i))*(nz+1)+p%zindex(i)+1 !J2=J1+1 !J3=(p%rindex(i)+1)*(nz+1)+p%zindex(i)+1 !J4=J3+1 ! ! Interpolation for magnetic field !BRZ=(1-p%BZ(i))*(1-p%BR(i))*Bz(J4) & !& +p%BZ(i)*(1-p%BR(i))*Bz(J3) & !& +(1-p%BZ(i))*p%BR(i)*Bz(J2) & !& +p%BZ(i)*p%BR(i)*Bz(J1) !BRR=(1-p%BZ(i))*(1-p%BR(i))*Br(J4) & !& +p%BZ(i)*(1-p%BR(i))*Br(J3) & !& +(1-p%BZ(i))*p%BR(i)*Br(J2) & !& +p%BZ(i)*p%BR(i)*Br(J1) ! Half inverse Rotation along magnetic field ZBZ=tau*p%BZ(i)/p%Gammaold(i) ZBR=tau*p%BR(i)/p%Gammaold(i) SQR=1+ZBZ*ZBZ+ZBR*ZBR ZPZ=(p%UZold(i)-ZBR*p%UTHETold(i))/SQR !u-_{z} ZPR=(p%URold(i)+ZBZ*p%UTHETold(i))/SQR !u-_{r} ZPTHET=p%UTHETold(i)+(ZBR*p%UZold(i)-ZBZ*p%URold(i))/SQR !u-_{theta} p%UZ(i)=ZPZ p%UR(i)=ZPR p%UTHET(i)=ZPTHET ! half of decceleration p%UZ(i)=p%UZ(i)+p%Ez(i)*tau p%UR(i)=p%UR(i)+p%Er(i)*tau IF(.not. nlclassical) THEN p%Gamma(i)=sqrt(1+p%UZ(i)**2+p%UR(i)**2+p%UTHET(i)**2) END IF END DO !$OMP END PARALLEL DO SIMD END IF END SUBROUTINE adapt_vinit !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief In the case of MPI parallelism, computes the indices of the particle assigned to the current process. !--------------------------------------------------------------------------- SUBROUTINE distribpartsonprocs(p) ! Computes the start and end indices for the Z boundaries on local processus ! Computes the particle indices from initial particle loading vector, that stay in current process USE basic, ONLY: nz, Zbounds TYPE(particles), INTENT(INOUT):: p INTEGER:: k, i, nbparts REAL(kind=db):: idealnbpartsperproc - INTEGER:: totparts ! Total number of particle from zgrid(0) to zgrid(k) - INTEGER:: partindexstart ! Starting index of local particle in the parts variable used later to move the local particles at the begining of the vector - INTEGER:: partindexlast ! Ending index of local particle in the parts variable used later to move the local particles at the begining of the vector INTEGER, DIMENSION(0:nz):: partspercol ! Vector containing the number of particles between zgrid(n) and zgrid(n+1) INTEGER:: Zmin, Zmax ! Minimum and maximum indices of particles in Z direction INTEGER:: Zperproc, zindex partspercol=0 idealnbpartsperproc = FLOOR(REAL(p%Nploc)/REAL(mpisize)) Zmin=MINVAL(p%Zindex(1:p%Nploc)) Zmax=MAXVAL(p%Zindex(1:p%Nploc)) Zperproc=(Zmax-Zmin)/mpisize IF(Zmax .eq. 0) Zmax=nz DO k=1,p%Nploc zindex=p%Zindex(k) partspercol(zindex)=partspercol(zindex)+1 END DO IF (Zperproc .eq. 0) THEN !! No particles are present initially Zperproc=nz/mpisize Zmin=0 DO k=1,mpisize-1 IF(k .lt. mpisize-1-MODULO(Zmax-Zmin,mpisize)) THEN Zbounds(k)=Zmin+k*Zperproc-1 ELSE Zbounds(k)=Zmin+k*Zperproc-1+k-mpisize+2+MODULO(Zmax-Zmin,mpisize) END IF END DO ELSE i=0 DO k=1,mpisize-1 nbparts=0 DO WHILE(nbparts<0.99*idealnbpartsperproc .and. i .lt. Zmax) nbparts=nbparts+partspercol(i) i=i+1 END DO Zbounds(k)=i END DO END IF - partindexstart=1 - totparts=0 - partindexlast=p%Nploc - DO k=0,mpisize-1 Nplocs_all(k)=SUM(partspercol(Zbounds(k):Zbounds(k+1)-1)) END DO WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " nptot", Nplocs_all(mpirank) END SUBROUTINE distribpartsonprocs SUBROUTINE keep_mpi_self_parts(p,Zbounds) TYPE(particles),INTENT(INOUT):: p INTEGER,INTENT(in)::Zbounds(0:) INTEGER :: i, partstart, old_sum,ierr partstart=1 p%Nploc=0 Do i=1,p%Nptot IF(p%Zindex(i).ge.Zbounds(mpirank).and.p%Zindex(i).lt.Zbounds(mpirank+1)) THEN p%Nploc=p%Nploc+1 CALL move_part(p,i,p%Nploc) END IF END DO old_sum=p%Nptot CALL MPI_REDUCE(p%Nploc, p%Nptot,1,MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr) IF(p%Nptot .ne. old_sum) WRITE(*,*) "Error in particle distribution kept: ", p%Nptot, "/",old_sum END SUBROUTINE keep_mpi_self_parts !_______________________________________________________________________________ !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Manage the particle communication between neighbours. !> This routine is responsible to receive the incoming particles from the MPI neighbours and to send its outgoing !> particles to these neighbours ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative !--------------------------------------------------------------------------- SUBROUTINE particlescommunication(p, lsendnbparts, rsendnbparts, sendholes, receivednbparts, procs) USE mpihelper, ONLY: particle_type #ifdef _DEBUG USE basic, ONLY: step #endif type(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(out) :: receivednbparts INTEGER, INTENT(in) :: sendholes(:) INTEGER, INTENT(in) :: procs(2) INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0 INTEGER, ASYNCHRONOUS :: sendrequest(2), recvrequest(2) INTEGER, ASYNCHRONOUS :: sendstatus(MPI_STATUS_SIZE,2), recvstatus(MPI_STATUS_SIZE,2) TYPE(particle), ALLOCATABLE :: rrecvpartbuff(:), lrecvpartbuff(:), rsendpartbuff(:), lsendpartbuff(:) ! buffers to send and receive particle from left and right processes INTEGER :: lsentnbparts, rsentnbparts INTEGER :: lreceivednbparts, rreceivednbparts, ierr lsentnbparts=lsendnbparts rsentnbparts=rsendnbparts sendrequest=MPI_REQUEST_NULL recvrequest=MPI_REQUEST_NULL lrecvnbparts=0 rrecvnbparts=0 ! Send and receive the number of particles to exchange CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, procs(1), 0, MPI_COMM_WORLD, recvrequest(1), ierr) CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, procs(2), 0, MPI_COMM_WORLD, recvrequest(2), ierr) CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, procs(1), 0, MPI_COMM_WORLD, sendrequest(1), ierr) CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, procs(2), 0, MPI_COMM_WORLD, sendrequest(2), ierr) CALL MPI_Waitall(2,recvrequest(1:2), recvstatus(:,1:2), ierr) recvrequest=MPI_REQUEST_NULL lreceivednbparts=lrecvnbparts rreceivednbparts=rrecvnbparts ! Re/allocate enough memory to store the incoming particles ALLOCATE(rrecvpartbuff(rreceivednbparts)) ALLOCATE(lrecvpartbuff(lreceivednbparts)) ! Receive particles from left and right processes to the corresponding buffers IF ( lrecvnbparts .gt. 0) THEN CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, procs(1), 1, MPI_COMM_WORLD, recvrequest(1), ierr) END IF IF( rrecvnbparts .gt. 0) THEN CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, procs(2), 1, MPI_COMM_WORLD, recvrequest(2), ierr) END IF ALLOCATE(rsendpartbuff(rsendnbparts)) ALLOCATE(lsendpartbuff(lsendnbparts)) ! Copy the leaving particles to the corresponding send buffers IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN CALL AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes, lsendpartbuff, rsendpartbuff) END IF CALL MPI_Waitall(2,sendrequest(1:2), sendstatus(:,1:2), ierr) ! Send the particles to the left and right neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, procs(1), 1, MPI_COMM_WORLD, sendrequest(1), ierr) #ifdef _DEBUG !WRITE(*,*)"snding ", lsendnbparts , " to left at step: ",step #endif END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, procs(2), 1, MPI_COMM_WORLD, sendrequest(2), ierr) #ifdef _DEBUG !WRITE(*,*)"snding ", rsendnbparts , " to right at step: ",step #endif END IF ! Receive the incoming parts in the receive buffers IF ( lreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(1), recvstatus(:,1), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,1) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF #ifdef _DEBUG !WRITE(*,*)"rcvd ", lreceivednbparts , " from left at step: ",step #endif END IF IF ( rreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(2), recvstatus(:,2), ierr) IF(ierr .ne. MPI_SUCCESS) THEN WRITE(*,*) "Error in particle reception on proc:", mpirank, " error code:", ierr, "status:", recvstatus(:,2) CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF #ifdef _DEBUG !WRITE(*,*)"rcvd ", rreceivednbparts , " from right at step: ",step #endif END IF receivednbparts=rreceivednbparts+lreceivednbparts IF(p%Nploc+receivednbparts-lsendnbparts-rsendnbparts .gt. size(p%R,1)) THEN CALL change_parts_allocation(p,receivednbparts) END IF ! Copy the incoming particles from the receive buffers to the simulation parts variable CALL Addincomingparts(p, rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, & & sendholes, lrecvpartbuff, rrecvpartbuff) ! Wait for the outgoing particles to be fully received by the neighbours IF( lsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(1), sendstatus(:,1), ierr) #ifdef _DEBUG !WRITE(*,*)"sent ", lsentnbparts , " to left at step: ",step #endif END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(2), sendstatus(:,2), ierr) #ifdef _DEBUG !WRITE(*,*)"sent ", rsentnbparts , " to right at step: ",step #endif END IF ! ! END SUBROUTINE particlescommunication !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the receive buffers to the local simulation variable parts. !> The incoming particles will first be stored in the holes left by the outgoing particles, then they !> will be added at the end of the parts variable ! !> @param [in] rrecvnbparts number of particles received from the right neighbour (mpirank+1) !> @param [in] lrecvnbparts number of particles received from the left neighbour (mpirank-1) !> @param [in] sendnbparts total number of particles having left the local domain !> @param [in] sendholes array containing the indices of the particle having left the local domain in ascending order. !--------------------------------------------------------------------------- SUBROUTINE Addincomingparts(p, rrecvnbparts, lrecvnbparts, sendnbparts, sendholes,lrecvpartbuff, rrecvpartbuff) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts INTEGER, INTENT(in) :: sendholes(:) TYPE(particle), INTENT(IN) :: rrecvpartbuff(:), lrecvpartbuff(:) INTEGER k,partpos ! First import the particles coming from the right IF(rrecvnbparts .gt. 0) THEN Do k=1,rrecvnbparts IF(k .le. sendnbparts) THEN ! Fill the holes left by sent parts partpos=abs(sendholes(k)) ELSE ! Add at the end of parts and keep track of number of parts p%Nploc=p%Nploc+1 partpos=p%Nploc END IF CALL Insertincomingpart(p, rrecvpartbuff, partpos, k) END DO END IF ! Then import the particles coming from the left IF(lrecvnbparts .gt. 0) THEN Do k=1,lrecvnbparts IF(k+rrecvnbparts .le. sendnbparts) THEN ! Fill the holes left by sent parts partpos=abs(sendholes(k+rrecvnbparts)) ELSE ! Add at the end of parts and keep track of number of parts p%Nploc=p%Nploc+1 partpos=p%Nploc END IF CALL Insertincomingpart(p, lrecvpartbuff, partpos, k) END DO END IF ! END SUBROUTINE Addincomingparts !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy one particle from the receive buffers to the local simulation variable parts. ! !> @param [in] buffer receive buffer containing the particles parameters to copy from !> @param [in] bufferindex particle index in the receive buffer !> @param [in] partsindex destination particle index in the local parts variable !--------------------------------------------------------------------------- SUBROUTINE Insertincomingpart(p, buffer, partsindex, bufferindex) USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), INTENT(in) :: buffer p%partindex(partsindex) = buffer(bufferindex)%partindex p%Rindex(partsindex) = buffer(bufferindex)%Rindex p%Zindex(partsindex) = buffer(bufferindex)%Zindex p%R(partsindex) = buffer(bufferindex)%R p%Z(partsindex) = buffer(bufferindex)%Z p%THET(partsindex) = buffer(bufferindex)%THET p%UZ(partsindex) = buffer(bufferindex)%UZ p%UR(partsindex) = buffer(bufferindex)%UR p%UTHET(partsindex) = buffer(bufferindex)%UTHET p%Gamma(partsindex) = buffer(bufferindex)%Gamma p%pot(partsindex) = buffer(bufferindex)%pot ! END SUBROUTINE Insertincomingpart !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy one particle from the local parts variable to the send buffer. ! !> @param [in] buffer send buffer to copy to !> @param [in] bufferindex particle index in the send buffer !> @param [in] partsindex origin particle index in the local parts variable !--------------------------------------------------------------------------- SUBROUTINE Insertsentpart(p, buffer, bufferindex, partsindex) USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), INTENT(inout) :: buffer buffer(bufferindex)%partindex = p%partindex(partsindex) buffer(bufferindex)%Rindex = p%Rindex(partsindex) buffer(bufferindex)%Zindex = p%Zindex(partsindex) buffer(bufferindex)%R = p%R(partsindex) buffer(bufferindex)%Z = p%Z(partsindex) buffer(bufferindex)%THET = p%THET(partsindex) buffer(bufferindex)%UZ = p%UZ(partsindex) buffer(bufferindex)%UR = p%UR(partsindex) buffer(bufferindex)%UTHET = p%UTHET(partsindex) buffer(bufferindex)%Gamma = p%Gamma(partsindex) buffer(bufferindex)%pot = p%pot(partsindex) ! END SUBROUTINE Insertsentpart !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy the particles from the local parts variable to the left and right send buffers. ! !> @param [in] lsendnbparts number of particles to send to the left neighbour (mpirank-1) !> @param [in] rsendnbparts number of particles to send to the right neighbour (mpirank+1) !> @param [in] sendholes array containing the indices of the particle leaving the local domain in ascending order. If the index is positive, the particle goes to the right neigbour, and to the left neighbour if the index is negative !--------------------------------------------------------------------------- SUBROUTINE AddPartSendBuffers(p, lsendnbparts, rsendnbparts, sendholes, lsendpartbuff, rsendpartbuff) ! USE mpihelper TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(:) TYPE(particle), INTENT(OUT) :: rsendpartbuff(:), lsendpartbuff(:) INTEGER:: partpos, k INTEGER:: lsendpos, rsendpos lsendpos=0 rsendpos=0 ! Loop over the outgoing particles and fill the correct send buffer Do k=lsendnbparts+rsendnbparts,1,-1 partpos=abs(sendholes(k)) IF(sendholes(k) .GT. 0) THEN rsendpos=rsendpos+1 CALL Insertsentpart(p, rsendpartbuff, rsendpos, partpos) ELSE IF(sendholes(k) .LT. 0) THEN lsendpos=lsendpos+1 CALL Insertsentpart(p, lsendpartbuff, lsendpos, partpos) END IF END DO ! ! END SUBROUTINE AddPartSendBuffers !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> @brief Exchange two particles in the parts variable. ! !> @param [in] index1 index in parts of the first particle to exchange. !> @param [in] index2 index in parts of the second particle to exchange. !--------------------------------------------------------------------------- SUBROUTINE exchange_parts(p, index1, index2) TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(IN) :: index1, index2 REAL(kind=db):: R, Z, THET, UR, UZ, UTHET, Gamma, geomweight(0:2),pot INTEGER :: Rindex, Zindex, partindex !! Exchange particle at index1 with particle at index2 ! Store part at index1 in temporary value partindex = p%partindex(index1) Gamma = p%Gamma(index1) pot = p%pot(index1) R = p%R(index1) Z = p%Z(index1) THET = p%THET(index1) UR = p%UR(index1) UTHET = p%UTHET(index1) UZ = p%UZ(index1) Rindex = p%Rindex(index1) Zindex = p%Zindex(index1) geomweight = p%geomweight(index1,:) ! Move part at index2 in part at index 1 p%partindex(index1) = p%partindex(index2) p%Gamma(index1) = p%Gamma(index2) p%pot(index1) = p%pot(index2) p%R(index1) = p%R(index2) p%Z(index1) = p%Z(index2) p%THET(index1) = p%THET(index2) p%UR(index1) = p%UR(index2) p%UTHET(index1) = p%UTHET(index2) p%UZ(index1) = p%UZ(index2) p%Rindex(index1) = p%Rindex(index2) p%Zindex(index1) = p%Zindex(index2) p%geomweight(index1,:) = p%geomweight(index2,:) ! Move temporary values from part(index1) to part(index2) p%partindex(index2) = partindex p%Gamma(index2) = Gamma p%pot(index2) = pot p%R(index2) = R p%Z(index2) = Z p%THET(index2) = THET p%UR(index2) = UR p%UTHET(index2) = UTHET p%UZ(index2) = UZ p%Rindex(index2) = Rindex p%Zindex(index2) = Zindex p%geomweight(index2,:) = geomweight END SUBROUTINE exchange_parts SUBROUTINE change_parts_allocation(p, sizedifference) implicit none TYPE(particles), INTENT(INOUT):: p INTEGER,INTENT(IN) :: sizedifference CALL change_array_size_int(p%Rindex, sizedifference) CALL change_array_size_int(p%Zindex, sizedifference) CALL change_array_size_int(p%partindex, sizedifference) CALL change_array_size_dp(p%ER,sizedifference) CALL change_array_size_dp(p%EZ,sizedifference) CALL change_array_size_dp(p%pot,sizedifference) CALL change_array_size_dp(p%potxt,sizedifference) CALL change_array_size_dp(p%R,sizedifference) CALL change_array_size_dp(p%Z,sizedifference) CALL change_array_size_dp(p%THET,sizedifference) CALL change_array_size_dp(p%BR,sizedifference) CALL change_array_size_dp(p%BZ,sizedifference) CALL change_array_size_dp2(p%geomweight,sizedifference) CALL change_array_size_dp_ptr(p%UR,sizedifference) CALL change_array_size_dp_ptr(p%URold,sizedifference) CALL change_array_size_dp_ptr(p%UZ,sizedifference) CALL change_array_size_dp_ptr(p%UZold,sizedifference) CALL change_array_size_dp_ptr(p%UTHET,sizedifference) CALL change_array_size_dp_ptr(p%UTHETold,sizedifference) CALL change_array_size_dp_ptr(p%Gamma,sizedifference) CALL change_array_size_dp_ptr(p%Gammaold,sizedifference) p%Nploc=MIN(p%Nploc,size(p%R)) END SUBROUTINE change_parts_allocation SUBROUTINE change_array_size_dp(arr, sizedifference) implicit none REAL(kind=db), ALLOCATABLE, INTENT(INOUT):: arr(:) INTEGER, INTENT(IN):: sizedifference REAL(kind=db), ALLOCATABLE:: temp(:) INTEGER:: current_size, new_size if(allocated(arr)) THEN current_size=size(arr) new_size=current_size+sizedifference ALLOCATE(temp(new_size)) temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size)) DEALLOCATE(arr) CALL move_alloc(temp, arr) END IF END SUBROUTINE change_array_size_dp SUBROUTINE change_array_size_dp2(arr, sizedifference) implicit none REAL(kind=db), ALLOCATABLE, INTENT(INOUT):: arr(:,:) INTEGER, INTENT(IN):: sizedifference REAL(kind=db), ALLOCATABLE:: temp(:,:) INTEGER:: current_size, new_size if(allocated(arr)) THEN current_size=size(arr,1) new_size=current_size+sizedifference ALLOCATE(temp(new_size,0:size(arr,2)-1)) temp(1:min(current_size,new_size),:)=arr(1:min(current_size,new_size),:) DEALLOCATE(arr) CALL move_alloc(temp, arr) END IF END SUBROUTINE change_array_size_dp2 SUBROUTINE change_array_size_dp_ptr(arr, sizedifference) implicit none REAL(kind=db), POINTER, INTENT(INOUT):: arr(:) INTEGER, INTENT(IN):: sizedifference REAL(kind=db), POINTER:: temp(:) INTEGER:: current_size, new_size if(associated(arr)) THEN current_size=size(arr) new_size=current_size+sizedifference ALLOCATE(temp(new_size)) temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size)) DEALLOCATE(arr) arr=> temp END IF END SUBROUTINE change_array_size_dp_ptr SUBROUTINE change_array_size_int(arr, sizedifference) implicit none INTEGER, ALLOCATABLE, INTENT(INOUT):: arr(:) INTEGER, INTENT(IN):: sizedifference INTEGER, ALLOCATABLE:: temp(:) INTEGER:: current_size, new_size if(allocated(arr)) THEN current_size=size(arr) new_size=current_size+sizedifference ALLOCATE(temp(new_size)) temp(1:min(current_size,new_size))=arr(1:min(current_size,new_size)) DEALLOCATE(arr) CALL move_alloc(temp,arr) END IF END SUBROUTINE change_array_size_int SUBROUTINE add_created_part(p, buffer) USE bsplines USE basic, ONLY: splrz, phinorm, nlclassical USE geometry IMPLICIT NONE TYPE(particles), INTENT(INOUT):: p TYPE(particle), ALLOCATABLE, INTENT(in) :: buffer(:) INTEGER:: i, nptotinit, parts_size_increase, nb_added, newindex real(kind=db), ALLOCATABLE::gtildeloc(:,:) nptotinit=p%Nploc+1 nb_added=size(buffer,1) IF(nb_added .le. 0) RETURN ! No particles to add ! if there is not enough space in the parts simulation buffer, increase the parst size IF(p%Nploc + nb_added .gt. size(p%Z,1)) THEN parts_size_increase=Max(floor(0.1*size(p%Z,1)),nb_added) CALL change_parts_allocation(p, parts_size_increase) END IF newindex=MAXVAL(p%partindex) DO i=1,nb_added p%Nploc=p%Nploc+1 newindex=newindex+1 CALL Insertincomingpart(p, buffer, p%Nploc, i) p%partindex(p%Nploc)=newindex CALL geom_weight(p%Z(p%Nploc),p%R(p%Nploc),p%geomweight(p%Nploc,:)) if( .not. is_inside(p,p%Nploc) ) then p%Nploc=p%Nploc-1 newindex=newindex-1 continue end if call p_calc_rzindex(p,p%Nploc) END DO nb_added=p%Nploc-nptotinit+1 if( .not. p%is_test) then IF(allocated(p%addedlist)) then call change_array_size_int(p%addedlist,2) else allocate(p%addedlist(2)) end if p%addedlist(size(p%addedlist)-1:size(p%addedlist))=(/nptotinit,nb_added/) end if END SUBROUTINE add_created_part function is_inside(p,id) Use basic, ONLY: rgrid,zgrid, nr, nz IMPLICIT NONE logical :: is_inside type(particles) :: p integer :: id is_inside=.true. if(p%geomweight(id,0).le.0)then is_inside=.false. return end if if(p%R(id).ge.rgrid(nr) .or. p%R(id) .le. rgrid(0))then is_inside=.false. return end if if(p%Z(id).ge.zgrid(nz) .or. p%Z(id) .le. zgrid(0))then is_inside=.false. return end if end function is_inside SUBROUTINE calc_newparts_energy(p) USE basic, ONLY: phinorm, nlclassical type(particles)::p integer::i,n,nptotinit,nbadded, nptotend if(p%is_test) return if( allocated(p%addedlist)) then n=size(p%addedlist) !write(*,*) n, "addedlist: ", p%addedlist Do i=1,n,2 nptotinit=p%addedlist(i) nbadded=p%addedlist(i+1) p%nbadded=p%nbadded+nbadded nptotend=nptotinit+nbadded-1 loc_etot0=loc_etot0+p%q*p%weight*sum(p%pot(nptotinit:nptotend))*phinorm IF(.not. nlclassical) THEN loc_etot0=loc_etot0+p%m*p%weight*vlight**2*sum(p%Gamma(nptotinit:nptotend)-1) ELSE loc_etot0=loc_etot0+0.5*p%m*p%weight*vlight**2*sum(p%UR(nptotinit:nptotend)**2 & & +p%UZ(nptotinit:nptotend)**2 & & +p%UTHET(nptotinit:nptotend)**2) END IF end do deallocate(p%addedlist) end if end subroutine calc_newparts_energy !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Delete particle at given index removing its energy from the system ! !> @param [in] index index of particle to be deleted !--------------------------------------------------------------------------- SUBROUTINE delete_part(p, index, replace) !! This will destroy particle at the given index USE constants, ONLY: vlight USE bsplines USE geometry USE basic, ONLY: phinorm, nlclassical,splrz TYPE(particles), INTENT(INOUT):: p INTEGER, INTENT(IN) :: index LOGICAL, OPTIONAL :: replace LOGICAL:: repl IF(present(replace)) THEN repl=replace ELSE repl=.true. END IF !Computes the potential at the particle position with phi_ext+phi_s IF(index .le. p%Nploc) THEN IF(.not. p%is_test) THEN loc_etot0=loc_etot0-p%q*p%weight*(p%pot(index))*phinorm IF(.not. nlclassical) THEN loc_etot0=loc_etot0-p%m*p%weight*vlight**2*(p%Gamma(index)-1) ELSE loc_etot0=loc_etot0-0.5*p%m*p%weight*vlight**2*(p%UR(index)**2+p%UZ(index)**2+p%UTHET(index)**2) END IF END IF IF(repl) THEN ! We fill the gap CALL move_part(p, p%Nploc, index) p%partindex(p%Nploc)=-1 ! Reduce the total number of simulated parts p%Nploc=p%Nploc-1 END IF END IF END SUBROUTINE delete_part !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Move particle with index sourceindex to particle with index destindex. !> !WARNING! This will overwrite particle at destindex. ! !> @param [in] sourceindex index in parts of the particle to move. !> @param [in] destindex index in parts of the moved particle destination. !--------------------------------------------------------------------------- SUBROUTINE move_part(p, sourceindex, destindex) !! This will destroy particle at destindex INTEGER, INTENT(IN) :: destindex, sourceindex TYPE(particles), INTENT(INOUT)::p IF(sourceindex .eq. destindex) RETURN IF(sourceindex .le. 0 .or. destindex .le. 0) RETURN ! Move part at sourceindex in part at destindex Call copy_part(p,sourceindex,destindex,p) END SUBROUTINE move_part !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Copy particle with index sourceindex in particles sourcep to particle with index destindex in particles destp. !> !WARNING! This will overwrite particle at destp(destindex). ! !> @param [inout] sourcep Structure of source particles. !> @param [in] sourceindex index in parts of the particle to move. !> @param [in] destindex index in parts of the moved particle destination. !> @param [inout] destp Structure of source particles. !--------------------------------------------------------------------------- SUBROUTINE copy_part(sourcep, sourceindex, destindex, destp) !! This will destroy particle at destindex INTEGER, INTENT(IN) :: destindex, sourceindex TYPE(particles), INTENT(IN)::sourcep TYPE(particles), INTENT(INOUT)::destp IF(sourceindex .le. 0 .or. destindex .le. 0) RETURN IF( destindex .gt. size(destp%R,1)) RETURN ! Move part at sourceindex in part at destindex destp%partindex(destindex) = sourcep%partindex(sourceindex) destp%Gamma(destindex) = sourcep%Gamma(sourceindex) destp%R(destindex) = sourcep%R(sourceindex) destp%Z(destindex) = sourcep%Z(sourceindex) destp%THET(destindex) = sourcep%THET(sourceindex) destp%UR(destindex) = sourcep%UR(sourceindex) destp%UTHET(destindex) = sourcep%UTHET(sourceindex) destp%UZ(destindex) = sourcep%UZ(sourceindex) destp%Rindex(destindex) = sourcep%Rindex(sourceindex) destp%Zindex(destindex) = sourcep%Zindex(sourceindex) destp%geomweight(destindex,:) = sourcep%geomweight(sourceindex,:) destp%pot(destindex) = sourcep%pot(sourceindex) END SUBROUTINE copy_part !________________________________________________________________________________ SUBROUTINE destroy_parts(p) TYPE(particles) :: p p%Nploc=0 IF(ALLOCATED(p%Z)) DEALLOCATE(p%Z) IF(ALLOCATED(p%R)) DEALLOCATE(p%R) IF(ALLOCATED(p%THET)) DEALLOCATE(p%THET) IF(ALLOCATED(p%BZ)) DEALLOCATE(p%BZ) IF(ALLOCATED(p%BR)) DEALLOCATE(p%BR) IF(ASSOCIATED(p%UR)) DEALLOCATE(p%UR) IF(Associated(p%URold)) DEALLOCATE(p%URold) IF(Associated(p%UZ)) DEALLOCATE(p%UZ) IF(Associated(p%UZold)) DEALLOCATE(p%UZold) IF(Associated(p%UTHET)) DEALLOCATE(p%UTHET) IF(Associated(p%UTHETold)) DEALLOCATE(p%UTHETold) IF(Associated(p%Gamma)) DEALLOCATE(p%Gamma) IF(Associated(p%Gammaold)) DEALLOCATE(p%Gammaold) IF(ALLOCATED(p%Rindex)) DEALLOCATE(p%Rindex) IF(ALLOCATED(p%Zindex)) DEALLOCATE(p%Zindex) IF(ALLOCATED(p%partindex)) DEALLOCATE(p%partindex) if(allocated(p%geomweight)) Deallocate(p%geomweight) END SUBROUTINE !________________________________________________________________________________ SUBROUTINE clean_beam ! INTEGER:: i Do i=1,size(partslist,1) CALL destroy_parts(partslist(i)) END DO ! END SUBROUTINE clean_beam !________________________________________________________________________________ SUBROUTINE swappointer( pointer1, pointer2) REAL(kind=db), DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2 REAL(kind=db), DIMENSION(:), POINTER:: temppointer temppointer=>pointer1 pointer1=>pointer2 pointer2=>temppointer END SUBROUTINE swappointer !_______________________________________________________________________________ SUBROUTINE loaduniformRZ(p, VR,VZ,VTHET) USE basic, ONLY: plasmadim, rnorm, temp, qsim, msim USE constants, ONLY: me, kb, elchar REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p CALL creat_parts(p, size(VR,1)) p%Nploc=size(VR,1) p%Nptot=size(VR,1) p%q=sign(elchar,qsim) p%weight=msim/me p%m=me p%qmRatio=qsim/msim ! Initial distribution in z with normalisation CALL loduni(1,p%Z(1:p%Nploc)) p%Z(1:p%Nploc)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*p%Z(1:p%Nploc))/rnorm ! Initial distribution in r with normalisation CALL loduni(2,p%R(1:p%Nploc)) p%R(1:p%Nploc)=(plasmadim(3)+p%R(1:p%Nploc)*(plasmadim(4)-plasmadim(3)))/rnorm ! Initial velocities distribution CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp) END SUBROUTINE loaduniformRZ !_______________________________________________________________________________ SUBROUTINE loadDavidson(p, VR,VZ,VTHET, lodr) USE constants, ONLY: me, kb, elchar USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, width, qsim, msim, & & omegac, zgrid, nz, rnorm, n0, nblock, temp interface subroutine lodr(nbase,y,rminus,rplus) USE constants REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rplus, rminus end subroutine end interface TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), DIMENSION(:), ALLOCATABLE::ra, rb, z REAL(kind=db) :: r0, deltar2, halfLz, Mirrorratio, Le, VOL INTEGER :: j, n, blockstart, blockend, addedpart, remainparts INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize CALL creat_parts(p, size(VR,1)) p%Nploc=size(VR,1) p%Nptot=p%Nploc Allocate(ra(nblock),rb(nblock), z(0:nblock)) !r0=(plasmadim(4)+plasmadim(3))/2 r0=sqrt(4*H0/(me*omegac**2)) halfLz=(zgrid(nz)+zgrid(0))/2 MirrorRatio=(Rcurv-1)/(Rcurv+1) z(0)=plasmadim(1) DO n=1,nblock ! Compute limits in radius and load radii for each part Le=(plasmadim(2)-plasmadim(1))/nblock*(n-0.5)-halfLz*rnorm+plasmadim(1) z(n)=z(0)+n*(plasmadim(2)-plasmadim(1))/nblock deltar2=1-MirrorRatio*cos(2*pi*Le/width) rb(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2+sqrt(1-P0*abs(omegac)/H0*deltar2)) ra(n)=r0/deltar2*sqrt(1-P0*abs(omegac)/2/H0*deltar2-sqrt(1-P0*abs(omegac)/H0*deltar2)) END DO VOL=SUM(2*pi*MINVAL(ra)*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock) qsim=VOL*n0*elchar/nplasma msim=abs(qsim)/elchar*me p%weight=abs(qsim)/elchar p%m=me p%q=sign(elchar,qsim) p%qmRatio=p%q/p%m blockstart=1 blockend=0 ALLOCATE(blocksize(nblock)) WRITE(*,*) "blocksize: ", size(blocksize), nblock DO n=1,nblock blocksize(n)=nplasma/VOL*2*pi*MINVAL(ra)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock END DO remainparts=p%Nploc-SUM(blocksize) addedpart=1 n=nblock/2 j=1 DO WHILE(remainparts .GT. 0) blocksize(n)=blocksize(n)+addedpart remainparts=remainparts-addedpart n=n+j j=-1*(j+SIGN(1,j)) END DO CALL loadPartSlices(p, lodr, ra, rb, z, blocksize) IF(distribtype .eq. 5) THEN CALL loadGaussianVelocities(p, VR, VZ, VTHET, temp) VZ=VZ/4 VR=VR*8 VTHET=VTHET*8 ELSE Call loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0) END IF END SUBROUTINE loadDavidson SUBROUTINE loadDavidsonVelocities(p, VR,VZ,VTHET, H0, P0) USE constants, ONLY: me, kb, elchar USE basic, ONLY: rnorm, Rcurv, B0, width, vnorm, zgrid, nz TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(INOUT)::VZ(:), VR(:), VTHET(:) REAL(kind=db), INTENT(IN):: H0, P0 REAL(kind=db) :: athetpos, rg, zg, halfLz, Mirrorratio, Pcomp, Acomp INTEGER :: i MirrorRatio=(Rcurv-1)/(Rcurv+1) halfLz=(zgrid(nz)+zgrid(0))/2 ! Load velocities theta velocity ! Loading of r and z velocity is done in adapt_vinit to have ! access to parts%pot DO i=1,p%Nploc ! Interpolation for Magnetic potential rg=p%R(i)*rnorm zg=(p%Z(i)-halfLz)*rnorm Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width)) Pcomp=P0/rg/p%m Acomp=-p%qmRatio*Athetpos VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/p%m)),Pcomp+Acomp) !VTHET(i)=Pcomp+Acomp END DO VTHET=VTHET/vnorm VZ=0._db VR=0._db p%Davidson=.true. p%H0=H0 p%P0=P0 END SUBROUTINE loadDavidsonvelocities SUBROUTINE loadGaussianVelocities(p, VR,VZ,VTHET, temperature) USE basic, ONLY: vnorm USE constants, ONLY: kb REAL(kind=db), INTENT(inout) ::VZ(:), VR(:), VTHET(:) TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(IN):: temperature REAL(kind=db):: vth ! Initial velocities distribution vth=sqrt(kb*temperature/p%m)/vnorm !thermal velocity CALL lodgaus(3,VZ) CALL lodgaus(5,VR) CALL lodgaus(7,VTHET) VZ=VZ*vth VR=VR*vth VTHET=VTHET*vth p%temperature=temperature p%Davidson=.false. END SUBROUTINE loadGaussianVelocities SUBROUTINE loadPartslices(p, lodr, ra, rb, z, blocksize) USE basic, ONLY: rnorm interface subroutine lodr(nbase,y,rminus,rplus) USE constants REAL(kind=db), INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase REAL(kind=db), INTENT(in) :: rplus, rminus end subroutine end interface TYPE(particles), INTENT(INOUT):: p REAL(kind=db), INTENT(IN)::ra(:), rb(:), z(0:) INTEGER, DIMENSION(:), INTENT(IN) :: blocksize INTEGER :: n, blockstart, blockend, nblock nblock=size(blocksize,1) blockstart=1 blockend=0 DO n=1,nblock blockstart=blockend+1 blockend=MIN(blockstart+blocksize(n)-1,p%Nploc) ! Initial distribution in z with normalisation between magnetic mirrors CALL loduni(1,p%Z(blockstart:blockend)) p%Z(blockstart:blockend)= (z(n-1)+p%Z(blockstart:blockend)*(z(n)-z(n-1)))/rnorm CALL lodr(2,p%R(blockstart:blockend),ra(n), rb(n)) p%R(blockstart:blockend)=p%R(blockstart:blockend)/rnorm END DO END SUBROUTINE loadPartslices SUBROUTINE loadPartFile(p, partfileindex, VR, VZ, VTHET) USE basic, ONLY: partfile, nplasma, lu_partfile implicit None TYPE(particles), INTENT(INOUT):: p REAL(kind=db), DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VR, VZ, VTHET INTEGER, INTENT(IN):: partfileindex INTEGER:: nblock = 0 REAL(kind=db), Dimension(:), ALLOCATABLE:: ra, rb, z INTEGER, Dimension(:), ALLOCATABLE:: npartsslice INTEGER:: velocitytype=1 !< 1) gaussian with temp 2) Davidson with H0, P0 INTEGER:: radialtype=1 !< 1) 1/R 2) uniform 3) 1/R^2 4) gauss INTEGER:: npartsalloc !< initial size of particles arrays REAL(kind=db):: mass=me REAL(kind=db):: charge=-elchar REAL(kind=db):: weight=1.0 CHARACTER(len=256) :: header=' ' !< header of csv file REAL(kind=db):: H0=3.2e-14 !< Total energy REAL(kind=db):: P0=8.66e-25 !< Canonical angular momentum REAL(kind=db):: temperature=10000 !< temperature in kelvins real(kind=db):: n0 !< density factor LOGICAL :: is_test = .false. !< Defines if particle are test particles or tracers INTEGER:: i, ierr, openerr NAMELIST /partsload/ nblock, mass, charge, weight, npartsalloc, velocitytype, & & radialtype, temperature, H0, P0, is_test, n0 OPEN(UNIT=lu_partfile,FILE=trim(partfile(partfileindex)),ACTION='READ',IOSTAT=openerr) header=' ' IF(openerr .ne. 0) THEN CLOSE(unit=lu_partfile) RETURN END IF READ(lu_partfile,partsload) IF(mpirank .eq.0) THEN WRITE(*,'(a,i2,a,a)')"reading partfile: ", partfileindex, " ", partfile(partfileindex) WRITE(*,partsload) END IF IF( nblock .ge. 1) THEN ALLOCATE(z(0:nblock),ra(nblock),rb(nblock), npartsslice(nblock)) DO WHILE(header(1:8) .ne. '//slices') READ(lu_partfile,'(a)') header END DO DO i=1,nblock READ(lu_partfile,*) z(i-1),ra(i),rb(i),npartsslice(i) END DO READ(lu_partfile,*) z(nblock) CALL creat_parts(p,max(npartsalloc,sum(npartsslice))) p%Nploc=sum(npartsslice) p%Nptot=p%Nploc IF(partfileindex .eq. 1) nplasma=p%Nploc IF( allocated(VR) ) THEN DEALLOCATE(VR,VZ,VTHET) end if if(.not. allocated(VR)) THEN ALLOCATE(VR(p%Nploc)) ALLOCATE(VZ(p%Nploc)) ALLOCATE(VTHET(p%Nploc)) END IF p%m=mass p%q=charge p%weight=weight p%qmRatio=charge/mass p%is_test=is_test SELECT CASE(radialtype) CASE(1) ! 1/R distribution in R CALL loadPartslices(p, lodunir, ra, rb, z, npartsslice) CASE(2) ! flat top distribution in R CALL loadPartslices(p, lodlinr, ra, rb, z, npartsslice) CASE(3) ! 1/R^2 distribution in R CALL loadPartslices(p, lodinvr, ra, rb, z, npartsslice) CASE(4) ! gaussian distribution in R CALL loadPartslices(p, lodgausr, ra, rb, z, npartsslice) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of radial distribution:", radialtype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT SELECT CASE(velocitytype) CASE(1) ! Gaussian with temperature CALL loadGaussianVelocities(p, VR, VZ, VTHET, temperature) CASE(2) CALL loadDavidsonVelocities(p, VR, VZ, VTHET, H0, P0) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of velocity distribution:", velocitytype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT END IF CLOSE(unit=lu_partfile) END SUBROUTINE !=============================================================================== SUBROUTINE Temp_rescale ! USE basic, ONLY: temprescale, nlclassical, nz ! INTEGER:: i, gridindex ! REAL(kind=db) :: vr, vz, vthet ! DO i=1,parts%Nptot ! gridindex=parts%zindex(i)+1+parts%rindex(i)*nz ! vr=parts%UR(i)/parts%Gamma(i)-fluidur(gridindex) ! vz=parts%UZ(i)/parts%Gamma(i)-fluiduz(gridindex) ! vthet=parts%UTHET(i)/parts%Gamma(i)-fluiduthet(gridindex) ! parts%UR(i)=vr*temprescale+fluidur(gridindex) ! parts%UTHET(i)=vthet*temprescale+fluiduthet(gridindex) ! parts%UZ(i)=vz*temprescale+fluiduz(gridindex) ! IF(nlclassical) THEN ! parts%Gamma(i)=1.0 ! ELSE ! parts%Gamma(i)=1/sqrt(1-(parts%UR(i)**2+parts%UTHET(i)**2+parts%UZ(i)**2)) ! END IF ! parts%UR(i)=parts%UR(i)*parts%Gamma(i) ! parts%UTHET(i)=parts%UTHET(i)*parts%Gamma(i) ! parts%UZ(i)=parts%UZ(i)*parts%Gamma(i) ! END DO END SUBROUTINE Temp_rescale !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief Increase the number of macroparticles by separating each previous macroparticles into !> samplefactor new macroparticles of equally divided weight. The new sub particles are distributed !> uniformly in space to maintain the density and other moments. ! !> @param [in] samplefactor multiplicator of the number of macroparticles. !> @param [in] p particles type to increase. !--------------------------------------------------------------------------- SUBROUTINE upsample(p, samplefactor) USE basic, ONLY : nplasma, dr, dz INTEGER, INTENT(IN) ::samplefactor TYPE(particles), INTENT(INOUT):: p INTEGER:: i, j, currentindex REAL(kind=db), DIMENSION(p%Nploc) :: spreaddir ! random direction for the spread of each initial macro particle REAL(kind=db) :: dir ! Direction in which the particle is moved REAL(kind=db) :: dl ! Particle displacement used for ! Load and scale the direction angle for spreading the new particles CALL loduni(2, spreaddir) spreaddir=spreaddir*2*pi/samplefactor dl=min(dz,minval(dr))/100 DO i=1,p%Nploc DO j=1,samplefactor-1 currentindex=p%Nploc+(i-1)*(samplefactor-1)+j CALL move_part(p,i,currentindex) p%partindex(currentindex)=currentindex dir = spreaddir(i)+2*pi*j/samplefactor p%R(currentindex)=p%R(currentindex) + dl*cos(dir) p%Z(currentindex)=p%Z(currentindex) + dl*sin(dir) END DO p%partindex(i)=i p%R(i)=p%R(i) + dl*cos(spreaddir(i)) p%Z(i)=p%Z(i) + dl*sin(spreaddir(i)) END DO nplasma=nplasma*samplefactor p%weight=p%weight/samplefactor p%Nploc=p%Nploc*samplefactor END SUBROUTINE upsample END MODULE beam diff --git a/src/diagnose.f90 b/src/diagnose.f90 index 53ded8f..c1a082b 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,369 +1,369 @@ SUBROUTINE diagnose(kstep) ! ! Diagnostics ! USE basic USE futils USE hashtable Use maxwsrce USE beam, ONLY : partslist, epot, ekin, etot, etot0, Nplocs_all, collectparts #if USE_X == 1 USE xg, ONLY : initw, updt_xg_var #endif USE fields, ONLY: phi_spline USE celldiag USE mpi Use geometry IMPLICIT NONE ! INTEGER, INTENT(in) :: kstep ! ! Local vars and arrays INTEGER, PARAMETER :: BUFSIZE = 20 CHARACTER(len=128) :: str, fname, grpname INTEGER:: Ntotal ! Total number of simulated particles remaining in the simulation INTEGER:: partsrank, partsdim(2) INTEGER:: Nplocs_all_save(64) INTEGER:: i !________________________________________________________________________________ ! 1. Initial diagnostics IF( kstep .EQ. 0 .and. mpirank .eq. 0) THEN ! Only process 0 should save on file ! WRITE(*,'(a)') ' Initial diagnostics' ! ! 1.1 Initial run or when NEWRES set to .TRUE. ! IF( .NOT. nlres .OR. newres) THEN CALL creatf(resfile, fidres, 'Espic2d result file', real_prec='d') WRITE(*,'(3x,a,a)') TRIM(resfile), ' created' ! ! Label the run IF( LEN_TRIM(label1).GT.0 ) CALL attach(fidres, "/", "label1", TRIM(label1)) IF( LEN_TRIM(label2).GT.0 ) CALL attach(fidres, "/", "label2", TRIM(label2)) IF( LEN_TRIM(label3).GT.0 ) CALL attach(fidres, "/", "label3", TRIM(label3)) IF( LEN_TRIM(label4).GT.0 ) CALL attach(fidres, "/", "label4", TRIM(label4)) ! ! Job number jobnum = 0 ! ! Data group CALL creatg(fidres, "/data", "data") CALL creatg(fidres, "/data/var0d", "0d history arrays") CALL creatg(fidres, "/data/var1d","1d history arrays") CALL creatg(fidres, "/data/part", "part phase space") CALL creatg(fidres, "/data/fields", "Electro static potential and Er Ez fields") ! ! File group CALL creatg(fidres, "/files", "files") CALL attach(fidres, "/files", "jobnum", jobnum) CALL putarr(fidres, "/data/var1d/rgrid", rgrid) CALL putarr(fidres, "/data/var1d/zgrid", zgrid) CALL creatd(fidres, 1, (/ 64 /), "/data/var0d/Nplocs_all", "Nplocs_all") ! Part and fields vectors ! Initialize time-dependent particle variables CALL creatd(fidres, 1, SHAPE(partslist(1)%R), "/data/part/R", "R") CALL creatd(fidres, 1, SHAPE(partslist(1)%Z), "/data/part/Z", "Z") CALL creatd(fidres, 1, SHAPE(partslist(1)%Z), "/data/part/THET", "THET") CALL creatd(fidres, 1, SHAPE(partslist(1)%UZ), "/data/part/UR", "UZ") CALL creatd(fidres, 1, SHAPE(partslist(1)%UR), "/data/part/UZ", "UR") CALL creatd(fidres, 1, SHAPE(partslist(1)%UTHET), "/data/part/UTHET", "UTHET") CALL creatd(fidres, 1, SHAPE(partslist(1)%Rindex), "/data/part/Rindex", "Rindex") CALL creatd(fidres, 1, SHAPE(partslist(1)%Zindex), "/data/part/Zindex", "Zindex") CALL creatd(fidres, 1, SHAPE(partslist(1)%partindex), "/data/part/partindex", "partindex") CALL creatd(fidres, 1, SHAPE(partslist(1)%pot), "/data/part/pot", "pot") CALL creatd(fidres, 0, SHAPE(time), "/data/part/time", "time" ) CALL creatd(fidres, 0, SHAPE(Ntotal), "/data/part/Nparts", "number of remaining parts in the simulation space") CALL attach(fidres,'/data/part', "q", partslist(1)%q) CALL attach(fidres,'/data/part', "m", partslist(1)%m) CALL attach(fidres,'/data/part', "weight", partslist(1)%weight) DO i=2,nbspecies IF ( .not. partslist(i)%is_test) CONTINUE WRITE(grpname,'(a,i2)')'/data/part/',i CALL creatg(fidres, grpname, "specific specie phase space") CALL creatd(fidres, 0, SHAPE(time), trim(grpname) // "/time", "time") CALL creatd(fidres, 0, SHAPE(time), trim(grpname) //"/Nparts", "number of remaining parts") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/R", "radial pos") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/Z", "axial pos") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/THET", "azimuthal pos") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/UZ", "axial beta*gamma") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/UR", "radial beta*gamma") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/UTHET", "azimuthal beta*gamma") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/pot", "electric potential") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/Rindex", "radial grid index") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/Zindex", "axial grid index") CALL creatd(fidres, 1, SHAPE(partslist(i)%R), trim(grpname) // "/partindex", "particle index") CALL attach(fidres,trim(grpname), "q", partslist(i)%q) CALL attach(fidres,trim(grpname), "m", partslist(i)%m) CALL attach(fidres,trim(grpname), "weight", partslist(i)%weight) END DO CALL creatd(fidres, 1, SHAPE(pot), "/data/fields/pot", "pot") CALL creatd(fidres, 1, SHAPE(Er), "/data/fields/Er", "Er") CALL creatd(fidres, 1, SHAPE(Ez), "/data/fields/Ez", "Ez") CALL creatd(fidres, 1, SHAPE(phi_spline), "/data/fields/phi", "spline form of Phi") CALL creatd(fidres, 2, SHAPE(moments), "/data/fields/moments", "moments",compress=.true.,chunking=(/1,nrank(1)*nrank(2)/)) !CALL creatd(fidres, 2, SHAPE(moments), "/data/fields/moments", "moments") CALL creatd(fidres, 0, SHAPE(time), "/data/fields/time", "time" ) CALL putarr(fidres, "/data/fields/Br", Br) CALL putarr(fidres, "/data/fields/Bz", Bz) CALL putarr(fidres, "/data/fields/Athet", Athet) CALL putarr(fidres, "/data/fields/volume", Volume) ! ! 1.2 Restart run ! ELSE CALL cp2bk(resfile) ! backup previous result file CALL openf(resfile, fidres, real_prec='d') WRITE(*,'(3x,a,a)') TRIM(resfile), ' open' CALL getatt(fidres, "/files", "jobnum", jobnum) jobnum = jobnum+1 WRITE(*,'(3x,a,i3)') "Current Job Number =", jobnum CALL attach(fidres, "/files", "jobnum", jobnum) END IF ! ! Add input namelist variables as attributes of /data/input WRITE(str,'(a,i2.2)') "/data/input.",jobnum CALL creatg(fidres, TRIM(str)) CALL attach(fidres, TRIM(str), "job_time", job_time) CALL attach(fidres, TRIM(str), "extra_time", extra_time) CALL attach(fidres, TRIM(str), "dt", dt) CALL attach(fidres, TRIM(str), "tmax", tmax) CALL attach(fidres, TRIM(str), "nrun", nrun) CALL attach(fidres, TRIM(str), "nlres", nlres) CALL attach(fidres, TRIM(str), "nlsave", nlsave) CALL attach(fidres, TRIM(str), "newres", newres) CALL attach(fidres, TRIM(str), "nz", nz) CALL putarr(fidres, TRIM(str)//"/lz", lz) CALL attach(fidres, TRIM(str), "nplasma", nplasma) CALL attach(fidres, TRIM(str), "potinn", potinn) CALL attach(fidres, TRIM(str), "potout", potout) CALL attach(fidres, TRIM(str), "B0", B0) CALL attach(fidres, TRIM(str), "Rcurv", Rcurv) CALL attach(fidres, TRIM(str), "width", width) CALL attach(fidres, TRIM(str), "n0", n0) CALL attach(fidres, TRIM(str), "temp", partslist(1)%temperature) CALL attach(fidres, TRIM(str), "it0d", it0d) CALL attach(fidres, TRIM(str), "it2d", it2d) CALL attach(fidres, TRIM(str), "itparts", itparts) CALL attach(fidres, TRIM(str), "nlclassical", nlclassical) CALL attach(fidres, TRIM(str), "nlPhis", nlPhis) CALL attach(fidres, TRIM(str), "qsim", partslist(1)%q*partslist(1)%weight) CALL attach(fidres, TRIM(str), "msim", partslist(1)%m*partslist(1)%weight) CALL attach(fidres, TRIM(str), "startstep", cstep) CALL attach(fidres, TRIM(str), "H0", partslist(1)%H0) CALL attach(fidres, TRIM(str), "P0", partslist(1)%P0) CALL putarr(fidres, TRIM(str)//"/femorder", femorder) CALL putarr(fidres, TRIM(str)//"/ngauss", ngauss) CALL putarr(fidres, TRIM(str)//"/nnr", nnr) CALL putarr(fidres, TRIM(str)//"/radii", radii) CALL putarr(fidres, TRIM(str)//"/plasmadim", plasmadim) CALL attach(fidres, TRIM(str), "rawparts", .true.) CALL attach(fidres, TRIM(str), "nbspecies", nbspecies) ! Save geometry parameters for non conforming boundary conditions Call geom_diag(fidres,str,rnorm) ! Save Maxwellsource parameters for the ad-hoc source Call maxwsrce_diag(fidres,str,vnorm) ! Save STDIN of this run WRITE(str,'(a,i2.2)') "/files/STDIN.",jobnum INQUIRE(unit=lu_in, name=fname) CALL putfile(fidres, TRIM(str), TRIM(fname)) ! ! Initialize buffers for 0d history arrays CALL htable_init(hbuf0, BUFSIZE) CALL set_htable_fileid(hbuf0, fidres, "/data/var0d") ! ! Initialize Xgrafix #if USE_X == 1 IF(nlxg) THEN CALL initw END IF #endif END IF IF(kstep .EQ. 0) THEN ! Initialize particle cell diagnostic CALL celldiag_init(lu_in,time, fidres) CLOSE(lu_in) END IF !________________________________________________________________________________ ! 2. Periodic diagnostics IF( kstep .NE. -1) THEN IF(modulo(step,ittracer) .eq. 0 .or. nlend) THEN ! We gather the traced particles on the mpi host DO i=1,nbspecies IF(partslist(i)%is_test) CALL collectparts(partslist(i)) END DO END IF - IF(modulo(step,itrestart) .eq. 0 .or. nlend) THEN + IF(modulo(step,itrestart) .eq. 0 .or. modulo(step,itparts) .eq. 0 .or. nlend) THEN ! We gather the traced particles on the mpi host DO i=1,nbspecies CALL collectparts(partslist(i)) END DO END IF IF(modulo(step,ittext) .eq. 0 .or. nlend) THEN ! We gather the number of gained and lost particles on the mpi host do i=1,nbspecies IF(mpirank .eq.0 ) THEN CALL MPI_REDUCE(MPI_IN_PLACE, partslist(i)%nbadded, 1, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) CALL MPI_REDUCE(MPI_IN_PLACE, partslist(i)%nblost, 4, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) ELSE CALL MPI_REDUCE(partslist(i)%nbadded, partslist(i)%nbadded, 1, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) CALL MPI_REDUCE(partslist(i)%nblost, partslist(i)%nblost, 4, MPI_INTEGER, MPI_SUM, & & 0, MPI_COMM_WORLD, ierr) partslist(i)%nbadded=0 partslist(i)%nblost=0 END IF end do end if ! Cell diag quantities IF(modulo(step,itcelldiag).eq. 0 .or. nlend) THEN ! CALL celldiag_save(time, fidres) END IF ! ! Only process 0 should save on file IF(mpirank .ne. 0) RETURN ! IF (mpisize .gt. 1) THEN partslist(1)%Nptot=sum(Nplocs_all) END IF ! IF(modulo(step,ittext).eq. 0 .or. nlend) THEN WRITE(*,'(a,1x,i8.8,a1,i8.8,20x,a,1pe10.3)') '*** Timestep (this run/total) =', & & step, '/', cstep, 'Time =', time if( abs(etot).gt. 0) then WRITE(*,'(a,5(1pe12.4),1x,i8.8,a1,i8.8)') 'Epot, Ekin, Etot, Eerr, Eerr rel, Nbparts/Ntotal', epot, ekin, etot, etot-etot0,(etot-etot0)/etot, partslist(1)%Nptot,'/',nplasma else WRITE(*,'(a,4(1pe12.4),1x,i8.8,a1,i8.8)') 'Epot, Ekin, Etot, Eerr, Nbparts/Ntotal', epot, ekin, etot, etot-etot0, partslist(1)%Nptot,'/',nplasma end if IF(mpisize .gt. 1 ) then WRITE(*,'(a,64i10.7)') 'Nbparts per proc', Nplocs_all end if Write(*,'(a)')"speci, added, zmloss, zploss, rmloss, rploss, tot var, tot" do i=1,nbspecies WRITE(*,'(i04,7i9.7)') i, partslist(i)%nbadded,-partslist(i)%nblost,partslist(i)%nbadded-sum(partslist(i)%nblost), partslist(i)%Nptot partslist(i)%nbadded=0 partslist(i)%nblost=0 end do END IF !________________________________________________________________________________ ! ! 2.1 0d history arrays ! IF(modulo(step,it0d).eq. 0 .or. nlend) THEN CALL add_record(hbuf0, "time", "simulation time", time) CALL add_record(hbuf0, "epot", "potential energy", epot) CALL add_record(hbuf0, "ekin", "kinetic energy", ekin) CALL add_record(hbuf0, "etot", "total energy", etot) CALL add_record(hbuf0, "etot0", "theoretical total energy", etot0) CALL add_record(hbuf0, "nbparts", "number of remaining parts in the simulation space", REAL(partslist(1)%Nptot,kind=db)) CALL htable_endstep(hbuf0) Nplocs_all_save=0 Nplocs_all_save(1:mpisize)=Nplocs_all(0:mpisize-1) CALL append(fidres, "/data/var0d/Nplocs_all", REAL(Nplocs_all_save,kind=db)) END IF ! ! 2.2 2d profiles IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL append(fidres, "/data/fields/time", time) CALL append(fidres, "/data/fields/pot", pot*phinorm) CALL append(fidres, "/data/fields/Er", Er*enorm) CALL append(fidres, "/data/fields/Ez", Ez*enorm) CALL append(fidres, "/data/fields/phi", phi_spline*phinorm) CALL append(fidres, "/data/fields/moments", moments) END IF ! ! 2.3 main specie quantities IF(modulo(step,itparts).eq. 0 .or. nlend) THEN !PRINT*, 'write particles to file_____________________' CALL append(fidres, "/data/part/time", time) CALL append(fidres, "/data/part/Nparts", REAL(partslist(1)%Nptot,kind=db)) IF ( isdataset(fidres,'/data/part/R') ) THEN CALL getdims(fidres, '/data/part/R', partsrank, partsdim) partsdim(1)=min(size(partslist(1)%R,1), partsdim(1)) CALL append(fidres, "/data/part/R", partslist(1)%R(1:partsdim(1))*rnorm) CALL append(fidres, "/data/part/Z", partslist(1)%Z(1:partsdim(1))*rnorm) CALL append(fidres, "/data/part/THET", partslist(1)%THET(1:partsdim(1))) CALL append(fidres, "/data/part/UZ", partslist(1)%UZ(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))) CALL append(fidres, "/data/part/UR", partslist(1)%UR(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))) CALL append(fidres, "/data/part/UTHET", partslist(1)%UTHET(1:partsdim(1))/partslist(1)%gamma(1:partsdim(1))) CALL append(fidres, "/data/part/pot", partslist(1)%pot(1:partsdim(1))*phinorm) CALL append(fidres, "/data/part/Rindex", REAL(partslist(1)%Rindex(1:partsdim(1)),kind=db)) CALL append(fidres, "/data/part/Zindex", REAL(partslist(1)%Zindex(1:partsdim(1)),kind=db)) CALL append(fidres, "/data/part/partindex", REAL(partslist(1)%partindex(1:partsdim(1)),kind=db)) END IF END IF ! ! 2.4 Tracer quantities IF(modulo(step,ittracer).eq. 0 .or. nlend) THEN !PRINT*, 'write particles to file_____________________' DO i=2,nbspecies IF ( .not. partslist(i)%is_test) CYCLE WRITE(grpname,'(a,i2,a)')'/data/part/',i,'/' CALL append(fidres, trim(grpname) // "time", time) CALL append(fidres, trim(grpname) //"Nparts", REAL(partslist(i)%Nptot,kind=db)) IF ( isdataset(fidres,trim(grpname)//'R') ) THEN CALL getdims(fidres, trim(grpname) // 'R', partsrank, partsdim) partsdim(1)=min(size(partslist(i)%R,1), partsdim(1)) CALL append(fidres, trim(grpname) // "R", partslist(i)%R(1:partsdim(1))*rnorm) CALL append(fidres, trim(grpname) // "Z", partslist(i)%Z(1:partsdim(1))*rnorm) CALL append(fidres, trim(grpname) // "THET", partslist(i)%THET(1:partsdim(1))) CALL append(fidres, trim(grpname) // "UZ", partslist(i)%UZ(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1))) CALL append(fidres, trim(grpname) // "UR", partslist(i)%UR(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1))) CALL append(fidres, trim(grpname) // "UTHET", partslist(i)%UTHET(1:partsdim(1))/partslist(i)%gamma(1:partsdim(1))) CALL append(fidres, trim(grpname) // "pot", partslist(i)%pot(1:partsdim(1))*phinorm) CALL append(fidres, trim(grpname) // "Rindex", REAL(partslist(i)%Rindex(1:partsdim(1)),kind=db)) CALL append(fidres, trim(grpname) // "Zindex", REAL(partslist(i)%Zindex(1:partsdim(1)),kind=db)) CALL append(fidres, trim(grpname) // "partindex", REAL(partslist(i)%partindex(1:partsdim(1)),kind=db)) END IF END DO ! END IF ! 2.5 3d profiles ! ! ! 2.6 Xgrafix ! #if USE_X == 1 IF(nlxg .AND. modulo(kstep,itgraph) .eq. 0) THEN call xgevent CALL updt_xg_var CALL xgupdate END IF #endif !________________________________________________________________________________ ! 3. Final diagnostics ELSE ! Only process 0 should save on file IF(mpirank .ne. 0) RETURN ! ! Flush 0d history array buffers CALL htable_hdf5_flush(hbuf0) ! ! Close all diagnostic files CALL closef(fidres) !________________________________________________________________________________ END IF ! END SUBROUTINE diagnose diff --git a/src/resume.f90 b/src/resume.f90 index f6c259d..d9015e8 100644 --- a/src/resume.f90 +++ b/src/resume.f90 @@ -1,56 +1,57 @@ SUBROUTINE resume ! ! Resume from previous run ! Use beam Use basic Use fields Use sort Use maxwsrce Use geometry IMPLICIT NONE ! ! Local vars and arrays INTEGER:: i !________________________________________________________________________________ WRITE(*,'(a)') ' Resume from previous run' !________________________________________________________________________________ ! ! Open and read initial conditions from restart file CALL chkrst(0) - - ! TODO: change subroutine to take into account multiple species IF (newres) THEN cstep=0 time=0 END IF call read_geom(lu_in, rnorm, rgrid, Potinn, Potout) ! Compute Self consistent electric field call init_mag CALL fields_init + call localisation(partslist(1)) IF(mpisize .gt. 1) THEN CALL distribpartsonprocs(partslist(1)) DO i=1,nbspecies CALL keep_mpi_self_parts(partslist(i),Zbounds) + call collectparts(partslist(i)) END DO END IF + CALL bound(partslist(1)) CALL localisation(partslist(1)) ! Init Electric and Magnetic Fields CALL fields_comm_init WRITE(*,*) "Fields initialized" CALL rhscon(partslist(1)) WRITE(*,*) "Initial rhs computed" CALL poisson(splrz) WRITE(*,*) "Initial field solved" CALL EFieldscompatparts(partslist(1)) WRITE(*,*) "Initial forces computed" IF (newres) THEN CALL diagnostics WRITE(*,*) "Initial beam%diagnostics done!" END IF IF (nlmaxwellsource) CALL maxwsrce_init(lu_in, time) ! END SUBROUTINE resume