diff --git a/.gitignore b/.gitignore index 085a0ab..6f37af2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ *.h5 +*.h5_old *_genmod.f90 *.mod *.o *.out *.m~ diff --git a/src/.vscode/launch.json b/.vscode/launch.json similarity index 63% rename from src/.vscode/launch.json rename to .vscode/launch.json index 664bc72..0005e2a 100644 --- a/src/.vscode/launch.json +++ b/.vscode/launch.json @@ -1,39 +1,47 @@ { // Use IntelliSense to learn about possible attributes. // Hover to view descriptions of existing attributes. // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 "version": "0.0.1", "configurations": [ + { + "name": "Python: Current File", + "type": "python", + "request": "launch", + "program": "${file}", + "args": [ "/misc/espic2dsimus/openmp/trajectories.h5" ], + "console": "integratedTerminal" + }, { "name": "Fortran Launch (GDB)", "type": "cppdbg", "request": "launch", - "program": "${workspaceFolder}/espic2d", - "args": ["<", "${workspaceFolder}/../wk/espic2d.in"], + "program": "${workspaceFolder}/srd/espic2d", + "args": ["<", "${workspaceFolder}/wk/espic2d.in"], "stopAtEntry": false, - "cwd": "${workspaceFolder}/../wk", + "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", }, { "name": "Fortran Launch stable.in", "type": "cppdbg", "request": "launch", - "program": "${workspaceFolder}/espic2d", - "args": ["${workspaceFolder}/../wk/stable.in"], + "program": "${workspaceFolder}/src/espic2d", + "args": ["${workspaceFolder}/wk/stable.in"], "stopAtEntry": false, - "cwd": "${workspaceFolder}/../wk", + "cwd": "${workspaceFolder}/wk", //"externalConsole": true, "preLaunchTask": "make debug", "miDebuggerPath": "gdb", }, { "name": "Intel Debug Attach", "type": "cppdbg", "request": "attach", "processId": "${command:pickProcess}", "program": "${workspaceFolder}/espic2d" } ] } diff --git a/src/.vscode/tasks.json b/.vscode/tasks.json similarity index 90% rename from src/.vscode/tasks.json rename to .vscode/tasks.json index 805b2e9..eaa071f 100644 --- a/src/.vscode/tasks.json +++ b/.vscode/tasks.json @@ -1,51 +1,51 @@ { "version": "2.0.0", "tasks": [ { "label": "make debug", "command": "bash", "type": "shell", "options": { - "cwd": "${workspaceFolder}/../src", + "cwd": "${workspaceFolder}/src", "env": { "PLATFORM": "intel", } }, "args": [ "-c", "make debug" ], "group": { "kind": "build", "isDefault": true } }, { "label": "make clean", "command": "bash", "type": "shell", "args": [ "-c", "cd ../src && make clean" ], "group": "build" }, { "label": "make release", "command": "bash", "type": "shell", "options": { - "cwd": "${workspaceFolder}/../src", + "cwd": "${workspaceFolder}/src", "env": { "PLATFORM": "intel", } }, "args": [ "-c", "make" ], "group": "build", "problemMatcher": [] } ] } diff --git a/matlab/PlotParticleTrajectory.m b/matlab/PlotParticleTrajectory.m index 422a2cc..6f69a27 100644 --- a/matlab/PlotParticleTrajectory.m +++ b/matlab/PlotParticleTrajectory.m @@ -1,28 +1,44 @@ -function PlotParticleTrajectory(M,partnumber) +function PlotParticleTrajectory(M,partnumber,t) %UNTITLED3 Summary of this function goes here % Detailed explanation goes here f=figure(); -sgtitle(sprintf('Particle %d',partnumber)); -ax=subplot(1,2,1); -plot(ax,M.Z(partnumber,:),M.R(partnumber,:)) -xlim(ax,[M.zgrid(1) M.zgrid(end)]) -ylim(ax,[M.rgrid(1) M.rgrid(end)]) -xlabel(ax,'Z') -ylabel(ax,'R') -title(ax,'Position') -ax.XTick=M.zgrid; -ax.YTick=M.rgrid; -xlim(ax,[M.zgrid(1) M.zgrid(end)]) -ylim(ax,[M.rgrid(1) M.rgrid(end)]) -grid on; -ax=subplot(1,2,2); -V=M.VR(partnumber,:); -t=[0:size(M.VR,2)-1]*M.dt*double(M.it2); -plot(ax,t,V','r-') -xlabel(ax,'t [s]') -ylabel(ax,'VR [m/s]') -title(ax,'Radial velocity') -grid on; +time=((t(1):t(end))-1)*M.dt*double(M.it2); +%sgtitle(sprintf('Particles %d',partnumber)); +for part=partnumber +ax1=subplot(1,2,1); +hold on +plot(ax1,M.Z(part,t),M.R(part,t),'x','Linewidth',1.2,'Displayname',sprintf('part=%d',part),'MarkerIndices',1:10:length(t)); +ax2=subplot(1,2,2); +hold on +plot(ax2,M.R(part,t).*cos(M.THET(part,t)),M.R(part,t).*sin(M.THET(part,t)),'x-.','Linewidth',1.2,'Displayname',sprintf('part=%d',part),'MarkerIndices',1:10:length(t)) + +% ax=subplot(1,2,2); +% hold on +% V=M.VR(part,t); +% plot(ax,time,V','Displayname',sprintf('part=%d',part)) +end +xlabel(ax1,'Z') +ylabel(ax1,'R') +sgtitle(sprintf('Position between t=[%3.2f-%3.2f] ns',time(1)*1e9,time(end)*1e9)) +xlim(ax1,[M.zgrid(1) M.zgrid(end)]) +ylim(ax1,[M.rgrid(1) M.rgrid(41)]) +legend(ax1) +grid(ax1,'on'); +xlabel(ax2,'X') +ylabel(ax2,'Y') +axis(ax2,'equal') +legend(ax2) +grid(ax2,'on'); + +% ax=subplot(1,2,2); +% xlabel(ax,'t [s]') +% ylabel(ax,'VR [m/s]') +% title(ax,'Radial velocity') +% grid on; + +f.PaperOrientation='landscape'; +[~, name, ~] = fileparts(M.file); +print(f,sprintf('%sTrajectories',name),'-dpdf','-fillpage') end diff --git a/matlab/analyse_espic2D.m b/matlab/analyse_espic2D.m index d8c36fd..1b962dd 100644 --- a/matlab/analyse_espic2D.m +++ b/matlab/analyse_espic2D.m @@ -1,150 +1,192 @@ -file='teststablefinefine.h5'; +file='teststable_13_fineinitDav3.h5'; %close all hidden; if (~exist('M','var')) M.file=file; 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-') +semilogy(M.t0d(tmin:tmax),abs(M.eerr(tmin:tmax)/M.etot(1)),'h-') xlabel('t [s]') -ylabel('E_{err} %') +ylabel('E_{err}/E_{tot}') 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'; f=figure(); zpos=65; tinit=1; tend=size(M.N,3); -deltat=200;%floor(tend/2); +deltat=1000;%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=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)); 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(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 [m^{-3}]') J=find(psi>0); 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(); +tstudied=0; 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)); +h2=histogram(P,10,'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) +%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)]) +xlim([0.7e-24 inf]) print(f,sprintf('%sParts_HP',name),'-dpdf','-fillpage') +f=figure(); +tstudied=0; +legtext=sprintf("t=%2.1f - %2.1f [ns]",M.tpart(end-tstudied)*1e9,M.tpart(end)*1e9); +subplot(1,3,1) +VR=M.VR(:,1); +h1=histogram(VR,20,'BinLimits',[min(VR(:)) max(VR(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +VR=mean(M.VR(:,end-tstudied:end),2); +h1=histogram(VR,'BinWidth',h1.BinWidth,'DisplayName',legtext); +ylabel('counts') +xlabel('V_r [m/s]') +legend + +subplot(1,3,2) +VTHET=M.VTHET(:,1); +h1=histogram(VTHET,20,'BinLimits',[min(VTHET(:)) max(VTHET(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +VTHET=mean(M.VTHET(:,end-tstudied:end),2); +h1=histogram(VTHET,'BinWidth',h1.BinWidth,'DisplayName',legtext); +ylabel('counts') +xlabel('V_\theta [m/s]') +legend + +subplot(1,3,3) +hold off +VZ=M.VZ(:,1); +h1=histogram(VZ,20,'BinLimits',[min(VZ(:)) max(VZ(:))],'DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +VZ=mean(M.VZ(:,end-tstudied:end),2); +h1=histogram(VZ,'BinWidth',h1.BinWidth,'DisplayName',legtext); +ylabel('counts') +xlabel('V_z [m/s]') +legend +f=gcf; +f.PaperOrientation='landscape'; +[~, name, ~] = fileparts(M.file); +print(f,sprintf('%sParts_V',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))) %dispespicFields(M) BrillouinRatio=2*M.omepe^2/M.omece^2 \ No newline at end of file diff --git a/matlab/distribution_function.m b/matlab/distribution_function.m new file mode 100644 index 0000000..d3a9db8 --- /dev/null +++ b/matlab/distribution_function.m @@ -0,0 +1,60 @@ +file='teststable_13_fineinithigh.h5'; +%close all hidden; +if (~exist('M','var')) + M.file=file; +end +M=load_espic2d(file,M,'all'); + +N=M.Npartscell(:,:,1); +[~,index]=max(N,[],'all','linear'); +[Rindex,Zindex]=ind2sub([70,128],index); +Rindex=16; Zindex=64; +Indices=M.Rindex(:,1)==Rindex & M.Zindex(:,1)==Zindex; + +Vr=M.VR(Indices,1)/M.vlight; +Vz=M.VZ(Indices,1)/M.vlight; +Vthet=M.VTHET(Indices,1)/M.vlight; + +Indicesend=M.Rindex(:,end)==Rindex & M.Zindex(:,end)==Zindex; + +Vrend=M.VR(Indicesend,end)/M.vlight; +Vzend=M.VZ(Indicesend,end)/M.vlight; +Vthetend=M.VTHET(Indicesend,end)/M.vlight; + +binwidth=0.1; +f=figure(); +tstudied=0; +legtext=sprintf("t=%2.1f [ns]",M.tpart(end)*1e9); +subplot(1,3,1) +h1=histogram(Vr,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +h1=histogram(Vrend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +ylabel('probability') +xlabel('\beta_r') +legend +grid on + +subplot(1,3,2) +binwidth=0.01; +h1=histogram(Vthet(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +h1=histogram(Vthetend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +ylabel('probability') +xlabel('\beta_\theta') +legend +grid on + +subplot(1,3,3) +binwidth=0.1; +h1=histogram(Vz(:,1),'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(1)*1e9)); +hold on +h1=histogram(Vzend,'Binwidth',binwidth,'Normalization','probability','DisplayName',sprintf("t=%2.3d [ns]",M.tpart(end)*1e9)); +ylabel('probability') +xlabel('\beta_z') +legend +grid on +f=gcf; +sgtitle(sprintf('R=%1.2e[m] Z=%1.2e[m] dt=%1.2e[ns]',M.rgrid(Rindex+1),M.zgrid(Zindex+1),M.dt*1e9)) +f.PaperOrientation='landscape'; +[~, name, ~] = fileparts(M.file); +print(f,sprintf('%sParts_V_RZ',name),'-dpdf','-fillpage') \ No newline at end of file diff --git a/matlab/load_espic2d.m b/matlab/load_espic2d.m index 9cacbad..5ffc965 100644 --- a/matlab/load_espic2d.m +++ b/matlab/load_espic2d.m @@ -1,164 +1,213 @@ 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(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.kb=1.38064852E-23; -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.job_time = h5readatt(file,'/data/input.00/','job_time'); +Data.extra_time = h5readatt(file,'/data/input.00/','extra_time'); +Data.dt = h5readatt(file,'/data/input.00/','dt'); +Data.tmax= h5readatt(file,'/data/input.00/','tmax'); +Data.nrun = h5readatt(file,'/data/input.00/','nrun'); 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'); +Data.nz = h5readatt(file,'/data/input.00/','nz'); +Data.nnr = h5read(file,'/data/input.00/nnr'); +Data.lz = h5read(file,'/data/input.00/lz'); +Data.qsim = h5readatt(file,'/data/input.00/','qsim'); +Data.msim = h5readatt(file,'/data/input.00/','msim'); +Data.nplasma = h5readatt(file,'/data/input.00/','nplasma'); +Data.potinn = h5readatt(file,'/data/input.00/','potinn'); +Data.potout = h5readatt(file,'/data/input.00/','potout'); +Data.B0 = h5readatt(file,'/data/input.00/','B0'); +Data.Rcurv = h5readatt(file,'/data/input.00/','Rcurv'); +Data.width = h5readatt(file,'/data/input.00/','width'); +Data.n0 = h5readatt(file,'/data/input.00/','n0'); +Data.temp = h5readatt(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'); + 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'); + 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.rgrid= h5read(file, '/data/var1d/rgrid')*Data.rnorm; +Data.zgrid= h5read(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; +Br = h5read(file,'/data/fields/Br')*Data.bnorm; Data.Br= reshape(Br,length(Data.zgrid),length(Data.rgrid)); -Bz = hdf5read(file,'/data/fields/Bz')*Data.bnorm; +Bz = h5read(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.plasmadim = hdf5read(file,'/data/input.00/plasmadim'); -Data.radii = hdf5read(file,'/data/input.00/radii'); +Data.femorder = h5read(file,'/data/input.00/femorder'); +Data.ngauss = h5read(file,'/data/input.00/ngauss'); +Data.plasmadim = h5read(file,'/data/input.00/plasmadim'); +Data.radii = h5read(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.epot = h5read(file,'/data/var0d/epot'); +Data.ekin = h5read(file,'/data/var0d/ekin'); +Data.etot = h5read(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; + pot = h5read(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; + Er = h5read(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; + Ez = h5read(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'); + N = h5read(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.Npartscell=Data.N; 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')) + 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.R = h5read(file,'/data/part/R')*Data.rnorm; + Data.Z = h5read(file,'/data/part/Z')*Data.rnorm; + try + Data.THET = h5read(file,'/data/part/THET'); + catch + clear Data.THET + end + try + Data.Rindex=h5read(file,'/data/part/Rindex'); + Data.Zindex=h5read(file,'/data/part/Zindex'); + catch + clear Data.Rindex Data.Zindex + end + Data.partepot = sign(Data.qsim)*h5read(file,'/data/part/pot')*Data.phinorm*Data.qe; + Data.UR = h5read(file,'/data/part/UR'); + Data.UZ = h5read(file,'/data/part/UZ'); + Data.UTHET= h5read(file,'/data/part/UTHET'); 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)); + try + Data.partindex = h5read(file,'/data/part/partindex'); + [~,Indices]=sort(Data.partindex); + for i=1:size(Indices,2) + Data.H(:,i)=Data.H(Indices(:,i),i); + Data.P(:,i)=Data.P(Indices(:,i),i); + Data.R(:,i)=Data.R(Indices(:,i),i); + Data.Z(:,i)=Data.Z(Indices(:,i),i); + Data.Rindex(:,i)=Data.Rindex(Indices(:,i),i); + Data.Zindex(:,i)=Data.Zindex(Indices(:,i),i); + Data.VR(:,i)=Data.VR(Indices(:,i),i); + Data.VZ(:,i)=Data.VZ(Indices(:,i),i); + Data.VTHET(:,i)=Data.VTHET(Indices(:,i),i); + Data.VPERP(:,i)=Data.VPERP(Indices(:,i),i); + Data.GAMM(:,i)=Data.GAMM(Indices(:,i),i); + Data.partepot(:,i)=Data.partepot(Indices(:,i),i); + Data.THET(:,i)=Data.THET(Indices(:,i),i); + end + catch + end + Data.ur=zeros(sum(Data.nnr),Data.nz,size(Data.R,2)); + Data.uz=zeros(sum(Data.nnr),Data.nz,size(Data.R,2)); + Data.uthet=zeros(sum(Data.nnr),Data.nz,size(Data.R,2)); + Data.Nlocal=zeros(sum(Data.nnr),Data.nz,size(Data.R,2)); + for i=1,size(Data.R,2) + for ipart=1,size(Data.R,1) + Data.ur(Data.Rindex(ipart,i)+1,Data.Zindex(ipart,i)+1,i)=Data.ur(Data.Rindex(ipart,i)+1,Data.Zindex(ipart,i)+1,i)+double(Data.VR(ipart,i)); + Data.uz(Data.Rindex(ipart,i)+1,Data.Zindex(ipart,i)+1,i)=Data.uz(Data.Rindex(ipart,i)+1,Data.Zindex(ipart,i)+1,i)+Data.VZ(ipart,i); + Data.uthet(Data.Rindex(ipart,i)+1,Data.Zindex(ipart,i)+1,i)=Data.uthet(Data.Rindex(ipart,i)+1,Data.Zindex(ipart,i)+1,i)+Data.VTHET(ipart,i); + Data.Nlocal(Data.Rindex(ipart,i)+1,Data.Zindex(ipart,i)+1,i)=Data.Nlocal(Data.Rindex(ipart,i)+1,Data.Zindex(ipart,i)+1,i)+1; + end + end + Data.ur=Data.ur./max(Data.Nlocal,1); + Data.uz=Data.uz./max(Data.Nlocal,1); + Data.uthet=Data.uthet./max(Data.Nlocal,1); + clear Data.partindex; + Data.tpart=Data.dt*(0:size(Data.R,2)-1)*double(Data.it2); - if(strcmp(Data.type,'fields')) + 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/python/plot_trajectories_espic2d.py b/python/plot_trajectories_espic2d.py new file mode 100644 index 0000000..aa11098 --- /dev/null +++ b/python/plot_trajectories_espic2d.py @@ -0,0 +1,76 @@ +import h5py +import numpy as numpy +import matplotlib as mplib +import matplotlib.pyplot as plt +import sys, getopt +import math + + +# read commandline arguments, first +fullCmdArguments = sys.argv + +# - further arguments +filename = fullCmdArguments[1] + +f=h5py.File(filename,'r') +r=f['/data/part/R'] +z=f['/data/part/Z'] +thet=f['/data/part/THET'] +ur=f['/data/part/UR'] +uz=f['/data/part/UZ'] +uthet=f['/data/part/UTHET'] +gamma=numpy.sqrt(1+ur[()]**2+uz[()]**2+uthet[()]**2) + +vlight=299792458 +qe=1.60217662E-19 +me=9.109383E-31 +eps_0=8.85418781762E-12 +kb=1.38064852E-23 +input00=f['/data/input.00'] +B0 =input00.attrs['B0'] +Rcurv =input00.attrs['Rcurv'] +n0 =input00.attrs['n0'] +dt =input00.attrs['dt'] +it0d =input00.attrs['it0d'] +it2d =input00.attrs['it2d'] +itparts = input00.attrs['itparts'] +nplasma = input00.attrs['nplasma'] +msim = input00.attrs['msim'] +qsim = input00.attrs['qsim'] +lz = f['/data/input.00/lz'] + + +omepe=math.sqrt(abs(n0)*qe**2/(me*eps_0)) +omece=qe*B0/me; + +# Normalizations +tnorm=1/omepe; +rnorm=vlight*tnorm; +bnorm=B0; +enorm=vlight*bnorm; +phinorm=enorm*rnorm; +vnorm=vlight; + +rgrid=f['/data/var1d/rgrid'][()]*rnorm +zgrid=f['/data/var1d/zgrid'][()]*rnorm + +x=r[()]*rnorm*numpy.cos(thet[()]) +y=r[()]*rnorm*numpy.sin(thet[()]) +fig, ax = plt.subplots(1,2) + +for i in range(132300,132310): + ax[0].plot(z[0:800,i]*rnorm,r[0:800,i]*rnorm,label='part %d' %i) + ax[1].plot(x[0:800,i],y[0:800,i],label='part %d' %i) + +ax[0].legend() +ax[0].grid() +ax[0].set_xlabel('Z [m]') +ax[0].set_ylabel('R [m]') +ax[1].legend() +ax[1].grid() +ax[1].set_xlabel('X [m]') +ax[1].set_ylabel('Y [m]') +plt.show() + +#plt.savefig(filename+'_trajectories.pdf',orientation='landscape',) +#print(gamma) diff --git a/src/auxval.f90 b/src/auxval.f90 index 6abb70a..2aadd21 100644 --- a/src/auxval.f90 +++ b/src/auxval.f90 @@ -1,101 +1,107 @@ SUBROUTINE auxval ! USE constants USE basic USE fields USE beam ! ! Set auxiliary values ! IMPLICIT NONE ! INTEGER:: i !________________________________________________________________________________ IF(mpirank .eq. 0) WRITE(*,'(a/)') '=== Set auxiliary values ===' !________________________________________________________________________________ ! ! Total number of intervals nr=sum(nnr) ! Normalisation constants qsim=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)*n0*elchar/nplasma msim=abs(qsim)/elchar*me vnorm=vlight tnorm=1/sqrt(abs(n0)*elchar*elchar/me/eps_0) rnorm=vnorm*tnorm bnorm=B0 enorm=vlight*bnorm phinorm=enorm*rnorm ! Normalised boundary conditions potinn=potinn/phinorm potout=potout/phinorm ! Characteristic frequencies and normalised volume omegac=qsim/msim*B0 omegap=sqrt(elchar**2*abs(n0)/me/eps_0) V=pi*(plasmadim(2)-plasmadim(1))*(plasmadim(4)**2-plasmadim(3)**2)/rnorm/rnorm/rnorm IF(mpirank .eq. 0) THEN IF(omegap.GT.omegac) THEN WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegap/omegac', omegap, omegac, omegap/omegac ELSE WRITE(*,'(a,3(1pe12.3))') 'omegap, omegac, omegac/omegap', omegap, omegac, omegac/omegap END IF END IF ! Construction of the mesh rgrid in r and zgrid in z and its normalisation CALL mesh rgrid=rgrid/rnorm zgrid=zgrid/rnorm dz=dz/rnorm dr=dr/rnorm ! Define Zindex bounds for the MPI process in the case of sequential run ALLOCATE( Zbounds(0:mpisize), Nplocs_all(0:mpisize-1)) Zbounds(0) = Zgrid(0) Zbounds(mpisize) = Zgrid(nz) - - leftproc=mpirank-1 - IF(mpirank .eq. 0) leftproc = mpisize-1 - rightproc = mpirank+1 - IF(mpirank .eq. mpisize-1) rightproc = 0 - + IF(mpisize.gt.1) THEN + leftproc=mpirank-1 + IF(mpirank .eq. 0) leftproc = mpisize-1 + rightproc = mpirank+1 + IF(mpirank .eq. mpisize-1) rightproc = 0 + ELSE + leftproc=-1 + rightproc=-1 + END IF WRITE(*,*) "mpirank, left, right", mpirank, leftproc, rightproc ! Auxiliary vectors ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1))) DO i=0,nr vec1(i*(nz+1)+1:(i+1)*(nz+1))=zgrid!(0:nz) vec2(i*(nz+1)+1:(i+1)*(nz+1))=rgrid(i) END DO ALLOCATE(partdensity(nz*nr)) + ALLOCATE(fluidur(nz*nr)) + ALLOCATE(fluiduz(nz*nr)) + ALLOCATE(fluiduthet(nz*nr)) !________________________________________________________________________________ CONTAINS !________________________________________________________________________________ SUBROUTINE mesh USE basic, ONLY: zgrid, rgrid, nr, nnr, nz, dr, dz, lz, radii INTEGER :: j ALLOCATE(zgrid(0:nz),rgrid(0:nr)) dz=(lz(2)-lz(1))/nz DO j=0,nz zgrid(j)=j*dz+lz(1) END DO dr(1)=(radii(2)-radii(1))/nnr(1) DO j=0,nnr(1) rgrid(j)=radii(1)+j*dr(1) END DO dr(2)=(radii(3)-radii(2))/nnr(2) DO j=1,nnr(2) rgrid(nnr(1)+j)=radii(2)+j*dr(2) END DO dr(3)=(radii(4)-radii(3))/nnr(3) DO j=1,nnr(3) rgrid(nnr(1)+nnr(2)+j)=radii(3)+j*dr(3) END DO END SUBROUTINE mesh !________________________________________________________________________________ END SUBROUTINE auxval diff --git a/src/basic_mod.f90 b/src/basic_mod.f90 index df6f469..e94da41 100644 --- a/src/basic_mod.f90 +++ b/src/basic_mod.f90 @@ -1,304 +1,305 @@ 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=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 + DOUBLE PRECISION, ALLOCATABLE:: fluidur(:), fluiduz(:), fluiduthet(:) ! Fluid velocities 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 3adf201..9105e03 100644 --- a/src/beam_mod.f90 +++ b/src/beam_mod.f90 @@ -1,1170 +1,1247 @@ 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 + INTEGER, DIMENSION(:), ALLOCATABLE :: Rindex, Zindex, partindex ! Indices of the particles on the grid DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: & - & R,Z, & + & R,Z,THET, & & 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%THET(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%partindex(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%THET=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 - + INTEGER:: i 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) - CASE(2) !Stable distribution from Davidson 4.95 p.119 - 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, lodinvr) - CASE DEFAULT - IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of distribution:", distribtype - CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) - - END SELECT + ! 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) + CASE(2) !Stable distribution from Davidson 4.95 p.119 + 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, lodinvr) + CASE(4) !Stable distribution from Davidson 4.95 p.119 but with distribution in R as + CALL loadDavidson(VR, VZ, VTHET, lodgausr) + CASE(5) !Stable distribution from Davidson 4.95 p.119 with gaussian in V around v_*0=sqrt(2*H0/(3*me)) + CALL loadDavidson(VR, VZ, VTHET, lodunir) + CASE DEFAULT + IF (mpirank .eq. 0) WRITE(*,*) "Unknown type of distribution:", distribtype + CALL MPI_Abort(MPI_COMM_WORLD, -1, ierr) + END SELECT + Do i=1,nplasma + parts%partindex(i)=i + END DO 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) + IF(mpisize .gt. 1) THEN + CALL particlescommunication(lsendnbparts, rsendnbparts, sendhole) + END IF 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 + DOUBLE PRECISION:: XP, YP, COSA, SINA, U1, U2, ALPHA INTEGER :: i IF (parts%Nptot .NE. 0) THEN - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i, XP, YP, COSA, SINA, U1, U2) + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i, XP, YP, COSA, SINA, U1, U2, ALPHA) 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 + ALPHA=0 ELSE COSA=XP/parts%R(i) SINA=YP/parts%R(i) + ALPHA=asin(SINA) END IF + parts%THET(i)=parts%THET(i)+ALPHA 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 + USE basic, ONLY: qsim, phinorm, msim, cstep, nlclassical, nz, partdensity, ierr, step, it2d, nlend, nr, itparts, nplasma, fluidur, fluiduz, fluiduthet DOUBLE PRECISION :: gammamid INTEGER :: i INTEGER, DIMENSION(:) :: stat(3*MPI_STATUS_SIZE) CALL MPI_WAIT(ireducerequest(3), stat(2*MPI_STATUS_SIZE+1:3*MPI_STATUS_SIZE), ierr) partdensity=0 + fluidur=0 + fluiduz=0 + fluiduthet=0 ekin=0 epot=0 etot=0 - !$OMP PARALLEL DO REDUCTION(+:epot, ekin, partdensity) DEFAULT(SHARED) PRIVATE(i,gammamid) + !$OMP PARALLEL DO REDUCTION(+:epot, ekin, partdensity,fluidur,fluiduz,fluiduthet) 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 + fluidur(parts%zindex(i)+1+parts%rindex(i)*nz)=fluidur(parts%zindex(i)+1+parts%rindex(i)*nz)+parts%UR/parts%Gamma + fluiduz(parts%zindex(i)+1+parts%rindex(i)*nz)=fluiduz(parts%zindex(i)+1+parts%rindex(i)*nz)+parts%UZ/parts%Gamma + fluiduthet(parts%zindex(i)+1+parts%rindex(i)*nz)=fluiduthet(parts%zindex(i)+1+parts%rindex(i)*nz)+parts%UTHET/parts%Gamma END IF END DO !$OMP END PARALLEL DO + fluiduz=fluiduz/partdensity + fluidur=fluidur/partdensity + fluiduthet=fluiduthet/partdensity 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) + CALL MPI_IREDUCE(MPI_IN_PLACE, fluidur, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & + & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) + CALL MPI_IREDUCE(MPI_IN_PLACE, fluiduz, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & + & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) + CALL MPI_IREDUCE(MPI_IN_PLACE, fluiduthet, 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) + CALL MPI_IREDUCE(fluidur, fluidur, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & + & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) + CALL MPI_IREDUCE(fluiduz, fluiduz, nz*nr, MPI_DOUBLE_PRECISION, MPI_SUM, & + & 0, MPI_COMM_WORLD, ireducerequest(3), ierr) + CALL MPI_IREDUCE(fluiduthet, fluiduthet, 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 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%THET, Nplocs_all(mpirank), MPI_DOUBLE_PRECISION, parts%THET, 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) + CALL MPI_Gatherv(parts%partindex, Nplocs_all(mpirank), MPI_INTEGER, parts%partindex, 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%THET, 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) + CALL MPI_Gatherv(MPI_IN_PLACE, Nplocs_all(mpirank), MPI_INTEGER, parts%partindex, Nplocs_all, displs, & + & MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) END IF 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 + IF(distribtype .EQ. 2 .OR. distribtype .EQ. 3 .OR. distribtype .EQ. 4) 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%partindex(partsindex) = buffer(bufferindex)%partindex 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%THET(partsindex) = buffer(bufferindex)%THET 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) + buffer(bufferindex)%partindex = parts%partindex(partsindex) + 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)%THET = parts%THET(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 + DOUBLE PRECISION:: R, Z, THET, UR, UZ, UTHET, Gamma + INTEGER :: Rindex, Zindex, partindex !! Exchange particle at index1 with particle at index2 ! Store part at index1 in temporary value + partindex= parts%partindex(index1) Gamma = parts%Gamma(index1) R = parts%R(index1) Z = parts%Z(index1) + THET = parts%THET(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%partindex(index1)= parts%partindex(index2) parts%Gamma(index1) = parts%Gamma(index2) parts%R(index1) = parts%R(index2) parts%Z(index1) = parts%Z(index2) + parts%THET(index1) = parts%THET(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%partindex(index2)= partindex parts%Gamma(index2) = Gamma parts%R(index2) = R parts%Z(index2) = Z + parts%THET(index2) = THET 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%partindex(destindex)= parts%partindex(sourceindex) parts%Gamma(destindex) = parts%Gamma(sourceindex) parts%R(destindex) = parts%R(sourceindex) parts%Z(destindex) = parts%Z(sourceindex) + parts%THET(destindex) = parts%THET(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 + !________________________________________________________________________________ + SUBROUTINE lodgausr(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 lodgaus(nbase, r) + y=0.5*(rb+ra)+ 0.1*(rb-ra)*r + END SUBROUTINE lodgausr !________________________________________________________________________________ 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%THET)) DEALLOCATE(p%THET) 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) + IF(ALLOCATED(p%partindex)) DEALLOCATE(p%partindex) 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 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,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 + & omegac, zgrid, nz, rnorm, V, n0, nblock, temp 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 + DOUBLE PRECISION :: r0, athetpos, deltar2, rg, zg, halfLz, Mirrorratio, Le, Pcomp, Acomp, gamma, VOL, vth 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)) 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) 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 - CALL lodr(3,parts%R(blockstart:blockend),ra(n), rb(n)) + CALL lodr(2,parts%R(blockstart:blockend),ra(n), rb(n)) parts%R(blockstart:blockend)=parts%R(blockstart:blockend)/rnorm END DO - + IF(distribtype .eq. 5) THEN + ! 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 + ELSE ! 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 IF END SUBROUTINE loadDavidson END MODULE beam diff --git a/src/chkrst.f90 b/src/chkrst.f90 index b0a2c87..1106147 100644 --- a/src/chkrst.f90 +++ b/src/chkrst.f90 @@ -1,104 +1,108 @@ 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/THET', parts%THET) 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 getarr(fidrst, '/Parts/partindex', parts%partindex) CALL closef(fidrst) WRITE(*,'(3x,a)') & & "Reading from restart file "//TRIM(rstfile)//" completed!" !________________________________________________________________________________ ! 2. Create and write to restart file (DP reals) ! CASE(1) IF(mpisize .gt. 1) THEN call collectparts END IF IF(mpirank .ne. 0) RETURN 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/THET', parts%THET) 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) + CALL putarr(fidrst, '/Parts/partindex', parts%partindex) ! ! 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 6b61a3c..dd1dd8c 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,213 +1,222 @@ 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%Z), "/data/part/THET", "THET") 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%Rindex), "/data/part/Rindex", "Rindex") CALL creatd(fidres, 1, SHAPE(parts%Zindex), "/data/part/Zindex", "Zindex") + CALL creatd(fidres, 1, SHAPE(parts%partindex), "/data/part/partindex", "partindex") 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") + CALL creatd(fidres, 1, SHAPE(fluidur), "/data/fields/fluidur", "fluidur") + CALL creatd(fidres, 1, SHAPE(fluiduz), "/data/fields/fluiduz", "fluiduz") + CALL creatd(fidres, 1, SHAPE(fluiduthet), "/data/fields/fluiduthet", "fluiduthet") 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,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) - + CALL append(fidres, "/data/fields/fluidur", fluidur) + CALL append(fidres, "/data/fields/fluiduz", fluiduz) + CALL append(fidres, "/data/fields/fluiduthet", fluiduthet) 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/THET", parts%THET) 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) CALL append(fidres, "/data/part/Rindex", REAL(parts%Rindex,kind=db)) CALL append(fidres, "/data/part/Zindex", REAL(parts%Zindex,kind=db)) + CALL append(fidres, "/data/part/partindex", REAL(parts%partindex,kind=db)) ! 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/mpihelper_mod.f90 b/src/mpihelper_mod.f90 index 44ba473..75f1e9d 100644 --- a/src/mpihelper_mod.f90 +++ b/src/mpihelper_mod.f90 @@ -1,138 +1,140 @@ MODULE mpihelper USE constants USE mpi IMPLICIT NONE TYPE BASICDATA LOGICAL :: nlres, nlsave, newres, nlxg, nlppform, nlPhis, nlclassical INTEGER :: nplasma, nz, it0d, it2d, itparts, & & itgraph, distribtype, nblock, nrun DOUBLE PRECISION :: job_time, extra_time, tmax, dt, & & potinn, potout, B0, n0, temp, & & Rcurv, width, H0, P0 INTEGER, DIMENSION(2) :: femorder, ngauss INTEGER, DIMENSION(3) :: nnr DOUBLE PRECISION, DIMENSION(2):: lz DOUBLE PRECISION, DIMENSION(4):: radii, plasmadim CHARACTER(len=64) :: resfile= "results.h5" END TYPE BASICDATA TYPE particle - INTEGER :: Rindex, Zindex - DOUBLE PRECISION :: R, Z, & + INTEGER :: Rindex, Zindex, partindex + DOUBLE PRECISION :: R, Z, THET,& & UZ, UR, UTHET, & & Gamma END TYPE particle INTEGER, SAVE :: basicdata_type INTEGER, SAVE :: particle_type INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: rhscol_type INTEGER, SAVE:: rhsoverlap_type INTEGER, DIMENSION(:), ALLOCATABLE, SAVE:: rhsrequests, rhsoverlaprequests INTEGER:: rhscommtag=10, rhsoverlapcommtag = 15 DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE, SAVE:: rhsoverlap CONTAINS !========================================================================= SUBROUTINE init_mpitypes CALL init_particlempi CALL init_basicdatampi END SUBROUTINE init_mpitypes !========================================================================= SUBROUTINE init_particlempi() - INTEGER :: nblock = 8 - INTEGER:: blocklength(8) - INTEGER(kind=MPI_ADDRESS_KIND):: displs(8), displ0 - INTEGER, DIMENSION(8) :: types + INTEGER :: nblock = 10 + INTEGER:: blocklength(10) + INTEGER(kind=MPI_ADDRESS_KIND):: displs(10), displ0 + INTEGER, DIMENSION(10) :: types TYPE(particle) :: part INTEGER:: ierr CALL mpi_get_address(part%Rindex, displs(1), ierr) CALL mpi_get_address(part%Zindex, displs(2), ierr) - types(1:2)=MPI_INTEGER - CALL mpi_get_address(part%R, displs(3), ierr) - CALL mpi_get_address(part%Z, displs(4), ierr) - CALL mpi_get_address(part%UZ, displs(5), ierr) - CALL mpi_get_address(part%UR, displs(6), ierr) - CALL mpi_get_address(part%UTHET, displs(7), ierr) - CALL mpi_get_address(part%GAMMA, displs(8), ierr) - types(3:8)=MPI_DOUBLE_PRECISION - blocklength(1:8) = 1 + CALL mpi_get_address(part%partindex, displs(3), ierr) + types(1:3)=MPI_INTEGER + CALL mpi_get_address(part%R, displs(4), ierr) + CALL mpi_get_address(part%Z, displs(5), ierr) + CALL mpi_get_address(part%THET, displs(6), ierr) + CALL mpi_get_address(part%UZ, displs(7), ierr) + CALL mpi_get_address(part%UR, displs(8), ierr) + CALL mpi_get_address(part%UTHET, displs(9), ierr) + CALL mpi_get_address(part%GAMMA, displs(10), ierr) + types(4:10)=MPI_DOUBLE_PRECISION + blocklength(1:10) = 1 CALL mpi_get_address(part, displ0, ierr) displs=displs-displ0 CALL MPI_Type_create_struct(nblock, blocklength, displs, types, particle_type, ierr) CALL MPI_Type_commit(particle_type,ierr) END SUBROUTINE init_particlempi !========================================================================= SUBROUTINE init_basicdatampi() INTEGER :: nblock = 36 INTEGER:: blocklength(36) INTEGER(kind=MPI_ADDRESS_KIND):: displs(36), displ0 INTEGER, DIMENSION(36) :: types TYPE(BASICDATA) :: bdata INTEGER:: ierr CALL mpi_get_address(bdata%nlres, displs(1), ierr) CALL mpi_get_address(bdata%nlsave, displs(2), ierr) CALL mpi_get_address(bdata%newres, displs(3), ierr) CALL mpi_get_address(bdata%nlxg, displs(4), ierr) CALL mpi_get_address(bdata%nlppform, displs(5), ierr) CALL mpi_get_address(bdata%nlPhis, displs(6), ierr) CALL mpi_get_address(bdata%nlclassical, displs(7), ierr) types(1:7)=MPI_LOGICAL CALL mpi_get_address(bdata%nplasma, displs(8), ierr) CALL mpi_get_address(bdata%nz, displs(9), ierr) CALL mpi_get_address(bdata%it0d, displs(10), ierr) CALL mpi_get_address(bdata%it2d, displs(11), ierr) CALL mpi_get_address(bdata%itparts, displs(12), ierr) CALL mpi_get_address(bdata%itgraph, displs(13), ierr) CALL mpi_get_address(bdata%distribtype, displs(14), ierr) CALL mpi_get_address(bdata%nblock, displs(15), ierr) CALL mpi_get_address(bdata%nrun, displs(16), ierr) types(8:16)=MPI_INTEGER CALL mpi_get_address(bdata%job_time, displs(17), ierr) CALL mpi_get_address(bdata%extra_time, displs(18), ierr) CALL mpi_get_address(bdata%tmax, displs(19), ierr) CALL mpi_get_address(bdata%dt, displs(20), ierr) CALL mpi_get_address(bdata%potinn, displs(21), ierr) CALL mpi_get_address(bdata%potout, displs(22), ierr) CALL mpi_get_address(bdata%B0, displs(23), ierr) CALL mpi_get_address(bdata%n0, displs(24), ierr) CALL mpi_get_address(bdata%temp, displs(25), ierr) CALL mpi_get_address(bdata%Rcurv, displs(26), ierr) CALL mpi_get_address(bdata%width, displs(27), ierr) CALL mpi_get_address(bdata%H0, displs(28), ierr) CALL mpi_get_address(bdata%P0, displs(29), ierr) types(17:29)=MPI_DOUBLE_PRECISION blocklength(1:29) = 1 CALL mpi_get_address(bdata%femorder, displs(30), ierr) CALL mpi_get_address(bdata%ngauss, displs(31), ierr) CALL mpi_get_address(bdata%nnr, displs(32), ierr) types(30:32)=MPI_INTEGER blocklength(30:31) = 2 blocklength(32) = 3 CALL mpi_get_address(bdata%lz, displs(33), ierr) blocklength(33) = 2 CALL mpi_get_address(bdata%radii, displs(34), ierr) CALL mpi_get_address(bdata%plasmadim, displs(35), ierr) types(33:35)=MPI_DOUBLE_PRECISION blocklength(34:35) = 4 CALL mpi_get_address(bdata%resfile, displs(36), ierr) types(36)=MPI_CHARACTER blocklength(36) = 64 CALL mpi_get_address(bdata, displ0, ierr) displs=displs-displ0 CALL MPI_Type_create_struct(nblock, blocklength, displs, types, basicdata_type, ierr) CALL MPI_Type_commit(basicdata_type,ierr) END SUBROUTINE init_basicdatampi !========================================================================= END MODULE mpihelper diff --git a/src/resume.f90 b/src/resume.f90 index 8b8d64a..76a3566 100644 --- a/src/resume.f90 +++ b/src/resume.f90 @@ -1,80 +1,84 @@ SUBROUTINE resume ! ! Resume from previous run ! Use beam Use basic Use fields IMPLICIT NONE ! ! Local vars and arrays INTEGER, DIMENSION(:), ALLOCATABLE:: displs,sendcounts INTEGER:: i !________________________________________________________________________________ WRITE(*,'(a)') ' Resume from previous run' !________________________________________________________________________________ ! ! Open and read initial conditions from restart file CALL chkrst(0) IF(mpisize .gt. 1) THEN ALLOCATE(sendcounts(0:mpisize-1), displs(0:mpisize-1)) IF(mpirank.eq.0) THEN CALL quicksortparts(1,nplasma) CALL distribpartsonprocs CALL MPI_Scatter(Nplocs_all,1,MPI_INTEGER,MPI_IN_PLACE, 1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) DO i=0,mpisize-1 displs(i)=i sendcounts(i)=2 END DO CALL MPI_Scatterv(Zbounds,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,2,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) displs(0)=0 sendcounts(0)=Nplocs_all(0) DO i=1,mpisize-1 displs(i)=displs(i-1)+Nplocs_all(i-1) sendcounts(i)=Nplocs_all(i) END DO WRITE(*,*) mpirank, "I am the root process" CALL MPI_BARRIER(MPI_COMM_WORLD, ierr) CALL MPI_Scatterv(parts%R,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Z,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%THET,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%UZ,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%UR,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%UTHET,sendcounts, displs,MPI_DOUBLE_PRECISION,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Rindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Zindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%partindex,sendcounts, displs,MPI_INTEGER,MPI_IN_PLACE,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) WRITE(*,*) mpirank, "I am the root process" ELSE WRITE(*,*) mpirank, "I am not the root process" CALL creat_parts(parts,nplasma) CALL MPI_Scatter(Nplocs_all,1,MPI_INTEGER, Nplocs_all(mpirank), 1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) sendcounts=2 CALL MPI_Scatterv(Zbounds,sendcounts,displs,MPI_DOUBLE_PRECISION,Zbounds(mpirank),2,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) parts%Nptot=Nplocs_all(mpirank) WRITE(*,*) mpirank, "received Nplocs_all : ", Nplocs_all WRITE(*,*) mpirank, "received Zbounds : ", Zbounds CALL MPI_BARRIER(MPI_COMM_WORLD, ierr) CALL MPI_Scatterv(parts%R,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%R,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Z,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%Z,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%THET,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%THET,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%UZ,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UZ,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%UR,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UR,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%UTHET,sendcounts, displs,MPI_DOUBLE_PRECISION,parts%UTHET,Nplocs_all(mpirank),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Rindex,sendcounts, displs,MPI_INTEGER,parts%Rindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) CALL MPI_Scatterv(parts%Zindex,sendcounts, displs,MPI_INTEGER,parts%Zindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) + CALL MPI_Scatterv(parts%partindex,sendcounts, displs,MPI_INTEGER,parts%partindex,Nplocs_all(mpirank),MPI_INTEGER,0,MPI_COMM_WORLD,ierr) WRITE(*,*) mpirank, "I am not the root process" END IF CALL MPI_Bcast(V,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(qsim,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) CALL MPI_Bcast(msim,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr) END IF ! Compute Self consistent electric field CALL bound CALL localisation ! Init Electric and Magnetic Fields CALL fields_init CALL rhscon CALL poisson ! END SUBROUTINE resume diff --git a/wk/stable.in b/wk/stable.in index b1bb615..a129529 100644 --- a/wk/stable.in +++ b/wk/stable.in @@ -1,45 +1,46 @@ A Skeleton for MPI Time-Dependent program ========================================= T.M. Tran SPC/EPFL - &BASIC job_time=144000.0, extra_time=360.0, - nrun=5, !# of steps - nlres=f, + nrun=600000, !# of steps + nlres=t, femorder=2,2, ngauss=4,4, nlppform=.TRUE. nlxg=f, ! Display graphical interface nlclassical=t, ! Solve classical equations of motion dt=8E-12, !timestep [s] nz=128, !# of intervals in z nnr=10,30,30 !# of intervals between radii(1) and radii(2), between radii(2) and radii(3) and between radii(3) and radii(4) - it0d=1, - it2d=10, - itparts=30000, + it0d=10, + it2d=100, + ittext=100, + itparts=1000, itgraph=10, - resfile='teststable126.h5' + resfile='Dav4.h5' radii=0.0,4e-3,25e-3,0.16 !radius of inner metallic wall, plasma boundaries and outer metallic wall plasmadim=-0.16, 0.16, 0.0070, 0.0074, !zmin zmax rmin rmax of initial plasma loading - distribtype=2 ! 1: uniform RZ gaussian in V, 2: stable eq 4.85 from Davidson + distribtype=4 ! 1: uniform RZ gaussian in V, 2: stable eq 4.85 from Davidson lz=-0.32,0.32, !Cylinder length nplasma=264600 !132300, ! 100000# of superparticles - nblock=60, - nlPhis=f, ! Deactivate calculation of self electric field + nblock=1000, + nlPhis=t, ! Deactivate calculation of self electric field ! nplasma=100 potinn=0, !potential at inner wall [V] potout=0, !30000, !potential at outer wall [V] B0=0.21, !Davidson chapter 4 formula 4.89 [T] Rcurv=1.5,!1.001, !Davidson chapter 4 formula 4.89 width=0.64, n0=-5e14, dt=8E-13 !Initial plasma density in plasmadim ! n0=-4E16, dt=8E-12 ! n0=-5.e15, dt=8E-12, ! n0=-3.e19, dt=8e-13, ! n0=-5.e8 temp=10000, !Initial temperature of plasma [K] (thermal velocity) H0=3.2e-14, - P0=8.78e-25, + P0=8.66e-25, /