diff --git a/matlab/Parallel.m b/matlab/Parallel.m index 330aefe..1c22402 100644 --- a/matlab/Parallel.m +++ b/matlab/Parallel.m @@ -1,29 +1,31 @@ clear all N=[1,2,4,8]; t=[469.8, 294.1,253.2, 158.5]; Speedup=t(1)./t; Efficiency=Speedup./N*100; figure plot(N,Speedup) xlabel('N') ylabel('Speed up') figure plot(N,Efficiency) xlabel('N') ylabel('Efficiency %') -N=[1,2,4]; -t=[182.6, 113.2,64.4]; +N=[ 1, 2, 4, 8, 12, 16, 24]; +t=[617.1, 389.4, 287.2, 211.9, 182.7, 142.4, 132.7]; Speedup=t(1)./t; Efficiency=Speedup./N*100; figure plot(N,Speedup) +hold on +plot(N,N) xlabel('N') ylabel('Speed up') figure plot(N,Efficiency) xlabel('N') ylabel('Efficiency %') \ No newline at end of file diff --git a/matlab/PlotEspic2dgriddata.m b/matlab/PlotEspic2dgriddata.m index 2aa63a5..1d13c5b 100644 --- a/matlab/PlotEspic2dgriddata.m +++ b/matlab/PlotEspic2dgriddata.m @@ -1,49 +1,49 @@ 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); pcolor(ax1,M.zgrid,M.rgrid,M.N(:,:,fieldstep)); 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'; view(ax1,2) -set(ax1,'colorscale','log') +%set(ax1,'colorscale','log') ax2=subplot(2,2,2,'Parent',fig); contourf(ax2,M.zgrid,M.rgrid,M.pot(:,:,fieldstep)'); %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') ax3=subplot(2,2,3,'Parent',fig); contourf(ax3,M.zgrid,M.rgrid,M.Er(:,:,fieldstep)') xlabel(ax3,'Z [m]') ylabel(ax3,'R [m]') c = colorbar(ax3); c.Label.String= 'E_r [V/m]'; title(ax3,'Radial Electric field') grid(ax3, 'on') ax4=subplot(2,2,4,'Parent',fig); contourf(ax4,M.zgrid,M.rgrid,M.Ez(:,:,fieldstep)') xlabel(ax4,'Z [m]') ylabel(ax4,'R [m]') c = colorbar(ax4); c.Label.String= 'E_z [V/m]'; title(ax4,'Axial Electric field') grid(ax4, 'on') end diff --git a/matlab/analyse_espic2D.m b/matlab/analyse_espic2D.m index 1899e50..d8c36fd 100644 --- a/matlab/analyse_espic2D.m +++ b/matlab/analyse_espic2D.m @@ -1,157 +1,150 @@ -file='teststable126.h5'; -close all hidden; +file='teststablefinefine.h5'; +%close all hidden; if (~exist('M','var')) M.file=file; - M.type='all'; end M=load_espic2d(file,M,'all'); tmin=2; tmax=length(M.ekin); figure plot(M.t0d(tmin:tmax),M.ekin(tmin:tmax),'o-',... M.t0d(tmin:tmax),M.epot(tmin:tmax),'d-',... M.t0d(tmin:tmax),M.etot(tmin:tmax),'h-',... M.t0d(tmin:tmax),M.eerr(tmin:tmax),'x--') legend('ekin', 'epot', 'etot','eerr') xlabel('Time [s]') ylabel('Energies [J]') grid on figure semilogy(M.t0d(tmin:tmax),abs(M.eerr(tmin:tmax)/M.etot(2)*100),'h-') xlabel('t [s]') ylabel('E_{err} %') grid on figure [Z,R] = meshgrid (M.zgrid,M.rgrid); subplot(2,1,1) contourf(Z,R,M.Br') xlabel('Z') ylabel('R') c = colorbar; c.Label.String= 'B_r [T]'; grid on subplot(2,1,2) contourf(Z,R,M.Bz') xlabel('Z [m]') ylabel('R [m]') c = colorbar; c.Label.String= 'B_z [T]'; grid on figure contourf(Z,R,M.B',100,'LineColor','none'); shading flat colormap('jet') hold on quiver(Z,R,M.Bz',M.Br','k') % h=streamslice(Z,R,M.Bz',M.Br'); % set(h,'Color','w') % figure % ax3=subplot(2,1,1); % contourf(ax3,M.zgrid,M.rgrid,mean(M.Er(:,:,end-2000:end),3)') % xlabel(ax3,'Z [m]') % ylabel(ax3,'R [m]') % c = colorbar(ax3); % c.Label.String= 'E_r [V/m]'; % title(ax3,'Average Radial Electric field') % grid on % ax3.Children(1).ZDataSource='data.Er'; % % ax4=subplot(2,1,2); % contourf(ax4,M.zgrid,M.rgrid,mean(M.Ez(:,:,end-2000:end),3)') % xlabel(ax4,'Z [m]') % ylabel(ax4,'R [m]') % c = colorbar(ax4); % c.Label.String= 'E_z [V/m]'; % title(ax4,'Average Axial Electric field') % grid on % ax4.Children(1).ZDataSource='data.Ez'; -figure -zpos=33; +f=figure(); +zpos=65; tinit=1; tend=size(M.N,3); -deltat=floor(tend/2); +deltat=200;%floor(tend/2); +tstudied=deltat/10; H0=3.2e-14; P0=8.78e-25; + +P0=double(mean(mean(M.P(:,end-tstudied:end)))); +H0=double(mean(mean(M.H(:,end-tstudied:end)))); Mirrorratio=(M.Rcurv-1)/(M.Rcurv+1); -pot=M.pot(zpos,:,tend); +pot=mean(M.pot(zpos,:,(tend-deltat):tend),3); psi=1+M.qe*pot(:)/H0-1/(2*M.me*H0)*(P0./M.rgrid+M.qe*0.5*M.B0.*(M.rgrid-M.width/pi*Mirrorratio*cos(2*pi*M.zgrid(zpos)/M.width)*besseli(1,2*pi*M.rgrid/M.width))).^2; [~,I]=max(M.N(:,zpos,tinit)); -plot(M.rgrid,M.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',M.t2d(tinit)*1e9)) +r=(M.rgrid(2:end)+M.rgrid(1:end-1))/2; +plot(r,M.N(1:end-1,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',M.t2d(tinit)*1e9)) hold on -plot(M.rgrid,mean(M.N(:,zpos,tend-deltat:tend),3),'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',M.t2d(tend-deltat)*1e9,M.t2d(tend)*1e9)) -%plot(M.rgrid(I-2:end),1./M.rgrid(I-2:end)*M.N(I,zpos,tinit)*M.rgrid(I),'DisplayName','N=c*1/r') +plot(r,mean(M.N(1:end-1,zpos,(tend-deltat):tend),3),'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',M.t2d(tend-deltat)*1e9,M.t2d(tend)*1e9)) +plot(r(I-2:end),1./r(I-2:end)*mean(M.N(I,zpos,(tend-deltat):tend),3)*r(I),'DisplayName','n=c*1/r') xlabel('R [m]') -ylabel('N') +ylabel('n [m^{-3}]') J=find(psi>0); -plot(M.rgrid(J),psi(J)) +Indices=(J(1)-1):(J(end)+1); +plot(M.rgrid(Indices),psi(Indices)/max(psi(Indices))*M.N(I,zpos,tinit),'+--','DisplayName','scaled \Psi(r) [a.u.]') title(sprintf('Number of particles at z=%1.2e[m]',M.zgrid(zpos))) legend +grid on +xlim([0 0.02]) +f.PaperOrientation='landscape'; +[~, name, ~] = fileparts(M.file); +print(f,sprintf('%s_Density_cut',name),'-dpdf','-fillpage') + +f=figure(); +legtext=sprintf("t=%2.1f - %2.1f [ns]",M.tpart(end-tstudied)*1e9,M.tpart(end)*1e9); +subplot(1,2,1) +H=M.H(:,1); +h1=histogram(H,20,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +H=mean(M.H(:,end-tstudied:end),2); +h1=histogram(H,'BinWidth',h1.BinWidth,'DisplayName',legtext); +ylabel('counts') +xlabel('H [J]') +legend + +subplot(1,2,2) +P=M.P(:,1); +h2=histogram(P,20,'BinLimits',[min(P(:)) max(P(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +P=mean(M.P(:,end-tstudied:end),2); +histogram(P,'BinWidth',h2.BinWidth,'DisplayName',legtext); +ylabel('counts') +xlabel('P [kg\cdotm^2\cdots^{-1}]') +xlim(h2.BinLimits) +%clear P +%clear H +legend +f.PaperOrientation='landscape'; +[~, name, ~] = fileparts(M.file); +xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)]) +print(f,sprintf('%sParts_HP',name),'-dpdf','-fillpage') figure N=M.N(I+2,zpos,:); plot(M.t2d,N(:)) xlabel('t [s]') ylabel('N') title(sprintf('Evolution of number of particles in cell at z=%1.2e[m] r=%1.2e[m]',M.zgrid(zpos),M.rgrid(I+2))) -fieldstep=1; - -data.pot=M.pot(:,:,fieldstep)'; -data.pos=[M.Z(:,fieldstep),M.R(:,fieldstep)]; -data.Ez=M.Ez(:,:,fieldstep)'; -data.Er=M.Er(:,:,fieldstep)'; - -Rlarmor=(M.me*sqrt(M.VR.^2+M.VTHET.^2)/(M.qe*M.B0)); -Rmax=M.R; -Rmin=M.R-Rlarmor; -[Isoutsideuprow,Isoutsideupcol]=find(Rmax-M.rgrid(end)>0); -[Isoutsidedownrow,Isoutsidedowncol]=find(M.rgrid(1)-Rmin>0); - -f=uifigure('Name','Grid data'); - -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))) -data.pos=[M.Z(:,fieldstep),M.R(:,fieldstep)]; -data.pot=M.pot(:,:,fieldstep)'; -data.Ez=M.Ez(:,:,fieldstep)'; -data.Er=M.Er(:,:,fieldstep)'; - -sld = uislider(m,'Position',[10 30 0.6*m.Position(3) 3]); -sld.Value=fieldstep; -sld.Limits=[1 size(M.pot,3)]; - -edt = uieditfield(m,'numeric','Limits',[1 size(M.pot,3)],'Value',1); -edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20]; -edt.RoundFractionalValues='on'; - - -Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save'); -%Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause'); - - -sld.ValueChangingFcn={@updatefigdata,edt,mf,M}; -%sld.ValueChangedFcn={@updatefigdata,edt,mf,M,data}; -edt.ValueChangedFcn={@updatefigdata,sld,mf,M}; -Printbt.ButtonPushedFcn={@plotGridButtonPushed,M,sld}; -%Playbt.ButtonPushedFcn={@playButtonPushed,M,data,sld,edt}; - -PlotEspic2dgriddata(mf,M,fieldstep); - - +%dispespicFields(M) BrillouinRatio=2*M.omepe^2/M.omece^2 \ No newline at end of file diff --git a/matlab/analyse_espicparts.m b/matlab/analyse_espicparts.m index b05c8ca..afb9d0e 100644 --- a/matlab/analyse_espicparts.m +++ b/matlab/analyse_espicparts.m @@ -1,65 +1,72 @@ -file='teststable10.h5'; +file='teststablefinefinenophis.h5'; +%close all hidden; if (~exist('M','var')) M.file=file; - M.type='all'; end -M=load_espic2d(file,M,'all'); -close all hidden; +M=load_espic2d(file,M,'parts'); + tmin=2; tmax=length(M.ekin); figure plot(M.t0d(tmin:tmax),M.ekin(tmin:tmax),'o-',... M.t0d(tmin:tmax),M.epot(tmin:tmax),'d-',... M.t0d(tmin:tmax),M.etot(tmin:tmax),'h-',... M.t0d(tmin:tmax),M.eerr(tmin:tmax),'x--') legend('ekin', 'epot', 'etot','eerr') xlabel('Time [s]') ylabel('Energies [J]') grid on figure semilogy(M.t0d(tmin:tmax),abs(M.eerr(tmin:tmax)/M.etot(2)*100),'h-') xlabel('t [s]') ylabel('E_{err} %') xlim([M.t0d(tmin),M.t0d(tmax)]) grid on +% Rlarmor=(M.me*sqrt(M.VR.^2+M.VTHET.^2)/(M.qe*M.B0)); +% Rmax=M.R; +% Rmin=M.R-Rlarmor; +% [Isoutsideuprow,Isoutsideupcol]=find(Rmax-M.rgrid(end)>0); +% [Isoutsidedownrow,Isoutsidedowncol]=find(M.rgrid(1)-Rmin>0); + +%dispespicParts(M) + + +f=figure(); +tstudied=20; +legtext=sprintf("t=%2.1f - %2.1f [ns]",M.tpart(end-tstudied)*1e9,M.tpart(end)*1e9); +subplot(1,2,1) +H=M.H(:,1); +h1=histogram(H,20,'BinLimits',[min(H(:)) max(H(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +H=mean(M.H(:,end-tstudied:end),2); +h1=histogram(H,'BinWidth',h1.BinWidth,'DisplayName',legtext); +ylabel('counts') +xlabel('H [J]') +legend + +subplot(1,2,2) +P=M.P(:,1); +h2=histogram(P,20,'BinLimits',[min(P(:)) max(P(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +P=mean(M.P(:,end-tstudied:end),2); +histogram(P,'BinWidth',h2.BinWidth,'DisplayName',legtext); +ylabel('counts') +xlabel('P [kg\cdotm^2\cdots^{-1}]') +xlim(h2.BinLimits) +%clear P +%clear H +legend +f.PaperOrientation='landscape'; +[~, name, ~] = fileparts(M.file); +xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)]) +print(f,sprintf('%sParts_HP',name),'-dpdf','-fillpage') -fieldstep=1; - -Rlarmor=(M.me*sqrt(M.VR.^2+M.VTHET.^2)/(M.qe*M.B0)); -Rmax=M.R; -Rmin=M.R-Rlarmor; -[Isoutsideuprow,Isoutsideupcol]=find(Rmax-M.rgrid(end)>0); -[Isoutsidedownrow,Isoutsidedowncol]=find(M.rgrid(1)-Rmin>0); - - -f2=uifigure('Name','Particles data'); - -mf2=uipanel(f2,'Position',[5 50 f2.Position(3)-10 f2.Position(4)-55]); -mf2.AutoResizeChildren='off'; -m2=uipanel(f2,'Position',[5 5 f2.Position(3)-10 40]); - -sgtitle(mf2,sprintf('step=%d t=%0.5e s',fieldstep*M.it2,M.tpart(fieldstep))) - -sld2 = uislider(m2,'Position',[10 30 0.6*m2.Position(3) 3]); -sld2.Value=fieldstep; -sld2.Limits=[1 size(M.VR,2)]; - -edt2 = uieditfield(m2,'numeric','Limits',[1 size(M.VR,2)],'Value',1); -edt2.Position=[sld2.Position(1)+sld2.Position(3)+25 5 40 20]; -edt2.RoundFractionalValues='on'; - -Printbt2=uibutton(m2,'Position',[edt2.Position(1)+edt2.Position(3)+10 5 40 20],'Text', 'Save'); - -sld2.ValueChangingFcn={@updatefigpartdata,edt2,mf2,M}; -edt2.ValueChangedFcn={@updatefigpartdata,sld2,mf2,M}; -Printbt2.ButtonPushedFcn={@plotPartButtonPushed,M,sld2}; -PlotEspic2dpartdata(mf2,M,fieldstep); BrillouinRatio=2*M.omepe^2/M.omece^2 \ No newline at end of file diff --git a/matlab/dispespicFields.m b/matlab/dispespicFields.m index 7bd5bfa..fbdfedd 100644 --- a/matlab/dispespicFields.m +++ b/matlab/dispespicFields.m @@ -1,124 +1,124 @@ function dispespicFields(M) %UNTITLED Summary of this function goes here % Detailed explanation goes here fieldstep=1; f=uifigure('Name','Grid data'); 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.pot,3)]; edt = uieditfield(m,'numeric','Limits',[1 size(M.pot,3)],'Value',1); edt.Position=[sld.Position(1)+sld.Position(3)+25 5 40 20]; edt.RoundFractionalValues='on'; Printbt=uibutton(m,'Position',[edt.Position(1)+edt.Position(3)+10 5 40 20],'Text', 'Save'); %Playbt=uibutton(m,'Position',[Printbt.Position(1)+Printbt.Position(3)+10 5 40 20],'Text', 'Play/Pause'); sld.ValueChangingFcn={@updatefigdata,edt,mf}; edt.ValueChangedFcn={@updatefigdata,sld,mf}; Printbt.ButtonPushedFcn={@plotGridButtonPushed}; PlotEspic2dgriddata(mf,M,fieldstep); 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); pcolor(ax1,M.zgrid,M.rgrid,M.N(:,:,fieldstep)); 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'; +c.Label.String= 'n[m^{-3}]'; view(ax1,2) %set(ax1,'colorscale','log') ax2=subplot(2,2,2,'Parent',fig); contourf(ax2,M.zgrid,M.rgrid,M.pot(:,:,fieldstep)'); %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') ax3=subplot(2,2,3,'Parent',fig); contourf(ax3,M.zgrid,M.rgrid,M.Er(:,:,fieldstep)') xlabel(ax3,'Z [m]') ylabel(ax3,'R [m]') c = colorbar(ax3); c.Label.String= 'E_r [V/m]'; title(ax3,'Radial Electric field') grid(ax3, 'on') ax4=subplot(2,2,4,'Parent',fig); contourf(ax4,M.zgrid,M.rgrid,M.Ez(:,:,fieldstep)') xlabel(ax4,'Z [m]') ylabel(ax4,'R [m]') c = colorbar(ax4); c.Label.String= 'E_z [V/m]'; title(ax4,'Axial Electric field') grid(ax4, 'on') 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; else fieldstep=floor(event.Value); end Othercontrol.Value=fieldstep; sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,double((fieldstep-1)*M.it1)*M.dt)) %% update Position histogram ax1=fig.Children(end); ax1.Children.CData=M.N(:,:,fieldstep); view(ax1,2) %% update ES Potential fig.Children(end-2).Children(1).ZData=M.pot(:,:,fieldstep)'; %% update Radial electric field fig.Children(end-4).Children(1).ZData=M.Er(:,:,fieldstep)'; %% update Axial electric field fig.Children(end-6).Children(1).ZData=M.Ez(:,:,fieldstep)'; drawnow limitrate end end diff --git a/matlab/dispespicParts.m b/matlab/dispespicParts.m index 83d7eeb..4b41db5 100644 --- a/matlab/dispespicParts.m +++ b/matlab/dispespicParts.m @@ -1,145 +1,138 @@ function dispespicParts(M) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here fieldstep=1; f2=uifigure('Name','Particles data'); mf2=uipanel(f2,'Position',[5 50 f2.Position(3)-10 f2.Position(4)-55]); mf2.AutoResizeChildren='off'; m2=uipanel(f2,'Position',[5 5 f2.Position(3)-10 40]); sgtitle(mf2,sprintf('step=%d t=%0.5e s',fieldstep*M.it2,M.tpart(fieldstep))) sld2 = uislider(m2,'Position',[10 30 0.6*m2.Position(3) 3]); sld2.Value=fieldstep; sld2.Limits=[1 size(M.VR,2)]; edt2 = uieditfield(m2,'numeric','Limits',[1 size(M.VR,2)],'Value',1); edt2.Position=[sld2.Position(1)+sld2.Position(3)+25 5 40 20]; edt2.RoundFractionalValues='on'; Printbt2=uibutton(m2,'Position',[edt2.Position(1)+edt2.Position(3)+10 5 40 20],'Text', 'Save'); sld2.ValueChangingFcn={@updatefigpartdata,edt2,mf2}; edt2.ValueChangedFcn={@updatefigpartdata,sld2,mf2}; Printbt2.ButtonPushedFcn={@plotPartButtonPushed}; PlotEspic2dpartdata(mf2,M,fieldstep); -function PlotEspic2dpartdata(fig,M,fieldstep) -%PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep -sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it1,M.t2d(fieldstep))) - -ax1=subplot(3,2,1,'Parent',fig); -plot(ax1,M.Z(:,fieldstep),M.R(:,fieldstep),'.'); -xlim(ax1,[M.zgrid(1) M.zgrid(end)]) -ylim(ax1,[M.rgrid(1) M.rgrid(end)]) -xlabel(ax1,'Z [m]') -ylabel(ax1,'R [m]') -title(ax1,'Position') -grid(ax1, 'on') - - -ax2=subplot(3,2,2,'Parent',fig); -plot(ax2,M.VZ(:,fieldstep),M.VR(:,fieldstep),'.'); -xlabel(ax2,'VZ [m/s]') -ylabel(ax2,'VR [m/s]') -axis(ax2, 'equal') -grid(ax2, 'on') - - -ax3=subplot(3,2,3,'Parent',fig); -plot(ax3,M.Z(:,fieldstep),M.VZ(:,fieldstep),'.'); -xlim(ax3,[M.zgrid(1) M.zgrid(end)]) -xlabel(ax3,'Z [m]') -ylabel(ax3,'VZ [m/s]') -grid(ax3, 'on') - -ax4=subplot(3,2,4,'Parent',fig); -plot(ax4,M.VZ(:,fieldstep),M.R(:,fieldstep),'.') -ylabel(ax4,'R [m]') -xlabel(ax4,'VZ [m/s]') -ylim(ax4,[M.rgrid(1) M.rgrid(end)]) -grid(ax4, 'on') - -ax5=subplot(3,2,5,'Parent',fig); -plot(ax5,M.Z(:,fieldstep),M.VR(:,fieldstep),'.'); -xlim(ax5,[M.zgrid(1) M.zgrid(end)]) -xlabel(ax5,'Z [m]') -ylabel(ax5,'VR [m/s]') -grid(ax5, 'on') - -ax6=subplot(3,2,6,'Parent',fig); -plot(ax6,M.VR(:,fieldstep),M.R(:,fieldstep),'.'); -ylim(ax6,[M.rgrid(1) M.rgrid(end)]) -ylabel(ax6,'R [m]') -xlabel(ax6,'VR [m/s]') -grid(ax6, 'on') -end + function PlotEspic2dpartdata(fig,M,fieldstep) + %PlotEspic2dgriddata Plot the 2d data of espic2d at time step fieldstep + sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it2,M.tpart(fieldstep))) + + ax1=subplot(3,2,1,'Parent',fig); + plot(ax1,M.Z(:,fieldstep),M.R(:,fieldstep),'.'); + xlim(ax1,[M.zgrid(1) M.zgrid(end)]) + ylim(ax1,[M.rgrid(1) M.rgrid(end)]) + xlabel(ax1,'Z [m]') + ylabel(ax1,'R [m]') + title(ax1,'Position') + grid(ax1, 'on') + + + ax2=subplot(3,2,2,'Parent',fig); + plot(ax2,M.VZ(:,fieldstep),M.VR(:,fieldstep),'.'); + xlabel(ax2,'VZ [m/s]') + ylabel(ax2,'VR [m/s]') + axis(ax2, 'equal') + grid(ax2, 'on') + + + ax3=subplot(3,2,3,'Parent',fig); + plot(ax3,M.Z(:,fieldstep),M.VZ(:,fieldstep),'.'); + xlim(ax3,[M.zgrid(1) M.zgrid(end)]) + xlabel(ax3,'Z [m]') + ylabel(ax3,'VZ [m/s]') + grid(ax3, 'on') + + ax4=subplot(3,2,4,'Parent',fig); + plot(ax4,M.VZ(:,fieldstep),M.R(:,fieldstep),'.') + ylabel(ax4,'R [m]') + xlabel(ax4,'VZ [m/s]') + ylim(ax4,[M.rgrid(1) M.rgrid(end)]) + grid(ax4, 'on') + + ax5=subplot(3,2,5,'Parent',fig); + plot(ax5,M.Z(:,fieldstep),M.VR(:,fieldstep),'.'); + xlim(ax5,[M.zgrid(1) M.zgrid(end)]) + xlabel(ax5,'Z [m]') + ylabel(ax5,'VR [m/s]') + grid(ax5, 'on') + + ax6=subplot(3,2,6,'Parent',fig); + plot(ax6,M.VR(:,fieldstep),M.R(:,fieldstep),'.'); + ylim(ax6,[M.rgrid(1) M.rgrid(end)]) + ylabel(ax6,'R [m]') + xlabel(ax6,'VR [m/s]') + grid(ax6, 'on') + end -function updatefigpartdata(control, event, Othercontrol, fig) - - - if strcmp(event.EventName,'ValueChanged') - fieldstep=floor(control.Value); - control.Value=fieldstep; - else - fieldstep=floor(event.Value); + function updatefigpartdata(control, event, Othercontrol, fig) + + + if strcmp(event.EventName,'ValueChanged') + fieldstep=floor(control.Value); + control.Value=fieldstep; + else + fieldstep=floor(event.Value); + end + Othercontrol.Value=fieldstep; + + sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it2,double((fieldstep-1)*M.it2)*M.dt)) + + %% update Position plot + ax1=fig.Children(end); + ax1.Children.XData=M.Z(:,fieldstep); + ax1.Children.YData=M.R(:,fieldstep); + + + view(ax1,2) + %% update VZ VPERP + fig.Children(end-1).Children(1).XData=M.VZ(:,fieldstep); + fig.Children(end-1).Children(1).YData=M.VR(:,fieldstep); + axis(fig.Children(end-1),'equal') + + %% update Z VZ + fig.Children(end-2).Children(1).XData=M.Z(:,fieldstep); + fig.Children(end-2).Children(1).YData=M.VZ(:,fieldstep); + %% update VZ VR + fig.Children(end-3).Children(1).XData=M.VZ(:,fieldstep); + fig.Children(end-3).Children(1).YData=M.R(:,fieldstep); + + %% update R VR + fig.Children(end-4).Children(1).XData=M.Z(:,fieldstep); + fig.Children(end-4).Children(1).YData=M.VR(:,fieldstep); + + %% update Z VR + fig.Children(end-5).Children(1).XData=M.VR(:,fieldstep); + fig.Children(end-5).Children(1).YData=M.R(:,fieldstep); + drawnow limitrate end - Othercontrol.Value=fieldstep; - - sgtitle(fig,sprintf('step=%d t=%0.5e s',(fieldstep-1)*M.it2,double((fieldstep-1)*M.it2)*M.dt)) - %data.pos=[M.Z(:,fieldstep),M.R(:,fieldstep)]; - data.pot=M.pot(:,:,fieldstep)'; - data.Ez=M.Ez(:,:,fieldstep)'; - data.Er=M.Er(:,:,fieldstep)'; - - %% update Position plot - ax1=fig.Children(end); - ax1.Children.XData=M.Z(:,fieldstep); - ax1.Children.YData=M.R(:,fieldstep); - - - view(ax1,2) - %% update VZ VPERP - fig.Children(end-1).Children(1).XData=M.VZ(:,fieldstep); - fig.Children(end-1).Children(1).YData=M.VR(:,fieldstep); - axis(fig.Children(end-1),'equal') - - %% update Z VZ - fig.Children(end-2).Children(1).XData=M.Z(:,fieldstep); - fig.Children(end-2).Children(1).YData=M.VZ(:,fieldstep); - %% update VZ VR - fig.Children(end-3).Children(1).XData=M.VZ(:,fieldstep); - fig.Children(end-3).Children(1).YData=M.R(:,fieldstep); - - %% update R VR - fig.Children(end-4).Children(1).XData=M.Z(:,fieldstep); - fig.Children(end-4).Children(1).YData=M.VR(:,fieldstep); - - %% update Z VR - fig.Children(end-5).Children(1).XData=M.VR(:,fieldstep); - fig.Children(end-5).Children(1).YData=M.R(:,fieldstep); - drawnow limitrate - - - -end -function plotPartButtonPushed(btn,ax) -%UNTITLED2 Summary of this function goes here -% Detailed explanation goes here -f=figure(); -PlotEspic2dpartdata(f,M,sld.Value); -f.PaperOrientation='portrait'; -[~, name, ~] = fileparts(M.file); -print(f,sprintf('%sParts%d',name,sld.Value),'-dpdf','-fillpage') -end + function plotPartButtonPushed(btn,ax) + %UNTITLED2 Summary of this function goes here + % Detailed explanation goes here + f=figure(); + PlotEspic2dpartdata(f,M,sld.Value); + f.PaperOrientation='portrait'; + [~, name, ~] = fileparts(M.file); + print(f,sprintf('%sParts%d',name,sld.Value),'-dpdf','-fillpage') + end end diff --git a/matlab/load_espic2d.m b/matlab/load_espic2d.m index f6df6be..9cacbad 100644 --- a/matlab/load_espic2d.m +++ b/matlab/load_espic2d.m @@ -1,136 +1,164 @@ -function Data = load_espic2d( file, M, type) +function Data = load_espic2d( file, Data, type) % Load Simulation data from espic2d simulation if(nargin==1) type='all'; else filedata=dir(file); + if (isempty(filedata)) + error("File: ""%s"" doesn't exist",file) + end timestamp=filedata.date; try - if (strcmp(M.file,file) && ~max(timestamp > M.timestamp) && strcmp(M.type,'all')) - Data=M; + if (strcmp(Data.file,file) && ~max(timestamp > Data.timestamp) && (strcmp(Data.type,type)||strcmp(Data.type,'all')) ) return end catch end if(nargin~=3) type='all'; end end - Data.type=type; Data.file=file; % Physical constants Data.vlight=299792458; Data.qe=1.60217662E-19; Data.me=9.109383E-31; Data.eps_0=8.85418781762E-12; Data.job_time = hdf5read(file,'/data/input.00/','job_time'); Data.extra_time = hdf5read(file,'/data/input.00/','extra_time'); Data.dt = hdf5read(file,'/data/input.00/','dt'); Data.tmax= hdf5read(file,'/data/input.00/','tmax'); Data.nrun = hdf5read(file,'/data/input.00/','nrun'); -Data.nlres = hdf5read(file,'/data/input.00/','nlres'); -Data.nlsave = hdf5read(file,'/data/input.00/','nlsave'); -Data.nlres = hdf5read(file,'/data/input.00/','nlres'); +Data.nlres = strcmp(h5readatt(file,'/data/input.00/','nlres'),'y'); +Data.nlsave = strcmp(h5readatt(file,'/data/input.00/','nlsave'),'y'); +Data.nlclassical =strcmp(h5readatt(file,'/data/input.00/','nlclassical'),'y'); +Data.nlPhis =strcmp(h5readatt(file,'/data/input.00/','nlPhis'),'y'); Data.nz = hdf5read(file,'/data/input.00/','nz'); +Data.nnr = hdf5read(file,'/data/input.00/nnr'); Data.lz = hdf5read(file,'/data/input.00/lz'); +Data.qsim = hdf5read(file,'/data/input.00/','qsim'); +Data.msim = hdf5read(file,'/data/input.00/','msim'); Data.nplasma = hdf5read(file,'/data/input.00/','nplasma'); Data.potinn = hdf5read(file,'/data/input.00/','potinn'); Data.potout = hdf5read(file,'/data/input.00/','potout'); Data.B0 = hdf5read(file,'/data/input.00/','B0'); Data.Rcurv = hdf5read(file,'/data/input.00/','Rcurv'); Data.width = hdf5read(file,'/data/input.00/','width'); Data.n0 = hdf5read(file,'/data/input.00/','n0'); Data.temp = hdf5read(file,'/data/input.00/','temp'); try Data.it0 = h5readatt(file,'/data/input.00/','it0d'); Data.it1 = h5readatt(file,'/data/input.00/','it2d'); Data.it2 = h5readatt(file,'/data/input.00/','itparts'); catch Data.it0 = h5readatt(file,'/data/input.00/','it0'); Data.it1 = h5readatt(file,'/data/input.00/','it1'); Data.it1 = h5readatt(file,'/data/input.00/','it2'); end Data.omepe=sqrt(abs(Data.n0)*Data.qe^2/(Data.me*Data.eps_0)); Data.omece=Data.qe*Data.B0/Data.me; % Normalizations Data.tnorm=1/Data.omepe; Data.rnorm=Data.vlight*Data.tnorm; Data.bnorm=Data.B0; Data.enorm=Data.vlight*Data.bnorm; Data.phinorm=Data.enorm*Data.rnorm; Data.vnorm=Data.vlight; % Grid data Data.rgrid= hdf5read(file, '/data/var1d/rgrid')*Data.rnorm; Data.zgrid= hdf5read(file, '/data/var1d/zgrid')*Data.rnorm; Data.dz=(Data.zgrid(end)-Data.zgrid(1))/double(Data.nz); Br = hdf5read(file,'/data/fields/Br')*Data.bnorm; Data.Br= reshape(Br,length(Data.zgrid),length(Data.rgrid)); Bz = hdf5read(file,'/data/fields/Bz')*Data.bnorm; Data.Bz= reshape(Bz,length(Data.zgrid),length(Data.rgrid)); Data.B=sqrt(Data.Bz.^2+Data.Br.^2); clear Br Bz Data.femorder = hdf5read(file,'/data/input.00/femorder'); Data.ngauss = hdf5read(file,'/data/input.00/ngauss'); -Data.nnr = hdf5read(file,'/data/input.00/femorder'); +Data.plasmadim = hdf5read(file,'/data/input.00/plasmadim'); Data.radii = hdf5read(file,'/data/input.00/radii'); Data.epot = hdf5read(file,'/data/var0d/epot'); Data.ekin = hdf5read(file,'/data/var0d/ekin'); Data.etot = hdf5read(file,'/data/var0d/etot'); Data.eerr = Data.etot-Data.etot(2); if(strcmp(type,'all') || strcmp(type,'fields')) pot = hdf5read(file,'/data/fields/pot')*Data.phinorm; Data.pot= reshape(pot,length(Data.zgrid),length(Data.rgrid),size(pot,2)); clear pot Er = hdf5read(file,'/data/fields/Er')*Data.enorm; Data.Er= reshape(Er,length(Data.zgrid),length(Data.rgrid),size(Er,2)); clear Er Ez = hdf5read(file,'/data/fields/Ez')*Data.enorm; Data.Ez= reshape(Ez,length(Data.zgrid),length(Data.rgrid),size(Ez,2)); clear Ez N = hdf5read(file,'/data/fields/partdensity'); Data.N=reshape(N,length(Data.zgrid)-1,length(Data.rgrid)-1,[]); Data.N=permute(Data.N,[2,1,3]); + Data.Vcell=(Data.zgrid(2:end)-Data.zgrid(1:end-1))*((Data.rgrid(2:end)-Data.rgrid(1:end-1)).*(Data.rgrid(2:end)+Data.rgrid(1:end-1))*pi)'; + for i=1:size(Data.N,3) + Data.N(:,:,i)=abs(Data.qsim/Data.qe)*Data.N(:,:,i)./Data.Vcell'; + end Data.N=padarray(Data.N,[1,1],0,'post'); + clear Data.Vcell clear N Data.t2d=Data.dt*(0:size(Data.N,3)-1)*double(Data.it1); + if(strcmp(Data.type,'parts')) + Data.type='all'; + end end if(strcmp(type,'all') || strcmp(type,'parts')) Data.R = hdf5read(file,'/data/part/R')*Data.rnorm; Data.Z = hdf5read(file,'/data/part/Z')*Data.rnorm; + Data.partepot = sign(Data.qsim)*hdf5read(file,'/data/part/pot')*Data.phinorm*Data.qe; Data.UR = hdf5read(file,'/data/part/UR'); Data.UZ = hdf5read(file,'/data/part/UZ'); Data.UTHET= hdf5read(file,'/data/part/UTHET'); - Data.GAMM = sqrt(1+Data.UR.^2+Data.UZ.^2+Data.UTHET.^2); + if(Data.nlclassical) + Data.GAMM=ones(size(Data.R)); + else + Data.GAMM = sqrt(1+Data.UR.^2+Data.UZ.^2+Data.UTHET.^2); + end Data.VR = Data.UR./Data.GAMM*Data.vnorm; clear Data.UR Data.VZ = Data.UZ./Data.GAMM*Data.vnorm; clear Data.UZ Data.VTHET = Data.UTHET./Data.GAMM*Data.vnorm; clear Data.UTHET Data.VPERP = sqrt(Data.VR.*Data.VR+Data.VTHET.*Data.VTHET); + Data.H=0.5*Data.me*(Data.VR.^2+Data.VTHET.^2+Data.VZ.^2)+Data.partepot; + Data.P=Data.R.*(Data.VTHET*Data.me+sign(Data.qsim)*Data.qe*Athet(Data.R,Data.Z)); Data.tpart=Data.dt*(0:size(Data.R,2)-1)*double(Data.it2); + if(strcmp(Data.type,'fields')) + Data.type='all'; + end end Data.t=Data.dt.*double(0:length(Data.epot)-1); Data.t0d=Data.dt.*double(0:length(Data.epot)-1); filedata=dir(file); Data.timestamp=filedata.date; disp('M reloaded') + function Atheta=Athet(R,Z) + Atheta=0.5*Data.B0*(R-Data.width/pi*(Data.Rcurv-1)/(Data.Rcurv+1)... + .*besseli(1,2*pi*R/Data.width).*cos(2*pi*Z/Data.width)); + end + end diff --git a/src/basic_mod.f90 b/src/basic_mod.f90 index afb8487..df6f469 100644 --- a/src/basic_mod.f90 +++ b/src/basic_mod.f90 @@ -1,303 +1,304 @@ MODULE basic ! USE hashtable USE constants USE bsplines USE mumps_bsplines USE futils USE mpihelper IMPLICIT NONE ! ! Basic module for time dependent problems ! CHARACTER(len=80) :: label1, label2, label3, label4 ! ! BASIC Namelist ! LOGICAL :: nlres = .FALSE. ! Restart flag LOGICAL :: nlsave = .TRUE. ! Checkpoint (save) flag LOGICAL :: newres=.FALSE. ! New result HDF5 file LOGICAL :: nlxg=.TRUE. ! Show graphical interface Xgrafix INTEGER :: nrun=1 ! Number of time steps to run DOUBLE PRECISION :: job_time=3600.0 ! Time allocated to this job in seconds DOUBLE PRECISION :: tmax=100000.0 ! Maximum simulation time DOUBLE PRECISION :: extra_time=60.0 ! Extra time allocated DOUBLE PRECISION :: dt=1 ! Time step DOUBLE PRECISION :: time=0 ! Current simulation time (Init from restart file) ! ! Other basic global vars and arrays ! INTEGER :: jobnum ! Job number INTEGER :: step ! Calculation step of this run INTEGER :: cstep=0 ! Current step number (Init from restart file) LOGICAL :: nlend ! Signal end of run INTEGER :: ierr ! Integer used for MPI - INTEGER :: it0d, it2d, itparts ! Number of iterations between 0d, 2d, particles, values writes to hdf5 - INTEGER :: ittext=1 ! Number of iterations between text outputs in the console + INTEGER :: it0d=1, it2d=1, itparts=1 ! Number of iterations between 0d, 2d, particles, values writes to hdf5 + INTEGER :: ittext=1 ! Number of iterations between text outputs in the console + INTEGER :: itrestart=10000 ! Number of iterations between save of restart.h5 file INTEGER :: itgraph ! Number of iterations between graphical interface updates INTEGER :: mpirank ! MPIrank of the current processus INTEGER :: mpisize ! Size of the MPI_COMM_WORLD communicator INTEGER :: rightproc, leftproc ! Rank of next and previous processor in the z decomposition ! ! List of logical file units INTEGER :: lu_in = 90 ! File duplicated from STDIN INTEGER :: lu_stop = 91 ! stop file, see subroutine TESEND ! ! HDF5 file CHARACTER(len=64) :: resfile = "results.h5" ! Main result file CHARACTER(len=64) :: rstfile = "restart.h5" ! Restart file INTEGER :: fidres ! FID for resfile INTEGER :: fidrst ! FID for restart file TYPE(BUFFER_TYPE) :: hbuf0 ! Hashtable for 0d var ! ! Plasma parameters LOGICAL :: nlPhis=.TRUE. ! Calculate self consistent electric field flag LOGICAL :: nlclassical=.FALSE. ! If true, solves the equation of motion classicaly LOGICAL :: nlperiod(2)=(/.false.,.false./) ! Set periodic splines on or off INTEGER :: nplasma ! Number of superparticles INTEGER :: nblock ! Number of intervals in Z for stable distribution initialisation DOUBLE PRECISION :: potinn,potout ! Potential at inner and outer metalic walls DOUBLE PRECISION :: B0 ! Max magnitude of magnetic field DOUBLE PRECISION, allocatable :: Bz(:), Br(:) ! Magnetic field components DOUBLE PRECISION, allocatable :: Athet(:) ! Magnetic potential TYPE(spline2d) :: splrz ! Spline at r and z DOUBLE PRECISION, allocatable :: Ez(:), Er(:) ! Electric field components DOUBLE PRECISION, allocatable :: pot(:) ! Electro static potential DOUBLE PRECISION :: radii(4) ! Inner and outer radius of cylinder and radii of fine mesh region DOUBLE PRECISION :: plasmadim(4) ! Zmin Zmax Rmin Rmax values for plasma particle loading INTEGER :: distribtype=1 ! Type of distribution function used to load the particles !1: gaussian, 2: Stable as defined in 4.95 of Davidson DOUBLE PRECISION :: H0=0, P0=0 ! Initial value of Hamiltonian and canonical angular momentum ! for distribution 2 DOUBLE PRECISION :: lz(2) ! Length of cylinder in z direction DOUBLE PRECISION :: n0 ! Average density of physical plasma DOUBLE PRECISION, ALLOCATABLE:: partdensity(:) ! Number of particles per gridcell computed every it2d INTEGER :: nz, nr, nnr(3) ! Number of intervals in z and r DOUBLE PRECISION :: dz, dr(3) ! Cell size in z and r for each region DOUBLE PRECISION, allocatable :: zgrid(:), rgrid(:) ! Nodes positions DOUBLE PRECISION :: bnorm,enorm,vnorm,tnorm,rnorm,phinorm ! Normalization constants DOUBLE PRECISION :: qsim ! Charge of superparticles DOUBLE PRECISION :: msim ! Mass of superparticles INTEGER :: femorder(2) ! FEM order INTEGER :: ngauss(2) ! Number of gauss points LOGICAL :: nlppform =.TRUE. ! Argument of set_spline INTEGER, SAVE :: nrank(2) ! Number of splines in both directions DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE:: rhs ! right hand side of the poisson equation solver DOUBLE PRECISION :: omegac ! Cyclotronic frequency DOUBLE PRECISION :: omegap ! Plasma frequency DOUBLE PRECISION :: V ! Normalised volume DOUBLE PRECISION :: temp ! Initial temperature of plasma DOUBLE PRECISION :: Rcurv ! Magnetic field curvature coefficient DOUBLE PRECISION :: Width ! Distance between two magnetic mirrors DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: Zbounds ! Index of bounds for local processus in Z direction TYPE(BASICDATA) :: bdata CONTAINS ! !================================================================================ SUBROUTINE basic_data ! ! Define basic data ! use mpihelper IMPLICIT NONE ! ! Local vars and arrays CHARACTER(len=256) :: line, inputfilename ! NAMELIST /BASIC/ job_time, extra_time, nrun, tmax, dt, nlres, nlsave, newres, nlxg, & & nplasma, potinn, potout, B0, lz, n0, nz, nnr, femorder, ngauss, & & nlppform, plasmadim, radii, temp, Rcurv, width, it0d, it2d, itparts, ittext, & & resfile, itgraph, nlPhis, distribtype, nblock, nlclassical, H0, P0 !________________________________________________________________________________ ! 1. Process Standard Input File ! IF(mpirank .eq. 0) THEN IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN WRITE(*,*)'ERROR, ONE COMMAND-LINE ARGUMENT REQUIRED, STOPPING' STOP ENDIF CALL GET_COMMAND_ARGUMENT(1,inputfilename) OPEN(UNIT=lu_in,FILE=trim(inputfilename)) !DO ! READ(*,'(a)', END=110) line ! WRITE(lu_in, '(a)') TRIM(line) !END DO !110 CONTINUE !REWIND(lu_in) !________________________________________________________________________________ ! 1. Label the run ! READ(lu_in,'(a)') label1 READ(lu_in,'(a)') label2 READ(lu_in,'(a)') label3 READ(lu_in,'(a)') label4 ! WRITE(*,'(12x,a/)') label1(1:len_trim(label1)) WRITE(*,'(12x,a/)') label2(1:len_trim(label2)) WRITE(*,'(12x,a/)') label3(1:len_trim(label3)) WRITE(*,'(12x,a/)') label4(1:len_trim(label4)) !________________________________________________________________________________ ! 2. Read in basic data specific to run ! READ(lu_in,basic) WRITE(*,basic) bdata=BASICDATA( nlres, nlsave, newres, nlxg, nlppform, nlPhis, nlclassical, & & nplasma, nz, it0d, it2d, itparts, itgraph, distribtype, nblock, & & nrun, job_time, extra_time, tmax, dt, potinn, potout, B0, n0, temp, & & Rcurv, width, H0, P0, femorder, ngauss, nnr, lz, radii, plasmadim, resfile ) END IF CALL init_mpitypes ! initialize all mpi types that will be needed in the simulation CALL MPI_Bcast(bdata, 1, basicdata_type, 0, MPI_COMM_WORLD, ierr) CALL readbdata ! END SUBROUTINE basic_data !================================================================================ SUBROUTINE daytim(str) ! ! Print date and time ! IMPLICIT NONE ! CHARACTER(len=*), INTENT(in) :: str ! ! Local vars and arrays CHARACTER(len=16) :: d, t, dat, functime !________________________________________________________________________________ ! CALL DATE_AND_TIME(d,t) dat=d(7:8) // '/' // d(5:6) // '/' // d(1:4) functime=t(1:2) // ':' // t(3:4) // ':' // t(5:10) WRITE(*,'(a,1x,a,1x,a)') str, dat(1:10), functime(1:12) ! END SUBROUTINE daytim !================================================================================ SUBROUTINE timera(cntrl, str, eltime) ! ! Timers (cntrl=0/1 to Init/Update) ! IMPLICIT NONE INTEGER, INTENT(in) :: cntrl CHARACTER(len=*), INTENT(in) :: str DOUBLE PRECISION, OPTIONAL, INTENT(out) :: eltime ! INTEGER, PARAMETER :: ncmax=128 INTEGER, SAVE :: icall=0, nc=0 DOUBLE PRECISION, SAVE :: startt0=0.0 CHARACTER(len=16), SAVE :: which(ncmax) INTEGER :: lstr, found, i DOUBLE PRECISION :: seconds DOUBLE PRECISION, DIMENSION(ncmax), SAVE :: startt = 0.0, endt = 0.0 !________________________________________________________________________________ IF( icall .EQ. 0 ) THEN icall = icall+1 startt0 = seconds() END IF lstr = LEN_TRIM(str) IF( lstr .GT. 0 ) found = loc(str) !________________________________________________________________________________ ! SELECT CASE (cntrl) ! CASE(-1) ! Current wall time IF( PRESENT(eltime) ) THEN eltime = seconds() - startt0 ELSE WRITE(*,'(/a,a,1pe10.3/)') "++ ", ' Wall time used so far = ', seconds() - startt0 END IF ! CASE(0) ! Init Timer IF( found .EQ. 0 ) THEN ! Called for the 1st time for 'str' nc = nc+1 which(nc) = str(1:lstr) found = nc END IF startt(found) = seconds() ! CASE(1) ! Update timer endt(found) = seconds() - startt(found) IF( PRESENT(eltime) ) THEN eltime = endt(found) ELSE WRITE(*,'(/a,a,1pe10.3/)') "++ "//str, ' wall clock time = ', endt(found) END IF ! CASE(2) ! Update and reset timer endt(found) = endt(found) + seconds() - startt(found) startt(found) = seconds() IF( PRESENT(eltime) ) THEN eltime = endt(found) END IF ! CASE(9) ! Display all timers IF( nc .GT. 0 ) THEN WRITE(*,'(a)') "Timer Summary" WRITE(*,'(a)') "=============" DO i=1,nc WRITE(*,'(a20,2x,2(1pe12.3))') TRIM(which(i))//":", endt(i) END DO END IF ! END SELECT ! CONTAINS INTEGER FUNCTION loc(str) CHARACTER(len=*), INTENT(in) :: str INTEGER :: i, ind loc = 0 DO i=1,nc ind = INDEX(which(i), str(1:lstr)) IF( ind .GT. 0 .AND. LEN_TRIM(which(i)) .EQ. lstr ) THEN loc = i EXIT END IF END DO END FUNCTION loc END SUBROUTINE timera !================================================================================ SUBROUTINE readbdata nlres = bdata%nlres nlsave = bdata%nlsave newres = bdata%newres nlxg = bdata%nlxg nlppform = bdata%nlppform nlPhis = bdata%nlPhis nlclassical = bdata%nlclassical nplasma = bdata%nplasma nz = bdata%nz it0d = bdata%it0d it2d = bdata%it2d itparts = bdata%itparts itgraph = bdata%itgraph distribtype = bdata%distribtype nblock = bdata%nblock nrun = bdata%nrun job_time = bdata%job_time extra_time = bdata%extra_time tmax = bdata%tmax dt = bdata%dt potinn = bdata%potinn potout = bdata%potout B0 = bdata%B0 n0 = bdata%n0 temp = bdata%temp Rcurv = bdata%Rcurv width = bdata%width H0 = bdata%H0 P0 = bdata%P0 femorder = bdata%femorder ngauss = bdata%ngauss nnr = bdata%nnr lz = bdata%lz radii = bdata%radii plasmadim = bdata%plasmadim resfile = bdata%resfile END SUBROUTINE readbdata !================================================================================ END MODULE basic diff --git a/src/beam_mod.f90 b/src/beam_mod.f90 index 448d3eb..65603a3 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,1136 +1,1171 @@ MODULE beam ! USE constants USE mpi USE mpihelper USE basic, ONLY: mpirank, mpisize IMPLICIT NONE TYPE particles INTEGER :: Nptot INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex, Zindex ! Indices of the particles on the grid DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: & & R,Z, & & WZ, WR, & & pot, Er, Ez DOUBLE PRECISION, DIMENSION(:),POINTER:: & & UR, URold, & & UTHET, UTHETold, & & UZ, UZold, & & Gamma, Gammaold END TYPE particles ! TYPE(particles) :: parts ! Storage for all the particles SAVE :: parts TYPE(particle), DIMENSION(:), ALLOCATABLE:: rrecvpartbuff, lrecvpartbuff, rsendpartbuff, lsendpartbuff ! buffer to send and receive particle from left and right processes ! Diagnostics (scalars) DOUBLE PRECISION :: ekin ! Total kinetic energz (J) DOUBLE PRECISION :: epot ! Total potential energy (J) DOUBLE PRECISION :: etot, etot0 ! Current and Initial total energy (J) - INTEGER, DIMENSION(3), SAVE :: ireducerequest=MPI_REQUEST_NULL ! INTEGER, DIMENSION(:), ALLOCATABLE :: Nplocs_all !F.B. get local numbers of particles + LOGICAL:: collected=.false. !Stores if the particles data have been collected to root process this timestep ! CONTAINS SUBROUTINE creat_parts(p, nparts) TYPE(particles) :: p INTEGER, INTENT(in) :: nparts IF (.NOT. ALLOCATED(p%Z) ) THEN p%Nptot = nparts ALLOCATE(p%Z(nparts)) ALLOCATE(p%R(nparts)) ALLOCATE(p%WZ(nparts)) ALLOCATE(p%WR(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%pot(nparts)) ALLOCATE(p%Er(nparts)) ALLOCATE(p%Ez(nparts)) ALLOCATE(p%GAMMAold(nparts)) parts%URold=0 parts%UZold=0 parts%UTHETold=0 parts%rindex=0 parts%zindex=0 parts%WR=0 parts%WZ=0 parts%UR=0 parts%UZ=0 parts%UTHET=0 parts%Z=0 parts%R=0 parts%Gamma=1 parts%Er=0 parts%Ez=0 parts%pot=0 parts%gammaold=1 END IF END SUBROUTINE creat_parts !______________________________________________________________________________ SUBROUTINE load_parts ! Loads the particles at the beginning of the simulation and create the parts variable if necessary USE basic, ONLY: nplasma, mpirank, ierr, distribtype, mpisize, nlclassical USE mpi DOUBLE PRECISION, 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 and uniform in R - CALL loaduniformRZ(VR,VZ,VTHET) + CALL loaduniformRZ(VR, VZ, VTHET) CASE(2) !Stable distribution from Davidson 4.95 p.119 - CALL loadDavidson(VR,VZ,VTHET) + CALL loadDavidson(VR, VZ, VTHET, lodunir) CASE(3) !Stable distribution from Davidson 4.95 p.119 but with distribution in R as 1/R**2 - CALL loadDavidson(VR,VZ,VTHET) + CALL loadDavidson(VR, VZ, VTHET, lodinvr) CASE DEFAULT IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of distribution:", distribtype CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END SELECT IF(nlclassical) THEN parts%Gamma=1 ELSE parts%Gamma=sqrt(1/(1-VR**2+VZ**2+VTHET**2)) END IF parts%UR=parts%Gamma*VR parts%UZ=parts%Gamma*VZ parts%UTHET=parts%Gamma*VTHET DEALLOCATE(VZ, VR, VTHET) CALL localisation IF(mpisize .gt. 1) THEN CALL distribpartsonprocs END IF END SUBROUTINE load_parts !________________________________________________________________________________ SUBROUTINE bound USE basic, ONLY: zgrid, nz, Zbounds, cstep, mpirank INTEGER :: i, rsendnbparts=0, lsendnbparts=0 INTEGER, DIMENSION(parts%Nptot) :: sendhole sendhole=0 lsendnbparts=0 rsendnbparts=0 IF (parts%Nptot .NE. 0) THEN ! Boundary condition at z direction !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot IF(cstep .ne. 0) THEN IF (parts%Z(i) .ge. Zbounds(mpirank+1)) THEN !$OMP CRITICAL (nbparts) rsendnbparts=rsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=i !$OMP END CRITICAL (nbparts) ELSE IF (parts%Z(i) .lt. Zbounds(mpirank)) THEN !$OMP CRITICAL (nbparts) lsendnbparts=lsendnbparts+1 sendhole(lsendnbparts+rsendnbparts)=-i !$OMP END CRITICAL (nbparts) END IF END IF DO WHILE (parts%Z(i) .GT. zgrid(nz)) parts%Z(i) = parts%Z(i) - zgrid(nz) + zgrid(0) END DO DO WHILE (parts%Z(i) .LT. zgrid(0)) parts%Z(i) = parts%Z(i) + zgrid(nz) - zgrid(0) END DO END DO !$OMP END PARALLEL DO CALL particlescommunication(lsendnbparts, rsendnbparts, sendhole) END IF END subroutine bound !________________________________________________________________________________ SUBROUTINE localisation USE basic, ONLY: zgrid, rgrid, dz, nr, nnr, dr, rnorm INTEGER :: j,i, nblostparts=0 INTEGER, DIMENSION(parts%Nptot) :: losthole i=1 losthole=0 IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot parts%zindex(i)=(parts%Z(i)-zgrid(0))/dz IF (parts%R(i) .GT. rgrid(0) .AND. parts%R(i) .LT. rgrid(nnr(1))) THEN parts%rindex(i)=(parts%R(i)-rgrid(0))/dr(1) ELSE IF(parts%R(i) .GT. rgrid(nnr(1)) .AND. parts%R(i) .LT. rgrid(nnr(1)+nnr(2))) THEN parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)))/dr(2)+nnr(1) ELSE IF(parts%R(i) .GT. rgrid(nnr(1)+nnr(2)) .AND. parts%R(i) .LT. rgrid(nr)) THEN parts%rindex(i)=(parts%R(i)-rgrid(nnr(1)+nnr(2)))/dr(3)+nnr(1)+nnr(2) ELSE WRITE(*,*) "Particle:",i , " of process:", mpirank, "is out of bound in r:", parts%R(i)*rnorm, rgrid(0)*rnorm, rgrid(nr)*rnorm !nlend=.TRUE. WRITE(*,*) "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Particle removed !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" !$OMP CRITICAL (lostparts) nblostparts=nblostparts+1 losthole(nblostparts)=i !$OMP END CRITICAL (lostparts) END IF END DO !$OMP END PARALLEL DO IF(nblostparts.gt.0) THEN DO j=1,nblostparts CALL move_part(parts%Nptot, losthole(j)) parts%Nptot=parts%Nptot-1 END DO END IF !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i) DO i=1,parts%Nptot parts%WZ(i)=(parts%Z(i)-zgrid(parts%zindex(i)))/(zgrid(parts%zindex(i)+1)-zgrid(parts%zindex(i))); parts%WR(i)=(parts%R(i)-rgrid(parts%rindex(i)))/(rgrid(parts%rindex(i)+1)-rgrid(parts%rindex(i))); END DO !$OMP END PARALLEL DO END IF END SUBROUTINE localisation !________________________________________________________________________________________________________ SUBROUTINE comp_velocity ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : nlclassical ! Store old Velocities CALL swappointer(parts%UZold, parts%UZ) CALL swappointer(parts%URold, parts%UR) CALL swappointer(parts%UTHETold, parts%UTHET) CALL swappointer(parts%Gammaold, parts%Gamma) IF (nlclassical) THEN CALL comp_velocity_fun(gamma_classical) ELSE CALL comp_velocity_fun(gamma_relativistic) END IF END SUBROUTINE comp_velocity !____________________________________________ SUBROUTINE comp_velocity_fun(gamma) ! ! Computes the new velocity of the particles due to Lorentz force ! USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nz interface subroutine gamma(gam, UZ, UR, UTHET) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, INTENT(OUT):: gam end subroutine end interface DOUBLE PRECISION :: tau DOUBLE PRECISION:: 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 IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2) DO i=1,parts%Nptot J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1 J2=J1+1 J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1 J4=J3+1 ! Interpolation for magnetic field BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Bz(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Bz(J2) & & +parts%WZ(i)*parts%WR(i)*Bz(J1) BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Br(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Br(J2) & & +parts%WZ(i)*parts%WR(i)*Br(J1) ! First half of electric pulse parts%UZ(i)=parts%UZold(i)+parts%Ez(i)*tau parts%UR(i)=parts%URold(i)+parts%ER(i)*tau CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i)) ! Rotation along magnetic field ZBZ=tau*BRZ/parts%Gamma(i) ZBR=tau*BRR/parts%Gamma(i) ZPZ=parts%UZ(i)-ZBR*parts%UTHETold(i) !u'_{z} ZPR=parts%UR(i)+ZBZ*parts%UTHETold(i) !u'_{r} ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZ(i)-ZBZ*parts%UR(i)) !u'_{theta} SQR=1+ZBZ*ZBZ+ZBR*ZBR ZBZ2=2*ZBZ/SQR ZBR2=2*ZBR/SQR parts%UZ(i)=parts%UZ(i)-ZBR2*ZPTHET !u+_{z} parts%UR(i)=parts%UR(i)+ZBZ2*ZPTHET !u+_{r} parts%UTHET(i)=parts%UTHETold(i)+(ZBR2*ZPZ-ZBZ2*ZPR) !u+_{theta} ! Second half of acceleration parts%UZ(i)=parts%UZ(i)+parts%EZ(i)*tau parts%UR(i)=parts%UR(i)+parts%ER(i)*tau CALL gamma(parts%Gamma(i), parts%UZ(i), parts%UR(i), parts%UTHETold(i)) END DO !$OMP END PARALLEL DO END IF + collected=.false. END SUBROUTINE comp_velocity_fun !______________________________________________________________________________ SUBROUTINE gamma_classical(gamma, UZ, UR, UTHET) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, INTENT(OUT):: gamma gamma=1.0 END SUBROUTINE gamma_classical !______________________________________________________________________________ SUBROUTINE gamma_relativistic(gamma, UZ, UR, UTHET) DOUBLE PRECISION, INTENT(IN):: UR,UZ,UTHET DOUBLE PRECISION, INTENT(OUT):: gamma gamma=sqrt(1+UZ**2+UR**2+UTHET**2) END SUBROUTINE !________________________________________________________________________________ SUBROUTINE push Use basic, ONLY: dt, tnorm DOUBLE PRECISION:: XP, YP, COSA, SINA, U1, U2 INTEGER :: i IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i, XP, YP, COSA, SINA, U1, U2) DO i=1,parts%Nptot ! Local Cartesian coordinates XP=parts%R(i)+dt/tnorm*parts%UR(i)/parts%Gamma(i) YP=dt/tnorm*parts%UTHET(i)/parts%Gamma(i) ! Conversion to cylindrical coordiantes parts%Z(i)=parts%Z(i)+dt/tnorm*parts%UZ(i)/parts%Gamma(i) parts%R(i)=sqrt(XP**2+YP**2) ! Velocity in rotated reference frame IF (parts%R(i) .EQ. 0) THEN COSA=1 SINA=0 ELSE COSA=XP/parts%R(i) SINA=YP/parts%R(i) END IF U1=COSA*parts%UR(i)+SINA*parts%UTHET(i) U2=-SINA*parts%UR(i)+COSA*parts%UTHET(i) parts%UR(i)=U1 parts%UTHET(i)=U2 END DO !$OMP END PARALLEL DO END IF + collected=.false. END SUBROUTINE push !_______________________________________________________________________________ !_______________________________________________________________________________ SUBROUTINE diagnostics ! ! Compute energies ! !!$ CALL energies USE constants, ONLY: vlight USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, partdensity, ierr, step, it2d, nlend, nr, itparts, nplasma DOUBLE PRECISION :: gammamid INTEGER :: i INTEGER, DIMENSION(:) :: stat(3*MPI_STATUS_SIZE) - INTEGER, DIMENSION(:), ALLOCATABLE :: displs CALL MPI_WAIT(ireducerequest(3), stat(2*MPI_STATUS_SIZE+1:3*MPI_STATUS_SIZE), ierr) partdensity=0 ekin=0 epot=0 etot=0 !$OMP PARALLEL DO REDUCTION(+:epot, ekin, partdensity) DEFAULT(SHARED) PRIVATE(i,gammamid) DO i=1,parts%Nptot ! epot=epot+qsim*parts%pot(i)*phinorm IF(.not. nlclassical) THEN gammamid=1/sqrt(1-(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i)))/(parts%Gammaold(i)*parts%Gamma(i))) ekin=ekin+msim*vlight**2*(gammamid-1) ELSE ekin=ekin+0.5*msim*vlight**2*(abs(parts%URold(i)*parts%UR(i))+abs(parts%UZold(i)*parts%UZ(i))+abs(parts%UTHETold(i)*parts%UTHET(i))) END IF IF(modulo(step,it2d).eq. 0 .or. nlend) THEN partdensity(parts%zindex(i)+1+parts%rindex(i)*nz)=partdensity(parts%zindex(i)+1+parts%rindex(i)*nz)+1 END IF END DO !$OMP END PARALLEL DO IF(mpirank .eq.0 ) THEN CALL MPI_IREDUCE(MPI_IN_PLACE, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(1), ierr) CALL MPI_IREDUCE(MPI_IN_PLACE, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(2), ierr) ELSE CALL MPI_IREDUCE(epot, epot, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(1), ierr) CALL MPI_IREDUCE(ekin, ekin, 1, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(2), ierr) END IF ! Stores the particle density on the grid IF(modulo(step,it2d).eq. 0 .or. nlend) THEN IF(mpirank .eq.0 ) THEN CALL MPI_IREDUCE(MPI_IN_PLACE, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) ELSE CALL MPI_IREDUCE(partdensity, partdensity, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) END IF END IF IF(mpisize .gt. 1) THEN Nplocs_all(mpirank)=parts%Nptot IF(mpirank .eq.0 ) THEN CALL MPI_gather(MPI_IN_PLACE, 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) ELSE CALL MPI_gather(Nplocs_all(mpirank), 1, MPI_INTEGER, Nplocs_all, 1, MPI_INTEGER,& & 0, MPI_COMM_WORLD, ierr) END IF IF(mpirank .eq.0 .and. INT(SUM(Nplocs_all)).ne. nplasma) THEN WRITE(*,*) "Error particle lost on some processus" !CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) END IF END IF CALL MPI_WAIT(ireducerequest(1), stat(1:MPI_STATUS_SIZE), ierr) CALL MPI_WAIT(ireducerequest(2), stat(MPI_STATUS_SIZE+1:2*MPI_STATUS_SIZE), ierr) etot=epot+ekin ! ! Shift to Etot at cstep=1 (not valable yet at cstep=0!) IF(cstep.EQ.1) THEN etot0 = etot END IF IF(mpisize .gt. 1 .and. (modulo(step,itparts) .eq. 0 .or. nlend)) THEN - IF(mpirank .ne. 0) THEN - ALLOCATE(displs(0:mpisize-1)) - displs=0 + CALL collectparts + END IF + end subroutine diagnostics + + SUBROUTINE collectparts + USE basic, ONLY: mpirank, mpisize, ierr + INTEGER, DIMENSION(:), ALLOCATABLE :: displs + INTEGER:: i + + IF(collected) RETURN ! exit subroutine if particles have already been collected this time step + + IF(mpirank .ne. 0) THEN + ALLOCATE(displs(0:mpisize-1)) + displs=0 ! Send Particles informations to root process - CALL MPI_Gatherv(parts%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(parts%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(parts%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(parts%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(parts%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(parts%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - ELSE + CALL MPI_Gatherv(parts%Z, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,& + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(parts%R, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,& + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(parts%UR, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,& + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(parts%UZ, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, & + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(parts%UTHET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, & + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(parts%pot, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, & + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(parts%Rindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, & + & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(parts%Zindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, & + & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + ELSE ALLOCATE(displs(0:mpisize-1)) displs(0)=0 Do i=1,mpisize-1 displs(i)=displs(i-1)+Nplocs_all(i-1) END DO ! Receive particle information from all processes - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,& - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, & - & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - END IF + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%Z, Nplocs_all, displs,& + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%R, Nplocs_all, displs,& + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UR, Nplocs_all, displs,& + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UZ, Nplocs_all, displs, & + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%UTHET, Nplocs_all, displs, & + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%pot, Nplocs_all, displs, & + & MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, & + & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, & + & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) END IF - end subroutine diagnostics + collected=.TRUE. + END SUBROUTINE collectparts !_______________________________________________________________________________ SUBROUTINE adapt_vinit ! ! Computes the velocity at time -dt/2 from velocities computed at time 0 ! USE basic, ONLY : omegac, omegap, dt, tnorm, BZ, BR, nlclassical, phinorm, nz, distribtype, H0, vnorm DOUBLE PRECISION :: tau, BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, & & SQR, ZBZ2, ZBR2, Vperp, v2 INTEGER :: J1, J2, J3, J4, i DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: VZ, VR, VTHET IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3) THEN ALLOCATE(VR(parts%Nptot),VZ(parts%Nptot),VTHET(parts%Nptot)) CALL loduni(7,VZ) VZ=VZ*2*pi VTHET=parts%UTHET/parts%Gamma*vnorm DO i=1,parts%Nptot Vperp=sqrt(MAX(2*H0/me+2*elchar/me*parts%pot(i)*phinorm-VTHET(i)**2,0.0_db)) VR(i)=Vperp*sin(VZ(i)) VZ(i)=Vperp*cos(VZ(i)) IF(nlclassical) THEN parts%Gamma(i)=1 ELSE v2=VR(i)**2+VZ(i)**2+VTHET(i)**2 parts%Gamma(i)=sqrt(1/(1-v2)) END IF parts%UR(i)=parts%Gamma(i)*VR(i)/vnorm parts%UZ(i)=parts%Gamma(i)*VZ(i)/vnorm parts%UTHET(i)=parts%Gamma(i)*VTHET(i)/vnorm END DO DEALLOCATE(VR,VZ,VTHET) END IF ! Normalised time increment tau=-omegac/2/omegap*dt/tnorm ! Store old Velocities CALL swappointer(parts%UZold, parts%UZ) CALL swappointer(parts%URold, parts%UR) CALL swappointer(parts%UTHETold, parts%UTHET) CALL swappointer(parts%Gammaold, parts%Gamma) IF (parts%Nptot .NE. 0) THEN !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,J1,J2,J3,J4,BRZ, BRR, ZBR, ZBZ, ZPR, ZPZ, ZPTHET, SQR, ZBZ2, ZBR2) DO i=1,parts%Nptot J1=(parts%rindex(i))*(nz+1)+parts%zindex(i)+1 J2=J1+1 J3=(parts%rindex(i)+1)*(nz+1)+parts%zindex(i)+1 J4=J3+1 ! Interpolation for magnetic field BRZ=(1-parts%WZ(i))*(1-parts%WR(i))*Bz(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Bz(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Bz(J2) & & +parts%WZ(i)*parts%WR(i)*Bz(J1) BRR=(1-parts%WZ(i))*(1-parts%WR(i))*Br(J4) & & +parts%WZ(i)*(1-parts%WR(i))*Br(J3) & & +(1-parts%WZ(i))*parts%WR(i)*Br(J2) & & +parts%WZ(i)*parts%WR(i)*Br(J1) ! Half inverse Rotation along magnetic field ZBZ=tau*BRZ/parts%Gammaold(i) ZBR=tau*BRR/parts%Gammaold(i) SQR=1+ZBZ*ZBZ+ZBR*ZBR ZPZ=(parts%UZold(i)-ZBR*parts%UTHETold(i))/SQR !u-_{z} ZPR=(parts%URold(i)+ZBZ*parts%UTHETold(i))/SQR !u-_{r} ZPTHET=parts%UTHETold(i)+(ZBR*parts%UZold(i)-ZBZ*parts%URold(i))/SQR !u-_{theta} parts%UZ(i)=ZPZ parts%UR(i)=ZPR parts%UTHET(i)=ZPTHET ! half of decceleration parts%UZ(i)=parts%UZ(i)+parts%Ez(i)*tau parts%UR(i)=parts%UR(i)+parts%Er(i)*tau IF(.not. nlclassical) THEN parts%Gamma(i)=sqrt(1+parts%UZ(i)**2+parts%UR(i)**2+parts%UTHET(i)**2) END IF END DO !$OMP END PARALLEL DO END IF END SUBROUTINE adapt_vinit !_______________________________________________________________________________ SUBROUTINE distribpartsonprocs ! 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, nplasma, ierr, Zbounds INTEGER:: k DOUBLE PRECISION:: 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 partspercol=0 idealnbpartsperproc = FLOOR(REAL(nplasma)/REAL(mpisize)) Zmin=MINVAL(parts%Zindex) Zmax=MAXVAL(parts%Zindex) Zperproc=(Zmax-Zmin)/mpisize partindexstart=1 totparts=0 partindexlast=parts%Nptot DO k=0,mpisize-2 Nplocs_all(k)=idealnbpartsperproc END DO NPlocs_all(mpisize-1)=idealnbpartsperproc+MODULO(nplasma,mpisize) DO k=1, mpisize-1 Zbounds(k)=parts%Z(INT(k*idealnbpartsperproc)) END DO ! Store local end of particle vector partindexlast=SUM(Nplocs_all(0:mpirank)) ! Store local start of particle vector IF(mpirank .ne. 0) partindexstart=SUM(Nplocs_all(0:mpirank-1))+1 IF(mpirank .ne. 0) THEN DO k= partindexstart, partindexlast CALL move_part(k,k-partindexstart+1) END DO END IF parts%Nptot=Nplocs_all(mpirank) WRITE(*,*) mpirank, " Zbounds: ", Zbounds(mpirank), Zbounds(mpirank+1), " indices ", partindexstart, partindexlast, " nptot", parts%Nptot IF(INT(SUM(Nplocs_all)) .ne. nplasma) THEN WRITE(*,*) "Error in particle counting on proc:", mpirank, " local sum:", INT(SUM(Nplocs_all)), " nplasma:", nplasma CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) END IF END SUBROUTINE distribpartsonprocs !_______________________________________________________________________________ SUBROUTINE particlescommunication(lsendnbparts, rsendnbparts, sendholes) USE mpihelper, ONLY: particle_type USE basic, ONLY: rightproc, leftproc, ierr INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) INTEGER, ASYNCHRONOUS :: rrecvnbparts=0, lrecvnbparts=0 INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendrequest, recvrequest INTEGER, DIMENSION(2), ASYNCHRONOUS :: sendstatus(2*MPI_STATUS_SIZE), recvstatus(2*MPI_STATUS_SIZE) INTEGER :: lsentnbparts, rsentnbparts INTEGER :: lreceivednbparts, rreceivednbparts lsentnbparts=lsendnbparts rsentnbparts=rsendnbparts sendrequest=MPI_REQUEST_NULL recvrequest=MPI_REQUEST_NULL CALL MPI_IRECV(lrecvnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, recvrequest(1), ierr) CALL MPI_IRECV(rrecvnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, recvrequest(2), ierr) CALL MPI_ISEND(lsentnbparts, 1, MPI_INTEGER, leftproc, 0, MPI_COMM_WORLD, sendrequest(1), ierr) CALL MPI_ISEND(rsentnbparts, 1, MPI_INTEGER, rightproc, 0, MPI_COMM_WORLD, sendrequest(2), ierr) CALL MPI_Wait(recvrequest(1), recvstatus(1), ierr) CALL MPI_Wait(recvrequest(2), recvstatus(2), ierr) recvrequest=MPI_REQUEST_NULL lreceivednbparts=lrecvnbparts rreceivednbparts=rrecvnbparts IF( ALLOCATED(rrecvpartbuff) .AND. (size(rrecvpartbuff) .lt. rrecvnbparts) ) DEALLOCATE(rrecvpartbuff) IF(.not. ALLOCATED(rrecvpartbuff) .and. rrecvnbparts .gt. 0) ALLOCATE(rrecvpartbuff(INT(rreceivednbparts*1.3))) IF( ALLOCATED(lrecvpartbuff) .AND. (size(lrecvpartbuff) .lt. lrecvnbparts) ) DEALLOCATE(lrecvpartbuff) IF((.not. ALLOCATED(lrecvpartbuff) ).and. lrecvnbparts .gt. 0) ALLOCATE(lrecvpartbuff(INT(lreceivednbparts*1.3))) ! Receive particles from left and right processes IF ( lrecvnbparts .gt. 0) THEN CALL MPI_IRECV(lrecvpartbuff, lreceivednbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, recvrequest(1), ierr) END IF IF( rrecvnbparts .gt. 0) THEN CALL MPI_IRECV(rrecvpartbuff, rreceivednbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, recvrequest(2), ierr) END IF IF ( (lsendnbparts + rsendnbparts) .gt. 0) THEN CALL AddPartSendBuffers(lsendnbparts, rsendnbparts, sendholes) END IF CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr) CALL MPI_Wait(sendrequest(2), sendstatus(MPI_STATUS_SIZE+1), ierr) IF( lsendnbparts .gt. 0) THEN CALL MPI_ISEND(lsendpartbuff, lsendnbparts, particle_type, leftproc, 1, MPI_COMM_WORLD, sendrequest(1), ierr) END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_ISEND(rsendpartbuff, rsendnbparts, particle_type, rightproc, 1, MPI_COMM_WORLD, sendrequest(2), ierr) END IF 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 END IF IF ( rreceivednbparts .gt. 0) THEN CALL MPI_Wait(recvrequest(2), recvstatus(MPI_STATUS_SIZE+1), 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 END IF CALL Addincomingparts(rreceivednbparts, lreceivednbparts, lsendnbparts+rsendnbparts, sendholes) IF( lsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(1), sendstatus(1), ierr) END IF IF( rsendnbparts .gt. 0) THEN CALL MPI_Wait(sendrequest(2), sendstatus(2), ierr) END IF ! ! END SUBROUTINE particlescommunication !_______________________________________________________________________________ SUBROUTINE Addincomingparts(rrecvnbparts, lrecvnbparts, sendnbparts, sendholes) ! USE mpihelper INTEGER, INTENT(in) :: rrecvnbparts, lrecvnbparts, sendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) INTEGER k,partpos, partdiff ! computes if we received less particles than we sent partdiff=sendnbparts-rrecvnbparts-lrecvnbparts IF(rrecvnbparts .gt. 0) THEN Do k=1,rrecvnbparts IF(k .le. sendnbparts) THEN partpos=abs(sendholes(k)) ELSE parts%Nptot=parts%Nptot+1 partpos=parts%Nptot END IF CALL Insertincomingpart(rrecvpartbuff, k, partpos) END DO END IF IF(lrecvnbparts .gt. 0) THEN Do k=1,lrecvnbparts IF(k+rrecvnbparts .le. sendnbparts) THEN partpos=abs(sendholes(k+rrecvnbparts)) ELSE parts%Nptot=parts%Nptot+1 partpos=parts%Nptot END IF CALL Insertincomingpart(lrecvpartbuff, k, partpos) END DO END IF IF(partdiff .gt. 0) THEN DO k=rrecvnbparts+lrecvnbparts+1, sendnbparts CALL move_part(parts%Nptot,abs(sendholes(k))) parts%Nptot = parts%Nptot-1 END DO END IF ! END SUBROUTINE Addincomingparts !_______________________________________________________________________________ SUBROUTINE Insertincomingpart(buffer, bufferindex, partsindex) USE mpihelper INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(in) :: buffer parts%Rindex(partsindex) = buffer(bufferindex)%Rindex parts%Zindex(partsindex) = buffer(bufferindex)%Zindex parts%R(partsindex) = buffer(bufferindex)%R parts%Z(partsindex) = buffer(bufferindex)%Z parts%UZ(partsindex) = buffer(bufferindex)%UZ parts%UR(partsindex) = buffer(bufferindex)%UR parts%UTHET(partsindex) = buffer(bufferindex)%UTHET parts%Gamma(partsindex) = buffer(bufferindex)%Gamma ! ! END SUBROUTINE Insertincomingpart !_______________________________________________________________________________ SUBROUTINE Insertsentpart(buffer, bufferindex, partsindex) USE mpihelper INTEGER, INTENT(in) :: bufferindex, partsindex TYPE(particle), DIMENSION(:), ALLOCATABLE, INTENT(inout) :: buffer buffer(bufferindex)%Rindex = parts%Rindex(partsindex) buffer(bufferindex)%Zindex = parts%Zindex(partsindex) buffer(bufferindex)%R = parts%R(partsindex) buffer(bufferindex)%Z = parts%Z(partsindex) buffer(bufferindex)%UZ = parts%UZ(partsindex) buffer(bufferindex)%UR = parts%UR(partsindex) buffer(bufferindex)%UTHET = parts%UTHET(partsindex) buffer(bufferindex)%Gamma = parts%Gamma(partsindex) ! ! END SUBROUTINE Insertsentpart !_______________________________________________________________________________ SUBROUTINE AddPartSendBuffers(lsendnbparts, rsendnbparts, sendholes) ! USE mpihelper INTEGER, INTENT(in) :: lsendnbparts, rsendnbparts INTEGER, INTENT(in) :: sendholes(parts%Nptot) INTEGER:: partpos, k INTEGER:: lsendpos, rsendpos lsendpos=0 rsendpos=0 IF( ALLOCATED(rsendpartbuff) .AND. size(rsendpartbuff) .lt. rsendnbparts ) DEALLOCATE(rsendpartbuff) IF( ALLOCATED(lsendpartbuff) .AND. size(lsendpartbuff) .lt. lsendnbparts ) DEALLOCATE(lsendpartbuff) IF(.not. ALLOCATED(rsendpartbuff) .and. rsendnbparts .gt. 0) ALLOCATE(rsendpartbuff(INT(rsendnbparts*1.3))) IF(.not. ALLOCATED(lsendpartbuff) .and. lsendnbparts .gt. 0) ALLOCATE(lsendpartbuff(INT(lsendnbparts*1.3))) Do k=lsendnbparts+rsendnbparts,1,-1 partpos=abs(sendholes(k)) IF(sendholes(k) .GT. 0) THEN rsendpos=rsendpos+1 CALL Insertsentpart(rsendpartbuff, rsendpos, partpos) ELSE IF(sendholes(k) .LT. 0) THEN lsendpos=lsendpos+1 CALL Insertsentpart(lsendpartbuff, lsendpos, partpos) END IF END DO ! ! END SUBROUTINE AddPartSendBuffers !_______________________________________________________________________________ SUBROUTINE exchange_parts(index1, index2) INTEGER, INTENT(IN) :: index1, index2 DOUBLE PRECISION:: R,Z,UR, UZ, UTHET, Gamma INTEGER :: Rindex, Zindex !! Exchange particle at index1 with particle at index2 ! Store part at index1 in temporary value Gamma = parts%Gamma(index1) R = parts%R(index1) Z = parts%Z(index1) UR = parts%UR(index1) UTHET = parts%UTHET(index1) UZ = parts%UZ(index1) Rindex = parts%Rindex(index1) Zindex = parts%Zindex(index1) ! Move part at index2 in part at index 1 parts%Gamma(index1) = parts%Gamma(index2) parts%R(index1) = parts%R(index2) parts%Z(index1) = parts%Z(index2) parts%UR(index1) = parts%UR(index2) parts%UTHET(index1) = parts%UTHET(index2) parts%UZ(index1) = parts%UZ(index2) parts%Rindex(index1) = parts%Rindex(index2) parts%Zindex(index1) = parts%Zindex(index2) ! Move temporary values from part(index1) to part(index2) parts%Gamma(index2) = Gamma parts%R(index2) = R parts%Z(index2) = Z parts%UR(index2) = UR parts%UTHET(index2) = UTHET parts%UZ(index2) = UZ parts%Rindex(index2) = Rindex parts%Zindex(index2) = Zindex END SUBROUTINE exchange_parts !_______________________________________________________________________________ SUBROUTINE move_part(sourceindex, destindex) !! Move particle with index sourceindex to particle with index destindex !! This will destroy particle at destindex INTEGER, INTENT(IN) :: destindex, sourceindex ! Move part at sourceindex in part at destindex parts%Gamma(destindex) = parts%Gamma(sourceindex) parts%R(destindex) = parts%R(sourceindex) parts%Z(destindex) = parts%Z(sourceindex) parts%UR(destindex) = parts%UR(sourceindex) parts%UTHET(destindex) = parts%UTHET(sourceindex) parts%UZ(destindex) = parts%UZ(sourceindex) parts%Rindex(destindex) = parts%Rindex(sourceindex) parts%Zindex(destindex) = parts%Zindex(sourceindex) END SUBROUTINE move_part !_______________________________________________________________________________ !________________________________________________________________________________ ! SUBROUTINE loduni(nbase, y) ! ! Load a uniform distribution using the Hammersley's sequence. ! (nbase = 0 => Random sampling) ! INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER :: n, i ! n = SIZE(y) ! SELECT CASE (nbase) CASE(0) CALL RANDOM_NUMBER(y) CASE(1) DO i=1,n y(i) = (i-0.5_db)/n END DO CASE(2:) DO i=1,n y(i) = rev(nbase,i) END DO CASE default WRITE(*,'(a,i5)') 'Invalid value of NBASE =', nbase END SELECT ! END SUBROUTINE loduni !________________________________________________________________________________ SUBROUTINE lodgaus(nbase, y) ! ! Sample y from the Gauss distributrion ! DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase ! INTEGER :: n, i, iflag DOUBLE PRECISION :: r(SIZE(y)) DOUBLE PRECISION :: t ! DOUBLE PRECISION :: c0, c1, c2 DOUBLE PRECISION :: d1, d2, d3 DATA c0, c1, c2/2.515517_db, 0.802853_db, 0.010328_db/ DATA d1, d2, d3/1.432788_db, 0.189269_db, 0.001308_db/ ! n = SIZE(y) CALL loduni(nbase, r) r = 1.0E-5_db + 0.99998_db*r ! DO i=1,n iflag = 1 IF (r(i) .GT. 0.5_db) THEN r(i) = 1.0_db - r(i) iflag = -1 END IF t = SQRT(LOG(1.0_db/(r(i)*r(i)))) y(i) = t - (c0+c1*t+c2*t**2) / (1.0_db+d1*t+d2*t**2+d3*t**3) y(i) = y(i) * iflag END DO y = y - SUM(y)/REAL(n,db) END SUBROUTINE lodgaus !________________________________________________________________________________ SUBROUTINE lodinvr(nbase, y, ra, rb) ! ! Sample y from the distribution f=1/(r)H(r-ra)H(rb-r) for where H is the heavyside function ! DOUBLE PRECISION, INTENT(out) :: y(:) INTEGER, INTENT(in) :: nbase DOUBLE PRECISION, INTENT(in) :: rb, ra ! DOUBLE PRECISION :: r(SIZE(y)) ! CALL loduni(nbase, r) r=r*log(rb/ra) y=ra*exp(r) END SUBROUTINE lodinvr + !________________________________________________________________________________ + SUBROUTINE lodunir(nbase, y, ra, rb) + DOUBLE PRECISION, INTENT(out) :: y(:) + INTEGER, INTENT(in) :: nbase + DOUBLE PRECISION, INTENT(in) :: rb, ra + ! + DOUBLE PRECISION :: r(SIZE(y)) + ! + CALL loduni(nbase, r) + y=ra+(rb-ra)*r + END SUBROUTINE lodunir !________________________________________________________________________________ REAL(db) FUNCTION rev(nbase,i) ! ! Return an element of the Hammersley's sequence ! INTEGER, INTENT(IN) :: nbase ! Base of the Hammersley's sequence (IN) INTEGER, INTENT(IN) :: i ! Index of the sequence (IN) ! ! Local vars INTEGER :: j1, j2 REAL(db) :: xs, xsi !----------------------------------------------------------------------- xs = 0._db xsi = 1.0_db j2 = i DO xsi = xsi/nbase j1 = j2/nbase xs = xs + (j2-nbase*j1)*xsi j2 = j1 IF( j2.LE.0 ) EXIT END DO rev = xs END FUNCTION rev !________________________________________________________________________________ !Modified Bessel functions of the first kind of the first order FUNCTION bessi1(x) DOUBLE PRECISION :: bessi1,x DOUBLE PRECISION :: ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END FUNCTION bessi1 !________________________________________________________________________________ SUBROUTINE destroy_parts(p) TYPE(particles) :: p p%nptot=0 IF(ALLOCATED(p%Z)) DEALLOCATE(p%Z) IF(ALLOCATED(p%R)) DEALLOCATE(p%R) IF(ALLOCATED(p%WZ)) DEALLOCATE(p%WZ) IF(ALLOCATED(p%WR)) DEALLOCATE(p%WR) 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) END SUBROUTINE !________________________________________________________________________________ SUBROUTINE clean_beam ! CALL destroy_parts(parts) ! ! END SUBROUTINE clean_beam !________________________________________________________________________________ RECURSIVE SUBROUTINE quicksortparts(leftlimit, rightlimit) ! Sorts the particle according to their Z position using quicksort algorithm INTEGER,INTENT(IN):: leftlimit, rightlimit DOUBLE PRECISION:: pivot INTEGER::i, cnt, mid IF(leftlimit .ge. rightlimit) RETURN mid=(leftlimit+rightlimit)/2 IF(parts%Z(mid).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,mid) IF(parts%Z(rightlimit).lt.parts%Z(leftlimit)) CALL exchange_parts(leftlimit,rightlimit) IF(parts%Z(mid).lt.parts%Z(rightlimit)) CALL exchange_parts(rightlimit,mid) ! Store the pivot point for comparison pivot=parts%Z(rightlimit) cnt=leftlimit ! Move all parts with Z smaller than pivot to the left of pivot DO i=leftlimit, rightlimit IF(parts%Z(i) .le. pivot) THEN CALL exchange_parts(i,cnt) cnt=cnt+1 END IF END DO WRITE(*,*)"sorted parts: ", leftlimit,"to", rightlimit CALL quicksortparts(leftlimit,cnt-2) CALL quicksortparts(cnt,rightlimit) END SUBROUTINE quicksortparts !________________________________________________________________________________ SUBROUTINE swappointer( pointer1, pointer2) DOUBLE PRECISION, DIMENSION(:), POINTER, INTENT(inout):: pointer1, pointer2 DOUBLE PRECISION, DIMENSION(:), POINTER:: temppointer temppointer=>pointer1 pointer1=>pointer2 pointer2=>temppointer END SUBROUTINE swappointer !_______________________________________________________________________________ SUBROUTINE loaduniformRZ(VR,VZ,VTHET) USE basic, ONLY: plasmadim, rnorm, temp, vnorm USE constants, ONLY: me, kb, elchar DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET DOUBLE PRECISION:: vth ! Initial distribution in z with normalisation CALL loduni(1,parts%Z(:)) parts%Z(:)=(plasmadim(1)+(plasmadim(2)-plasmadim(1))*parts%Z(:))/rnorm ! Initial distribution in r with normalisation CALL loduni(2,parts%R(:)) parts%R=(plasmadim(3)+parts%R(:)*(plasmadim(4)-plasmadim(3)))/rnorm ! Initial velocities distribution vth=sqrt(3*kb*temp/me)/vnorm !thermal velocity CALL lodgaus(3,VZ(:)) CALL lodgaus(5,VR(:)) CALL lodgaus(7,VTHET(:)) VZ=VZ*vth VR=VR*vth VTHET=VTHET*vth END SUBROUTINE loaduniformRZ !_______________________________________________________________________________ - SUBROUTINE loadDavidson(VR,VZ,VTHET) + SUBROUTINE loadDavidson(VR,VZ,VTHET,lodr) USE constants, ONLY: me, kb, elchar USE basic, ONLY: nplasma, rnorm, plasmadim, distribtype, H0, P0, Rcurv, B0, width, qsim, msim, vnorm, & & omegac, zgrid, nz, rnorm, V, n0, nblock + interface + subroutine lodr(nbase,y,ra,rb) + DOUBLE PRECISION, INTENT(out) :: y(:) + INTEGER, INTENT(in) :: nbase + DOUBLE PRECISION, INTENT(in) :: rb, ra + end subroutine + end interface DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, INTENT(INOUT)::VZ, VR, VTHET DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE::ra, rb DOUBLE PRECISION :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL INTEGER :: i, j, n, blockstart, blockend, addedpart, remainparts INTEGER, DIMENSION(:), ALLOCATABLE :: blocksize Allocate(ra(nblock),rb(nblock)) - r0=(plasmadim(4)+plasmadim(3))/2 - !r0=sqrt(2*H0/(me*omegac**2)) + !r0=(plasmadim(4)+plasmadim(3))/2 + r0=sqrt(2*H0/(me*omegac**2)) halfLz=(zgrid(nz)+zgrid(0))/2 + MirrorRatio=(Rcurv-1)/(Rcurv+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) - MirrorRatio=(Rcurv-1)/(Rcurv+1) - deltar2=MirrorRatio*cos(2*pi*Le/width) - rb(n)=r0*(1+sqrt(deltar2))/(1-deltar2) - ra(n)=r0*(1-sqrt(deltar2))/(1-deltar2) + deltar2=1-MirrorRatio*cos(2*pi*Le/width) + rb(n)=r0*(1+sqrt(1-abs(omegac)*P0*deltar2/H0))/deltar2 + ra(n)=r0*(1-sqrt(1-abs(omegac)*P0*deltar2/H0))/deltar2 + WRITE(*,'(a,7(1pe12.4))') "ra,rb: ", omegac, H0, P0, r0, deltar2, ra(n), rb(n) END DO + VOL=SUM(2*pi*ra*(rb-ra)*(plasmadim(2)-plasmadim(1))/nblock) qsim=VOL*n0*elchar/nplasma msim=abs(qsim)/elchar*me V=VOL/rnorm**3 !H0=me*omegac**2*r0**2*0.5 gamma=1/sqrt(1-omegac**2*r0**2/vlight**2) !P0=H0/omegac blockstart=1 blockend=0 ALLOCATE(blocksize(nblock)) DO n=1,nblock blocksize(n)=nplasma/VOL*2*pi*ra(n)*(rb(n)-ra(n))*(plasmadim(2)-plasmadim(1))/nblock END DO remainparts=parts%Nptot-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 DO n=1,nblock blockstart=blockend+1 blockend=MIN(blockstart+blocksize(n)-1,parts%Nptot) ! Initial distribution in z with normalisation between magnetic mirrors CALL loduni(1,parts%Z(blockstart:blockend)) parts%Z(blockstart:blockend)=(plasmadim(2)-plasmadim(1))/rnorm/nblock*((n-1)+parts%Z(blockstart:blockend))+plasmadim(1)/rnorm - IF(distribtype.eq. 2) THEN - CALL loduni(3,parts%R(blockstart:blockend)) - parts%R(blockstart:blockend)=ra(n)+(rb(n)-ra(n))*parts%R(blockstart:blockend) - ELSE IF(distribtype.eq.3) THEN - CALL lodinvr(3,parts%R(blockstart:blockend),ra(n),rb(n)) - END IF + CALL lodr(3,parts%R(blockstart:blockend),ra(n), rb(n)) parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm END DO ! Load velocities theta velocity ! Loading of r and z velocity is done in adapt_vinit to have ! access to parts%pot DO i=1,parts%Nptot ! Interpolation for Magnetic potential rg=parts%R(i)*rnorm zg=(parts%Z(i)-halfLz)*rnorm Athetpos=0.5*B0*(rg - width/pi*MirrorRatio*bessi1(2*pi*rg/width)*COS(2*pi*zg/width)) Pcomp=abs(P0/rg/me) Acomp=-SIGN(elchar/me*Athetpos,qsim) VTHET(i)=SIGN(MIN(abs(Pcomp+Acomp),sqrt(2*H0/me)),Pcomp+Acomp) END DO VTHET=VTHET/vnorm VZ=0._db VR=0._db END SUBROUTINE loadDavidson END MODULE beam diff --git a/src/chkrst.f90 b/src/chkrst.f90 index fcd4a69..91cd720 100644 --- a/src/chkrst.f90 +++ b/src/chkrst.f90 @@ -1,122 +1,103 @@ SUBROUTINE chkrst(flag) ! ! Process checkpoint/restart file ! USE basic USE futils USE beam USE fields IMPLICIT NONE INTEGER, INTENT(in) :: flag INTEGER, DIMENSION(:), ALLOCATABLE:: displs INTEGER :: i ! Only process 0 should save on file ! ! Local vars and arrays !________________________________________________________________________________ ! SELECT CASE(flag) !________________________________________________________________________________ ! 1. Open and read restart file ! CASE(0) IF(mpirank .ne. 0) RETURN CALL openf(rstfile, fidrst) CALL getatt(fidrst, '/Basic', 'cstep', cstep) CALL getatt(fidrst, '/Basic', 'time', time) CALL getatt(fidrst, '/Basic', 'msim', msim) CALL getatt(fidrst, '/Basic', 'qsim', qsim) CALL getatt(fidrst, '/Basic', 'V', V) CALL getatt(fidrst, '/var0d', 'etot0', etot0) CALL getatt(fidrst, '/var0d', 'epot', epot) CALL getatt(fidrst, '/var0d', 'ekin', ekin) CALL getatt(fidrst, '/var0d', 'etot', etot) CALL getatt(fidrst, '/Parts', 'nplasma', nplasma) CALL creat_parts(parts,nplasma) CALL getarr(fidrst, '/Parts/Z', parts%Z) CALL getarr(fidrst, '/Parts/R', parts%R) CALL getarr(fidrst, '/Parts/UZ', parts%UZ) CALL getarr(fidrst, '/Parts/UR', parts%UR) CALL getarr(fidrst, '/Parts/UTHET', parts%UTHET) CALL getarr(fidrst, '/Parts/Zindex', parts%Zindex) CALL getarr(fidrst, '/Parts/Rindex', parts%Rindex) CALL closef(fidrst) WRITE(*,'(3x,a)') & & "Reading from restart file "//TRIM(rstfile)//" completed!" !________________________________________________________________________________ ! 2. Create and write to restart file (DP reals) ! CASE(1) - ALLOCATE(displs(0:mpisize-1)) - displs(0)=1 - DO i=1, mpisize-1 - displs(i)=displs(i-1)+Nplocs_all(i-1)-1 - END DO - IF(mpisize .gt. 1) THEN - IF(mpirank .ne. 0) THEN - ! Send Particles informations to root process - CALL MPI_Gatherv(parts%Zindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, & - & MPI_INTEGER,0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(parts%Rindex, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, & - & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - RETURN - ELSE - ! Receive particle information from all processes - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Zindex, Nplocs_all, displs, & - & MPI_INTEGER,0, MPI_COMM_WORLD, ierr) - CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%Rindex, Nplocs_all, displs, & - & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - END IF - parts%Nptot=sum(Nplocs_all) - END IF + IF(mpisize .gt. 1) THEN + call collectparts + END IF IF( .NOT. nlsave ) RETURN CALL mv2bk(rstfile) CALL creatf(rstfile, fidrst, real_prec='d', desc='Restart file') CALL creatg(fidrst, '/Basic', 'Basic data') CALL attach(fidrst, '/Basic', 'cstep', cstep) CALL attach(fidrst, '/Basic', 'time', time) CALL attach(fidrst, '/Basic', 'jobnum', jobnum) CALL attach(fidrst, '/Basic', 'qsim', qsim) CALL attach(fidrst, '/Basic', 'msim', msim) CALL attach(fidrst, '/Basic', 'V', V) ! ! 0D variables ! CALL creatg(fidrst, '/var0d', '0D variables') CALL attach(fidrst, '/var0d','etot0', etot0) CALL attach(fidrst, '/var0d','epot', epot) CALL attach(fidrst, '/var0d','ekin', ekin) CALL attach(fidrst, '/var0d','etot', etot) ! ! Parts ! CALL creatg(fidrst, '/Parts', 'Particles data') CALL attach(fidrst, '/Parts', 'nplasma', nplasma) CALL putarr(fidrst, '/Parts/Z', parts%Z) CALL putarr(fidrst, '/Parts/R', parts%R) CALL putarr(fidrst, '/Parts/UZ', parts%UZ) CALL putarr(fidrst, '/Parts/UR', parts%UR) CALL putarr(fidrst, '/Parts/UTHET', parts%UTHET) CALL putarr(fidrst, '/Parts/Zindex', parts%Zindex) CALL putarr(fidrst, '/Parts/Rindex', parts%Rindex) ! ! Fields ! CALL closef(fidrst) WRITE(*,'(3x,a)') & & "Writing to restart file "//TRIM(rstfile)//" completed!" ! END SELECT ! END SUBROUTINE chkrst diff --git a/src/diagnose.f90 b/src/diagnose.f90 index b4dd903..48aaf46 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,214 +1,209 @@ SUBROUTINE diagnose(kstep) ! ! Diagnostics ! USE basic USE futils USE hashtable USE beam, ONLY : parts, epot, ekin, etot, etot0, ireducerequest, Nplocs_all USE xg, ONLY : initw, updt_xg_var USE mpi IMPLICIT NONE ! INTEGER, INTENT(in) :: kstep ! ! Local vars and arrays INTEGER, PARAMETER :: BUFSIZE = 20 CHARACTER(len=128) :: str, fname INTEGER :: partdensitystatus(MPI_STATUS_SIZE) ! Only process 0 should save on file IF(mpirank .ne. 0) RETURN !________________________________________________________________________________ ! 1. Initial diagnostics IF( kstep .EQ. 0) THEN ! WRITE(*,'(a)') ' Initial diagnostics' ! ! 1.1 Initial run or when NEWRES set to .TRUE. ! IF( .NOT. nlres .OR. newres) THEN CALL creatf(resfile, fidres) 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/input.", "input") 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) ! Part and fields vectors ! Initialize time-dependent particle variables CALL creatd(fidres, 1, SHAPE(parts%R), "/data/part/R", "R") CALL creatd(fidres, 1, SHAPE(parts%Z), "/data/part/Z", "Z") CALL creatd(fidres, 1, SHAPE(parts%UZ), "/data/part/UR", "UZ") CALL creatd(fidres, 1, SHAPE(parts%UR), "/data/part/UZ", "UR") CALL creatd(fidres, 1, SHAPE(parts%UTHET), "/data/part/UTHET", "UTHET") CALL creatd(fidres, 1, SHAPE(parts%pot), "/data/part/pot", "pot") 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(partdensity), "/data/fields/partdensity", "partdensity") IF(mpisize .gt. 1) CALL creatd(fidres, 1, SHAPE(Nplocs_all), "/data/part/Nplocs_all", "Nplocs_all") CALL putarr(fidres, "/data/fields/Br", Br) CALL putarr(fidres, "/data/fields/Bz", Bz) ! ! 1.2 Restart run ! ELSE CALL cp2bk(resfile) ! backup previous result file CALL openf(resfile, fidres) 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", temp) 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", qsim) CALL attach(fidres, TRIM(str), "msim", msim) CALL attach(fidres, TRIM(str), "V", V) 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) ! 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)) CLOSE(lu_in) ! ! Initialize buffers for 0d history arrays CALL htable_init(hbuf0, BUFSIZE) CALL set_htable_fileid(hbuf0, fidres, "/data/var0d") ! ! ! Initialize Xgrafix ! IF(nlxg) THEN CALL initw END IF END IF !________________________________________________________________________________ ! 2. Periodic diagnostics IF( kstep .NE. -1) THEN ! IF(modulo(step,ittext).eq. 0 .or. nlend) THEN - WRITE(*,'(a,1x,i6.6,a1,i6.6,20x,a,1pe10.3)') & + WRITE(*,'(a,1x,i8.8,a1,i8.8,20x,a,1pe10.3)') & & '*** Timestep (this run/total) =', step, '/', cstep, 'Time =', time WRITE(*,'(a,4(1pe12.4))') 'Epot, Ekin, Etot, Eerr', epot, ekin, etot, etot-etot0 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 htable_endstep(hbuf0) IF(mpisize .gt. 1) CALL append(fidres, "/data/part/Nplocs_all", REAL(Nplocs_all,kind=db)) END IF ! ! 2.2 1d profiles IF(modulo(step,it2d).eq. 0 .or. nlend) THEN CALL append(fidres, "/data/fields/pot", pot) CALL append(fidres, "/data/fields/Er", Er) CALL append(fidres, "/data/fields/Ez", Ez) CALL MPI_WAIT(ireducerequest(3), partdensitystatus, ierr) CALL append(fidres, "/data/fields/partdensity", partdensity) END IF ! ! 2.3 2d profiles IF(modulo(step,itparts).eq. 0 .or. nlend) THEN !PRINT*, 'write particles to file_____________________' CALL append(fidres, "/data/part/R", parts%R) CALL append(fidres, "/data/part/Z", parts%Z) CALL append(fidres, "/data/part/UZ", parts%UZ) CALL append(fidres, "/data/part/UR", parts%UR) CALL append(fidres, "/data/part/UTHET", parts%UTHET) CALL append(fidres, "/data/part/pot", parts%pot) - - - !WRITE(*,'(a,6(1pe12.3))') 'R-Z bounds of space = ', rgrid(0), rgrid(nnr(1)), rgrid(nnr(2)+nnr(1)), rgrid(nr), zgrid(0), zgrid(nz) - !WRITE(*,'(a,4(1pe12.3))') 'MinMax part. R, Z =', & - ! & MINVAL(parts%R), MAXVAL(parts%R), MINVAL(parts%Z), MAXVAL(parts%Z) ! END IF ! ! 2.4 3d profiles ! ! ! 2.5 Xgrafix ! IF(nlxg .AND. modulo(kstep,itgraph) .eq. 0 .AND. mpisize .eq. 1) THEN call xgevent CALL updt_xg_var CALL xgupdate END IF !________________________________________________________________________________ ! 3. Final diagnostics ELSE ! ! 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/fields_mod.f90 b/src/fields_mod.f90 index f06f7c5..6516d92 100644 --- a/src/fields_mod.f90 +++ b/src/fields_mod.f90 @@ -1,320 +1,324 @@ MODULE fields ! ! Module containing all variables and routines used to obtain |E(r,z) and |B(r,z) ! USE constants USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, splrz, nlperiod, phinorm, nlPhis, nrank, rhs, mpirank, mpisize USE beam, ONLY: parts USE bsplines USE mumps_bsplines USE mpi IMPLICIT NONE DOUBLE PRECISION, allocatable, SAVE :: matcoef(:,:), sol(:), vec1(:), vec2(:) -DOUBLE PRECISION, ALLOCATABLE, SAVE :: fun(:,:), fun2(:,:) TYPE(mumps_mat), SAVE :: femat ! FEM matrix CONTAINS !________________________________________________________________________________ SUBROUTINE fields_init USE basic, ONLY: pot, nlperiod, nrank, rhs USE bsplines USE mumps_bsplines INTEGER :: nrz(2) ! Set up 2d spline splrz CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz, nlppform=nlppform, period=nlperiod) ! Calculate dimension of splines nrz(1)=nz nrz(2)=nr CALL get_dim(splrz, nrank, nrz, femorder) ALLOCATE(matcoef(nrank(1),nrank(2))) ALLOCATE(pot((nr+1)*(nz+1))) ALLOCATE(rhs(nrank(1)*nrank(2))) - ALLOCATE(sol(nrank(1)*nrank(2))) - ALLOCATE(fun(1:femorder(1)+1,0:1),fun2(1:femorder(2)+1,0:1))!Arrays keeping values of b-splines at gauss node - + ALLOCATE(sol(nrank(1)*nrank(2))) ! Calculate magnetic field components in grid points (Davidson analitical formula employed) ALLOCATE(Br((nr+1)*(nz+1)),Bz((nr+1)*(nz+1))) ALLOCATE(Athet((nr+1)*(nz+1))) ALLOCATE(Er((nr+1)*(nz+1)),Ez((nr+1)*(nz+1))) CALL magnet ! Calculate and factorise FEM matrix (depended only on mesh) CALL fematrix CALL factor(femat) END SUBROUTINE fields_init !________________________________________________________________________________ SUBROUTINE rhscon USE bsplines USE mumps_bsplines USE mpi USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, nplasma USE beam, ONLY: parts USE mpihelper - INTEGER :: zleft, rleft, irow, jcol, it, jw, mu, i + INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch + INTEGER,DIMENSION(:),ALLOCATABLE::zleft, rleft + DOUBLE PRECISION, ALLOCATABLE :: fun(:,:,:), fun2(:,:,:) + nbunch=20 + ALLOCATE(zleft(nbunch),rleft(nbunch)) + ALLOCATE(fun(1:femorder(1)+1,0:0,nbunch),fun2(1:femorder(2)+1,0:0,nbunch))!Arrays keeping values of b-splines at gauss node rhs=0 !Initialise rhs ! Assemble rhs IF (parts%Nptot .ne. 0) THEN -!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(zleft, rleft, jw, it, i, irow, jcol, mu) PRIVATE(fun, fun2) REDUCTION(+:rhs) - DO i=1,parts%Nptot - CALL locintv(splrz%sp1, parts%Z(i), zleft) - CALL locintv(splrz%sp2, parts%R(i), rleft) - CALL basfun(parts%Z(i), splrz%sp1, fun, zleft+1) - CALL basfun(parts%R(i), splrz%sp2, fun2, rleft+1) +!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(zleft, rleft, jw, it, i, iend, irow, jcol, mu, k) PRIVATE(fun, fun2) REDUCTION(+:rhs) + DO i=1,parts%Nptot,nbunch + iend=min(i+nbunch-1,parts%Nptot) + CALL locintv(splrz%sp1, parts%Z(i:iend), zleft) + CALL locintv(splrz%sp2, parts%R(i:iend), rleft) + CALL basfun(parts%Z(i:iend), splrz%sp1, fun, zleft+1) + CALL basfun(parts%R(i:iend), splrz%sp2, fun2, rleft+1) - DO jw=1,(femorder(2)+1) + DO k=1,nbunch + DO jw=1,(femorder(2)+1) DO it=1,(femorder(1)+1) - irow=zleft+it - jcol=rleft+jw + irow=zleft(k)+it + jcol=rleft(k)+jw mu=irow+(jcol-1)*(nrank(1)) - rhs(mu)=rhs(mu)+fun(it,0)*fun2(jw,0)/2/pi*omegap/omegac*V/nplasma + rhs(mu)=rhs(mu)+fun(it,0,k)*fun2(jw,0,k)/2/pi*omegap/omegac*V/nplasma END DO + END DO END DO END DO !$OMP END PARALLEL DO END IF IF(mpisize .gt. 1) THEN IF (mpirank.eq.0) THEN CALL MPI_REDUCE(MPI_IN_PLACE, rhs, nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, & & MPI_SUM, 0, MPI_COMM_WORLD, ierr) ELSE CALL MPI_REDUCE(rhs, rhs, nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, & & MPI_SUM, 0, MPI_COMM_WORLD, ierr) END IF END IF ! Dirichlet BC on RHS IF(rgrid(0).ne.0.0_db) rhs(1:nrank(1))=potinn rhs((nrank(2)-1)*nrank(1)+1:nrank(2)*nrank(1))=potout - + DEALLOCATE(fun,fun2, zleft, rleft) END SUBROUTINE rhscon !========================================================================! SUBROUTINE poisson USE bsplines USE mumps_bsplines USE futils INTEGER:: partend, i, ierr, jder(2) jder(1)=0 jder(2)=0 partend=parts%Nptot IF(mpirank.eq.0) CALL bsolve(femat,rhs,sol) CALL MPI_Bcast(sol(1), nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) matcoef= reshape(sol, (/nrank(1),nrank(2)/)) CALL gridval(splrz, vec1, vec2, pot, jder, matcoef) CALL getgrad(splrz, parts%Z(1:partend), parts%R(1:partend), parts%pot(1:partend), parts%EZ(1:partend), parts%ER(1:partend)) IF (nlphis) THEN Do i=1,partend parts%EZ(i)=-parts%Ez(i) parts%ER(i)=-parts%ER(i) END DO ELSE Do i=1,partend parts%EZ(i)=0 parts%ER(i)=0 END DO END IF IF( mpirank .eq. 0) THEN CALL gridval(splrz, vec1, vec2, Ez, (/1,0/)) CALL gridval(splrz, vec1, vec2, Er, (/0,1/)) Ez=-Ez Er=-Er END IF END SUBROUTINE poisson !========================================================================! SUBROUTINE fematrix USE bsplines USE mumps_bsplines DOUBLE PRECISION, ALLOCATABLE :: xgauss(:,:), wgauss(:), zg(:), rg(:), wzg(:), wrg(:) DOUBLE PRECISION, ALLOCATABLE :: f(:,:), aux(:) DOUBLE PRECISION, ALLOCATABLE :: coefs(:) + DOUBLE PRECISION, ALLOCATABLE :: fun(:,:), fun2(:,:) DOUBLE PRECISION :: contrib INTEGER, ALLOCATABLE :: idert(:), iderw(:) INTEGER :: i, j, k, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms kterms=2 !Number of terms in weak form - + ALLOCATE(fun(1:femorder(1)+1,0:1),fun2(1:femorder(2)+1,0:1))!Arrays keeping values of b-splines at gauss node ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE(f((femorder(1)+1)*(femorder(2)+1),2),aux(femorder(1)+1)) !Auxiliary arrays ordering bsplines ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms)) !Pointers on the order of derivatives ! Constuction of auxiliary array ordering bsplines in given interval DO i=1,(femorder(1)+1) aux(i)=i END DO DO i=1,(femorder(2)+1) f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),1)=aux f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),2)=i END DO ! Initialisation of FEM matrix CALL init(nrank(1)*nrank(2),2,femat) ! Assemble FEM matrix DO j=1,nr ! Loop on r position DO i=1,nz ! Loop on z position ! Computation of gauss weight and position in r and z direction for gaussian integration CALL get_gauss(splrz%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(splrz%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration DO k=1,ngauss(2) xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k) wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg END DO DO igauss=1,ngauss(1)*ngauss(2) ! Loop on gaussian weights and positions CALL basfun(xgauss(igauss,1), splrz%sp1, fun, i) CALL basfun(xgauss(igauss,2), splrz%sp2, fun2, j) CALL coefeq(xgauss(igauss,:), idert, iderw, coefs) DO iterm=1,kterms ! Loop on the two integration dimensions DO jt=1,(1+femorder(1))*(femorder(2)+1) DO iw=1,(1+femorder(1))*(femorder(2)+1) irow=i+f(jt,1)-1; jcol=j+f(jt,2)-1 irow2=i+f(iw,1)-1; jcol2=j+f(iw,2)-1 mu=irow+(jcol-1)*nrank(1) mu2=irow2+(jcol2-1)*nrank(1) contrib=fun(f(jt,1),idert(iterm))* fun(f(iw,1),idert(iterm))* & & fun2(f(jt,2),iderw(iterm))*fun2(f(iw,2),iderw(iterm))* & & wgauss(igauss)*xgauss(igauss,2) CALL updtmat(femat, mu, mu2, contrib) END DO END DO END DO END DO END DO END DO ! Impose Dirichlet boundary conditions on FEM matrix CALL fe_dirichlet DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE(f, aux) DEALLOCATE(idert, iderw, coefs, fun,fun2) - ALLOCATE(fun(1:femorder(1)+1,0:0),fun2(1:femorder(2)+1,0:0))!Arrays keeping values of b-splines at gauss node END SUBROUTINE fematrix SUBROUTINE fe_dirichlet DOUBLE PRECISION, ALLOCATABLE :: arr(:) INTEGER :: i ALLOCATE(arr(nrank(1)*nrank(2))) DO i=1,nrank(1)*nrank(2) CALL getrow(femat,i,arr) END DO DO i=1,nrank(1) IF(rgrid(0).ne.0.0_db) THEN arr=0; arr(i)=1; CALL putrow(femat,i,arr) END IF arr=0; arr(nrank(1)*nrank(2)+1-i)=1 ; CALL putrow(femat,nrank(1)*nrank(2)+1-i,arr) END DO DEALLOCATE(arr) END SUBROUTINE fe_dirichlet !________________________________________________________________________________ SUBROUTINE coefeq(x, idt, idw, c) DOUBLE PRECISION, INTENT(in) :: x(2) INTEGER, INTENT(out) :: idt(:), idw(:) DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt)) c(1) = x(2) idt(1) = 0 idw(1) = 1 c(2) = x(2) idt(2) = 1 idw(2) = 0 END SUBROUTINE coefeq !________________________________________________________________________________ SUBROUTINE magnet USE basic, ONLY: B0, Rcurv, rgrid, zgrid, width, rnorm, nr, nz, bnorm USE constants, ONLY: Pi DOUBLE PRECISION :: rg, zg,halfLz, MirrorRatio INTEGER :: i, rindex halfLz=(zgrid(nz)+zgrid(0))/2 MirrorRatio=(Rcurv-1)/(Rcurv+1) DO i=1,(nr+1)*(nz+1) rindex=(i-1)/(nz+1) rg=rgrid(rindex) zg=zgrid(i-rindex*(nz+1)-1)-halfLz Br(i)=-B0*MirrorRatio*SIN(2*pi*zg/width*rnorm)*bessi1(2*pi*rg/width*rnorm)/bnorm Bz(i)=B0*(1-MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm))/bnorm Athet(i)=0.5*B0*(rg*rnorm - width/pi*MirrorRatio*bessi1(2*pi*rg/width*rnorm)*COS(2*pi*zg/width*rnorm)) END DO END SUBROUTINE magnet !________________________________________________________________________________ !Modified Bessel functions of the first kind of the zero order FUNCTION bessi0(x) DOUBLE PRECISION :: bessi0,x DOUBLE PRECISION :: ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, 0.1328592d-1, 0.225319d-2, -0.157565d-2 ,0.916281d-2, & & -0.2057706d-1, 0.2635537d-1,-0.1647633d-1, 0.392377d-2 / if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))) else ax=abs(x) y=3.75/ax bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) endif return END FUNCTION bessi0 !________________________________________________________________________________ !Modified Bessel functions of the first kind of the first order FUNCTION bessi1(x) DOUBLE PRECISION :: bessi1,x DOUBLE PRECISION :: ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/ax bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endif return END FUNCTION bessi1 !________________________________________________________________________________ SUBROUTINE clean_fields Use bsplines DEALLOCATE(matcoef) DEALLOCATE(pot) DEALLOCATE(rhs) DEALLOCATE(sol) - DEALLOCATE(fun,fun2) DEALLOCATE(Br,Bz) DEALLOCATE(Er,Ez) DEALLOCATE(vec1,vec2) Call DESTROY_SP(splrz) END SUBROUTINE clean_fields END MODULE fields diff --git a/src/main.f90 b/src/main.f90 index a1fa92d..51d19e1 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -1,84 +1,85 @@ PROGRAM main ! ! Skeleton for a time dependent program ! Note: Even in this sequential version, MPI is required ! because of FUTILS (more specifcally because ! of the HASTABLE module)! ! USE basic USE mpi USE bsplines USE mumps_bsplines USE futils IMPLICIT NONE INTEGER:: required, provided ! ! required=MPI_THREAD_FUNNELED CALL mpi_init_thread(required,provided,ierr) CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) CALL MPI_COMM_SIZE(MPI_COMM_WORLD, mpisize, ierr) !-------------------------------------------------------------------------------- ! 1. Prologue CALL timera(0, 'Prologue') CALL daytim('Start at') ! ! Define data specific to run ! CALL basic_data !Definition of global variables and input paramaters loading step=0 ! IF( .NOT. nlres ) THEN CALL newrun !not implemented yet ELSE CALL restart !not implemented yet END IF ! ! Compute auxilliary values ! CALL auxval !time independent values ! ! Initial conditions ! IF( .NOT. nlres ) THEN CALL inital !plasma initialisation ELSE CALL resume !loads restart.h5 file END IF ! ! Start or restart the run ! CALL start !not implemented yet ! ! Initial diagnostocs ! CALL diagnose(0) CALL timera(1, 'Prologue') !-------------------------------------------------------------------------------- ! 2. Time stepping CALL timera(0, 'Main loop') ! DO step = step+1 cstep = cstep+1 time = time+dt CALL tesend CALL stepon CALL diagnose(step) + IF(modulo(step,itrestart) .eq. 0) CALL chkrst(1) IF( nlend ) EXIT END DO CALL timera(1, 'Main loop') !-------------------------------------------------------------------------------- ! 9. Epilogue CALL timera(0, 'Epilogue') ! CALL diagnose(-1) CALL endrun IF(mpirank .eq. 0) THEN CALL timera(1, 'Epilogue') CALL timera(9, '') CALL timera(-1, '') CALL daytim('Done at ') END IF CALL mpi_finalize(ierr) END PROGRAM main