diff --git a/Matlab_Analysis/KORC_postproc/plot_GC.m b/Matlab_Analysis/KORC_postproc/plot_GC.m new file mode 100644 index 00000000..54d4b4fb --- /dev/null +++ b/Matlab_Analysis/KORC_postproc/plot_GC.m @@ -0,0 +1,827 @@ +format longE + +%% load data from file + +join=0; +if(join==1) + + load('../TEST29d/29d/ST.mat') + ST29d=ST; + load('../TEST29d/29d1/ST.mat') + ST29d1=ST; + + ST.time=cat(2,ST29d.time,ST29d1.time); + + ST.data.sp1.Y=cat(3,ST29d.data.sp1.Y,ST29d1.data.sp1.Y); + ST.data.sp1.V=cat(3,ST29d.data.sp1.V,ST29d1.data.sp1.V); + ST.data.sp1.B=cat(3,ST29d.data.sp1.B,ST29d1.data.sp1.B); + ST.data.sp1.E=cat(3,ST29d.data.sp1.E,ST29d1.data.sp1.E); + ST.data.sp1.g=cat(2,ST29d.data.sp1.g,ST29d1.data.sp1.g); + ST.data.sp1.eta=cat(2,ST29d.data.sp1.eta,ST29d1.data.sp1.eta); + ST.data.sp1.flag=cat(2,ST29d.data.sp1.flag,ST29d1.data.sp1.flag); + ST.data.sp1.PSIp=cat(2,ST29d.data.sp1.PSIp,ST29d1.data.sp1.PSIp); + ST.data.sp1.ne=cat(2,ST29d.data.sp1.ne,ST29d1.data.sp1.ne); + +end + +%% extract data from ST structure + +c=ST.params.scales.v; +me=ST.params.scales.m; +qe=ST.params.scales.q; +e0=8.8542e-12; + +time=ST.time; + +R=squeeze(ST.data.sp1.Y(:,1,:)); +PHI=squeeze(ST.data.sp1.Y(:,2,:)); +Z=squeeze(ST.data.sp1.Y(:,3,:)); + +%R0=ST.params.fields.Ro; +%Z0=ST.params.fields.Zo; +R0=ST.params.species.Ro; +Z0=ST.params.species.Zo; +rm=sqrt((R-R0).^2+(Z-Z0).^2); +theta=atan2((Z-Z0),(R-R0)); + +BR=squeeze(ST.data.sp1.B(:,1,:)); +BPHI=squeeze(ST.data.sp1.B(:,2,:)); +BZ=squeeze(ST.data.sp1.B(:,3,:)); + +EPHI=squeeze(ST.data.sp1.E(:,2,:)); + +%PSIp=ST.data.sp1.PSIp; + +%lam=ST.params.fields.lambda; +%q0=ST.params.fields.qo; +%B0=ST.params.fields.Bo; +%PSIanaly=R*lam^2*B0./(2*q0*(R0+rm.*cos(theta))).*log(1+(rm/lam).^2); + +BMAG=sqrt(BR.^2+BPHI.^2+BZ.^2); + +ppll=squeeze(ST.data.sp1.V(:,1,:)); +mu=squeeze(ST.data.sp1.V(:,2,:)); + +%momgcv=ppll.*R.*BPHI./BMAG; +%momgcb=qe*PSIp; +%momgc=momgcv+momgcb; + +pmag=sqrt(ppll.^2+mu.*BMAG*2*me); + +gam=ST.data.sp1.g(:,:); +eta=ST.data.sp1.eta(:,:); +xi=cos(deg2rad(eta)); + +vpll=ppll./(gam*me); + +flagCon=ST.data.sp1.flagCon; +flagCol=ST.data.sp1.flagCol; +%flagRE=ST.data.sp1.flagRE; + +Ipart=qe*vpll.*BPHI./BMAG; +Ipart(flagCon==0)=0; +Ipart(flagCol==0)=0; +I=sum(Ipart,1); + +gam(flagCon==0)=0; +gam(flagCol==0)=0; +E=me*c^2/qe*sum(gam); +KE=me*c^2/qe*sum((gam-1).*flagCon.*flagCol); + +AveE=me*c^2/qe*sum(gam.*flagCon.*flagCol)/size(R,1); +AveE(isnan(AveE))=0; +AveKE=me*c^2/qe*sum((gam-1).*flagCon.*flagCol)/size(R,1); +AveKE(isnan(AveKE))=0; + +dgamdt=-qe^4*BMAG.^2/(6*pi*e0*(me*c)^3).*(gam.^2-1).*(1-xi.^2); +dgam=dgamdt*ST.params.simulation.dt; + +gamanaly=zeros(size(time)); +gamanaly(1)=gam(1,1); +for i=2:size(time,2) + gamanaly(i)=gamanaly(i-1)+dgam(1,i-1); +end + +dpplldt=-ppll.*(1-xi.^2).*(gam-1./gam)./(6*pi*e0*(me*c)^3./(qe^4*BMAG.^2)); +dmudt=-2*mu./(6*pi*e0*(me*c)^3./(qe^4*BMAG.^2)).*(gam.*(1-xi.^2)+xi.^2./gam); + +%ne=ST.data.sp1.ne; + +%% construct background toroidal magnetic vector potential + +if (strcmp(ST.params.simulation.field_model,'ANALYTICAL')) + lam=ST.params.fields.lambda; + R0=ST.params.fields.Ro; + Z0=ST.params.fields.Zo; + a=ST.params.fields.a; + B0=ST.params.fields.Bo; + q0=ST.params.fields.qo; + + RF=linspace(R0-a,R0+a,100); + ZF=linspace(Z0-a,Z0+a,100); + + rm=zeros(100); + PSIP=zeros(100); + FLAG=ones(100); + for ii=1:100 + for jj=1:100 + rm(ii,jj)=sqrt((RF(ii)-R0)^2+(ZF(jj)-Z0)^2); + PSIP(ii,jj)=lam^2*B0./(2*q0)*log(1+(rm(ii,jj)/lam)^2); + if rm(ii,jj)>a + FLAG(ii,jj)=0; + end + end + end + +elseif (strcmp(ST.params.simulation.field_model{1},'EXTERNAL') &&... + ST.params.fields.Axisymmetric==0) + PSIP=ST.params.fields.psi_p; + FLAG=ST.params.fields.Flag; + + RF=ST.params.fields.R; + ZF=ST.params.fields.Z; + + RF1=RF(1); + ZF1=ZF(1); + + skipres=2^0; + + for i=skipres:skipres:size(RF) + RF1=cat(1,RF1,RF(i)); + end + for i=skipres:skipres:size(ZF) + ZF1=cat(1,ZF1,ZF(i)); + end + +elseif (strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') &&... + ST.params.fields.Axisymmetric==0) + + RF=ST.params.fields.R; + %PHIF=ST.params.fields.PHI; + ZF=ST.params.fields.Z; + + %PSIP=zeros(size(RF,1),size(PHIF,1),size(ZF,1)); + + filename='/Users/21b/Desktop/M3D-C1/it1ip1_p32.rs16/I_psi_0.h5'; + PSI=h5read(filename,'/PSI'); + + PSIP1=squeeze(PSI(1,:,:)); + FLAG=squeeze(ST.params.fields.Flag(:,1,:)); + + PHI=mod(PHI,2*pi); + + %[RR,PP,ZZ]=meshgrid(RF,PHIF,ZF); + [RR,ZZ]=meshgrid(RF,ZF); + %PSIp_part=interp3(RR,PP,ZZ,PSIP1,R,PHI,Z); + PSIp_part=interp2(RR,ZZ,PSIP1',R,Z); + + momgcb=-qe*PSIp_part; + momgc=momgcv+momgcb; + + PSIP=PSIP1; + +elseif (strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') &&... + ST.params.fields.Axisymmetric==1 && isfield(ST.params.fields,'psi_p')) + + PSIP=ST.params.fields.psi_p; + FLAG=ST.params.fields.Flag2D; + + RF=ST.params.fields.R; + ZF=ST.params.fields.Z; + +elseif (strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') &&... + ST.params.fields.Axisymmetric==1 && isfield(ST.params.fields,'psi_p3D')) + + PSIP=ST.params.fields.psi_p3D; + FLAG=ST.params.fields.Flag3D; + + RF=ST.params.fields.R; + TF=ST.params.fields.PHI; + ZF=ST.params.fields.Z; + +end + + + +%% plot data +%close all + +histpeta=0; +histreta=0; +histp=0; +histpmulti=0; +histxi=0; +histRZ=0; +evoEave=0; +evoI=0; +evomomgc=0; +momgccomp=0; +peta_RE_combo=1; +multiRZ=0; +NREevo=0; +fRE=0; + +timeslice=size(time,2); + +%figure; +%hold on +%scatter(R(:,timeslice),Z(:,timeslice),'o') +%hold off +%daspect([1,1,1]) + +if fRE==1 + + fig=figure; + + KEp=me*c^2*(gam-1)/qe; + + KEbins=10.^(linspace(log10(min(KEp/(10^6))),log10(max(KEp/(10^6))),50)); + + [histKE,edgesKE]=histcounts(KEp(flagCon(:,timeslice)>0,timeslice)./(10^6),... + KEbins); + + hist_binsKE=zeros(size(histKE)); + for i=1:size(histKE,2) + hist_binsKE(i)=(edgesKE(i)+edgesKE(i+1))/2; + end + + binwidth=diff(edgesKE); + normKE=trapz(hist_binsKE,histKE); + + + loglog(hist_binsKE,histKE/normKE./binwidth,'Linewidth',2) + xlabel('E [MeV]') + ylabel('f_E') + + fig=figure; + loglog(hist_binsKE,histKE,'Linewidth',2) + xlabel('E [MeV]') + ylabel('N_{RE}') +end + +if NREevo==1 + + E0=ST.params.fields.Eo; + %E_CH=ST.params.collisions_ss.Ec; + %tau_c=ST.params.collisions_ss.Tau; + %Clog0=ST.params.collisions_ss.Clogee; + %Zeff_c=ST.params.collisions_ss.Zeff; + %dt_c=ST.params.simulation.dt*ST.params.collisions_ss.subcycling_iterations; + %ne_c=ST.params.collisions_ss.ne; + + + RPgrowthrate=((E0/E_CH)-1)/(2*tau_c*Clog0); + + NREtot=sum(flagRE); + NREact=sum(flagRE.*flagCol); + + startind=int64(.5*size(time,2)); + %startind=1; + + myExp='a*exp(b*x)'; + startPointsE=[NREtot(startind),RPgrowthrate]; + startPointsA=[NREact(startind),RPgrowthrate]; + + + f1=fit(time(startind:end)',NREtot(startind:end)',myExp,'Start',startPointsE); + f2=fit(time(startind:end)',NREact(startind:end)',myExp,'Start',startPointsA); + + fig=figure; + p1=semilogy(time,NREtot,'linewidth',3,'color',[0,0,1]); + hold on + p1a=semilogy(time,f1.a*exp(time*f1.b),'--','linewidth',3,'color',[0,0,1]); + p2=semilogy(time,NREact,'linewidth',3,'color',[1,0,0]); + p2a=semilogy(time,f2.a*exp(time*f2.b),'--','linewidth',3,'color',[1,0,0]); + + legend([p1,p1a,p2,p2a],{'All RE','All RE fit','Active RE','Active RE fit'},... + 'numcolumns',2,'Location','southeast') + + xlabel('Time(s)') + ylabel('N_{\rm RE}') + + %txt={strcat('Zeff: ',num2str(Zeff_c)),... + % strcat('dt=',num2str(dt_c),'(s)'),... + % strcat('n_e= ',num2str(ne_c),'m^{-3}'),... + % strcat('E/E_{CH}: ',num2str(E0/E_CH)),... + % strcat('\Gamma_{\rm RP}: ',num2str(RPgrowthrate)),... + % strcat('All RE fit \gamma: ',num2str(f1.b)),... + % strcat('Active RE fit \gamma: ',num2str(f2.b))}; + %text([.025],[.7],txt,'FontSize',14,'Units','normalized') + + set(gca,'Fontsize',20) + + print(fig,'NREevo','-depsc') + +end + +if multiRZ==1 + + + + fig=figure('Renderer', 'painters', 'Position', [10 10 700 400]); + + + + subplot('Position',[.05, .1, .2, .65]) + box on + hold on + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,PSIP',20,'k','LineWidth',1) + end + + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,FLAG',1,'k','LineWidth',2) + end + + hist0=histogram2(Rj(flag(:,timeslice)>0,timeslice),Zj(flag(:,timeslice)>0,timeslice),30,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + hold off + + set(gca,'Fontsize',12) + + cb=colorbar(); + cb.Label.String='$N_{\rm RE}$'; + cb.Label.Interpreter='latex'; + cb.Label.FontSize=14; + + xlabel('$R\,{\rm (m)}$','Interpreter','latex','FontSize',14) + ylabel('$Z\,{\rm (m)}$','Interpreter','latex','FontSize',14) + title('$\psi_{N,{\rm max}}=0.8446$','interpreter','latex') + daspect([1,1,1]) + yticks([-1,0,1]); + + colormap(gca,'jet') + + subplot('Position',[.275, .1, .2, .65]) + box on + hold on + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,PSIP',20,'k','LineWidth',1) + end + + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,FLAG',1,'k','LineWidth',2) + end + + hist2=histogram2(Rj2(flag(:,timeslice)>0,timeslice),Zj2(flag(:,timeslice)>0,timeslice),... + 'XBinEdges',hist0.XBinEdges,'YBinEdges',hist0.YBinEdges,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + hold off + + set(gca,'Fontsize',12) + + + cb=colorbar(); + cb.Label.String='$N_{\rm RE}$'; + cb.Label.Interpreter='latex'; + cb.Label.FontSize=14; + + xlabel('$R\,{\rm (m)}$','Interpreter','latex','FontSize',14) + %ylabel('$Z\,{\rm (m)}$','Interpreter','latex','FontSize',14) + title('$\psi_{N,{\rm max}}=0.63345$','interpreter','latex') + daspect([1,1,1]) + %yticks([-1,0,1]); + yticklabels([]); + + colormap(gca,'jet') + + subplot('Position',[.5, .1, .2, .65]) + box on + hold on + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,PSIP',20,'k','LineWidth',1) + end + + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,FLAG',1,'k','LineWidth',2) + end + + hist1=histogram2(Rj1(flag(:,timeslice)>0,timeslice),Zj1(flag(:,timeslice)>0,timeslice),... + 'XBinEdges',hist0.XBinEdges,'YBinEdges',hist0.YBinEdges,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + hold off + + set(gca,'Fontsize',12) + + + cb=colorbar(); + cb.Label.String='$N_{\rm RE}$'; + cb.Label.Interpreter='latex'; + cb.Label.FontSize=14; + + xlabel('$R\,{\rm (m)}$','Interpreter','latex','FontSize',14) + %ylabel('$Z\,{\rm (m)}$','Interpreter','latex','FontSize',14) + title('$\psi_{N,{\rm max}}=0.4223$','interpreter','latex') + daspect([1,1,1]) + %yticks([-1,0,1]); + yticklabels([]); + + colormap(gca,'jet') + + subplot('Position',[.725, .1, .2, .65]) + box on + hold on + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,PSIP',20,'k','LineWidth',1) + end + + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,FLAG',1,'k','LineWidth',2) + end + + hist3=histogram2(Rj3(flag(:,timeslice)>0,timeslice),Zj3(flag(:,timeslice)>0,timeslice),... + 'XBinEdges',hist0.XBinEdges,'YBinEdges',hist0.YBinEdges,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + hold off + + set(gca,'Fontsize',12) + + + cb=colorbar(); + cb.Label.String='$N_{\rm RE}$'; + cb.Label.Interpreter='latex'; + cb.Label.FontSize=14; + + xlabel('$R\,{\rm (m)}$','Interpreter','latex','FontSize',14) + %ylabel('$Z\,{\rm (m)}$','Interpreter','latex','FontSize',14) + title('$\psi_{N,{\rm max}}=0.21115$','interpreter','latex') + daspect([1,1,1]) + %yticks([-1,0,1]); + yticklabels([]); + + colormap(gca,'jet') + + orient(fig,'landscape') + print(fig,'D3D164409_initspatRZ','-dpdf') + +end + +if histpmulti==1 + fig=figure; + + pmax=20; + + timeslices=[1 20 40 60 80]; + pedges=linspace(min(pmag1(flag(:,timeslices(end))>0,... + timeslices(end))/(me*c)),pmax,100); + + hold on + + cm=colormap(jet(size(timeslices,2))); + legends={size(timeslices,2)}; + + for ii=1:size(timeslices,2) + + [histp,edgesp]=histcounts(pmag(flag(:,timeslices(ii))>0,timeslices(ii))/(me*c)... + ,pedges,'Normalization','count'); + + hist_binsp=zeros(size(histp)); + for i=1:size(histp,2) + hist_binsp(i)=(edgesp(i)+edgesp(i+1))/2; + end + + legends{ii}=strcat('t=',num2str(time(ii)),'s'); + + plot(hist_binsp,histp,'Linewidth',2,'color',cm(ii,:)) + end + + + set(gca,'Fontsize',12) + + %text(-.25,.98,'a)','FontSize',14,'Units','normalized'); + + ylabel('$N_{\rm RE}$','Interpreter','latex','FontSize',16) + xlabel('$p$ ($m_ec$)','Interpreter','latex','FontSize',16) +% axis([0,55,hist.YBinLimits(1),hist.YBinLimits(2)]) + hold off + + legend(legends,'location','northwest','fontsize',14) + + print(fig,'TEST25e_histp_evo','-dpdf') +end + +if peta_RE_combo==1 + fig=figure('Renderer', 'painters', 'Position', [10 10 750 400]); + fig.PaperUnits='points'; + fig.PaperPosition = [0 50 650 350]; + + + pfilter=pmag(:,timeslice); + etafilter=eta(:,timeslice); + Rfilter=R(:,timeslice); + Zfilter=Z(:,timeslice); + + zmax=0.1; + zmin=-0.1; + rmax=1.1; + rmin=1.0; + + Zfilter(Rfilter>rmax)=0; + Zfilter(Rfilterzmax)=0; + Rfilter(Zfilterrmax)=0; + Rfilter(Rfilterzmax)=0; + Zfilter(Zfilter0,timeslice),... + me*c^2*(gam(flagCon(:,timeslice)>0,timeslice)-1)/(10^6*qe),50,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + colormap jet + cb=colorbar(); + cb.Label.String='$N_{\rm RE}$'; + cb.Label.Interpreter='latex'; + cb.Label.FontSize=14; + + set(gca,'Fontsize',12) + + text(-.2,.98,'a)','FontSize',14,'Units','normalized'); + + xlabel('$\eta\,({\rm rad})$','Interpreter','latex','FontSize',16) + ylabel('$E\,({\rm MeV})$','Interpreter','latex','FontSize',16) +% axis([0,55,hist.YBinLimits(1),hist.YBinLimits(2)]) + hold off + + subplot(1,2,2) + + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,squeeze(PSIP(:,11,:))',50,'k','LineWidth',1) + end + hold on + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,squeeze(FLAG(:,11,:))',1,'k','LineWidth',2) + end + + +% hist=histogram2(Rfilter,Zfilter,30,... +% 'DisplayStyle','tile','LineStyle','none',... +% 'Normalization','count'); + hist=histogram2(R(flagCon(:,timeslice)>0,timeslice),Z(flagCon(:,timeslice)>0,timeslice),30,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + colormap jet + cb=colorbar(); + cb.Label.String='$N_{\rm RE}$'; + cb.Label.Interpreter='latex'; + cb.Label.FontSize=14; + + set(gca,'Fontsize',12) + + + text(-.3,.98,'b)','FontSize',14,'Units','normalized'); + + + xlabel('$R\,{\rm (m)}$','Interpreter','latex','FontSize',16) + ylabel('$Z\,{\rm (m)}$','Interpreter','latex','FontSize',16) + %axis([rmin,rmax,zmin,zmax]) + daspect([1,1,1]) + hold off + + print(fig,'D3D164409_initial_dist','-djpeg') +end + +if histpeta==1 + fig=figure; + hist=histogram2(eta(:,timeslice),pmag(:,timeslice)/(me*c),50,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + colormap jet + cb=colorbar(); + cb.Label.String='$N_{\rm RE}$'; + cb.Label.Interpreter='latex'; + cb.Label.FontSize=14; + + set(gca,'Fontsize',12) + + xlabel('Pitch Angle $\eta$','Interpreter','latex','FontSize',16) + ylabel('$p$ ($m_ec$)','Interpreter','latex','FontSize',16) +% axis([0,55,hist.YBinLimits(1),hist.YBinLimits(2)]) + hold off + +% txt={strcat('E:',num2str(ST.params.fields.Eo)),... +% strcat('Impurity:',num2str(ST.params.collisions_ms.Zj)),... +% strcat('Z0= ',num2str(ST.params.collisions_ms.Zj)),... +% strcat('dt= ',num2str(ST.params.simulation.dt)),... +% strcat('ne= ',num2str(ST.params.profiles.neo))}; +% text(.25,.9,txt,'FontSize',14,'Units','normalized'); + + print(fig,'Initial_peta_distribution','-dpdf') +end + +if histreta==1 + fig=figure; + hist=histogram2(R(:,timeslice),eta(:,timeslice),50,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + colormap jet + colorbar() + xlabel('$R (m)$','Interpreter','latex','FontSize',14) + ylabel('Pitch Angle $\theta$','Interpreter','latex','FontSize',14) +% axis([0,55,hist.YBinLimits(1),hist.YBinLimits(2)]) + hold off + + print(fig,'TEST13_V0','-dpdf') +end + +if histp==1 + fig=figure; + hist=histogram(pmag(:,timeslice)/(me*c),50,'Normalization','count'); + hold on + +% line([pmag(1,1)/(me*c),pmag(1,1)/(me*c)],[0,max(hist.Values)],... +% 'LineStyle','--','Linewidth',2.,'Color','red') + + xlabel('$p$ ($m_ec$)','Interpreter','latex','FontSize',14) + ylabel('$f(p)$','Interpreter','latex','FontSize',14) +% axis([0,21,0,max(hist.Values)]) + hold off + set(gca,'Yscale','log') + + + txt={strcat('E:',num2str(ST.params.fields.Eo)),... +% strcat('Impurity:',num2str(ST.params.collisions_ms.Zj)),... +% strcat('Z0= ',num2str(ST.params.collisions_ms.Zj)),... +% strcat('dt= ',num2str(ST.params.simulation.dt)),... + strcat('ne= ',num2str(ST.params.profiles.neo))}; + text(.25,.9,txt,'FontSize',14,'Units','normalized'); + + print(fig,'TEST13_V0','-dpdf') +end + +if histxi==1 + fig=figure; + hist=histogram(xi(flag(:,timeslice)>0,timeslice),50,'Normalization','count'); + hold on + set(gca,'Yscale','log') + +% line([xi(1,1),xi(1,1)],[0,max(hist.Values)],... +% 'LineStyle','--','Linewidth',2.,'Color','red') + + xlabel('$\xi$','Interpreter','latex','FontSize',14) + ylabel('$f(\xi)$','Interpreter','latex','FontSize',14) +% axis([-1.1,1.1,0,inf]) + hold off + + txt={strcat('E:',num2str(ST.params.fields.Eo)),... +% strcat('Impurity:',num2str(ST.params.collisions_ms.Zj)),... +% strcat('Z0= ',num2str(ST.params.collisions_ms.Zj)),... +% strcat('dt= ',num2str(ST.params.simulation.dt)),... + strcat('ne= ',num2str(ST.params.profiles.neo))}; + text(.25,.9,txt,'FontSize',14,'Units','normalized'); + + print(fig,'TEST13_V0','-dpdf') +end + +if histRZ==1 + fig=figure; + + if isfield(ST.params.fields,'psi_p3D') + t0=1.405; + tt=t0+time(timeslice); + + [TT,RR,ZZ]=meshgrid(TF,RF,ZF); + + RRF=squeeze(RR(:,1,:)); + ZZF=squeeze(ZZ(:,1,:)); + TTF=tt*ones(size(RRF)); + + PSIPtt=interp3(TT,RR,ZZ,PSIP,TTF,RRF,ZZF); + FLAGtt=interp3(TT,RR,ZZ,FLAG,TTF,RRF,ZZF); + else + PSIPtt=PSIP; + FLAGtt=FLAG; + end + + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,PSIPtt',50,'k','LineWidth',1,... + 'LineStyle','--') + end + hold on + if strcmp(ST.params.simulation.field_model{1}(1:8),'EXTERNAL') + contour(RF,ZF,FLAGtt',1,'k','LineWidth',2) + end + hist=histogram2(R(flagCon(:,timeslice)>0 & flagCol(:,timeslice)>0,timeslice),... + Z(flagCon(:,timeslice)>0 & flagCol(:,timeslice)>0,timeslice),50,... + 'DisplayStyle','tile','LineStyle','none',... + 'Normalization','count'); + colormap jet + cb=colorbar(); + cb.Label.String='$N_{\rm RE}$'; + cb.Label.Interpreter='latex'; + cb.Label.FontSize=14; + + set(gca,'Fontsize',12) + + xlabel('$R\,{\rm (m)}$','Interpreter','latex','FontSize',16) + ylabel('$Z\,{\rm (m)}$','Interpreter','latex','FontSize',16) +% axis([1.4,1.85,-.175,.175]) + daspect([1,1,1]) + + print(fig,'Initial_RZ_distribution','-dpdf') +end + +if evoEave==1 + fig=figure; + plot(time,AveE) + xlabel('${\rm t\,(s)}$','Interpreter','latex','FontSize',14) + ylabel('${\rm E_{\rm ave}\,(eV)}$','Interpreter','latex','FontSize',14) + + txt={strcat('E:',num2str(ST.params.fields.Eo)),... +% strcat('Impurity:',num2str(ST.params.collisions_ms.Zj)),... +% strcat('Z0= ',num2str(ST.params.collisions_ms.Zj)),... +% strcat('dt= ',num2str(ST.params.simulation.dt)),... + strcat('ne= ',num2str(ST.params.profiles.neo))}; + text(.25,.9,txt,'FontSize',14,'Units','normalized'); + + print(fig,'TEST13_evoEave','-dpdf') + +end + +if evoI==1 + fig=figure; + plot(time,I/I(1)) + xlabel('${\rm t\,(s)}$','Interpreter','latex','FontSize',14) + ylabel('${\rm I/I(0)}$','Interpreter','latex','FontSize',14) + + txt={strcat('E:',num2str(ST.params.fields.Eo)),... +% strcat('Impurity:',num2str(ST.params.collisions_ms.Zj)),... +% strcat('Z0= ',num2str(ST.params.collisions_ms.Zj)),... +% strcat('dt= ',num2str(ST.params.simulation.dt)),... + strcat('ne= ',num2str(ST.params.profiles.neo))}; + text(.25,.9,txt,'FontSize',14,'Units','normalized'); + + print(fig,'evoI','-dpdf') + +end + +if evomomgc==1 + fig=figure; + + momnorm=zeros(size(momgc)); + for i=1:size(momgc,1) + momnorm(i,:)=(momgc(i,:)-momgc(i,1))/momgc(i,1); + end + + avemomnorm=sum(momnorm.*flag)./sum(flag); + + plot(time,momnorm(1,:)) + hold on + line([0,time(end)],[0,0],'Color','k','LineStyle',':') + xlabel('${\rm t\,(s)}$','Interpreter','latex','FontSize',14) + ylabel('${\rm \Delta p_\zeta/p_\zeta(0)}$','Interpreter','latex','FontSize',14) + hold off + + txt={strcat('dt= ',num2str(ST.params.simulation.dt))}; + text(.25,.9,txt,'FontSize',14,'Units','normalized'); + + print(fig,'evoI','-dpdf') + +end + +if momgccomp==1 + fig=figure; + + + momnorm0=(momgc0(1,:)-momgc0(1,1))/momgc0(1,1); + momnorm1=(momgc1(1,:)-momgc1(1,1))/momgc1(1,1); + momnorm2=(momgc2(1,:)-momgc2(1,1))/momgc2(1,1); + momnorm3=(momgc3(1,:)-momgc3(1,1))/momgc3(1,1); + momnorm4=(momgc4(1,:)-momgc4(1,1))/momgc4(1,1); + momnorm5=(momgc5(1,:)-momgc5(1,1))/momgc5(1,1); + momnorm6=(momgc6(1,:)-momgc6(1,1))/momgc6(1,1); + momnorm7=(momgc7(1,:)-momgc7(1,1))/momgc7(1,1); + + semilogy(time0,abs(momnorm0)) + hold on + semilogy(time1,abs(momnorm1)) + semilogy(time2,abs(momnorm2)) + semilogy(time3,abs(momnorm3)) + semilogy(time4,abs(momnorm4)) + semilogy(time5,abs(momnorm5)) + semilogy(time6,abs(momnorm6)) + semilogy(time7,abs(momnorm7)) + + xlabel('${\rm t\,(s)}$','Interpreter','latex','FontSize',14) + ylabel('${\rm p_\zeta/p_\zeta(0)}$','Interpreter','latex','FontSize',14) + hold off + legend({'0','1','2','3','4','5','6','7'},'Location','Northwest') + + + print(fig,'compare_momgct','-dpdf') + +end diff --git a/ci_stub/onyx_job_test_cpu.sh b/ci_stub/onyx_job_test_cpu.sh index edeed8bd..d854cf52 100755 --- a/ci_stub/onyx_job_test_cpu.sh +++ b/ci_stub/onyx_job_test_cpu.sh @@ -4,5 +4,5 @@ env | grep JOB_COUNT JOB_COUNT=1 -cd ./KORC/build_cpu && ctest -j --output-on-failure +cd ./KORC/build_cpu && ctest --output-on-failure diff --git a/ci_stub/onyx_job_test_gpu.sh b/ci_stub/onyx_job_test_gpu.sh index ba57b941..0b267778 100755 --- a/ci_stub/onyx_job_test_gpu.sh +++ b/ci_stub/onyx_job_test_gpu.sh @@ -6,7 +6,8 @@ JOB_COUNT=1 pushd ./KORC/build_gpu >/dev/null -ctest -j --output-on-failure +ctest --output-on-failure +testexit=$? pushd ./bin >/dev/null @@ -19,6 +20,12 @@ for j in "mars_test" "egyro_test"; do done done - popd >/dev/null popd >/dev/null + +if [ $testexit -eq 0 ] +then + exit 0 +else + exit 1 +fi \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 11cfb18b..d6773d30 100755 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -13,7 +13,6 @@ target_sources(korc $ $ $ - $ $ $ $ diff --git a/src/korc_collisions.f90 b/src/korc_collisions.f90 index 05ba3fd7..cde3e84f 100755 --- a/src/korc_collisions.f90 +++ b/src/korc_collisions.f90 @@ -2443,6 +2443,13 @@ subroutine include_CoulombCollisions_GC_p(tt,params,random,Y_R,Y_PHI,Y_Z, & CB_ei_SD(params,v(cc),ne(cc),Te(cc),Zeff(cc),P,Y_R(cc),Y_Z(cc))) endif + !write(6,*) 'v,xi,Te,ne',v*C_C,xi,Te*params%cpp%temperature/C_E, & + ! ne*params%cpp%density + !write(6,*) 'CAL',CAL/(params%cpp%time/(C_C*C_ME)**2) + !write(6,*) 'dCAL',dCAL/(sqrt(params%cpp%time)/(C_C*C_ME)) + !write(6,*) 'CFL',CFL/(params%cpp%time/(C_C*C_ME)) + !write(6,*) 'CBL',CBL/(params%cpp%time/(C_C*C_ME)**2) + if (.not.cparams_ss%slowing_down) CFL(cc)=0._rp if (.not.cparams_ss%pitch_diffusion) CBL(cc)=0._rp if (.not.cparams_ss%energy_diffusion) THEN diff --git a/src/korc_fields.f90 b/src/korc_fields.f90 index 69e17d79..bc4cec50 100755 --- a/src/korc_fields.f90 +++ b/src/korc_fields.f90 @@ -911,6 +911,48 @@ subroutine uniform_fields_p(pchunk,F,B_X,B_Y,B_Z,E_X,E_Y,E_Z) end subroutine uniform_fields_p +subroutine uniform_fields_GC_p(pchunk,F,B_X,B_Y,B_Z,E_X,E_Y,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,PSIp) + INTEGER, INTENT(IN) :: pchunk + !! @note Subroutine that returns the value of a uniform magnetic + !! field. @endnote + !! This subroutine is used only when the simulation is ran for a + !! 'UNIFORM' plasma. As a convention, in a uniform plasma we + !! set \(\mathbf{B} = B_0 \hat{x}\). + TYPE(FIELDS), INTENT(IN) :: F + !! An instance of the KORC derived type FIELDS. + REAL(rp),DIMENSION(pchunk), INTENT(OUT) :: B_X,B_Y,B_Z + REAL(rp),DIMENSION(pchunk), INTENT(OUT) :: E_X,E_Y,E_Z + REAL(rp),DIMENSION(pchunk), INTENT(OUT) :: curlb_R,curlb_PHI,curlb_Z + REAL(rp),DIMENSION(pchunk), INTENT(OUT) :: gradB_R,gradB_PHI,gradB_Z,PSIp + !! Magnetic field components in Cartesian coordinates; + !! B(1,:) = \(B_x\), B(2,:) = \(B_y\), B(3,:) = \(B_z\) + integer(ip) :: cc + + !$OMP SIMD + do cc=1_idef,pchunk + B_X(cc) = F%Bo + B_Y(cc) = 0._rp + B_Z(cc) = 0._rp + + E_X(cc) = F%Eo + E_Y(cc) = 0._rp + E_Z(cc) = 0._rp + + curlb_R(cc) = 0._rp + curlb_PHI(cc) = 0._rp + curlb_Z(cc) = 0._rp + + gradB_R(cc) = 0._rp + gradB_PHI(cc) = 0._rp + gradB_Z(cc) = 0._rp + + PSIp(cc) = 0._rp + end do + !$OMP END SIMD + +end subroutine uniform_fields_GC_p + subroutine uniform_electric_field(F,E) !! @note Subroutine that returns the value of a uniform electric !! field. @endnote @@ -1130,7 +1172,8 @@ subroutine unitVectors(params,Xo,F,b1,b2,b3,flag,cart,hint,Bo) b2(ii,:)=b2tmp b3(ii,:)=b3tmp - !write(6,*) 'unitvectors',ii,'B',vars%B(ii,:),'psi',vars%PSI_P(ii) + !write(6,*) 'unitvectors',ii,'B',vars%B(ii,:)*params%cpp%Bo,'psi',vars%PSI_P(ii)* & + ! params%cpp%Bo*params%cpp%length**2 end do diff --git a/src/korc_finalize.f90 b/src/korc_finalize.f90 index 34588306..96e06b99 100755 --- a/src/korc_finalize.f90 +++ b/src/korc_finalize.f90 @@ -74,6 +74,11 @@ subroutine deallocate_variables(params,F,P,spp) DEALLOCATE(spp(ii)%vars%k3) DEALLOCATE(spp(ii)%vars%k4) DEALLOCATE(spp(ii)%vars%RHS) + DEALLOCATE(spp(ii)%vars%gradB) + DEALLOCATE(spp(ii)%vars%curlb) + DEALLOCATE(spp(ii)%vars%BR) + DEALLOCATE(spp(ii)%vars%BPHI) + DEALLOCATE(spp(ii)%vars%BZ) end if if (ALLOCATED(spp(ii)%BMC_ra)) DEALLOCATE(spp(ii)%BMC_ra) diff --git a/src/korc_hammersley_generator.f90 b/src/korc_hammersley_generator.f90 deleted file mode 100755 index 76dce6bc..00000000 --- a/src/korc_hammersley_generator.f90 +++ /dev/null @@ -1,826 +0,0 @@ -!> @brief Module containing subroutines for generating 1-D and 2-D Hammersley quasi-Monte Carlo sequences. -!! @details The algorithm and code for generating the 1-D Hammersley sequence was developed by John Burkardt at the Florida State University. -!! Visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html" for more information. -!! The algorithm and code for generating the 2-D Hammersley sequence was developed by L. Carbajal at the Oak Ridge National Lab. -MODULE korc_hammersley_generator - use korc_types - - IMPLICIT NONE - - PUBLIC :: hammersley,& - generate_2D_hammersley_sequence - PRIVATE :: hammersley_inverse,& - hammersley_sequence,& - r8mat_print,& - r8mat_print_some,& - timestamp,& - prime - -CONTAINS - - - !> @brief Subroutine for generating a 2-D Hammersley sequence. - !! @details This subroutine uses the algorithm for generating a 1-D Hammersley sequence. - !! Each MPI process in KORC generates a (different) subset of pairs (X,Y) of a 2-D Hammersley sequence. The total number of pairs (X,Y) - !! is NMPIS*N, where NMPIS is the number of MPI processes in the simulation and N is the number of particles followed by each MPI process. - !! Each subset of pairs (X,Y) has N elements. - !! - !! @param[in,out] X 1-D array with elements of a 2-D Hammersley sequence. - !! @param[in,out] Y 1-D array with elements of a 2-D Hammersley sequence. - !! @param[in] ID MPI rank of MPI process. - !! @param[in] NMPIS Total number of MPI processes in the simulation. - !! @param N Number of particles per MPI process. - !! @param offset An offset to indicate the subroutine what subset of the 2-D Harmmersley sequence will be generated. - !! @param ii Particle iterator. - subroutine generate_2D_hammersley_sequence(ID,NMPIS,X,Y) - REAL(rp), DIMENSION(:), INTENT(INOUT) :: X - REAL(rp), DIMENSION(:), INTENT(INOUT) :: Y - INTEGER, INTENT(IN) :: ID - INTEGER, INTENT(IN) :: NMPIS - INTEGER(4) :: N - INTEGER(4) :: offset - REAL(8), DIMENSION(2) :: R - INTEGER(4) :: ii - - N = INT(SIZE(X),4) - offset = (INT(ID+1_idef,4) - 1_4)*N - - ! write(output_unit_write,'("MPI process: ",I5," offset: ",I5)') ID, offset - - do ii=1_4,N - call hammersley(ii+INT(offset,4),2_4,N*INT(NMPIS,4),R) - - X(ii) = REAL(R(1),rp) - Y(ii) = REAL(R(2),rp) - end do - end subroutine generate_2D_hammersley_sequence - - - !> @brief For more info please visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html". - subroutine hammersley ( i, m, n, r ) - !*****************************************************************************80 - ! - !! HAMMERSLEY computes an element of a Hammersley sequence. - ! - ! Licensing: - ! - ! This code is distributed under the GNU LGPL license. - ! - ! Modified: - ! - ! 20 August 2016 - ! - ! Author: - ! - ! John Burkardt - ! - ! Reference: - ! - ! John Hammersley, - ! Monte Carlo methods for solving multivariable problems, - ! Proceedings of the New York Academy of Science, - ! Volume 86, 1960, pages 844-874. - ! - ! Parameters: - ! - ! Input, integer ( kind = 4 ) I, the index of the element of the sequence. - ! 0 <= I. - ! - ! Input, integer ( kind = 4 ) M, the spatial dimension. - ! - ! Input, integer ( kind = 4 ) N, the "base" for the first component. - ! 1 <= N. - ! - ! Output, real ( kind = 8 ) R(M), the element of the sequence with index I. - ! - implicit none - - integer ( kind = 4 ) m - - integer ( kind = 4 ) d - integer ( kind = 4 ) i - integer ( kind = 4 ) i1 - integer ( kind = 4 ) j - integer ( kind = 4 ) n - ! integer ( kind = 4 ) prime - real ( kind = 8 ) prime_inv(2:m) - real ( kind = 8 ) r(m) - integer ( kind = 4 ) t(2:m) - - t(2:m) = abs ( i ) - ! - ! Carry out the computation. - ! - do i1 = 2, m - prime_inv(i1) = 1.0D+00 / real ( prime ( i1 - 1 ), kind = 8 ) - end do - - r(1) = real ( mod ( i, n + 1 ), kind = 8 ) / real ( n, kind = 8 ) - - r(2:m) = 0.0D+00 - - do while ( any ( t(2:m) /= 0 ) ) - do j = 2, m - d = mod ( t(j), prime ( j - 1 ) ) - r(j) = r(j) + real ( d, kind = 8 ) * prime_inv(j) - prime_inv(j) = prime_inv(j) / real ( prime ( j - 1 ), kind = 8 ) - t(j) = ( t(j) / prime ( j - 1 ) ) - end do - end do - - return - end subroutine hammersley - - - !> @brief For more info please visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html". - subroutine hammersley_inverse ( r, m, n, i ) - - !*****************************************************************************80 - ! - !! HAMMERSLEY_INVERSE inverts an element of the Hammersley sequence. - ! - ! Licensing: - ! - ! This code is distributed under the GNU LGPL license. - ! - ! Modified: - ! - ! 20 August 2016 - ! - ! Author: - ! - ! John Burkardt - ! - ! Parameters: - ! - ! Input, real ( kind = 8 ) R(M), the I-th element of the Hammersley sequence. - ! 0 <= R < 1.0 - ! - ! Input, integer ( kind = 4 ) M, the spatial dimension. - ! - ! Input, integer ( kind = 4 ) N, the "base" for the first component. - ! 1 <= N. - ! - ! Output, integer ( kind = 4 ) I, the index of the element of the sequence. - ! - implicit none - - integer ( kind = 4 ) m - - integer ( kind = 4 ) d - integer ( kind = 4 ) i - integer ( kind = 4 ) n - integer ( kind = 4 ) p - real ( kind = 8 ) r(m) - real ( kind = 8 ) t - - if ( any ( r(1:m) < 0.0D+00 ) .or. any ( 1.0D+00 < r(1:m) ) ) then - write ( *, '(a)' ) '' - write ( *, '(a)' ) 'HAMMERSLEY_INVERSE - Fatal error!' - write ( *, '(a)' ) ' 0 <= R <= 1.0 is required.' - stop 1 - end if - - if ( m < 1 ) then - write ( *, '(a)' ) '' - write ( *, '(a)' ) 'HAMMERSLEY_INVERSE - Fatal error!' - write ( *, '(a)' ) ' 1 <= M is required.' - stop 1 - end if - - if ( n < 1 ) then - write ( *, '(a)' ) '' - write ( *, '(a)' ) 'HAMMERSLEY_INVERSE - Fatal error!' - write ( *, '(a)' ) ' 1 <= N is required.' - stop 1 - end if - ! - ! Invert using the second component only, because working with base - ! 2 is accurate. - ! - if ( 2 <= m ) then - - i = 0 - t = r(2) - p = 1 - - do while ( t /= 0.0D+00 ) - t = t * 2.0D+00 - d = int ( t ) - i = i + d * p - p = p * 2 - t = mod ( t, 1.0D+00 ) - end do - - else - - i = nint ( real ( n, kind = 8 ) * r(1) ) - - end if - - return - end subroutine hammersley_inverse - - - !> @brief For more info please visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html". - function prime ( n ) - !*****************************************************************************80 - ! - !! PRIME returns any of the first PRIME_MAX prime numbers. - ! - ! Discussion: - ! - ! PRIME_MAX is 1600, and the largest prime stored is 13499. - ! - ! Thanks to Bart Vandewoestyne for pointing out a typo, 18 February 2005. - ! - ! Licensing: - ! - ! This code is distributed under the GNU LGPL license. - ! - ! Modified: - ! - ! 18 February 2005 - ! - ! Author: - ! - ! John Burkardt - ! - ! Reference: - ! - ! Milton Abramowitz, Irene Stegun, - ! Handbook of Mathematical Functions, - ! US Department of Commerce, 1964, pages 870-873. - ! - ! Daniel Zwillinger, - ! CRC Standard Mathematical Tables and Formulae, - ! 30th Edition, - ! CRC Press, 1996, pages 95-98. - ! - ! Parameters: - ! - ! Input, integer ( kind = 4 ) N, the index of the desired prime number. - ! In general, is should be true that 0 <= N <= PRIME_MAX. - ! N = -1 returns PRIME_MAX, the index of the largest prime available. - ! N = 0 is legal, returning PRIME = 1. - ! - ! Output, integer ( kind = 4 ) PRIME, the N-th prime. If N is out of range, - ! PRIME is returned as -1. - ! - implicit none - - integer ( kind = 4 ), parameter :: prime_max = 1600 - - integer ( kind = 4 ), save :: icall = 0 - integer ( kind = 4 ) n - integer ( kind = 4 ), save, dimension ( prime_max ) :: npvec - integer ( kind = 4 ) prime - - if ( icall == 0 ) then - - icall = 1 - - npvec(1:100) = (/ & - 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, & - 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, & - 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, & - 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, & - 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, & - 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, & - 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, & - 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, & - 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, & - 467, 479, 487, 491, 499, 503, 509, 521, 523, 541 /) - - npvec(101:200) = (/ & - 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, & - 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, & - 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, & - 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, & - 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, & - 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, & - 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, & - 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, & - 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, & - 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223 /) - - npvec(201:300) = (/ & - 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, & - 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, & - 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, & - 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, & - 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, & - 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, & - 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, & - 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, & - 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, & - 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987 /) - - npvec(301:400) = (/ & - 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, & - 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, & - 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, & - 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, & - 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, & - 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, & - 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, & - 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, & - 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, & - 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741 /) - - npvec(401:500) = (/ & - 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, & - 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, & - 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, & - 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, & - 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, & - 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, & - 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, & - 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, & - 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, & - 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571 /) - - npvec(501:600) = (/ & - 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, & - 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, & - 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, & - 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, & - 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, & - 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, & - 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, & - 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, & - 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, & - 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409 /) - - npvec(601:700) = (/ & - 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, & - 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, & - 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, & - 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, & - 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, & - 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, & - 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, & - 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, & - 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, & - 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279 /) - - npvec(701:800) = (/ & - 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, & - 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, & - 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, & - 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, & - 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, & - 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, & - 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, & - 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, & - 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, & - 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133 /) - - npvec(801:900) = (/ & - 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, & - 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, & - 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, & - 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, & - 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, & - 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, & - 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, & - 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, & - 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, & - 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997 /) - - npvec(901:1000) = (/ & - 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, & - 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, & - 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, & - 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, & - 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, & - 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, & - 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, & - 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, & - 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, & - 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919 /) - - npvec(1001:1100) = (/ & - 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, & - 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, & - 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, & - 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, & - 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, & - 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, & - 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, & - 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, & - 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, & - 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831 /) - - npvec(1101:1200) = (/ & - 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, & - 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, & - 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, & - 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, & - 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, & - 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, & - 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, & - 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, & - 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, & - 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733 /) - - npvec(1201:1300) = (/ & - 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, & - 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, & - 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973,10007, & - 10009,10037,10039,10061,10067,10069,10079,10091,10093,10099, & - 10103,10111,10133,10139,10141,10151,10159,10163,10169,10177, & - 10181,10193,10211,10223,10243,10247,10253,10259,10267,10271, & - 10273,10289,10301,10303,10313,10321,10331,10333,10337,10343, & - 10357,10369,10391,10399,10427,10429,10433,10453,10457,10459, & - 10463,10477,10487,10499,10501,10513,10529,10531,10559,10567, & - 10589,10597,10601,10607,10613,10627,10631,10639,10651,10657 /) - - npvec(1301:1400) = (/ & - 10663,10667,10687,10691,10709,10711,10723,10729,10733,10739, & - 10753,10771,10781,10789,10799,10831,10837,10847,10853,10859, & - 10861,10867,10883,10889,10891,10903,10909,10937,10939,10949, & - 10957,10973,10979,10987,10993,11003,11027,11047,11057,11059, & - 11069,11071,11083,11087,11093,11113,11117,11119,11131,11149, & - 11159,11161,11171,11173,11177,11197,11213,11239,11243,11251, & - 11257,11261,11273,11279,11287,11299,11311,11317,11321,11329, & - 11351,11353,11369,11383,11393,11399,11411,11423,11437,11443, & - 11447,11467,11471,11483,11489,11491,11497,11503,11519,11527, & - 11549,11551,11579,11587,11593,11597,11617,11621,11633,11657 /) - - npvec(1401:1500) = (/ & - 11677,11681,11689,11699,11701,11717,11719,11731,11743,11777, & - 11779,11783,11789,11801,11807,11813,11821,11827,11831,11833, & - 11839,11863,11867,11887,11897,11903,11909,11923,11927,11933, & - 11939,11941,11953,11959,11969,11971,11981,11987,12007,12011, & - 12037,12041,12043,12049,12071,12073,12097,12101,12107,12109, & - 12113,12119,12143,12149,12157,12161,12163,12197,12203,12211, & - 12227,12239,12241,12251,12253,12263,12269,12277,12281,12289, & - 12301,12323,12329,12343,12347,12373,12377,12379,12391,12401, & - 12409,12413,12421,12433,12437,12451,12457,12473,12479,12487, & - 12491,12497,12503,12511,12517,12527,12539,12541,12547,12553 /) - - npvec(1501:1600) = (/ & - 12569,12577,12583,12589,12601,12611,12613,12619,12637,12641, & - 12647,12653,12659,12671,12689,12697,12703,12713,12721,12739, & - 12743,12757,12763,12781,12791,12799,12809,12821,12823,12829, & - 12841,12853,12889,12893,12899,12907,12911,12917,12919,12923, & - 12941,12953,12959,12967,12973,12979,12983,13001,13003,13007, & - 13009,13033,13037,13043,13049,13063,13093,13099,13103,13109, & - 13121,13127,13147,13151,13159,13163,13171,13177,13183,13187, & - 13217,13219,13229,13241,13249,13259,13267,13291,13297,13309, & - 13313,13327,13331,13337,13339,13367,13381,13397,13399,13411, & - 13417,13421,13441,13451,13457,13463,13469,13477,13487,13499 /) - - end if - - if ( n == -1 ) then - prime = prime_max - else if ( n == 0 ) then - prime = 1 - else if ( n <= prime_max ) then - prime = npvec(n) - else - prime = -1 - write ( *, '(a)' ) ' ' - write ( *, '(a)' ) 'PRIME - Fatal error!' - write ( *, '(a,i8)' ) ' Illegal prime index N = ', n - write ( *, '(a,i8)' ) ' N should be between 1 and PRIME_MAX =', prime_max - stop 1 - end if - - return - end function prime - - - !> @brief For more info please visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html". - subroutine hammersley_sequence ( i1, i2, m, n, r ) - !*****************************************************************************80 - ! - !! HAMMERSLEY_SEQUENCE computes elements I1 through I2 of a Hammersley sequence. - ! - ! Licensing: - ! - ! This code is distributed under the GNU LGPL license. - ! - ! Modified: - ! - ! 20 August 2016 - ! - ! Author: - ! - ! John Burkardt - ! - ! Reference: - ! - ! John Hammersley, - ! Monte Carlo methods for solving multivariable problems, - ! Proceedings of the New York Academy of Science, - ! Volume 86, 1960, pages 844-874. - ! - ! Parameters: - ! - ! Input, integer ( kind = 4 ) I1, I2, the indices of the first and last - ! elements of the sequence. 0 <= I1, I2. - ! - ! Input, integer ( kind = 4 ) M, the spatial dimension. - ! - ! Input, integer ( kind = 4 ) N, the "base" for the first component. - ! 1 <= N. - ! - ! Output, real ( kind = 8 ) R(M,abs(I1-I2)+1), the elements of the sequence - ! with indices I1 through I2. - ! - implicit none - - integer ( kind = 4 ) m - - integer ( kind = 4 ) d - integer ( kind = 4 ) i - integer ( kind = 4 ) i1 - integer ( kind = 4 ) i2 - integer ( kind = 4 ) i3 - integer ( kind = 4 ) j - integer ( kind = 4 ) k - integer ( kind = 4 ) l - integer ( kind = 4 ) n - ! integer ( kind = 4 ) prime - real ( kind = 8 ) prime_inv(2:m) - real ( kind = 8 ) r(m,abs(i1-i2)+1) - integer ( kind = 4 ) t(2:m) - - if ( i1 <= i2 ) then - i3 = +1 - else - i3 = -1 - end if - - l = abs ( i2 - i1 ) + 1 - r(1:m,1:l) = 0.0D+00 - k = 0 - - do i = i1, i2, i3 - - t(2:m) = i - ! - ! Carry out the computation. - ! - do j = 2, m - prime_inv(j) = 1.0D+00 / real ( prime ( j - 1 ), kind = 8 ) - end do - - k = k + 1; - - r(1,k) = real ( mod ( i, n + 1 ), kind = 8 ) / real ( n, kind = 8 ) - - do while ( any ( t(2:m) /= 0 ) ) - do j = 2, m - d = mod ( t(j), prime ( j - 1 ) ) - r(j,k) = r(j,k) + real ( d, kind = 8 ) * prime_inv(j) - prime_inv(j) = prime_inv(j) / real ( prime ( j - 1 ), kind = 8 ) - t(j) = ( t(j) / prime ( j - 1 ) ) - end do - end do - - end do - - return - end subroutine hammersley_sequence - - - !> @brief For more info please visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html". - subroutine r8mat_print ( m, n, a, title ) - !*****************************************************************************80 - ! - !! R8MAT_PRINT prints an R8MAT. - ! - ! Discussion: - ! - ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. - ! - ! Licensing: - ! - ! This code is distributed under the GNU LGPL license. - ! - ! Modified: - ! - ! 12 September 2004 - ! - ! Author: - ! - ! John Burkardt - ! - ! Parameters: - ! - ! Input, integer ( kind = 4 ) M, the number of rows in A. - ! - ! Input, integer ( kind = 4 ) N, the number of columns in A. - ! - ! Input, real ( kind = 8 ) A(M,N), the matrix. - ! - ! Input, character ( len = * ) TITLE, a title. - ! - implicit none - - integer ( kind = 4 ) m - integer ( kind = 4 ) n - - real ( kind = 8 ) a(m,n) - character ( len = * ) title - - call r8mat_print_some ( m, n, a, 1, 1, m, n, title ) - - return - end subroutine r8mat_print - - - !> @brief For more info please visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html". - subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) - !*****************************************************************************80 - ! - !! R8MAT_PRINT_SOME prints some of an R8MAT. - ! - ! Discussion: - ! - ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. - ! - ! Licensing: - ! - ! This code is distributed under the GNU LGPL license. - ! - ! Modified: - ! - ! 10 September 2009 - ! - ! Author: - ! - ! John Burkardt - ! - ! Parameters: - ! - ! Input, integer ( kind = 4 ) M, N, the number of rows and columns. - ! - ! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed. - ! - ! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print. - ! - ! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print. - ! - ! Input, character ( len = * ) TITLE, a title. - ! - implicit none - - integer ( kind = 4 ), parameter :: incx = 5 - integer ( kind = 4 ) m - integer ( kind = 4 ) n - - real ( kind = 8 ) a(m,n) - character ( len = 14 ) ctemp(incx) - integer ( kind = 4 ) i - integer ( kind = 4 ) i2hi - integer ( kind = 4 ) i2lo - integer ( kind = 4 ) ihi - integer ( kind = 4 ) ilo - integer ( kind = 4 ) inc - integer ( kind = 4 ) j - integer ( kind = 4 ) j2 - integer ( kind = 4 ) j2hi - integer ( kind = 4 ) j2lo - integer ( kind = 4 ) jhi - integer ( kind = 4 ) jlo - character ( len = * ) title - - write ( *, '(a)' ) ' ' - write ( *, '(a)' ) trim ( title ) - - if ( m <= 0 .or. n <= 0 ) then - write ( *, '(a)' ) ' ' - write ( *, '(a)' ) ' (None)' - return - end if - - do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx - - j2hi = j2lo + incx - 1 - j2hi = min ( j2hi, n ) - j2hi = min ( j2hi, jhi ) - - inc = j2hi + 1 - j2lo - - write ( *, '(a)' ) ' ' - - do j = j2lo, j2hi - j2 = j + 1 - j2lo - write ( ctemp(j2), '(i8,6x)' ) j - end do - - write ( *, '('' Col '',5a14)' ) ctemp(1:inc) - write ( *, '(a)' ) ' Row' - write ( *, '(a)' ) ' ' - - i2lo = max ( ilo, 1 ) - i2hi = min ( ihi, m ) - - do i = i2lo, i2hi - - do j2 = 1, inc - - j = j2lo - 1 + j2 - - if ( a(i,j) == real ( int ( a(i,j) ), kind = 8 ) ) then - write ( ctemp(j2), '(f8.0,6x)' ) a(i,j) - else - write ( ctemp(j2), '(g14.6)' ) a(i,j) - end if - - end do - - write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc ) - - end do - - end do - - return - end subroutine r8mat_print_some - - - !> @brief For more info please visit "https://people.sc.fsu.edu/~jburkardt/f_src/hammersley/hammersley.html". - subroutine timestamp ( ) - !*****************************************************************************80 - ! - !! TIMESTAMP prints the current YMDHMS date as a time stamp. - ! - ! Example: - ! - ! 31 May 2001 9:45:54.872 AM - ! - ! Licensing: - ! - ! This code is distributed under the GNU LGPL license. - ! - ! Modified: - ! - ! 18 May 2013 - ! - ! Author: - ! - ! John Burkardt - ! - ! Parameters: - ! - ! None - ! - implicit none - - character ( len = 8 ) ampm - integer ( kind = 4 ) d - integer ( kind = 4 ) h - integer ( kind = 4 ) m - integer ( kind = 4 ) mm - character ( len = 9 ), parameter, dimension(12) :: month = (/ & - 'January ', 'February ', 'March ', 'April ', & - 'May ', 'June ', 'July ', 'August ', & - 'September', 'October ', 'November ', 'December ' /) - integer ( kind = 4 ) n - integer ( kind = 4 ) s - integer ( kind = 4 ) values(8) - integer ( kind = 4 ) y - - call date_and_time ( values = values ) - - y = values(1) - m = values(2) - d = values(3) - h = values(5) - n = values(6) - s = values(7) - mm = values(8) - - if ( h < 12 ) then - ampm = 'AM' - else if ( h == 12 ) then - if ( n == 0 .and. s == 0 ) then - ampm = 'Noon' - else - ampm = 'PM' - end if - else - h = h - 12 - if ( h < 12 ) then - ampm = 'PM' - else if ( h == 12 ) then - if ( n == 0 .and. s == 0 ) then - ampm = 'Midnight' - else - ampm = 'AM' - end if - end if - end if - - write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & - d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) - - return - end subroutine timestamp -END MODULE korc_hammersley_generator diff --git a/src/korc_initialize.f90 b/src/korc_initialize.f90 index e0624c8a..1de0ebb8 100755 --- a/src/korc_initialize.f90 +++ b/src/korc_initialize.f90 @@ -276,8 +276,10 @@ subroutine define_time_step(params,F) if (params%output_cadence.EQ.0_ip) params%output_cadence = 1_ip - params%num_snapshots = params%t_steps/params%output_cadence - + !params%num_snapshots = params%t_steps/params%output_cadence + params%num_snapshots = ceiling(params%simulation_time/ & + params%snapshot_frequency) + if (params%t_steps.gt.params%output_cadence) then #ifdef __NVCOMPILER params%dt=params%snapshot_frequency/real(params%output_cadence) @@ -305,7 +307,9 @@ subroutine define_time_step(params,F) if (params%output_cadence.EQ.0_ip) params%output_cadence = 1_ip - params%num_snapshots = params%t_steps/params%output_cadence + !params%num_snapshots = params%t_steps/params%output_cadence + params%num_snapshots = ceiling(params%simulation_time/ & + params%snapshot_frequency) if (params%t_steps.gt.params%output_cadence) then if (params%snapshot_frequency.lt.params%dt) then @@ -317,10 +321,13 @@ subroutine define_time_step(params,F) #else params%dt=params%snapshot_frequency/float(params%output_cadence) #endif + + params%t_steps = CEILING(params%simulation_time/params%dt,ip) + endif - params%restart_output_cadence = CEILING(params%restart_overwrite_frequency/ & - params%dt,ip) + params%restart_output_cadence = & + CEILING(params%restart_overwrite_frequency/params%dt,ip) params%t_skip=min(params%t_steps,params%output_cadence) params%t_skip=max(1_ip,params%t_skip) @@ -596,25 +603,18 @@ subroutine initialize_particles(params,random,F,P,spp) ALLOCATE( spp(ii)%vars%k4(spp(ii)%ppp,4) ) ALLOCATE( spp(ii)%vars%k5(spp(ii)%ppp,4) ) ALLOCATE( spp(ii)%vars%k6(spp(ii)%ppp,4) ) - if (params%orbit_model(3:5)=='pre'.or. & - TRIM(params%field_model)=='M3D_C1'.or. & - TRIM(params%field_model)=='NIMROD') then - ALLOCATE( spp(ii)%vars%gradB(spp(ii)%ppp,3) ) - ALLOCATE( spp(ii)%vars%curlb(spp(ii)%ppp,3) ) - - spp(ii)%vars%gradB = 0.0_rp - spp(ii)%vars%curlb = 0.0_rp - else if (params%orbit_model(3:6)=='grad') then - ALLOCATE( spp(ii)%vars%BR(spp(ii)%ppp,3) ) - ALLOCATE( spp(ii)%vars%BPHI(spp(ii)%ppp,3) ) - ALLOCATE( spp(ii)%vars%BZ(spp(ii)%ppp,3) ) - - spp(ii)%vars%BR = 0.0_rp - spp(ii)%vars%BPHI = 0.0_rp - spp(ii)%vars%BZ = 0.0_rp - end if + ALLOCATE( spp(ii)%vars%gradB(spp(ii)%ppp,3) ) + ALLOCATE( spp(ii)%vars%curlb(spp(ii)%ppp,3) ) + ALLOCATE( spp(ii)%vars%BR(spp(ii)%ppp,3) ) + ALLOCATE( spp(ii)%vars%BPHI(spp(ii)%ppp,3) ) + ALLOCATE( spp(ii)%vars%BZ(spp(ii)%ppp,3) ) ALLOCATE( spp(ii)%vars%RHS(spp(ii)%ppp,5) ) + spp(ii)%vars%gradB = 0.0_rp + spp(ii)%vars%curlb = 0.0_rp + spp(ii)%vars%BR = 0.0_rp + spp(ii)%vars%BPHI = 0.0_rp + spp(ii)%vars%BZ = 0.0_rp spp(ii)%vars%Y0 = 0.0_rp spp(ii)%vars%Y1 = 0.0_rp spp(ii)%vars%V0 = 0.0_rp diff --git a/src/korc_input.f90 b/src/korc_input.f90 index f01646e2..ba871d1b 100755 --- a/src/korc_input.f90 +++ b/src/korc_input.f90 @@ -214,6 +214,8 @@ module korc_input REAL(rp) :: MARS_phase=0.0 REAL(rp) :: AORSA_AMP_Scale=1.0 REAL(rp) :: AORSA_freq=0.0 + REAL(rp) :: psir=0.0 + REAL(rp) :: width=0.0 REAL(rp) :: AORSA_nmode=0.0 REAL(rp) :: AORSA_mmode=0.0 LOGICAL :: Analytic_D3D_IWL=.FALSE. @@ -465,7 +467,7 @@ subroutine read_namelist(params,infile,echo_in,outdir) axisymmetric_fields, Eo,E_dyn,E_pulse,E_width,res_double, & dim_1D,dt_E_SC,Ip_exp,PSIp_lim,Dim2x1t,t0_2x1t,E_2x1t,ReInterp_2x1t, & ind0_2x1t,PSIp_0,B1field,psip_conv,MARS_AMP_Scale,Analytic_D3D_IWL, & - ntiles,circumradius,AORSA_AMP_Scale,AORSA_freq,AORSA_nmode,AORSA_mmode,E1field, & + ntiles,circumradius,AORSA_AMP_Scale,AORSA_freq,AORSA_nmode,AORSA_mmode,width,psir,E1field, & useLCFS,useDiMES,DiMESloc,DiMESdims,MARS_phase NAMELIST /plasmaProfiles/ radius_profile,ne_profile,neo,n_ne,a_ne, & Te_profile,Teo,n_Te,a_Te,n_REr0,n_tauion,n_lamfront,n_lamback, & diff --git a/src/korc_interp.f90 b/src/korc_interp.f90 index 7dd08f9f..f43b1de2 100755 --- a/src/korc_interp.f90 +++ b/src/korc_interp.f90 @@ -2575,8 +2575,9 @@ subroutine interp_FOfields_aorsa(prtcls, F, params) TYPE(FIELDS), INTENT(IN) :: F REAL(rp) :: Y_R,Y_PHI,Y_Z REAL(rp) :: B0_R,B0_PHI,B0_Z - REAL(rp) :: B0_X,B0_Y + REAL(rp) :: B0_X,B0_Y,theta REAL(rp) :: B1_X,B1_Y,B1_Z + REAL(rp) :: E1_X,E1_Y,E1_Z REAL(rp) :: B1Re_X,B1Re_Y,B1Re_Z REAL(rp) :: B1Im_X,B1Im_Y,B1Im_Z REAL(rp) :: E1Re_X,E1Re_Y,E1Re_Z @@ -2586,8 +2587,8 @@ subroutine interp_FOfields_aorsa(prtcls, F, params) ! INTEGER(ip) :: ezerr INTEGER :: pp,ss !! Particle chunk iterator. - REAL(rp) :: psip_conv - REAL(rp) :: amp,nmode + REAL(rp) :: psip_conv,psirr,widthh + REAL(rp) :: amp,nmode,omega,BXavg,BYavg,BZavg,EXavg,EYavg,EZavg,mmode if (size(prtcls%Y,1).eq.1) then ss = size(prtcls%Y,1) @@ -2601,6 +2602,7 @@ subroutine interp_FOfields_aorsa(prtcls, F, params) psip_conv=F%psip_conv amp=F%AMP + mmode=F%AORSA_mmode nmode=F%AORSA_nmode do pp = 1,ss @@ -2608,6 +2610,7 @@ subroutine interp_FOfields_aorsa(prtcls, F, params) Y_R=prtcls%Y(pp,1) Y_PHI=prtcls%Y(pp,2) Y_Z=prtcls%Y(pp,3) + theta=atan2(Y_Z,(Y_R-1.7)) call EZspline_interp2_FOaorsa(bfield_2d%A,& b1Refield_2dx%X,b1Refield_2dx%Y,b1Refield_2dx%Z, & @@ -2619,6 +2622,8 @@ subroutine interp_FOfields_aorsa(prtcls, F, params) call EZspline_error(ezerr) prtcls%PSI_P=A(1) + widthh=F%width + psirr=F%psir B0_R = psip_conv*A(3)/Y_R B0_PHI = -F%Bo*F%Ro/Y_R @@ -2627,13 +2632,31 @@ subroutine interp_FOfields_aorsa(prtcls, F, params) cnP=cos(nmode*Y_PHI) snP=sin(nmode*Y_PHI) + !cnP=cos(mmode*theta-nmode*Y_PHI) + !snP=sin(mmode*theta-nmode*Y_PHI) + cP=cos(Y_PHI) sP=sin(Y_PHI) + !BXavg=6.95E-07/params%cpp%Bo + !BYavg=1.184E-06/params%cpp%Bo + !BZavg=1.118E-06/params%cpp%Bo + + !B1_X = amp*(BXavg*((1/cosh((prtcls%PSI_P-psirr)/widthh))**2)*cnp) + !B1_Y = amp*(BYavg*((1/cosh((prtcls%PSI_P-psirr)/widthh))**2)*snp) + !B1_Z =-amp*(BZavg*((1/cosh((prtcls%PSI_P-psirr)/widthh))**2)*cnp) + B1_X = amp*(B1Re_X*cnP-B1Im_X*snP) B1_Y = amp*(B1Re_Y*cnP-B1Im_Y*snP) B1_Z = amp*(B1Re_Z*cnP-B1Im_Z*snP) + !EXavg=76.2/params%cpp%Eo + !EYavg=57.7/params%cpp%Eo + !EZavg=9.25/params%cpp%Eo + !prtcls%E(pp,1)= amp*((1/cosh((prtcls%PSI_P-psirr)/widthh))**2)*EXavg*snp + !prtcls%E(pp,2)= amp*((1/cosh((prtcls%PSI_P-psirr)/widthh))**2)*EYavg*cnp + !prtcls%E(pp,3)=-amp*((1/cosh((prtcls%PSI_P-psirr)/widthh))**2)*EZavg*snp + prtcls%E(pp,1) = amp*(E1Re_X*cnP-E1Im_X*snP) prtcls%E(pp,2) = amp*(E1Re_Y*cnP-E1Im_Y*snP) prtcls%E(pp,3) = amp*(E1Re_Z*cnP-E1Im_Z*snP) @@ -2648,15 +2671,15 @@ subroutine interp_FOfields_aorsa(prtcls, F, params) !write(6,*) '(R,PHI,Z)',Y_R*params%cpp%length,Y_PHI, & ! Y_Z*params%cpp%length !write(6,*) 'amp',amp,'cP,sP',cP,sP,'cnP,snP',cnP,snP - !write(6,*) 'psi',PSIp*params%cpp%Bo*params%cpp%length**2 - !write(6,*) 'dpsidR',A(:,2)*params%cpp%Bo*params%cpp%length - !write(6,*) 'dpsidZ',A(:,3)*params%cpp%Bo*params%cpp%length + !write(6,*) 'psi',A(1)*params%cpp%Bo*params%cpp%length**2 + !write(6,*) 'dpsidR',A(2)*params%cpp%Bo*params%cpp%length + !write(6,*) 'dpsidZ',A(3)*params%cpp%Bo*params%cpp%length !write(6,*) 'B0',B0_R*params%cpp%Bo,B0_PHI*params%cpp%Bo,B0_Z*params%cpp%Bo !write(6,*) 'AMP',amp !write(6,*) 'B1Re',B1Re_X*params%cpp%Bo,B1Re_Y*params%cpp%Bo,B1Re_Z*params%cpp%Bo !write(6,*) 'B1Im',B1Im_X*params%cpp%Bo,B1Im_Y*params%cpp%Bo,B1Im_Z*params%cpp%Bo !write(6,*) 'B1',B1_X*params%cpp%Bo,B1_Y*params%cpp%Bo,B1_Z*params%cpp%Bo - !write(6,*) 'B',B_X*params%cpp%Bo,B_Y*params%cpp%Bo,B_Z*params%cpp%Bo + !write(6,*) 'B',(B0_X+B1_X)*params%cpp%Bo,(B0_Y+B1_Y)*params%cpp%Bo,(B0_Z+B1_Z)*params%cpp%Bo !write(6,*) 'unitvectors',pp,prtcls%B(pp,:) @@ -2702,6 +2725,13 @@ subroutine interp_FOfields_aorsa_p_ACC(time,bfield_2d_local,b1Refield_2dx_local, E1Re_X,E1Re_Y,E1Re_Z,E1Im_X,E1Im_Y,E1Im_Z,ezerr_local) call EZspline_error(ezerr_local) + !EXavg=76.2/params%cpp%Eo + !EYavg=57.7/params%cpp%Eo + !EZavg=9.25/params%cpp%Eo + !E_X(cc) = -amp*(EXavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*snp(cc) + !E_Y(cc) = amp*(EYavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*cnp(cc) + !E_Z(cc) = amp*(EZavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*snp(cc) + PSIp=A(1) @@ -2748,6 +2778,7 @@ subroutine interp_FOfields_aorsa_p(time,params,pchunk,F,Y_R,Y_PHI,Y_Z, & REAL(rp) :: B1Im_X,B1Im_Y,B1Im_Z REAL(rp) :: E1Re_X,E1Re_Y,E1Re_Z REAL(rp) :: E1Im_X,E1Im_Y,E1Im_Z + REAL(rp) :: radius,theta REAL(rp),DIMENSION(pchunk),INTENT(OUT) :: PSIp REAL(rp) :: cP,sP,cnP,snP REAL(rp), DIMENSION(3) :: A @@ -2755,13 +2786,14 @@ subroutine interp_FOfields_aorsa_p(time,params,pchunk,F,Y_R,Y_PHI,Y_Z, & INTEGER :: cc !! Particle chunk iterator. INTEGER(is),DIMENSION(pchunk),INTENT(INOUT) :: flag_cache - REAL(rp) :: psip_conv - REAL(rp) :: amp,nmode,omega + REAL(rp) :: psip_conv,psirr,widthh + REAL(rp) :: amp,nmode,omega,BXavg,BYavg,BZavg,EXavg,EYavg,EZavg,mmode psip_conv=F%psip_conv amp=F%AMP nmode=F%AORSA_nmode omega=2*C_PI*F%AORSA_freq + mmode=F%AORSA_mmode call check_if_in_fields_domain_p(pchunk,F,Y_R,Y_PHI,Y_Z,flag_cache) @@ -2777,6 +2809,10 @@ subroutine interp_FOfields_aorsa_p(time,params,pchunk,F,Y_R,Y_PHI,Y_Z, & PSIp(cc)=A(1) + radius=sqrt((Y_R(cc)-1.7)**2-Y_Z(cc)**2) + theta=atan2(Y_Z(cc),(Y_R(cc)-F%AB%Ro)) + psirr=F%psir + widthh=F%width B0_R = psip_conv*A(3)/Y_R(cc) B0_PHI = -F%Bo*F%Ro/Y_R(cc) @@ -2789,6 +2825,9 @@ subroutine interp_FOfields_aorsa_p(time,params,pchunk,F,Y_R,Y_PHI,Y_Z, & B0_Y = B0_R*sP + B0_PHI*cP + !cnP=cos(mmode*theta-omega*time-nmode*Y_PHI(cc)) + !snP=sin(mmode*theta-omega*time-nmode*Y_PHI(cc)) + cnP=cos(omega*time+nmode*Y_PHI(cc)) snP=sin(omega*time+nmode*Y_PHI(cc)) @@ -2796,10 +2835,25 @@ subroutine interp_FOfields_aorsa_p(time,params,pchunk,F,Y_R,Y_PHI,Y_Z, & B1_Y = amp*(B1Re_Y*cnP-B1Im_Y*snP) B1_Z = amp*(B1Re_Z*cnP-B1Im_Z*snP) + !BXavg=6.95E-07/params%cpp%Bo + !BYavg=1.184E-06/params%cpp%Bo + !BZavg=1.118E-06/params%cpp%Bo + + !B1_X = amp*(BXavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*cnp + !B1_Y = amp*(BYavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*snp + !B1_Z = -amp*(BZavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*cnp + B_X(cc) = B0_X+B1_X B_Y(cc) = B0_Y+B1_Y B_Z(cc) = B0_Z+B1_Z + !EXavg=76.2/params%cpp%Eo + !EYavg=57.7/params%cpp%Eo + !EZavg=9.25/params%cpp%Eo + !E_X(cc) = -amp*(EXavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*snp + !E_Y(cc) = amp*(EYavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*cnp + !E_Z(cc) = amp*(EZavg*((1/cosh((PSIp(cc)-psirr)/widthh))**2))*snp + E_X(cc) = amp*(E1Re_X*cnP-E1Im_X*snP) E_Y(cc) = amp*(E1Re_Y*cnP-E1Im_Y*snP) E_Z(cc) = amp*(E1Re_Z*cnP-E1Im_Z*snP) @@ -2808,19 +2862,18 @@ subroutine interp_FOfields_aorsa_p(time,params,pchunk,F,Y_R,Y_PHI,Y_Z, & end do !$OMP END SIMD -#if DBG_CHECK !write(6,*) '(R,PHI,Z,time)',Y_R*params%cpp%length,Y_PHI, & ! Y_Z*params%cpp%length,time !write(6,*) 'psi',PSIp*params%cpp%Bo*params%cpp%length**2 - !write(6,*) 'dpsidR',A(:,2)*params%cpp%Bo*params%cpp%length - !write(6,*) 'dpsidZ',A(:,3)*params%cpp%Bo*params%cpp%length + !write(6,*) 'dpsidR',A(2)*params%cpp%Bo*params%cpp%length + !write(6,*) 'dpsidZ',A(3)*params%cpp%Bo*params%cpp%length !write(6,*) 'B0',B0_R*params%cpp%Bo,B0_PHI*params%cpp%Bo,B0_Z*params%cpp%Bo !write(6,*) 'AMP',amp !write(6,*) 'B1Re',B1Re_X*params%cpp%Bo,B1Re_Y*params%cpp%Bo,B1Re_Z*params%cpp%Bo !write(6,*) 'B1Im',B1Im_X*params%cpp%Bo,B1Im_Y*params%cpp%Bo,B1Im_Z*params%cpp%Bo !write(6,*) 'B1',B1_X*params%cpp%Bo,B1_Y*params%cpp%Bo,B1_Z*params%cpp%Bo !write(6,*) 'B',B_X*params%cpp%Bo,B_Y*params%cpp%Bo,B_Z*params%cpp%Bo -#endif + end subroutine interp_FOfields_aorsa_p subroutine interp_FOcollision_p(pchunk,Y_R,Y_PHI,Y_Z,ne,Te,Zeff,flag_cache) @@ -3601,6 +3654,7 @@ subroutine interp_fields(params,prtcls,F) if (((ALLOCATED(F%PSIp).and.F%Bflux).or. & F%ReInterp_2x1t).and.(.not.((TRIM(params%field_model).eq.'M3D_C1').or. & (params%field_model(10:13).eq.'MARS').or. & + (params%field_model(10:14).eq.'AORSA').or. & (TRIM(params%field_model).eq.'NIMROD')))) then ! write(output_unit_write,'("3 size of PSI_P: ",I16)') size(prtcls%PSI_P) diff --git a/src/korc_ppusher.f90 b/src/korc_ppusher.f90 index 3eb33840..1903a09d 100755 --- a/src/korc_ppusher.f90 +++ b/src/korc_ppusher.f90 @@ -2486,17 +2486,18 @@ subroutine adv_FOinterp_top(params,random,F,P,spp) call cart_to_cyl_p(pchunk,X_X,X_Y,X_Z,Y_R,Y_PHI,Y_Z) - if (F%axisymmetric_fields.and. & - (params%orbit_model(3:3).eq.'B')) then + if (F%axisymmetric_fields.and.F%Bfield) then call interp_FOfields_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z, & E_X,E_Y,E_Z,PSIp,flagCon) - else if ((.not.F%axisymmetric_fields).and. & - (params%orbit_model(3:3).eq.'B')) then + else if ((.not.F%axisymmetric_fields).and.F%Bfield) then call interp_FO3Dfields_p(pchunk,F,Y_R,Y_PHI,Y_Z, & B_X,B_Y,B_Z,E_X,E_Y,E_Z,flagCon) - else if (params%orbit_model(3:5).eq.'psi') then + else if (F%axisymmetric_fields.and.F%Bflux) then call interp_FOfields1_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_X,B_Y,B_Z, & E_X,E_Y,E_Z,PSIp,flagCon) + else + write(6,*) 'No fields interpolated!!!' + call KORC_ABORT(25) end if @@ -4402,7 +4403,7 @@ subroutine GC_init(params,F,spp) Bmag1 = SQRT( DOT_PRODUCT(spp(ii)%vars%B(pp,:), & spp(ii)%vars%B(pp,:))) - pmag=sqrt(spp(ii)%vars%g(pp)**2-1) + pmag=spp(ii)%m*sqrt(spp(ii)%vars%g(pp)**2-1) spp(ii)%vars%V(pp,1)=pmag*cos(deg2rad(spp(ii)%vars%eta(pp))) @@ -4489,7 +4490,10 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !! and simulation variables of the different species in the simulation. REAL(rp), DIMENSION(params%pchunk) :: Bmag REAL(rp),DIMENSION(params%pchunk) :: Y_R,Y_PHI,Y_Z - REAL(rp),DIMENSION(params%pchunk) :: B_R,B_PHI,B_Z,E_PHI + REAL(rp),DIMENSION(params%pchunk) :: B_R,B_PHI,B_Z + REAL(rp),DIMENSION(params%pchunk) :: E_R,E_PHI,E_Z + REAL(rp),DIMENSION(params%pchunk) :: gradB_R,gradB_PHI,gradB_Z + REAL(rp),DIMENSION(params%pchunk) :: curlb_R,curlb_PHI,curlb_Z REAL(rp),DIMENSION(params%pchunk) :: PSIp,ne,Te REAL(rp),DIMENSION(params%pchunk) :: V_PLL,V_MU REAL(rp) :: B0,EF0,R0,q0,lam,ar,m_cache,q_cache,ne0,Te0,Zeff0 @@ -4529,8 +4533,9 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !$OMP& FIRSTPRIVATE(E0,q_cache,m_cache,pchunk) & !$OMP& shared(F,P,params,ii,spp,random) & !$OMP& PRIVATE(pp,tt,ttt,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & - !$OMP& flagCon,flagCol,B_R,B_PHI,B_Z,E_PHI,PSIp,ne, & - !$OMP& Vden,Vdenave) & + !$OMP& flagCon,flagCol,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,PSIp,ne, & + !$OMP& Vden,Vdenave,gradB_R,gradB_PHI,gradB_Z, & + !$OMP& curlb_R,curlb_PHI,curlb_Z) & !$OMP& REDUCTION(+:VdenOMP) do pp=1_idef,spp(ii)%ppp,pchunk @@ -4561,7 +4566,7 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) call advance_GCeqn_vars(spp(ii)%vars,pp, & tt+params%t_skip*(ttt-1),params, & Y_R,Y_PHI, Y_Z,V_PLL,V_MU,flagCon,flagCol,q_cache,m_cache, & - B_R,B_PHI,B_Z,F,P,PSIp,E_PHI) + B_R,B_PHI,B_Z,F,P,PSIp,E_R,E_PHI,E_Z) ! write(output_unit_write,*) params%mpi_params%rank,'Y_R',Y_R @@ -4595,7 +4600,9 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc) spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc) + spp(ii)%vars%E(pp-1+cc,1) = E_R(cc) spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) + spp(ii)%vars%E(pp-1+cc,3) = E_Z(cc) end do !$OMP END SIMD @@ -4631,17 +4638,24 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) end if - call analytical_fields_Bmag_p(pchunk,F,Y_R,Y_PHI,Y_Z, & - Bmag,E_PHI) + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_Bmag_p(pchunk,F,Y_R,Y_PHI,Y_Z, & + Bmag,E_PHI) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z, & + E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) + + Bmag=sqrt(B_R*B_R+B_PHI*B_PHI+B_Z*B_Z) + end if !$OMP SIMD do cc=1_idef,pchunk - spp(ii)%vars%g(pp-1+cc)=sqrt(1+V_PLL(cc)**2+ & - 2*V_MU(cc)*Bmag(cc)*m_cache) + spp(ii)%vars%g(pp-1+cc)=sqrt(1+(V_PLL(cc)/m_cache)**2+ & + 2*V_MU(cc)*Bmag(cc)/m_cache) spp(ii)%vars%eta(pp-1+cc) = rad2deg(atan2(sqrt(2*m_cache* & - Bmag(cc)*spp(ii)%vars%V(pp-1+cc,2)), & - spp(ii)%vars%V(pp-1+cc,1))) + Bmag(cc)*V_MU(cc)),V_PLL(cc))) end do !$OMP END SIMD @@ -4715,8 +4729,10 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !$OMP PARALLEL DO default(none) & !$OMP& FIRSTPRIVATE(m_cache,pchunk) & - !$OMP& shared(F,ii,spp) & - !$OMP& PRIVATE(pp,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,achunk,E_PHI) + !$OMP& shared(F,ii,spp,params) & + !$OMP& PRIVATE(pp,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,achunk, & + !$OMP& E_R,E_PHI,E_Z,gradB_R,gradB_PHI,gradB_Z, & + !$OMP& curlb_R,curlb_PHI,curlb_Z,PSIp,B_R,B_PHI,B_Z) do pp=1_idef,spp(ii)%pRE,pchunk if ((spp(ii)%pRE-pp).lt.pchunk) then @@ -4734,21 +4750,29 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !$OMP END SIMD !write(6,*) 'Y_R',Y_R + + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_Bmag_p(pchunk,F,Y_R,Y_PHI,Y_Z, & + Bmag,E_PHI) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z, & + E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) - call analytical_fields_Bmag_p(achunk,F,Y_R,Y_PHI,Y_Z, & - Bmag,E_PHI) + Bmag=sqrt(B_R*B_R+B_PHI*B_PHI+B_Z*B_Z) + end if !$OMP SIMD do cc=1_idef,achunk V_PLL(cc)=spp(ii)%vars%V(pp-1+cc,1) V_MU(cc)=spp(ii)%vars%V(pp-1+cc,2) - spp(ii)%vars%g(pp-1+cc)=sqrt(1+V_PLL(cc)**2+ & - 2*V_MU(cc)*Bmag(cc)*m_cache) + + spp(ii)%vars%g(pp-1+cc)=sqrt(1+(V_PLL(cc)/m_cache)**2+ & + 2*V_MU(cc)*Bmag(cc)/m_cache) spp(ii)%vars%eta(pp-1+cc) = rad2deg(atan2(sqrt(2*m_cache* & - Bmag(cc)*spp(ii)%vars%V(pp-1+cc,2)), & - spp(ii)%vars%V(pp-1+cc,1))) + Bmag(cc)*V_MU(cc)),V_PLL(cc))) end do !$OMP END SIMD @@ -4763,7 +4787,7 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !$OMP& FIRSTPRIVATE(m_cache,pchunk,q_cache) & !$OMP& shared(F,P,params,ii,spp,tt,random) & !$OMP& PRIVATE(pp,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & - !$OMP& flagCon,flagCol,B_R,B_PHI,B_Z,E_PHI,PSIp,ne, & + !$OMP& flagCon,flagCol,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,PSIp,ne, & !$OMP& achunk,tttt,Te) do pp=1_idef,spp(ii)%pRE,pchunk @@ -4798,10 +4822,9 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) if (.not.params%FokPlan) then do tttt=1_ip,params%orbits_per_coll - call advance_GCeqn_vars(spp(ii)%vars,pp, & - tttt,params, & + call advance_GCeqn_vars(spp(ii)%vars,pp,tttt,params, & Y_R,Y_PHI,Y_Z,V_PLL,V_MU,flagCon,flagCol,q_cache,m_cache, & - B_R,B_PHI,B_Z,F,P,PSIp,E_PHI) + B_R,B_PHI,B_Z,F,P,PSIp,E_R,E_PHI,E_Z) end do endif @@ -4829,7 +4852,9 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc) spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc) + spp(ii)%vars%E(pp-1+cc,1) = E_R(cc) spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) + spp(ii)%vars%E(pp-1+cc,3) = E_Z(cc) end do !$OMP END SIMD @@ -4840,8 +4865,10 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !$OMP PARALLEL DO default(none) & !$OMP& FIRSTPRIVATE(m_cache,pchunk) & - !$OMP& shared(F,ii,spp) & - !$OMP& PRIVATE(pp,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,achunk,E_PHI) + !$OMP& shared(F,ii,spp,params) & + !$OMP& PRIVATE(pp,Bmag,cc,Y_R,Y_PHI,Y_Z,V_PLL,V_MU,achunk, & + !$OMP& E_R,E_PHI,E_Z,gradB_R,gradB_PHI,gradB_Z, & + !$OMP& curlb_R,curlb_PHI,curlb_Z,PSIp,B_R,B_PHI,B_Z) do pp=1_idef,spp(ii)%pRE,pchunk if ((spp(ii)%pRE-pp).lt.pchunk) then @@ -4859,21 +4886,28 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) !$OMP END SIMD !write(6,*) 'Y_R',Y_R + + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_Bmag_p(pchunk,F,Y_R,Y_PHI,Y_Z, & + Bmag,E_PHI) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z, & + E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) - call analytical_fields_Bmag_p(achunk,F,Y_R,Y_PHI,Y_Z, & - Bmag,E_PHI) - + Bmag=sqrt(B_R*B_R+B_PHI*B_PHI+B_Z*B_Z) + end if !$OMP SIMD do cc=1_idef,achunk V_PLL(cc)=spp(ii)%vars%V(pp-1+cc,1) V_MU(cc)=spp(ii)%vars%V(pp-1+cc,2) - spp(ii)%vars%g(pp-1+cc)=sqrt(1+V_PLL(cc)**2+ & - 2*V_MU(cc)*Bmag(cc)*m_cache) + + spp(ii)%vars%g(pp-1+cc)=sqrt(1+(V_PLL(cc)/m_cache)**2+ & + 2*V_MU(cc)*Bmag(cc)/m_cache) spp(ii)%vars%eta(pp-1+cc) = rad2deg(atan2(sqrt(2*m_cache* & - Bmag(cc)*spp(ii)%vars%V(pp-1+cc,2)), & - spp(ii)%vars%V(pp-1+cc,1))) + Bmag(cc)*V_MU(cc)),V_PLL(cc))) end do !$OMP END SIMD @@ -4890,8 +4924,8 @@ subroutine adv_GCeqn_top(params,random,F,P,spp) end subroutine adv_GCeqn_top -subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & - flagCon,flagCol,q_cache,m_cache,B_R,B_PHI,B_Z,F,P,PSIp,E_PHI) + subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & + flagCon,flagCol,q_cache,m_cache,B_R,B_PHI,B_Z,F,P,PSIp,E_R,E_PHI,E_Z) !! @note Subroutine to advance GC variables \(({\bf X},p_\parallel)\) !! @endnote !! Comment this section further with evolution equations, numerical @@ -4930,8 +4964,8 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: PSIp REAL(rp),DIMENSION(params%pchunk) :: curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z REAL(rp),DIMENSION(params%pchunk),INTENT(INOUT) :: V_PLL,V_MU - REAL(rp),DIMENSION(params%pchunk) :: RHS_R,RHS_PHI,RHS_Z,RHS_PLL,V0,E_Z,E_R - REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: E_PHI + REAL(rp),DIMENSION(params%pchunk) :: RHS_R,RHS_PHI,RHS_Z,RHS_PLL,V0 + REAL(rp),DIMENSION(params%pchunk),INTENT(OUT) :: E_PHI,E_Z,E_R REAL(rp),DIMENSION(params%pchunk) :: Bmag,ne,Te,Zeff INTEGER(is),dimension(params%pchunk), intent(inout) :: flagCon,flagCol @@ -4957,10 +4991,14 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & V0(cc)=V_PLL(cc) end do !$OMP END SIMD - - call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & - Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & - gradB_R,gradB_PHI,gradB_Z,PSIp) + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & + Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,PSIp) + end if ! write(output_unit_write,'("ER:",E17.10)') E_R ! write(output_unit_write,'("EPHI:",E17.10)') E_PHI @@ -5002,9 +5040,14 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & ! write(output_unit_write,'("k1Z: ",E17.10)') k1_Z(1) ! write(output_unit_write,'("k1PLL: ",E17.10)') k1_PLL(1) - call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & - Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & - gradB_R,gradB_PHI,gradB_Z,PSIp) + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & + Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,PSIp) + end if call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, & B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & @@ -5031,9 +5074,14 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & ! write(output_unit_write,'("Y_PHI 2: ",E17.10)') Y_PHI(1) ! write(output_unit_write,'("Y_Z 2: ",E17.10)') Y_Z(1) - call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & - Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & - gradB_R,gradB_PHI,gradB_Z,PSIp) + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & + Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,PSIp) + end if call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, & B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & @@ -5060,9 +5108,14 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & ! write(output_unit_write,'("Y_PHI 3: ",E17.10)') Y_PHI(1) ! write(output_unit_write,'("Y_Z 3: ",E17.10)') Y_Z(1) - call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & - Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & - gradB_R,gradB_PHI,gradB_Z,PSIp) + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & + Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,PSIp) + end if call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, & B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & @@ -5092,9 +5145,14 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & ! write(output_unit_write,'("Y_PHI 4: ",E17.10)') Y_PHI(1) ! write(output_unit_write,'("Y_Z 4: ",E17.10)') Y_Z(1) - call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & - Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & - gradB_R,gradB_PHI,gradB_Z,PSIp) + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & + Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,PSIp) + end if call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, & B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & @@ -5124,9 +5182,14 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & ! write(output_unit_write,'("Y_PHI 5: ",E17.10)') Y_PHI(1) ! write(output_unit_write,'("Y_Z 5: ",E17.10)') Y_Z(1) - call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & - Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & - gradB_R,gradB_PHI,gradB_Z,PSIp) + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & + Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,PSIp) + end if call GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, & B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z,gradB_R, & @@ -5156,7 +5219,9 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & ! write(output_unit_write,'("Y_PHI 6: ",E17.10)') Y_PHI(1) ! write(output_unit_write,'("Y_Z 6: ",E17.10)') Y_Z(1) - call cyl_check_if_confined_p(pchunk,ar,R0,Y_R,Y_Z,flagCon) + if (params%field_model(1:3).eq.'ANA') then + call cyl_check_if_confined_p(pchunk,ar,R0,Y_R,Y_Z,flagCon) + end if !$OMP SIMD ! !$OMP& aligned(Y_R,Y_PHI,Y_Z,V_PLL,Y0_R,Y0_PHI,Y0_Z,V0) @@ -5172,9 +5237,14 @@ subroutine advance_GCeqn_vars(vars,pp,tt,params,Y_R,Y_PHI,Y_Z,V_PLL,V_MU, & end do !$OMP END SIMD - call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & - Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & - gradB_R,gradB_PHI,gradB_Z,PSIp) + if (params%field_model(1:3).eq.'ANA') then + call analytical_fields_GC_p(pchunk,F,Y_R,Y_PHI, & + Y_Z,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,PSIp) + else if (params%field_model(1:3).eq.'UNI') then + call uniform_fields_GC_p(pchunk,F,B_R,B_PHI,B_Z,E_R,E_PHI,E_Z, & + curlb_R,curlb_PHI,curlb_Z,gradB_R,gradB_PHI,gradB_Z,PSIp) + end if end subroutine advance_GCeqn_vars @@ -5208,7 +5278,7 @@ end subroutine advance_FPeqn_vars #ifdef PSPLINE -subroutine adv_GCinterp_psi_top(params,random,spp,P,F) + subroutine adv_GCinterp_psi_top(params,random,spp,P,F) TYPE(KORC_PARAMS), INTENT(INOUT) :: params !! Core KORC simulation parameters. @@ -5691,12 +5761,13 @@ subroutine adv_GCinterp_fio_top(params,spp,P,F) Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ & B_Z(cc)*B_Z(cc)) - spp(ii)%vars%g(pp-1+cc)=sqrt(1+V_PLL(cc)**2+ & - 2*V_MU(cc)*Bmag(cc)) - spp(ii)%vars%eta(pp-1+cc) = atan2(sqrt(2*m_cache*Bmag(cc)* & - spp(ii)%vars%V(pp-1+cc,2)),spp(ii)%vars%V(pp-1+cc,1))* & - 180.0_rp/C_PI + spp(ii)%vars%g(pp-1+cc)=sqrt(1+(V_PLL(cc)/m_cache)**2+ & + 2*V_MU(cc)*Bmag(cc)/m_cache) + + spp(ii)%vars%eta(pp-1+cc) = rad2deg(atan2(sqrt(2*m_cache* & + Bmag(cc)*V_MU(cc)),V_PLL(cc))) + end do !$OMP END SIMD @@ -5916,12 +5987,16 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ & B_Z(cc)*B_Z(cc)) - spp(ii)%vars%g(pp-1+cc)=sqrt(1+V_PLL(cc)**2+ & - 2*V_MU(cc)*Bmag(cc)) + V_PLL(cc)=spp(ii)%vars%V(pp-1+cc,1) + V_MU(cc)=spp(ii)%vars%V(pp-1+cc,2) + + + spp(ii)%vars%g(pp-1+cc)=sqrt(1+(V_PLL(cc)/m_cache)**2+ & + 2*V_MU(cc)*Bmag(cc)/m_cache) + + spp(ii)%vars%eta(pp-1+cc) = rad2deg(atan2(sqrt(2*m_cache* & + Bmag(cc)*V_MU(cc)),V_PLL(cc))) - spp(ii)%vars%eta(pp-1+cc) = atan2(sqrt(2*m_cache*Bmag(cc)* & - spp(ii)%vars%V(pp-1+cc,2)),spp(ii)%vars%V(pp-1+cc,1))* & - 180.0_rp/C_PI end do !$OMP END SIMD @@ -6127,45 +6202,43 @@ subroutine adv_GCinterp_psiwE_top(params,random,spp,P,F) end do !timestep iterator - !$OMP PARALLEL DO default(none) & - !$OMP& FIRSTPRIVATE(m_cache,pchunk) & - !$OMP& SHARED(ii,spp) & - !$OMP& PRIVATE(pp,Bmag,cc, & - !$OMP& B_R,B_PHI,B_Z,achunk) - - do pp=1_idef,spp(ii)%pRE,pchunk - - if ((spp(ii)%pRE-pp).lt.pchunk) then - achunk=spp(ii)%pRE-pp+1 - else - achunk=pchunk - end if + !$OMP PARALLEL DO default(none) & + !$OMP& FIRSTPRIVATE(m_cache,pchunk) & + !$OMP& SHARED(ii,spp) & + !$OMP& PRIVATE(pp,Bmag,cc, & + !$OMP& B_R,B_PHI,B_Z,achunk,V_PLL,V_MU) + do pp=1_idef,spp(ii)%pRE,pchunk + + if ((spp(ii)%pRE-pp).lt.pchunk) then + achunk=spp(ii)%pRE-pp+1 + else + achunk=pchunk + end if - !$OMP SIMD - do cc=1_idef,achunk - B_R(cc)=spp(ii)%vars%B(pp-1+cc,1) - B_PHI(cc)=spp(ii)%vars%B(pp-1+cc,2) - B_Z(cc)=spp(ii)%vars%B(pp-1+cc,3) + !$OMP SIMD + do cc=1_idef,achunk + B_R(cc)=spp(ii)%vars%B(pp-1+cc,1) + B_PHI(cc)=spp(ii)%vars%B(pp-1+cc,2) + B_Z(cc)=spp(ii)%vars%B(pp-1+cc,3) - Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ & - B_Z(cc)*B_Z(cc)) + Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ & + B_Z(cc)*B_Z(cc)) - spp(ii)%vars%g(pp-1+cc)=sqrt(1+spp(ii)%vars%V(pp-1+cc,1)**2 & - +2*spp(ii)%vars%V(pp-1+cc,2)*Bmag(cc)) + spp(ii)%vars%V(pp-1+cc,1)=V_PLL(cc) + spp(ii)%vars%V(pp-1+cc,2)=V_MU(cc) - spp(ii)%vars%eta(pp-1+cc) = atan2(sqrt(2*m_cache*Bmag(cc)* & - spp(ii)%vars%V(pp-1+cc,2)),spp(ii)%vars%V(pp-1+cc,1))* & - 180.0_rp/C_PI - end do - !$OMP END SIMD + spp(ii)%vars%g(pp-1+cc)=sqrt(1+(V_PLL(cc)/m_cache)**2+ & + 2*V_MU(cc)*Bmag(cc)/m_cache) - end do !particle chunk iterator - !$OMP END PARALLEL DO + spp(ii)%vars%eta(pp-1+cc) = rad2deg(atan2(sqrt(2*m_cache* & + Bmag(cc)*V_MU(cc)),V_PLL(cc))) + end do + !$OMP END SIMD - endif + end do !particle chunk iterator + !$OMP END PARALLEL DO - !write(6,*) 'Out Y',spp(ii)%vars%Y(1,1)*params%cpp%length,spp(ii)%vars%Y(1,2),spp(ii)%vars%Y(1,3)*params%cpp%length - !write(6,*) 'Out Y0',spp(ii)%vars%Y0(1,1)*params%cpp%length,spp(ii)%vars%Y0(1,2),spp(ii)%vars%Y0(1,3)*params%cpp%length + endif end do !species iterator @@ -6284,11 +6357,11 @@ subroutine adv_GCinterp_psi2x1t_top(params,random,spp,P,F) else if (params%FokPlan.and.params%collisions) then - call advance_FPinterp_vars(params,random,Y_R,Y_PHI, & - Y_Z,V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,PSIp) + call advance_FPinterp_vars(params,random,Y_R,Y_PHI, & + Y_Z,V_PLL,V_MU,m_cache,flagCon,flagCol,F,P,E_PHI,ne,PSIp) - !$OMP SIMD - do cc=1_idef,pchunk + !$OMP SIMD + do cc=1_idef,pchunk spp(ii)%vars%V(pp-1+cc,1)=V_PLL(cc) spp(ii)%vars%V(pp-1+cc,2)=V_MU(cc) @@ -6297,61 +6370,58 @@ subroutine adv_GCinterp_psi2x1t_top(params,random,spp,P,F) spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) spp(ii)%vars%ne(pp-1+cc) = ne(cc) - end do - !$OMP END SIMD + end do + !$OMP END SIMD - else - do tt=1_ip,params%t_skip + else + do tt=1_ip,params%t_skip time=params%init_time+(params%it-1+tt)* & - params%dt + params%dt call calculate_GCfields_2x1t_p(pchunk,F,Y_R,Y_PHI,Y_Z,B_R,B_PHI,B_Z, & - E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & - gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp,time) - end do + E_R,E_PHI,E_Z,curlb_R,curlb_PHI,curlb_Z, & + gradB_R,gradB_PHI,gradB_Z,flagCon,PSIp,time) + end do - !$OMP SIMD - do cc=1_idef,pchunk + !$OMP SIMD + do cc=1_idef,pchunk spp(ii)%vars%B(pp-1+cc,1) = B_R(cc) spp(ii)%vars%B(pp-1+cc,2) = B_PHI(cc) spp(ii)%vars%B(pp-1+cc,3) = B_Z(cc) spp(ii)%vars%E(pp-1+cc,2) = E_PHI(cc) spp(ii)%vars%PSI_P(pp-1+cc) = PSIp(cc) - end do - !$OMP END SIMD - - end if - + end do + !$OMP END SIMD - !$OMP SIMD - do cc=1_idef,pchunk - B_R(cc)=spp(ii)%vars%B(pp-1+cc,1) - B_PHI(cc)=spp(ii)%vars%B(pp-1+cc,2) - B_Z(cc)=spp(ii)%vars%B(pp-1+cc,3) + end if - Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ & - B_Z(cc)*B_Z(cc)) - spp(ii)%vars%g(pp-1+cc)=sqrt(1+V_PLL(cc)**2+ & - 2*V_MU(cc)*Bmag(cc)) + !$OMP SIMD + do cc=1_idef,pchunk + B_R(cc)=spp(ii)%vars%B(pp-1+cc,1) + B_PHI(cc)=spp(ii)%vars%B(pp-1+cc,2) + B_Z(cc)=spp(ii)%vars%B(pp-1+cc,3) - spp(ii)%vars%eta(pp-1+cc) = atan2(sqrt(2*m_cache*Bmag(cc)* & - spp(ii)%vars%V(pp-1+cc,2)),spp(ii)%vars%V(pp-1+cc,1))* & - 180.0_rp/C_PI - end do - !$OMP END SIMD + Bmag(cc)=sqrt(B_R(cc)*B_R(cc)+B_PHI(cc)*B_PHI(cc)+ & + B_Z(cc)*B_Z(cc)) - end do !particle chunk iterator - !$OMP END PARALLEL DO + spp(ii)%vars%g(pp-1+cc)=sqrt(1+V_PLL(cc)**2+ & + 2*V_MU(cc)*Bmag(cc)) + spp(ii)%vars%eta(pp-1+cc) = atan2(sqrt(2*m_cache*Bmag(cc)* & + spp(ii)%vars%V(pp-1+cc,2)),spp(ii)%vars%V(pp-1+cc,1))* & + 180.0_rp/C_PI + end do + !$OMP END SIMD + end do !particle chunk iterator + !$OMP END PARALLEL DO end do !species iterator end subroutine adv_GCinterp_psi2x1t_top - subroutine advance_GCinterp_psi_vars(pchunk,spp,pp,tt,params,random,Y_R,Y_PHI,Y_Z, & V_PLL,V_MU,q_cache,m_cache,flagCon,flagCol,F,P,B_R,B_PHI,B_Z,E_PHI,PSIp, & @@ -8124,9 +8194,9 @@ subroutine GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, & bhat_PHI(cc) = B_PHI(cc)/Bmag(cc) bhat_Z(cc) = B_Z(cc)/Bmag(cc) - Bst_R(cc)=q_cache*B_R(cc)+V_PLL(cc)*curlb_R(cc) - Bst_PHI(cc)=q_cache*B_PHI(cc)+V_PLL(cc)*curlb_PHI(cc) - Bst_Z(cc)=q_cache*B_Z(cc)+V_PLL(cc)*curlb_Z(cc) + Bst_R(cc)=q_cache*(B_R(cc)+V_PLL(cc)*curlb_R(cc)/m_cache) + Bst_PHI(cc)=q_cache*(B_PHI(cc)+V_PLL(cc)*curlb_PHI(cc)/m_cache) + Bst_Z(cc)=q_cache*(B_Z(cc)+V_PLL(cc)*curlb_Z(cc)/m_cache) bdotBst(cc)=bhat_R(cc)*Bst_R(cc)+bhat_PHI(cc)*Bst_PHI(cc)+ & bhat_Z(cc)*Bst_Z(cc) @@ -8143,21 +8213,21 @@ subroutine GCEoM_p(params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,B_R,B_PHI, & bcrossgradB_PHI(cc)=bhat_Z(cc)*gradB_R(cc)-bhat_R(cc)*gradB_Z(cc) bcrossgradB_Z(cc)=bhat_R(cc)*gradB_PHI(cc)-bhat_PHI(cc)*gradB_R(cc) - gamgc(cc)=sqrt(1+V_PLL(cc)*V_PLL(cc)+2*V_MU(cc)*Bmag(cc)) + gamgc(cc)=sqrt(1+(V_PLL(cc)/m_cache)**2+2*V_MU(cc)*Bmag(cc)/m_cache) - pm(cc)=sqrt(gamgc(cc)**2-1) + pm(cc)=m_cache*sqrt(gamgc(cc)**2-1) xi(cc)=V_PLL(cc)/pm(cc) - RHS_R(cc)=(q_cache*Ecrossb_R(cc)+(m_cache*V_MU(cc)* & + RHS_R(cc)=(q_cache*Ecrossb_R(cc)+(q_cache*V_MU(cc)* & bcrossgradB_R(cc)+V_PLL(cc)*Bst_R(cc))/(m_cache*gamgc(cc)))/ & bdotBst(cc) - RHS_PHI(cc)=(q_cache*Ecrossb_PHI(cc)+(m_cache*V_MU(cc)* & + RHS_PHI(cc)=(q_cache*Ecrossb_PHI(cc)+(q_cache*V_MU(cc)* & bcrossgradB_PHI(cc)+V_PLL(cc)*Bst_PHI(cc))/(m_cache*gamgc(cc)))/ & (Y_R(cc)*bdotBst(cc)) - RHS_Z(cc)=(q_cache*Ecrossb_Z(cc)+(m_cache*V_MU(cc)* & + RHS_Z(cc)=(q_cache*Ecrossb_Z(cc)+(q_cache*V_MU(cc)* & bcrossgradB_Z(cc)+V_PLL(cc)*Bst_Z(cc))/(m_cache*gamgc(cc)))/ & bdotBst(cc) - RHS_PLL(cc)=(q_cache*BstdotE(cc)-V_MU(cc)*BstdotgradB(cc)/gamgc(cc))/ & + RHS_PLL(cc)=(sign(m_cache,q_cache)*BstdotE(cc)-V_MU(cc)*BstdotgradB(cc)/gamgc(cc))/ & bdotBst(cc) end do @@ -8221,9 +8291,9 @@ subroutine GCEoM1_p(pchunk,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, & bhat_PHI(cc) = B_PHI(cc)/Bmag(cc) bhat_Z(cc) = B_Z(cc)/Bmag(cc) - Bst_R(cc)=q_cache*B_R(cc)+V_PLL(cc)*curlb_R(cc) - Bst_PHI(cc)=q_cache*B_PHI(cc)+V_PLL(cc)*curlb_PHI(cc) - Bst_Z(cc)=q_cache*B_Z(cc)+V_PLL(cc)*curlb_Z(cc) + Bst_R(cc)=q_cache*(B_R(cc)+V_PLL(cc)*curlb_R(cc)/m_cache) + Bst_PHI(cc)=q_cache*(B_PHI(cc)+V_PLL(cc)*curlb_PHI(cc)/m_cache) + Bst_Z(cc)=q_cache*(B_Z(cc)+V_PLL(cc)*curlb_Z(cc)/m_cache) ! write(output_unit_write,*) 'bmag',Bmag(cc),'bhat',bhat_R(cc),bhat_PHI(cc),bhat_Z(cc),'Bst',Bst_R(cc),Bst_PHI(cc),Bst_Z(cc) @@ -8246,22 +8316,22 @@ subroutine GCEoM1_p(pchunk,tt,P,F,params,RHS_R,RHS_PHI,RHS_Z,RHS_PLL,RHS_MU, & bcrossgradB_Z(cc)=bhat_R(cc)*gradB_PHI(cc)-bhat_PHI(cc)*gradB_R(cc) ! write(output_unit_write,*) 'bcrossgradB',bcrossgradB_R(cc),bcrossgradB_PHI(cc),bcrossgradB_Z(cc) - - gamgc(cc)=sqrt(1+V_PLL(cc)*V_PLL(cc)+2*V_MU(cc)*Bmag(cc)) - - pm(cc)=sqrt(gamgc(cc)**2-1) + + gamgc(cc)=sqrt(1+(V_PLL(cc)/m_cache)**2+2*V_MU(cc)*Bmag(cc)/m_cache) + + pm(cc)=m_cache*sqrt(gamgc(cc)**2-1) xi(cc)=V_PLL(cc)/pm(cc) - RHS_R(cc)=(q_cache*Ecrossb_R(cc)+(m_cache*V_MU(cc)* & + RHS_R(cc)=(q_cache*Ecrossb_R(cc)+(q_cache*V_MU(cc)* & bcrossgradB_R(cc)+V_PLL(cc)*Bst_R(cc))/(m_cache*gamgc(cc)))/ & bdotBst(cc) - RHS_PHI(cc)=(q_cache*Ecrossb_PHI(cc)+(m_cache*V_MU(cc)* & + RHS_PHI(cc)=(q_cache*Ecrossb_PHI(cc)+(q_cache*V_MU(cc)* & bcrossgradB_PHI(cc)+V_PLL(cc)*Bst_PHI(cc))/(m_cache*gamgc(cc)))/ & (Y_R(cc)*bdotBst(cc)) - RHS_Z(cc)=(q_cache*Ecrossb_Z(cc)+(m_cache*V_MU(cc)* & + RHS_Z(cc)=(q_cache*Ecrossb_Z(cc)+(q_cache*V_MU(cc)* & bcrossgradB_Z(cc)+V_PLL(cc)*Bst_Z(cc))/(m_cache*gamgc(cc)))/ & bdotBst(cc) - RHS_PLL(cc)=(q_cache*BstdotE(cc)-V_MU(cc)*BstdotgradB(cc)/gamgc(cc))/ & + RHS_PLL(cc)=(m_cache*BstdotE(cc)-V_MU(cc)*BstdotgradB(cc)/gamgc(cc))/ & bdotBst(cc) RHS_MU(cc)=0._rp diff --git a/src/korc_spatial_distribution.f90 b/src/korc_spatial_distribution.f90 index 76c69212..4b935385 100755 --- a/src/korc_spatial_distribution.f90 +++ b/src/korc_spatial_distribution.f90 @@ -8,7 +8,6 @@ MODULE korc_spatial_distribution use korc_fields use korc_profiles use korc_random - use korc_hammersley_generator use korc_avalanche use korc_experimental_pdf USE korc_input @@ -30,6 +29,7 @@ MODULE korc_spatial_distribution fzero,& MH_gaussian_elliptic_torus,& indicator,& + indicatortwo,& PSI_ROT,& Spong_3D,& Spong_2D, & @@ -711,6 +711,31 @@ FUNCTION indicator(psi,psi_max) END FUNCTION indicator + +FUNCTION indicatortwo(psi,PSIp_min) + + REAL(rp), INTENT(IN) :: psi + REAL(rp), INTENT(IN) :: PSIp_min + REAL(rp) :: indicatortwo + + IF (psi.GT.PSIp_min) THEN + indicatortwo=1 + ELSE + indicatortwo=0 + END IF + END FUNCTION indicatortwo +FUNCTION random_norm(mean,sigma) + REAL(rp), INTENT(IN) :: mean + REAL(rp), INTENT(IN) :: sigma + REAL(rp) :: random_norm + REAL(rp) :: rand1, rand2 + + call RANDOM_NUMBER(rand1) + call RANDOM_NUMBER(rand2) + + random_norm = mean+sigma*SQRT(-2.0_rp*LOG(rand1))*COS(2.0_rp*C_PI*rand2); +END FUNCTION random_norm + subroutine MH_gaussian_elliptic_torus(params,random,spp) !! @note Subroutine that generates a 2D Gaussian distribution in an !! elliptic torus as the initial spatial condition of a given particle @@ -1301,7 +1326,7 @@ subroutine MH_psi(params,random,spp,F) params%GC_coords=.TRUE. PSIp_lim=F%PSIp_lim - + !changed_YG if (params%field_model.eq.'M3D_C1') then min_R=params%rmin/params%cpp%length max_R=params%rmax/params%cpp%length @@ -1322,11 +1347,15 @@ subroutine MH_psi(params,random,spp,F) psi_max = spp%psi_max psi_min = spp%psi_min psi_max_buff = spp%psi_max*2._rp + !PSIp_min=spp%PSIp_min/(params%cpp%Bo*params%cpp%length**2) + ! psi_max=(psi_max-PSIp0)/(PSIp_lim-PSIp0) + ! PSIp_min=(PSIp_min-PSIp0)/(PSIp_lim-PSIp0) +! write(6,*)"psip0",PSIp0,"psi_max",psi_max,"PSIp_min",PSIp_min,"PSIp_lim", PSIp_lim end if sigma=spp%sigmaR*params%cpp%length - !write(output_unit_write,*) min_R,max_R + !write(6,*) "sigma",sigma !write(output_unit_write,*) min_Z,max_Z @@ -1457,7 +1486,6 @@ subroutine MH_psi(params,random,spp,F) indicator(PSIN1,psi_max)*indicator(psi_min,PSIN1)* & R_test*EXP(-PSIN1/sigma)/ & (R_buffer*EXP(-PSIN0/sigma)) - ! ratio = f1*sin(deg2rad(eta_test))/(f0*sin(deg2rad(eta_buffer))) accepted=.false. diff --git a/src/korc_types.f90 b/src/korc_types.f90 index 41dfa7f1..17ae17d6 100755 --- a/src/korc_types.f90 +++ b/src/korc_types.f90 @@ -623,7 +623,10 @@ module korc_types REAL(rp) :: MARS_phase REAL(rp) :: AORSA_AMP_Scale REAL(rp) :: AORSA_freq - REAL(rp) :: AORSA_nmode,AORSA_mmode + REAL(rp) :: psir + REAL(rp) :: width + REAL(rp) :: AORSA_nmode + REAL(rp) :: AORSA_mmode !! interpolated E field LOGICAL :: Analytic_D3D_IWL INTEGER :: ntiles diff --git a/src/korc_units.f90 b/src/korc_units.f90 index 578f14cd..896403b2 100755 --- a/src/korc_units.f90 +++ b/src/korc_units.f90 @@ -151,7 +151,8 @@ subroutine normalize_variables(params,spp,F,P) F%Bo = F%Bo/params%cpp%Bo F%Eo = F%Eo/params%cpp%Eo F%AB%Ero = F%AB%Ero/params%cpp%Eo - F%AB%rmn=F%AB%rmn/params%cpp%length + !F%width=F%width/(params%cpp%length) + F%AB%rmn=F%AB%rmn/(params%cpp%Bo*params%cpp%length**2) F%AB%sigmamn=F%AB%sigmamn/params%cpp%length F%Ro = F%Ro/params%cpp%length F%Zo = F%Zo/params%cpp%length @@ -407,7 +408,9 @@ subroutine normalize_variables(params,spp,F,P) params%cpp%Bo else if (params%field_model(10:14).eq.'AORSA') then - + F%AORSA_freq=F%AORSA_freq*params%cpp%time + F%psir=F%psir/(params%cpp%Bo*params%cpp%length**2) + F%width=F%width/(params%cpp%Bo*params%cpp%length**2) if (ALLOCATED(F%B1Re_2DX%X)) F%B1Re_2DX%X = F%B1Re_2DX%X/ & params%cpp%Bo if (ALLOCATED(F%B1Re_2DX%Y)) F%B1Re_2DX%Y = F%B1Re_2DX%Y/ & diff --git a/src/korc_velocity_distribution.f90 b/src/korc_velocity_distribution.f90 index aa43eb18..4685242e 100755 --- a/src/korc_velocity_distribution.f90 +++ b/src/korc_velocity_distribution.f90 @@ -7,7 +7,6 @@ MODULE korc_velocity_distribution USE korc_hpc use korc_fields use korc_random - use korc_hammersley_generator use korc_avalanche use korc_experimental_pdf @@ -229,8 +228,10 @@ subroutine initial_energy_pitch_dist(params,random,spp) spp(ii)%Eo = spp(ii)%Eo_lims(1) spp(ii)%go = (spp(ii)%Eo + spp(ii)%m*C_C**2)/(spp(ii)%m*C_C**2) - call generate_2D_hammersley_sequence(params%mpi_params%rank, & - params%mpi_params%nmpi,spp(ii)%vars%g,spp(ii)%vars%eta) + CALL random%uniform%set(0.0_rp,1.0_rp) + do pp=1_idef,spp(ii)%ppp + spp(ii)%vars%g(pp)=random%uniform%get() + enddo spp(ii)%vars%g = (spp(ii)%Eo_lims(2) - & spp(ii)%Eo_lims(1))*spp(ii)%vars%g/(spp(ii)%m*C_C**2) + & @@ -281,6 +282,11 @@ subroutine initial_energy_pitch_dist(params,random,spp) CASE ('UNIFORM') spp(ii)%etao = spp(ii)%etao_lims(1) + CALL random%uniform%set(0.0_rp,1.0_rp) + do pp=1_idef,spp(ii)%ppp + spp(ii)%vars%eta(pp)=random%uniform%get() + enddo + spp(ii)%vars%eta = (spp(ii)%etao_lims(2) - & spp(ii)%etao_lims(1))*spp(ii)%vars%eta + spp(ii)%etao_lims(1) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d72a8852..873d5e40 100755 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -22,6 +22,7 @@ set_property(TARGET xtest PROPERTY LINKER_LANGUAGE Fortran) configure_file(${CMAKE_SOURCE_DIR}/test/egyro/korc_egyro.sh.in ${CMAKE_BINARY_DIR}/egyro_test/korc_egyro.sh) configure_file(${CMAKE_SOURCE_DIR}/test/mars/korc_mars.sh.in ${CMAKE_BINARY_DIR}/mars_test/korc_mars.sh) +configure_file(${CMAKE_SOURCE_DIR}/test/aorsa/korc_aorsa.sh.in ${CMAKE_BINARY_DIR}/aorsa_test/korc_aorsa.sh) #if(USE_PSPLINE) # add_test (NAME mars_test_1 @@ -32,6 +33,7 @@ configure_file(${CMAKE_SOURCE_DIR}/test/mars/korc_mars.sh.in ${CMAKE_BINARY_DIR} file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/egyro_test) file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/mars_test) +file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/bin/aorsa_test) foreach(RANK IN ITEMS 1 2 4 8) #16) if(${RANK} LESS_EQUAL ${MPIEXEC_MAX_NUMPROCS}) @@ -45,6 +47,11 @@ foreach(RANK IN ITEMS 1 2 4 8) #16) COMMAND ${CMAKE_BINARY_DIR}/mars_test/korc_mars.sh ${RANK} WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bin) set_tests_properties (mars_test_${RANK} PROPERTIES PROCESSORS ${RANK} ENVIRONMENT OMP_NUM_THREADS=1) + + add_test (NAME aorsa_test_${RANK} + COMMAND ${CMAKE_BINARY_DIR}/aorsa_test/korc_aorsa.sh ${RANK} + WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bin) + set_tests_properties (aorsa_test_${RANK} PROPERTIES PROCESSORS ${RANK} ENVIRONMENT OMP_NUM_THREADS=1) endif() add_test (NAME unit_testing_${RANK} diff --git a/test/aorsa/AORSA_D3D_171089_200MHz_EFITgrid.h5 b/test/aorsa/AORSA_D3D_171089_200MHz_EFITgrid.h5 new file mode 100644 index 00000000..c0a17817 Binary files /dev/null and b/test/aorsa/AORSA_D3D_171089_200MHz_EFITgrid.h5 differ diff --git a/test/aorsa/file_0.h5 b/test/aorsa/file_0.h5 new file mode 100644 index 00000000..4ced9f84 Binary files /dev/null and b/test/aorsa/file_0.h5 differ diff --git a/test/aorsa/input_file_D3D_171089_AORSA.korc b/test/aorsa/input_file_D3D_171089_AORSA.korc new file mode 100644 index 00000000..b0d931ec --- /dev/null +++ b/test/aorsa/input_file_D3D_171089_AORSA.korc @@ -0,0 +1,139 @@ +&input_parameters + restart = F + ! Restart simulation that exited before simulation_time reached + proceed = F + ! Append simulation results after previous simulation_time reached + reinit = F + ! Begin a new simulation, reinitializing from restart file state + simulation_time = 5.0E-9 + ! Total aimed simulation time in seconds + ! Run 10 mu s If transients exist put 5 extra mu s. + snapshot_frequency = 5.0E-10 + ! Time between snapshots in seconds + restart_overwrite_frequency = 1.E-1 + ! Time between overwritting of restart file in seconds + dt = 1.E-2 + ! Time step as fraction of relativistic gyro-period + num_species = 1 + minimum_particle_energy = 1.0E2 + ! Minimum energy of simulated particles in eV + radiation = F + GC_rad_model='SDE' + collisions = F + collisions_model = 'SINGLE_SPECIES' + ! Options are: 'NONE','SINGLE_SPECIES' and 'MULTIPLE_SPECIES' + bound_electron_model = 'HESSLOW' + ! Options are: 'NO_BOUND', 'HESSLOW', and 'ROSENBLUTH' + field_model = 'EXTERNAL-AORSA' + profile_model = 'ANALYTIC' + ! The two options for this parameter are 'ANALYTICAL' or 'EXTERNAL'. + ! For 'ANALYTICAL', the magnetic field is calculated based on + ! the parameters given in the "analytic_mag_field_params" section. + ! For 'EXTERNAL', the magnetic field is loaded from the file + ! specified in "filename". + ! 'UNIFORM' A uniform magnetic field used to advance only electrons' + ! velocity. + magnetic_field_filename = 'AORSA_D3D_171089_200MHz_EFITgrid.h5' + time_slice = 0 + rmax = 2.1 + rmin = 0.9 + zmax = 0.6 + zmin = -0.6 + outputs_list = '{X,Y,V,B,E,g,eta,PSIp,flagCon,flagCol}' + ! List of outputs + !'{X,Y,V,E,B,g,mu,eta,Prad,Pin,flagCon,flagCol,gradB, + ! curlb,ne,Te,Zeff,PSIp,nimp}' + HDF5_error_handling = T + FO_GC_compare = F + orbit_model = 'FO' + ! 'FO' for full orbit, 'GCpre' for guiding center with pre-computed + ! auxiliary fields, 'GCgrad' for guiding center with auxiliary + ! fields computed with PSPLINE. + field_eval = 'interp' + ! Set for plasma_model='ANALYTICAL'. Can be 'interp' or 'eqn', + ! where 'eqn' evaluates particle fields at particle positions and + ! 'interp' interpolates precomputed fields. + FokPlan = F + SameRandSeed = T + SC_E = F + SC_E_add = F + pchunk = 1 +/ + +&plasma_species + runaway = T + ppp = 1000 + ! Number of particles per process (mpi) +! pinit = 1000 + q = -1.0 + ! Electric charge + m = 1.0 + ! In units of electron mass + spatial_distribution = 'MH_psi' + !spatial_distribution = 'TRACER' + ! Options are: 'UNIFORM', 'DISK', 'TORUS', 'EXPONENTIAL-TORUS', + ! 'GAUSSIAN-TORUS', 'ELLIPTIC-TORUS', 'EXPONENTIAL-ELLIPTIC-TORUS', + ! 'GAUSSIAN-ELLIPTICAL-TORUS', '2D-GAUSSIAN-ELLIPTIC-TORUS-MH', + ! 'AVALANCHE-4D','TRACER','SPONG-3D','HOLLMANN-3D' + Ro = 1.7 + PHIo = 0.0 + ! In degrees + Zo = 0.0 + r_inner = 0.0 + r_outter = 0.49 + shear_factor = 0.35 + sigmaR = 1.e6 + sigmaZ = 0.2 + theta_gauss = 0.0 + psi_max=1 + ! goes as R^2 for HOLLMANN-3D, is psiN_max for HOLLMANN-3D-PSI +! PSIp_min=0.01 + falloff_rate = 0.0 + energy_distribution = 'UNIFORM' + ! Options are: 'MONOENERGETIC', 'THERMAL', 'AVALANCHE', + ! 'EXPERIMENTAL', and 'MONOPITCH' + pitch_distribution = 'UNIFORM' + ! Options are: 'MONOPITCH', 'AVALANCHE', 'EXPERIMENTAL', and 'UNIFORM'. + Eno = 1.0E0 + ! Initial energy in eV + etao = 2.0 + ! Initial pitch angle + Eo_lims = 3.5E6,5.0E6 + ! Lower and upper limit of simulated energy range, in eV. + etao_lims = 0.0,20.0 + ! Lower and upper limit of simulated pitch-angle range, in degrees. + Xtrace = 1.7,0.0,0.0 + ! Initial position of tracer particle for debugging with + ! spatial_distribution='TRACER' + Spong_b = 0.2 + Spong_w = 0.1 + Spong_dlam = 0.1 + dth = 3. + ! Variance of sampling normal variate for pitch angle + dgam = 3. + ! Variance of sampling normal variate for pitch angle + dR = 0.1 + ! Variance of sampling normal variate for R location + dZ = 0.1 + ! Variance of sampling normal variate for Z location +/ + +&externalPlasmaModel + Bfield = F + B1field=T + E1field=T + dBfield = F + axisymmetric_fields = F + Bflux = T + Bflux3D = F + Efield = F + PSIp_lim=2.037510076503873e-01 + PSIp_0=-4.932491356369065e-02 + psip_conv=-1 + ! sign appended to \nabla\phi\times\nabla\psi_p for definition of + ! poloidal magnetic field. +1 for DIII-D, -1 for JET + AORSA_AMP_Scale=1. + AORSA_freq=200e6 + AORSA_nmode=35 + useLCFS = T +/ \ No newline at end of file diff --git a/test/aorsa/korc_aorsa.sh.in b/test/aorsa/korc_aorsa.sh.in new file mode 100755 index 00000000..8dd02641 --- /dev/null +++ b/test/aorsa/korc_aorsa.sh.in @@ -0,0 +1,23 @@ +#!/bin/bash + +set -ex + +#define input file +INPUT_FILE="${CMAKE_SOURCE_DIR}/test/aorsa/input_file_D3D_171089_AORSA.korc" +#define output directory +OUT_DIR="aorsa_test/rank_$1" + +#check that output directory doesn't exist so bash doesn't complain +if [ ! -d $OUT_DIR ]; then + mkdir -p $OUT_DIR +fi +if [ ! -f AORSA_D3D_171089_200MHz_EFITgrid.h5 ]; then + ln -s ${CMAKE_SOURCE_DIR}/test/aorsa/AORSA_D3D_171089_200MHz_EFITgrid.h5 AORSA_D3D_171089_200MHz_EFITgrid.h5 +fi + +${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} $1 ./xkorc $INPUT_FILE $OUT_DIR/ + +for i in $(seq 0 $(($1-1))); +do + h5diff -r -p 0.000008 $OUT_DIR/file_$i.h5 ${CMAKE_SOURCE_DIR}/test/aorsa/rank_$1/file_$i.h5 +done \ No newline at end of file diff --git a/test/aorsa/rank_1/file_0.h5 b/test/aorsa/rank_1/file_0.h5 new file mode 100644 index 00000000..cd0780fd Binary files /dev/null and b/test/aorsa/rank_1/file_0.h5 differ diff --git a/test/aorsa/rank_2/file_0.h5 b/test/aorsa/rank_2/file_0.h5 new file mode 100644 index 00000000..934ceb0f Binary files /dev/null and b/test/aorsa/rank_2/file_0.h5 differ diff --git a/test/aorsa/rank_2/file_1.h5 b/test/aorsa/rank_2/file_1.h5 new file mode 100644 index 00000000..d2b29b7d Binary files /dev/null and b/test/aorsa/rank_2/file_1.h5 differ diff --git a/test/aorsa/rank_4/file_0.h5 b/test/aorsa/rank_4/file_0.h5 new file mode 100644 index 00000000..28a20a06 Binary files /dev/null and b/test/aorsa/rank_4/file_0.h5 differ diff --git a/test/aorsa/rank_4/file_1.h5 b/test/aorsa/rank_4/file_1.h5 new file mode 100644 index 00000000..0429dd33 Binary files /dev/null and b/test/aorsa/rank_4/file_1.h5 differ diff --git a/test/aorsa/rank_4/file_2.h5 b/test/aorsa/rank_4/file_2.h5 new file mode 100644 index 00000000..8449191b Binary files /dev/null and b/test/aorsa/rank_4/file_2.h5 differ diff --git a/test/aorsa/rank_4/file_3.h5 b/test/aorsa/rank_4/file_3.h5 new file mode 100644 index 00000000..6fb73eb6 Binary files /dev/null and b/test/aorsa/rank_4/file_3.h5 differ diff --git a/test/aorsa/rank_8/file_0.h5 b/test/aorsa/rank_8/file_0.h5 new file mode 100644 index 00000000..12adebed Binary files /dev/null and b/test/aorsa/rank_8/file_0.h5 differ diff --git a/test/aorsa/rank_8/file_1.h5 b/test/aorsa/rank_8/file_1.h5 new file mode 100644 index 00000000..d2766320 Binary files /dev/null and b/test/aorsa/rank_8/file_1.h5 differ diff --git a/test/aorsa/rank_8/file_2.h5 b/test/aorsa/rank_8/file_2.h5 new file mode 100644 index 00000000..c6257246 Binary files /dev/null and b/test/aorsa/rank_8/file_2.h5 differ diff --git a/test/aorsa/rank_8/file_3.h5 b/test/aorsa/rank_8/file_3.h5 new file mode 100644 index 00000000..7e24e4cf Binary files /dev/null and b/test/aorsa/rank_8/file_3.h5 differ diff --git a/test/aorsa/rank_8/file_4.h5 b/test/aorsa/rank_8/file_4.h5 new file mode 100644 index 00000000..dfeffc75 Binary files /dev/null and b/test/aorsa/rank_8/file_4.h5 differ diff --git a/test/aorsa/rank_8/file_5.h5 b/test/aorsa/rank_8/file_5.h5 new file mode 100644 index 00000000..45e8cb1f Binary files /dev/null and b/test/aorsa/rank_8/file_5.h5 differ diff --git a/test/aorsa/rank_8/file_6.h5 b/test/aorsa/rank_8/file_6.h5 new file mode 100644 index 00000000..1c3fe516 Binary files /dev/null and b/test/aorsa/rank_8/file_6.h5 differ diff --git a/test/aorsa/rank_8/file_7.h5 b/test/aorsa/rank_8/file_7.h5 new file mode 100644 index 00000000..95934df1 Binary files /dev/null and b/test/aorsa/rank_8/file_7.h5 differ diff --git a/test/mars/file_0.h5 b/test/mars/file_0.h5 deleted file mode 100644 index 89e39280..00000000 Binary files a/test/mars/file_0.h5 and /dev/null differ diff --git a/test/mars/korc_mars.sh.in b/test/mars/korc_mars.sh.in index 5cf6a7ad..d30cdd19 100755 --- a/test/mars/korc_mars.sh.in +++ b/test/mars/korc_mars.sh.in @@ -17,4 +17,8 @@ fi ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} $1 ./xkorc $INPUT_FILE $OUT_DIR/ -h5diff -r -p 0.000008 $OUT_DIR/file_0.h5 ${CMAKE_SOURCE_DIR}/test/mars/file_0_$1.h5 +for i in $(seq 0 $(($1-1))); +do + h5diff -r -p 0.000008 $OUT_DIR/file_$i.h5 ${CMAKE_SOURCE_DIR}/test/mars/rank_$1/file_$i.h5 +done + diff --git a/test/mars/file_0_1.h5 b/test/mars/rank_1/file_0.h5 similarity index 99% rename from test/mars/file_0_1.h5 rename to test/mars/rank_1/file_0.h5 index 291d9cac..d483c9bf 100644 Binary files a/test/mars/file_0_1.h5 and b/test/mars/rank_1/file_0.h5 differ diff --git a/test/mars/rank_16/file_0.h5 b/test/mars/rank_16/file_0.h5 new file mode 100644 index 00000000..4d5019f3 Binary files /dev/null and b/test/mars/rank_16/file_0.h5 differ diff --git a/test/mars/rank_16/file_1.h5 b/test/mars/rank_16/file_1.h5 new file mode 100644 index 00000000..5acb1189 Binary files /dev/null and b/test/mars/rank_16/file_1.h5 differ diff --git a/test/mars/rank_16/file_10.h5 b/test/mars/rank_16/file_10.h5 new file mode 100644 index 00000000..8ba9819c Binary files /dev/null and b/test/mars/rank_16/file_10.h5 differ diff --git a/test/mars/rank_16/file_11.h5 b/test/mars/rank_16/file_11.h5 new file mode 100644 index 00000000..b8e34866 Binary files /dev/null and b/test/mars/rank_16/file_11.h5 differ diff --git a/test/mars/rank_16/file_12.h5 b/test/mars/rank_16/file_12.h5 new file mode 100644 index 00000000..98b98b1c Binary files /dev/null and b/test/mars/rank_16/file_12.h5 differ diff --git a/test/mars/rank_16/file_13.h5 b/test/mars/rank_16/file_13.h5 new file mode 100644 index 00000000..b14ce9ec Binary files /dev/null and b/test/mars/rank_16/file_13.h5 differ diff --git a/test/mars/rank_16/file_14.h5 b/test/mars/rank_16/file_14.h5 new file mode 100644 index 00000000..845b7aed Binary files /dev/null and b/test/mars/rank_16/file_14.h5 differ diff --git a/test/mars/rank_16/file_15.h5 b/test/mars/rank_16/file_15.h5 new file mode 100644 index 00000000..3f9bb28b Binary files /dev/null and b/test/mars/rank_16/file_15.h5 differ diff --git a/test/mars/rank_16/file_2.h5 b/test/mars/rank_16/file_2.h5 new file mode 100644 index 00000000..79e6f3f5 Binary files /dev/null and b/test/mars/rank_16/file_2.h5 differ diff --git a/test/mars/rank_16/file_3.h5 b/test/mars/rank_16/file_3.h5 new file mode 100644 index 00000000..4e2aff96 Binary files /dev/null and b/test/mars/rank_16/file_3.h5 differ diff --git a/test/mars/rank_16/file_4.h5 b/test/mars/rank_16/file_4.h5 new file mode 100644 index 00000000..8d796f3f Binary files /dev/null and b/test/mars/rank_16/file_4.h5 differ diff --git a/test/mars/rank_16/file_5.h5 b/test/mars/rank_16/file_5.h5 new file mode 100644 index 00000000..b109aabd Binary files /dev/null and b/test/mars/rank_16/file_5.h5 differ diff --git a/test/mars/rank_16/file_6.h5 b/test/mars/rank_16/file_6.h5 new file mode 100644 index 00000000..4493f672 Binary files /dev/null and b/test/mars/rank_16/file_6.h5 differ diff --git a/test/mars/rank_16/file_7.h5 b/test/mars/rank_16/file_7.h5 new file mode 100644 index 00000000..833fb228 Binary files /dev/null and b/test/mars/rank_16/file_7.h5 differ diff --git a/test/mars/rank_16/file_8.h5 b/test/mars/rank_16/file_8.h5 new file mode 100644 index 00000000..fee134e7 Binary files /dev/null and b/test/mars/rank_16/file_8.h5 differ diff --git a/test/mars/rank_16/file_9.h5 b/test/mars/rank_16/file_9.h5 new file mode 100644 index 00000000..493390c1 Binary files /dev/null and b/test/mars/rank_16/file_9.h5 differ diff --git a/test/mars/file_0_2.h5 b/test/mars/rank_2/file_0.h5 similarity index 99% rename from test/mars/file_0_2.h5 rename to test/mars/rank_2/file_0.h5 index 7a11abed..b51d46e7 100644 Binary files a/test/mars/file_0_2.h5 and b/test/mars/rank_2/file_0.h5 differ diff --git a/test/mars/rank_2/file_1.h5 b/test/mars/rank_2/file_1.h5 new file mode 100644 index 00000000..e499b66b Binary files /dev/null and b/test/mars/rank_2/file_1.h5 differ diff --git a/test/mars/file_0_4.h5 b/test/mars/rank_4/file_0.h5 similarity index 99% rename from test/mars/file_0_4.h5 rename to test/mars/rank_4/file_0.h5 index 659136ba..ca34f8d4 100644 Binary files a/test/mars/file_0_4.h5 and b/test/mars/rank_4/file_0.h5 differ diff --git a/test/mars/rank_4/file_1.h5 b/test/mars/rank_4/file_1.h5 new file mode 100644 index 00000000..fd9d94d3 Binary files /dev/null and b/test/mars/rank_4/file_1.h5 differ diff --git a/test/mars/rank_4/file_2.h5 b/test/mars/rank_4/file_2.h5 new file mode 100644 index 00000000..12a13691 Binary files /dev/null and b/test/mars/rank_4/file_2.h5 differ diff --git a/test/mars/rank_4/file_3.h5 b/test/mars/rank_4/file_3.h5 new file mode 100644 index 00000000..69363333 Binary files /dev/null and b/test/mars/rank_4/file_3.h5 differ diff --git a/test/mars/file_0_8.h5 b/test/mars/rank_8/file_0.h5 similarity index 99% rename from test/mars/file_0_8.h5 rename to test/mars/rank_8/file_0.h5 index 75a86a14..7d44f54a 100644 Binary files a/test/mars/file_0_8.h5 and b/test/mars/rank_8/file_0.h5 differ diff --git a/test/mars/rank_8/file_1.h5 b/test/mars/rank_8/file_1.h5 new file mode 100644 index 00000000..89966f2e Binary files /dev/null and b/test/mars/rank_8/file_1.h5 differ diff --git a/test/mars/rank_8/file_2.h5 b/test/mars/rank_8/file_2.h5 new file mode 100644 index 00000000..99bacc7c Binary files /dev/null and b/test/mars/rank_8/file_2.h5 differ diff --git a/test/mars/rank_8/file_3.h5 b/test/mars/rank_8/file_3.h5 new file mode 100644 index 00000000..60ade3af Binary files /dev/null and b/test/mars/rank_8/file_3.h5 differ diff --git a/test/mars/rank_8/file_4.h5 b/test/mars/rank_8/file_4.h5 new file mode 100644 index 00000000..9ac75902 Binary files /dev/null and b/test/mars/rank_8/file_4.h5 differ diff --git a/test/mars/rank_8/file_5.h5 b/test/mars/rank_8/file_5.h5 new file mode 100644 index 00000000..9a9178a2 Binary files /dev/null and b/test/mars/rank_8/file_5.h5 differ diff --git a/test/mars/rank_8/file_6.h5 b/test/mars/rank_8/file_6.h5 new file mode 100644 index 00000000..bf67d7ec Binary files /dev/null and b/test/mars/rank_8/file_6.h5 differ diff --git a/test/mars/rank_8/file_7.h5 b/test/mars/rank_8/file_7.h5 new file mode 100644 index 00000000..b10522a9 Binary files /dev/null and b/test/mars/rank_8/file_7.h5 differ