diff --git a/matlab/fft_espic2d.m b/matlab/fft_espic2d.m index bc2d44d..bfd4922 100644 --- a/matlab/fft_espic2d.m +++ b/matlab/fft_espic2d.m @@ -1,139 +1,144 @@ -file='teststablegauss_13.h5'; +file='teststable_13_1058400.h5'; %close all hidden; if (~exist('M','var')) M.file=file; end M=load_espic2d(file,M,'fields'); %close all hidden; figure; ax1=gca; pcolor(ax1,M.zgrid,M.rgrid,M.N(:,:,end)); hold on [r,z]=find(M.N(:,:,end)~=0); xlim(ax1,[M.zgrid(min(z)) M.zgrid(max(z))]) ylim(ax1,[M.rgrid(min(r)) M.rgrid(max(r))]) xlabel(ax1,'Z [m]') ylabel(ax1,'R [m]') c = colorbar(ax1); c.Label.String= 'n [m^{-3}]'; view(ax1,2) %set(ax1,'colorscale','log') -stepmin=900; +stepmin=700; -[x,y]=ginput(1); -nz=find(x>M.zgrid,1,'last'); -nr=find(y>M.rgrid,1,'last'); +% [x,y]=ginput(1); +% nz=find(x>M.zgrid,1,'last'); +% nr=find(y>M.rgrid,1,'last'); -% nr=35; -% nz=65; +%nr=20;nz=129; +nr=19; nz=76; z=(M.zgrid(nz)+M.zgrid(nz+1))/2; r=(M.rgrid(nr)+M.rgrid(nr+1))/2; plot(z,r,'rx','Markersize',12); dt=M.dt*double(M.it1); FS=1/dt; L=size(M.N,3)-stepmin; f=FS*double(0:L/2)/L/1e9; %S=M.N(nr,nz,stepmin:end); S=M.pot(nr,nz,stepmin:end); S=detrend(S(:),0); S=S/max(S); S=lowpass(S(:),.499*FS,FS); P2=abs(fft(S(:)))/L; P1=P2(1:floor(L/2)+1); P1(2:end-1) = 2*P1(2:end-1); tmin=stepmin*dt; tmax=dt*size(M.N,3); fce=M.qe/M.me/2/pi*M.B(nz,nr)/1e9; N=M.N(:,:,end); -max(max(N)) -mean(mean(N(N>0))) -fpemin=sqrt(min(min(N(N>0)))*M.qe^2/(M.eps_0*M.me))/(2*pi*1e9); +fpe=sqrt(abs(M.n0)*M.qe^2/(M.eps_0*M.me))/(2*pi*1e9); +fpemin=sqrt(min(min(N(N>0.1e14)))*M.qe^2/(M.eps_0*M.me))/(2*pi*1e9); fpemax=sqrt(max(max(N))*M.qe^2/(M.eps_0*M.me))/(2*pi*1e9); +freplus=-sign(M.qsim)*fce/2*(1+sqrt(1-2*fpe^2/fce^2)); +freminus=-sign(M.qsim)*fce/2*(1-sqrt(1-2*fpe^2/fce^2)); fig=figure(); subplot(2,1,1) pl=plot(f(1:end),P1(1:end),'Displayname','Spectrum'); title(sprintf('r=%3.2e[m],z=%3.2e[m], t=[%3.2e,%3.2e][s]',r,z,tmin,tmax)) xlabel('f [GHz]','Interpreter','latex','fontsize',14) ylabel('$\left|\hat{\Phi}(f)\right|$ [a.u.]','Interpreter','latex','fontsize',14) grid on %xlim([0 10]) ylimits=ylim(); hold on plot(fce*[1 1],ylimits,'r--','Displayname','f_{ce}','linewidth',1.2); -fill([fpemin fpemax fpemax fpemin],[ylimits(1) ylimits(1) ylimits(2) ylimits(2)],'r','FaceColor',[1 1 1]*0.5,'FaceAlpha',.3,'Displayname','f_{pe}'); +plot(freminus*[1 1],ylimits,'k-.','Displayname','f_{re}^-','linewidth',1.2); +plot(freplus*[1 1],ylimits,'k-','Displayname','f_{re}^+','linewidth',1.2); +fill([fpemin fpemax fpemax fpemin],[ylimits(1) ylimits(1) ylimits(2) ylimits(2)],'r','FaceColor',[1 0 0]*0.5,'FaceAlpha',.3,'Displayname','f_{pe}'); legend uistack(pl,'top'); subplot(2,1,2) Signal=M.pot(nr,nz,1:end); %Signal=M.N(nr,nz,stepmin:end); pl=plot(dt*((1:length(Signal))-1),Signal(:)); xlabel('t [s]','Interpreter','latex','fontsize',14) ylabel('$\Phi(t)$ [V]','Interpreter','latex','fontsize',14) ylimits=ylim(); hold on x=[stepmin-1 length(Signal)-1 length(Signal)-1 stepmin-1]*dt; y=[ylimits(1) ylimits(1) ylimits(2) ylimits(2)]; fill(x,y,'r','FaceColor',[1 1 1]*0.5,'FaceAlpha',.3) uistack(pl,'top'); fig.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); fig.PaperUnits='centimeters'; fig.PaperSize=[16 10]; ylim(ylimits) print(fig,sprintf('%s_fftPhi_R%dZ%d',name,nr,nz),'-dpdf','-fillpage') -fce=M.qe/M.me/2/pi*M.B(nz,nr)/1e9; -fpe=sqrt(M.N(nr,nz,end)*M.qe^2/(M.eps_0*M.me))/(2*pi*1e9); + + S=M.N(nr,nz,stepmin:end); S=detrend(S(:),0); S=S/max(S); S=lowpass(S(:),.499*FS,FS); P2=abs(fft(S(:)))/L; P1=P2(1:floor(L/2)+1); P1(2:end-1) = 2*P1(2:end-1); fig=figure(); subplot(2,1,1) pl=plot(f(1:end),P1(1:end),'Displayname','Spectrum'); title(sprintf('r=%3.2e[m],z=%3.2e[m], t=[%3.2e,%3.2e][s]',r,z,tmin,tmax)) xlabel('f [GHz]','Interpreter','latex','fontsize',14) ylabel('$\left|\hat{n}(f)\right|$ [a.u.]','Interpreter','latex','fontsize',14) grid on %xlim([0 10]) ylimits=ylim(); hold on plot(fce*[1 1],ylimits,'r--','Displayname','f_{ce}','linewidth',1.2); -fill([fpemin fpemax fpemax fpemin],[ylimits(1) ylimits(1) ylimits(2) ylimits(2)],'r','FaceColor',[1 1 1]*0.5,'FaceAlpha',.3,'Displayname','f_{pe}'); +plot(freminus*[1 1],ylimits,'k-.','Displayname','f_{re}^-','linewidth',1.2); +plot(freplus*[1 1],ylimits,'k-','Displayname','f_{re}^+','linewidth',1.2); +fill([fpemin fpemax fpemax fpemin],[ylimits(1) ylimits(1) ylimits(2) ylimits(2)],'r','FaceColor',[1 0 0]*0.5,'FaceAlpha',.3,'Displayname','f_{pe}'); legend uistack(pl,'top'); subplot(2,1,2) Signal=M.N(nr,nz,1:end); %Signal=M.N(nr,nz,stepmin:end); pl=plot(dt*((1:length(Signal))-1),Signal(:)); xlabel('t [s]','Interpreter','latex','fontsize',14) ylabel('$n(t)$ [m$^{-3}$]','Interpreter','latex','fontsize',14) ylimits=ylim(); hold on x=[stepmin-1 length(Signal)-1 length(Signal)-1 stepmin-1]*dt; y=[ylimits(1) ylimits(1) ylimits(2) ylimits(2)]; fill(x,y,'r','FaceColor',[1 1 1]*0.5,'FaceAlpha',.3) uistack(pl,'top'); fig.PaperOrientation='landscape'; [~, name, ~] = fileparts(M.file); fig.PaperUnits='centimeters'; fig.PaperSize=[16 10]; ylim(ylimits) print(fig,sprintf('%s_fftdens_R%dZ%d',name,nr,nz),'-dpdf','-fillpage') \ No newline at end of file diff --git a/matlab/load_espic2d.m b/matlab/load_espic2d.m index 82a6bac..7fc049f 100644 --- a/matlab/load_espic2d.m +++ b/matlab/load_espic2d.m @@ -1,229 +1,237 @@ function Data = load_espic2d( file, Data, type, partstime) % 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 if(nargin<4) partstime=[0 1]; 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 = 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 = 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'); catch Data.it0 = h5readatt(file,'/data/input.00/','it0'); Data.it1 = h5readatt(file,'/data/input.00/','it1'); Data.it1 = h5readatt(file,'/data/input.00/','it2'); end Data.omepe=sqrt(abs(Data.n0)*Data.qe^2/(Data.me*Data.eps_0)); Data.omece=Data.qe*Data.B0/Data.me; % Normalizations Data.tnorm=1/Data.omepe; Data.rnorm=Data.vlight*Data.tnorm; Data.bnorm=Data.B0; Data.enorm=Data.vlight*Data.bnorm; Data.phinorm=Data.enorm*Data.rnorm; Data.vnorm=Data.vlight; % Grid data Data.rgrid= 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 = h5read(file,'/data/fields/Br')*Data.bnorm; Data.Br= reshape(Br,length(Data.zgrid),length(Data.rgrid)); Bz = h5read(file,'/data/fields/Bz')*Data.bnorm; Data.Bz= reshape(Bz,length(Data.zgrid),length(Data.rgrid)); try Atheta = h5read(file,'/data/fields/Athet')*Data.bnorm; Data.Athet= reshape(Atheta,length(Data.zgrid),length(Data.rgrid)); catch end Data.B=sqrt(Data.Bz.^2+Data.Br.^2); clear Br Bz 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 = 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 = h5read(file,'/data/fields/pot')*Data.phinorm; Data.pot= reshape(pot,length(Data.zgrid),length(Data.rgrid),size(pot,2)); clear pot Data.pot=permute(Data.pot,[2,1,3]); Er = h5read(file,'/data/fields/Er')*Data.enorm; Data.Er= reshape(Er,length(Data.zgrid),length(Data.rgrid),size(Er,2)); clear Er Data.Er=permute(Data.Er,[2,1,3]); Ez = h5read(file,'/data/fields/Ez')*Data.enorm; Data.Ez= reshape(Ez,length(Data.zgrid),length(Data.rgrid),size(Ez,2)); clear Ez Data.Ez=permute(Data.Ez,[2,1,3]); 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; Vcell=(Data.zgrid(2:end)-Data.zgrid(1:end-1))*((Data.rgrid(2:end).^2-Data.rgrid(1:end-1).^2)*pi)'; for i=1:size(Data.N,3) Data.N(:,:,i)=abs(Data.qsim/Data.qe)*Data.N(:,:,i)./Vcell'; end Data.N=padarray(Data.N,[1,1],0,'post'); clear Vcell clear N Data.fluidUR=reshape(h5read(file,'/data/fields/fluidur'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]); Data.fluidUR=permute(Data.fluidUR,[2,1,3]); Data.fluidUR=padarray(Data.fluidUR,[1,1],0,'post'); Data.fluidUZ=reshape(h5read(file,'/data/fields/fluiduz'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]); Data.fluidUZ=permute(Data.fluidUZ,[2,1,3]); Data.fluidUZ=padarray(Data.fluidUZ,[1,1],0,'post'); Data.fluidUTHET=reshape(h5read(file,'/data/fields/fluiduthet'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]); Data.fluidUTHET=permute(Data.fluidUTHET,[2,1,3]); Data.fluidUTHET=padarray(Data.fluidUTHET,[1,1],0,'post'); Data.Presstens=reshape(h5read(file,'/data/fields/presstens'),6,length(Data.zgrid)-1,length(Data.rgrid)-1,[]); Data.Presstens=permute(Data.Presstens,[1,3,2,4]); - Data.t2d=Data.dt*(0:size(Data.N,3)-1)*double(Data.it1); + try + Data.t2d = h5read(file,'/data/fields/time'); + catch + Data.t2d=Data.dt*(0:size(Data.N,3)-1)*double(Data.it1); + end if(strcmp(Data.type,'parts')) Data.type='all'; end end if(strcmp(type,'all') || strcmp(type,'parts')) Data.R = h5read(file,'/data/part/R')*Data.rnorm; partsbegin=floor(partstime(1)*size(Data.R,2))+1; partsend=floor(partstime(2)*size(Data.R,2)); Data.R=Data.R(:,1:partsend); Data.Z = h5read(file,'/data/part/Z',[1,partsbegin],[size(Data.R,1) partsend])*Data.rnorm; try Data.THET = h5read(file,'/data/part/THET',[1,partsbegin],[size(Data.R,1) partsend]); catch clear Data.THET end try Data.Rindex=h5read(file,'/data/part/Rindex',[1,partsbegin],[size(Data.R,1) partsend]); Data.Zindex=h5read(file,'/data/part/Zindex',[1,partsbegin],[size(Data.R,1) partsend]); catch clear Data.Rindex Data.Zindex end UR = h5read(file,'/data/part/UR',[1,partsbegin],[size(Data.R,1) partsend]); UZ = h5read(file,'/data/part/UZ',[1,partsbegin],[size(Data.R,1) partsend]); UTHET= h5read(file,'/data/part/UTHET',[1,partsbegin],[size(Data.R,1) partsend]); if(Data.nlclassical) Data.VR = UR*Data.vnorm; clear UR Data.VZ = UZ*Data.vnorm; clear UZ Data.VTHET = UTHET*Data.vnorm; clear UTHET else GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2); Data.VR = UR./GAMM*Data.vnorm; clear UR Data.VZ = UZ./GAMM*Data.vnorm; clear UZ Data.VTHET = UTHET./GAMM*Data.vnorm; clear UTHET clear GAMM end %Data.VPERP = sqrt(Data.VR.*Data.VR+Data.VTHET.*Data.VTHET); partepot = sign(Data.qsim)*h5read(file,'/data/part/pot',[1,partsbegin],[size(Data.R,1) partsend])*Data.phinorm*Data.qe; Data.H=0.5*Data.me*(Data.VR.^2+Data.VTHET.^2+Data.VZ.^2)+partepot; clear partepot Data.P=Data.R.*(Data.VTHET*Data.me+sign(Data.qsim)*Data.qe*Athet(Data.R,Data.Z)); try partindex = h5read(file,'/data/part/partindex',[1,partsbegin],[size(Data.R,1) partsend]); [~,Indices]=sort(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.THET(:,i)=Data.THET(Indices(:,i),i); end catch end clear partindex; - Data.tpart=Data.dt*(0:size(Data.R,2)-1)*double(Data.it2); + try + Data.tpart = h5read(file,'/data/part/time'); + catch + Data.tpart=Data.dt*(0:size(Data.R,2)-1)*double(Data.it2); + end if(strcmp(Data.type,'fields')) Data.type='all'; end end Data.t0d=Data.dt.*double(0:length(Data.epot)-1); filedata=dir(file); Data.timestamp=filedata.date; disp('M reloaded') function Atheta=Athet(R,Z) Atheta=0.5*Data.B0*(R-Data.width/pi*(Data.Rcurv-1)/(Data.Rcurv+1)... .*besseli(1,2*pi*R/Data.width).*cos(2*pi*Z/Data.width)); end end diff --git a/src/basic_mod.f90 b/src/basic_mod.f90 index 6702376..ed1b91b 100644 --- a/src/basic_mod.f90 +++ b/src/basic_mod.f90 @@ -1,306 +1,306 @@ 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=.FALSE. ! 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, DIMENSION(:), ALLOCATABLE:: fluidur, fluiduz, fluiduthet ! Fluid velocities per gridcell computed every it2d DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE:: Presstens ! Pressure tensor 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 + & resfile, rstfile, 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/diagnose.f90 b/src/diagnose.f90 index a08c23a..11d9ebe 100644 --- a/src/diagnose.f90 +++ b/src/diagnose.f90 @@ -1,224 +1,228 @@ 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) + INTEGER :: partdensitystatus(5*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/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, 0, SHAPE(time), "/data/part/time", "time" ) + 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") CALL creatd(fidres, 2, SHAPE(Presstens), "/data/fields/presstens", "presstens") + CALL creatd(fidres, 0, SHAPE(time), "/data/fields/time", "time" ) 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) CALL putarr(fidres, "/data/fields/Athet", Athet) ! ! 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/time", time) 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 MPI_WAITALL(5, ireducerequest(3:7), 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) CALL append(fidres, "/data/fields/presstens", Presstens) 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/time", time) 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