diff --git a/matlab/@espic2dhdf5/display2Dpotentialwell.m b/matlab/@espic2dhdf5/display2Dpotentialwell.m index 1053f62..4e5399d 100644 --- a/matlab/@espic2dhdf5/display2Dpotentialwell.m +++ b/matlab/@espic2dhdf5/display2Dpotentialwell.m @@ -1,176 +1,197 @@ function f=display2Dpotentialwell(obj,timestep,rcoord,clims,rescale, mag,nblv) % Display the 2D potential well at time obj.t2d(timestep) % if rcoord is true, the potential is evaluated at grid points in r,z coordinates % if false, the potential is evaluated at grid points in magnetic field line coordinates % clims are the values limits in eV for the color coding of the % potential well if iscell(timestep) timestep=cell2mat(timestep); end if nargin <3 rcoord=true; end if nargin<4 clims=[-inf inf]; end if nargin<5 rescale=false; end if nargin<6 mag=[]; end if nargin<7 - nblv=400; + nblv=500; end f=figure('Name',sprintf('%s Potential well',obj.name)); ax1=gca; model=obj.potentialwellmodel(timestep,false,mag,nblv); z=model.z; r=model.r; Pot=model.pot; rathet=model.rathet; N0=obj.N(:,:,1); id=find(timestep==0); timestep(id)=[]; Nend=obj.N(:,:,timestep); if(~isempty(id)) N0=zeros(obj.N.nr,obj.N.nz); Nend=cat(3,Nend(:,:,1:id-1),N0,Nend(:,:,id:end)); end Nend=mean(Nend,3); geomw=obj.geomweight(:,:,1); %z(isnan(Pot))=[]; %r(isnan(Pot))=[]; %Pot(isnan(Pot))=[]; if rescale if obj.spl_bound.nbsplines >0 Pot=Pot/(obj.spl_bound.boundary(2).Dval-obj.spl_bound.boundary(1).Dval); else Pot=Pot/(obj.potout-obj.potinn)/obj.phinorm; end end if rcoord - [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); - Pot=griddata(z,r,Pot,Zmesh,Rmesh,'natural'); + [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); %Pot(obj.geomweight(:,:,1)<=0)=NaN; %Pot=imgaussfilt(Pot,'FilterSize',11); - Pot(Pot<=0)=NaN; - Pot(obj.geomweight(:,:,1)<0)=NaN; - contourf(obj.zgrid(1:end)*1e3,obj.rgrid(1:end)*1e3,Pot(1:end,1:end),50,'edgecolor','none','Displayname','Well') + Pot(Pot<0)=NaN; + geow=griddedInterpolant(Zmesh',Rmesh',obj.geomweight(:,:,1)'); + [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); + finezgrid=obj.zgrid; + finergrid=obj.rgrid; +% finergrid=[obj.rgrid';[0.5*(obj.rgrid(1:end-1)+obj.rgrid(2:end)); 0]']; +% finergrid=finergrid(1:end-1)'; +% finezgrid=[obj.zgrid';[0.5*(obj.zgrid(1:end-1)+obj.zgrid(2:end)); 0]']; +% finezgrid=finezgrid(1:end-1)'; +% [Zmesh,Rmesh]=meshgrid(finezgrid,finergrid); + Pot=griddata(z,r,Pot,Zmesh,Rmesh,'natural'); + boundaries=geow(Zmesh',Rmesh')'; + Pot(boundaries<0)=NaN; + contourf(finezgrid*1e3,finergrid*1e3,Pot(1:end,1:end),50,'edgecolor','none','Displayname','Well') xlabel('z [mm]') ylabel('r [mm]') xlim([obj.zgrid(1) obj.zgrid(end)]*1e3) ylim([obj.rgrid(1) obj.rgrid(end)]*1e3) hold(gca, 'on') rdisp=obj.rgrid; %% Magnetic field lines Blines=obj.rAthet; levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),20); Blines(obj.geomweight(:,:,1)<0)=NaN; [~,h1]=contour(obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'k-.','linewidth',1.2,'Displayname','Magnetic field lines'); - + geomw=obj.geomweight(:,:,1); else - rdisp=sort(obj.rAthet(:,end)); + lvls=linspace(min(obj.rAthet(:)),max(obj.rAthet(:)),400); + rdisp=lvls; [Zmesh,Rmesh]=meshgrid(obj.zgrid,rdisp); - Pot=griddata(z,rathet,Pot,Zmesh,Rmesh); + finergrid=[rdisp';[0.5*(drisp(1:end-1)+rdisp(2:end)); 0]']; + finergrid=finergrid(1:end-1)'; + finezgrid=[obj.zgrid';[0.5*(obj.zgrid(1:end-1)+obj.zgrid(2:end)); 0]']; + finezgrid=finezgrid(1:end-1)'; + [Zmesh,Rmesh]=meshgrid(finezgrid,finergrid); + Pot=griddata(z,rathet,Pot,Zmesh,Rmesh,'natural'); + boundaries=geow(Zmesh',Rmesh')'; + Pot(boundaries<0)=NaN; + + %Pot=griddata(z,rathet,Pot,Zmesh,Rmesh); [Zinit,~]=meshgrid(obj.zgrid,obj.rAthet(:,1)); % if isempty(obj.maxwellsrce) % end N0=griddata(Zinit,obj.rAthet,N0,Zmesh,Rmesh); Nend=griddata(Zinit,obj.rAthet,Nend,Zmesh,Rmesh); geomw=griddata(Zinit,obj.rAthet,geomw,Zmesh,Rmesh); - Pot(geomw<=0)=0; - surface(obj.zgrid(1:end),rdisp,Pot(1:end,1:end),'edgecolor','none') + Pot(geomw<0)=NaN; + %geomw(geomw2>0)=-1; + contourf(finezgrid,finergrid*1000,Pot(1:end,1:end),'edgecolor','none') ylabel('rA_\theta [Tm^2]') xlabel('z [m]') - xlim([obj.zgrid(1) obj.zgrid(end)]) - ylim([min(rdisp) max(rdisp)]) + xlim([obj.zgrid(1) obj.zgrid(end)]*1e3) + ylim([min(rdisp) max(rdisp)]*1e3) hold(gca, 'on') end if (timestep==0) title(sprintf('Potential well Vacuum')) else title(sprintf('Potential well t=%1.2f [ns]',obj.t2d(timestep)*1e9)) end maxdensend=max(Nend(:)); contourscale=0.1; Nend=(Nend-contourscale*maxdensend)/maxdensend; maxdens0=max(N0(:)); contourscale=0.1; N0=(N0-contourscale*maxdens0)/maxdens0; contour(obj.zgrid*1e3,rdisp*1e3,Nend,linspace(0,1-contourscale,5),'k--','linewidth',1.5,'Displayname','Cloud Boundaries'); contour(obj.zgrid*1e3,rdisp*1e3,N0,[0 0],'k-.','linewidth',1.5,'Displayname','Source boundaries'); %contour(obj.zgrid*1e3,rdisp*1e3,geomw,[0 0],'-','linecolor',[0.5 0.5 0.5],'linewidth',1.5,'Displayname','Vessel Boundaries'); c=colorbar; colormap('jet'); % Grey outline showing the metalic walls - geomw(obj.geomweight(:,:,1)>0)=-1; - geomw(obj.geomweight(:,:,1)<=0)=1; - [c1,hContour]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,geomw, [0 0]); + + [c1,hContour]=contourf(ax1,obj.zgrid*1000,rdisp*1000,-geomw, [0 0]); drawnow; xlim(ax1,[obj.zgrid(1)*1000 obj.zgrid(end)*1000]) if(obj.conformgeom) - ylim([ax1 ],[obj.rgrid(1)*1000 obj.rgrid(rgridend)*1000]) + %ylim([ax1 ],[obj.rgrid(1)*1000 obj.rgrid(rgridend)*1000]) else - ylim([ax1],[obj.rgrid(1)*1000 obj.rgrid(end)*1000]) + % ylim([ax1],[obj.rgrid(1)*1000 obj.rgrid(end)*1000]) end %ylim(ax1,[0.05*1000 obj.rgrid(end)*1000]) %xlim([obj.zgrid(1) 0.185]*1e3) xlabel(ax1,'z [mm]') ylabel(ax1,'r [mm]') view(ax1,2) if rescale c.Label.String= 'Normalized well depth (U_{well}/\Delta\phi) [eV/V]'; else c.Label.String= 'well depth [eV]'; end f.PaperUnits='centimeters'; caxis(clims) grid(ax1, 'on'); hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' try drawnow + hFills(1).ColorData = uint8([150;150;150;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end catch end % add central and external metallic walls if we have a coaxial % configuration if( obj.walltype >=2 && obj.walltype<=4) rectangle('Position',[obj.zgrid(1) obj.r_b obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none') ylimits=ylim; ylim([ylimits(1),ylimits(2)+1]) end if sum(obj.geomweight(:,1,1))==0 rectangle('Position',[obj.zgrid(1) obj.r_a-0.001 obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none') ylimits=ylim; ylim([ylimits(1)-1,ylimits(2)]) end %axis equal %xlim([-100 200]) [max_depth,id]=max(abs(Pot(:))); [idr,idz]=ind2sub(size(Pot),id); fprintf('Maximum potential wel depth: %f eV\n',max_depth) - fprintf('at location r=%f z=%f [mm]\n',obj.rgrid(idr)*1e3, obj.zgrid(idz)*1e3) + fprintf('at location r=%f z=%f [mm]\n',finergrid(idr)*1e3, finezgrid(idz)*1e3) papsize=[14 8]; if rcoord obj.savegraph(f,sprintf('%s/%s_wellr%i',obj.folder,obj.name,floor(mean(timestep))),papsize); else obj.savegraph(f,sprintf('%s/%s_wellpsi%i',obj.folder,obj.name,floor(mean(timestep))),papsize); end end \ No newline at end of file diff --git a/matlab/@espic2dhdf5/displaybiasevol.m b/matlab/@espic2dhdf5/displaybiasevol.m index 745b89d..4048ba1 100644 --- a/matlab/@espic2dhdf5/displaybiasevol.m +++ b/matlab/@espic2dhdf5/displaybiasevol.m @@ -1,20 +1,72 @@ -function displaybiasevol(obj) +function displaybiasevol(obj,scalet) %% Plot the time evolution of the non-ideal power supply bias + if nargin<2 + scalet=true; + end tmin=1; tmax=length(obj.ekin); f=figure('Name', sprintf('%s Non ideal bias',obj.name)); + subplot(2,1,1) time=obj.t0d(tmin:tmax); if ~isfield(obj.psupply,'biases') - plot(time*1e9,(obj.potout-obj.potinn)*obj.phinorm*ones(size(time))) + time=time*1e9; + biases=(obj.potout-obj.potinn)*obj.phinorm*ones(size(time)); else - plot(time*(obj.neutcol.neutdens/obj.psupply.expdens),obj.psupply.biases(:,tmin:tmax),'linewidth',2) + if scalet && obj.neutcol.present + time=time*(obj.neutcol.neutdens/obj.psupply.expdens); + else + time=time*1e9; + end + biases=obj.psupply.biases(:,tmin:tmax); + end + for i=1:size(biases,1)-1 + if obj.spl_bound.nbsplines>0 && isfield(obj.spl_bound.boundary(i),'type') && obj.spl_bound.boundary(i).type~=0 + continue + end + p(i)=plot(time,biases(i,:),'linewidth',2); + hold on end set(gca,'fontsize',12); %obj.t0d(tmin:tmax),obj.ekin(tmin:tmax)-obj.epot(tmin:tmax),'--', xlabel('Time [s]') ylabel('Bias [V]') grid on + + subplot(2,1,2) + + geomw=obj.geomweight(:,:,1); + geomw(geomw<0)=0; + geomw(geomw>0)=NaN; + [c1,hContour]=contourf(obj.zgrid*1000,obj.rgrid*1000,geomw, [0 0]); + hold on + drawnow; + grid on; + if obj.spl_bound.nbsplines>0 + for i=1:length(obj.spl_bound.boundary)-1 + if isfield(obj.spl_bound.boundary(i),'type') && obj.spl_bound.boundary(i).type~=0 + continue + end + plot(obj.spl_bound.boundary(i).coefs(:,1)*1000,obj.spl_bound.boundary(i).coefs(:,2)*1000,'linestyle',p(i).LineStyle,... + 'color',p(i).Color,'marker',p(i).Marker,... + 'displayname',sprintf('border %i',i),'linewidth',1.8) + hold on + end + end + title('Domain') + xlabel('z [mm]') + ylabel('r [mm]') + grid on + set(gca,'fontsize',12) + hFills=hContour.FacePrims; + [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' + try + hFills(1).ColorData = uint8([150;150;150;255]); + for idx = 2 : numel(hFills) + hFills(idx).ColorData(4) = 0; % default=255 + end + catch + end obj.savegraph(f,sprintf('%s/%sAppliedBias',obj.folder,obj.name)); end \ No newline at end of file diff --git a/matlab/@espic2dhdf5/displayconfiguration.m b/matlab/@espic2dhdf5/displayconfiguration.m index fdaa083..eb71e31 100644 --- a/matlab/@espic2dhdf5/displayconfiguration.m +++ b/matlab/@espic2dhdf5/displayconfiguration.m @@ -1,113 +1,113 @@ function f=displayconfiguration(obj,fieldsteps,xlimits,ylimits) %displayconfiguration plot the configuration of the simulation % domain withe boundaries, the magnetic field lines the % electric equipotential lines and the electron density % averaged in time between t2d(fieldsteps(1)) and t2d(fieldsteps(end)) fieldstart=fieldsteps(1); fieldend=fieldsteps(end); if nargin<3 xlimits=inf*[-1 1]; end if nargin<4 ylimits=inf*[-1 1]; end dens=mean(obj.N(:,:,fieldstart:fieldend),3); geomw=obj.geomweight(:,:,1); maxdens=max(dens(:)); geomw(geomw<0)=0; geomw(geomw>0)=maxdens; dens(geomw<=0)=0; geomw(geomw>0)=NaN; f=figure('Name', sprintf('%s fields',obj.name)); ax1=gca; title(sprintf('Configuration')) %dens(dens<=1e13)=NaN; %% electron density h=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,dens,50,'Displayname','n_e [m^{-3}]', 'linestyle','none'); hold on; colormap(flipud(hot)); %% Magnetic field lines Blines=obj.rAthet; levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),50); [~,h1]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'-.','color','k','linewidth',1.5,'Displayname','Magnetic field lines'); %% Equipotential lines Pot=mean(obj.pot(:,:,fieldstart:fieldend),3); Pot(obj.geomweight(:,:,1)<0)=NaN; %levels=8;%[-3.4 -5 -10 -15 -20 -25];%7; potcolor='b';%[0.3660 0.6740 0.1880]; - [c1,h2]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Pot/1e3,50,'--','color',potcolor,'ShowText','on','linewidth',1.2,'Displayname','Equipotentials [kV]'); + [c1,h2]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Pot/1e3,'--','color',potcolor,'ShowText','on','linewidth',1.2,'Displayname','Equipotentials [kV]'); clabel(c1,h2,'Color',potcolor) clabel(c1,h2, 'labelspacing', 200); % Grey outline shows metallic parts [c1,hContour]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,geomw, [0 0]); drawnow; % set the axia limits xlim(ax1,[obj.zgrid(1)*1000 obj.zgrid(end)*1000]) if(obj.conformgeom) ylim([ax1 ],[obj.rgrid(1)*1000 obj.rgrid(rgridend)*1000]) else ylim([ax1],[obj.rgrid(1)*1000 obj.rgrid(end)*1000]) end legend([h1,h2],{'Magnetic field lines','Equipotentials [kV]'},'location','northeast') xlabel(ax1,'z [mm]') ylabel(ax1,'r [mm]') c = colorbar(ax1); c.Label.String= 'n[m^{-3}]'; view(ax1,2) grid on; hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' try hFills(1).ColorData = uint8([150;150;150;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end catch end [~, name, ~] = fileparts(obj.file); % with this you could show the outline of the maxwellian source % if obj.maxwellsrce.present % rlen=diff(obj.maxwellsrce.rlim); % zlen=diff(obj.maxwellsrce.zlim); % rectangle('Position',[obj.maxwellsrce.zlim(1) obj.maxwellsrce.rlim(1) zlen rlen]*1000,'Edgecolor','g','Linewidth',2,'Linestyle','--') % end % in case of coaxial configuration, extend the display domain % and add grey rectangles to show metallic boundaries if( obj.walltype >=2 && obj.walltype<=4) rectangle('Position',[obj.zgrid(1) obj.r_b obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none') ylimits=ylim; ylim([ylimits(1),ylimits(2)+1]) end if sum(obj.geomweight(:,1,1))==0 rectangle('Position',[obj.zgrid(1) obj.r_a-0.001 obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none') ylimits=ylim; ylim([ylimits(1)-1,ylimits(2)]) end if nargin>2 xlim(xlimits) end if nargin>3 ylim(ylimits) end f.PaperUnits='centimeters'; %axis equal papsize=[14 5 ]; obj.savegraph(f,sprintf('%s/%sFields',obj.folder,obj.name),papsize); end diff --git a/matlab/@espic2dhdf5/displaytotcurrevol_geom.m b/matlab/@espic2dhdf5/displaytotcurrevol_geom.m index 28a0573..a55f4b1 100644 --- a/matlab/@espic2dhdf5/displaytotcurrevol_geom.m +++ b/matlab/@espic2dhdf5/displaytotcurrevol_geom.m @@ -1,186 +1,193 @@ function displaytotcurrevol_geom(obj,timesteps,toptitle,scalet,dens,subdiv,nmean) % Computes and display the time evolution of the outgoing % currents at time obj.t2d(timesteps) %scalet=true scales the time by the ellastic collision %frequency %dens = true plot the time evolution of the maximum electron %density in the simulation domain otherwise plot the total %number of electrons in the domain % also plot in a subplot the color coded boundary corresponding % to each current if nargin<2 timesteps=1:length(obj.t2d); end if nargin<3 scalet=true; end if nargin <4 dens=true; end if nargin<5 subdiv=1; end if nargin<6 nmean=1; end if scalet if obj.neutcol.present vexb0=(obj.Ezxt(:,:,1).*obj.Br'-obj.Erxt(:,:,1).*obj.Bz')./(obj.B'.^2); vexb0(obj.geomweight(:,:,1)<=0)=0; potwell=obj.PotentialWell(0); vexb0(isnan(potwell))=NaN; vexb0=mean(abs(vexb0(:)),'omitnan'); E=0.5*obj.msim/obj.weight*vexb0^2/obj.qe; taucol=1/(obj.neutcol.neutdens*vexb0*(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E))); try Sio_S=1e17*(obj.neutcol.neutdens*vexb0*obj.sigio(E))/(obj.maxwellsrce.frequency*obj.weight/(pi*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)*diff(obj.maxwellsrce.zlim))) catch end tlabel='t/\tau_d [-]'; else taucol=2*pi/obj.omece; tlabel='t/\tau_ce [-]'; end else taucol=1e-9; tlabel='t [ns]'; end if dens N=obj.N(:,:,timesteps); geomw=obj.geomweight(:,:,1); geomw(geomw<0)=0; geomw(geomw>0)=1; N=N.*geomw; %[~,idl]=max(N(:,:,end),[],'all','linear'); %[ir,iz]=ind2sub(size(geomw),idl); %nmax=squeeze(max(max(N,[],1),[],2)); tn=(obj.t2d(timesteps)); nrhalf=find(obj.rgrid>0.07,1,'first'); if(nrhalf=obj.t2d(timesteps(1)),1,'first'); t0dlst=find(obj.t0d<=obj.t2d(timesteps(end)),1,'last'); tn=obj.t0d(t0dst:t0dlst); nmax=obj.npart(t0dst:t0dlst)*obj.weight; nlabel='Nb e^-'; ndlabel='Nb e^-'; end [currents,pos]=obj.OutCurrents(timesteps,subdiv); - P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar + P=1; + if obj.neutcol.present + P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar + end currents=currents/P; f=figure('Name',sprintf('%s Charges',obj.name)); tiledlayout(2,1) nexttile % Plot the evolution of nb of particles yyaxis right if size(nmax,2)>1 nmax=nmax'; end for i=1:size(nmax,2) p=plot(tn/taucol,nmax(:,i),'b','linewidth',1.8,'Displayname',sprintf('%s, %d',ndlabel,i)); hold on end ylabel(nlabel) axl=gca; axl.YAxis(2).Color=p.Color; ylim([0 inf]) if(obj.B(1,1)>obj.B(end,1)) lname='HFS'; rname='LFS'; else lname='LFS'; rname='HFS'; end yyaxis( 'left'); map=colormap(lines); set(axl,'linestyleorder',{'-',':','--','*','+'},... 'ColorOrder',map(2:7,:), 'NextPlot','replacechildren') p(1)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(1,:),nmean),'Displayname',lname,'linewidth',1.8); hold on p(2)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(2,:),nmean),'Displayname',rname,'linewidth',1.8); % Plot the currents for i=3:size(currents,1) p(i)=plot(axl,obj.t2d(timesteps)/taucol,movmean(currents(i,:),nmean),'Displayname',sprintf('border %i',i-2),'linewidth',1.8); end plot(axl,obj.t2d(timesteps)/taucol,movmean(sum(currents(:,:),1,'omitnan'),nmean),'k-','Displayname','total','linewidth',1.8); xlabel(tlabel) - ylabel('I/p_n [A/mbar]') + if obj.neutcol.present + ylabel('I/p_n [A/mbar]') + else + ylabel('I [A]') + end grid on set(gca,'fontsize',12) ax.YAxis(1).Color='black'; %legend('Orientation','horizontal','location','north','numcolumns',3) if ~isempty(toptitle) title(toptitle) end ax2=nexttile; geomw=obj.geomweight(:,:,1); geomw(geomw<0)=0; geomw(geomw>0)=NaN; [c1,hContour]=contourf(ax2,obj.zgrid*1000,obj.rgrid*1000,geomw, [0 0]); hold on drawnow; grid on; for i=1:length(pos) plot(ax2,pos{i}(1,:)*1000,pos{i}(2,:)*1000,'linestyle',p(i+2).LineStyle,... 'color',p(i+2).Color,'marker',p(i+2).Marker,... 'displayname',sprintf('border %i',i),'linewidth',1.8) hold on end title('Domain') plot(ax2,ones(size(obj.rgrid))*obj.zgrid(1)*1000,obj.rgrid*1000,'linestyle',p(1).LineStyle,... 'color',p(1).Color,'marker',p(1).Marker,... 'displayname',lname,'linewidth',1.8) plot(ax2,ones(size(obj.rgrid))*obj.zgrid(end)*1000,obj.rgrid*1000,'linestyle',p(2).LineStyle,... 'color',p(2).Color,'marker',p(2).Marker,... 'displayname',rname,'linewidth',1.8) xlabel('z [mm]') ylabel('r [mm]') grid on set(gca,'fontsize',12) hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' try hFills(1).ColorData = uint8([150;150;150;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end catch end %legend('Orientation','horizontal','location','north','numcolumns',4) fprintf('mean total current: %f [A/mbar]\n',mean(sum(currents(:,max(1,size(currents,2)-30):end),1,'omitnan'))); %if nargin <3 % sgtitle(sprintf('\\phi_b-\\phi_a=%.2g kV, B=%f T',(obj.potout-obj.potinn)*obj.phinorm/1e3,mean(obj.B(:)))) %elseif ~isempty(toptitle) % sgtitle(toptitle) %end if length(subdiv)>1 obj.savegraph(f,sprintf('%s/%s_totIEvol%i%i_div%i',obj.folder,obj.name,scalet,dens,nmean),[16 14]); else obj.savegraph(f,sprintf('%s/%s_totIEvol%i%i_%i',obj.folder,obj.name,scalet,dens,nmean),[16 14]); end end \ No newline at end of file diff --git a/matlab/@espic2dhdf5/espic2dhdf5.m b/matlab/@espic2dhdf5/espic2dhdf5.m index 793e40a..d65e177 100644 --- a/matlab/@espic2dhdf5/espic2dhdf5.m +++ b/matlab/@espic2dhdf5/espic2dhdf5.m @@ -1,2994 +1,3008 @@ classdef espic2dhdf5 %espic2dhdf5 General class used to treat hdf5 result files of espic2d code % A result file is loaded with a call to M=espic2dhdf5(filename) where filename is the relative or absolute file path % after loading, several quantities and composite diagnostics such as moments of the distribution function or individual particles % quantities can be accessed. properties filename name folder fullpath timestamp info t0d t1d t2d tpart it0 it1 it2 restartsteps restarttimes %% Physical constants vlight=299792458; qe=1.60217662E-19; me=9.109383E-31; eps_0=8.85418781762E-12; kb=1.38064852E-23; %% Run parameters dt % simulation time step nrun % number of time steps simulated nlres nlsave nlclassical % Was the equation of motion solved in the classical framework nlPhis % Was the self-consistent electric field computed nz % number of intervals in the z direction for the grid nnr % number of intervals in the r direction for the grid for each of the 3 mesh regions lz % physical axial dimension of the simulation space nplasma % Number of initial macro particles potinn % Normalized electric potential at the coaxial insert potout % Normalized electric potential at the cylinder surface B0 % Normalization for the magnetic field Rcurv % Magnetic mirror ratio width % Magnetic mirror length n0 % Initial particle density in case of old particle loading temp % Initial particle temperature in case of old particle loading femorder % finite element method order in z and r direction ngauss % Order of the Gauss integration method for the FEM plasmadim % initial dimensions of the plasma for the old particle loading system radii % Radial limits of the three mesh regions coarse,fine,coarse H0 % Initial particle Energy for Davidsons distribution function P0 % Initial particle Angular momentum for Davidsons distribution function normalized % Are the parts quantities normalized in the h5 file nbspecies % Number of species simulated %% Frequencies omepe % Reference plasma frequency used for normalization omece % Reference cyclotronic frequency for normalization %% Normalizations tnorm % Time normalization rnorm % Dimension normalization bnorm % Magnetic field normalization enorm % Electric field normalization phinorm % Electric potential normalization vnorm % Velocity normalization %% Grid data rgrid % Radial grid position points zgrid % Axial grid position points dz % Axial grid step dr % Radial grid step for the three mesh regions CellVol % Volume of the cell used for density calculation celltype % type of cell -1 outside 1 inside 0 border linked_s % location of linked spline bsplinetype %% Magnetic field Br % Radial magnetic field Bz % Axial magnetic field Athet % Azimuthal component of the Magnetic potential vector rAthet % r*Athet used for the representation of magnetic field lines B % Magnetic field amplitude sinthet % ratio to project quantities along the magnetic field lines costhet % ratio to project quantities along the magnetic field lines %% Energies epot % Time evolution of the particles potential energy ekin % Time evolution of the particles kinetic energy etot % Time evolution of the particles total energy etot0 % Time evolution of the reference particle total energy eerr % Time evolution of the error on the energy conservation npart % Time evolution of the number of simulated %% 2D time data evaluated on grid points N % main specie Density fluidUR % main specie radial fluid velocity fluidUZ % main specie axial fluid velocity fluidUTHET % main specie azimuthal fluid velocity pot % Electric potential evaluated at grid points potxt % External Electric potential evaluated at grid points phi % Electric potential in spline form Er % Radial electric field Ez % Axial electric field Erxt % External Radial electric field Ezxt % External Axial electric field Presstens % Pressure tensor fluidEkin % average kinetic energy in each direction %% Splines knotsr % Spline radial knots knotsz % Spline axial knots %% Particle parameters weight % Macro particle numerical weight of the main specie qsim % Macro particle charge msim % Macro particle mass nbparts % Time evolution of the number of simulated particles partepot % Electric potential at the particles positions R % Particles radial position Z % Particles axial position Rindex % Particles radial grid index Zindex % Particles axial grid index partindex % Particles unique id for tracing trajectories VR % Particles radial velocity VZ % Particles axial velocity VTHET % Particles azimuthal velocity THET % Particles azimuthal position species % Array containing the other simulated species %% Celldiag celldiag % Array containing the cell diagnostic data nbcelldiag % Total number of cell diagnostics %% Curvilinear geometry conformgeom % stores if we use the conforming or nonconforming boundary conditions r_a r_b z_r z_0 r_0 r_r L_r L_z Interior above1 above2 interior walltype geomweight dirichletweight gtilde spl_bound %% Maxwell source parameters maxwellsrce %% Collision with neutral parameters neutcol nudcol % effective momentum collision frequency %% Non ideal power supply psupply end methods function file=file(obj) % returns the h5 file name file=obj.filename; end function obj = espic2dhdf5(filename,readparts,old) % Reads the new result file filename and read the parts data if readparts==true % adds the helper_classes folder to the path matlabfuncpath = dir([mfilename('fullpath'),'.m']); addpath(sprintf('%s/../helper_classes',matlabfuncpath.folder)); addpath(sprintf('%s/../extrema',matlabfuncpath.folder)); addpath(sprintf('%s/../export_fig',matlabfuncpath.folder)); rehash path % Try catch are there for compatibility with older simulation files filedata=dir(filename); if (isempty(filedata)) error("File: ""%s"" doesn't exist",filename) end obj.folder=filedata.folder; obj.filename=filename; [~, obj.name, ext] = fileparts(obj.filename); obj.filename=[obj.name,ext]; obj.fullpath=[obj.folder,'/',obj.filename]; obj.timestamp=filedata.date; if nargin==1 readparts=true; end if nargin<3 old=false; end %obj.info=h5info(filename); %% Read the run parameters obj.dt = h5readatt(obj.fullpath,'/data/input.00/','dt'); obj.nrun = h5readatt(obj.fullpath,'/data/input.00/','nrun'); obj.nlres = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlres'),'y'); obj.nlsave = strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlsave'),'y'); obj.nlclassical =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlclassical'),'y'); obj.nlPhis =strcmp(h5readatt(obj.fullpath,'/data/input.00/','nlPhis'),'y'); obj.nz = h5readatt(obj.fullpath,'/data/input.00/','nz'); obj.nnr = h5read(obj.fullpath,'/data/input.00/nnr'); obj.lz = h5read(obj.fullpath,'/data/input.00/lz'); obj.qsim = h5readatt(obj.fullpath,'/data/input.00/','qsim'); obj.msim = h5readatt(obj.fullpath,'/data/input.00/','msim'); try obj.r_a=h5readatt(obj.fullpath,'/data/input.00/geometry','r_a'); obj.r_b=h5readatt(obj.fullpath,'/data/input.00/geometry','r_b'); obj.z_r=h5readatt(obj.fullpath,'/data/input.00/geometry','z_r'); obj.r_r=h5readatt(obj.fullpath,'/data/input.00/geometry','r_r'); obj.r_0=h5readatt(obj.fullpath,'/data/input.00/geometry','r_0'); obj.z_0=h5readatt(obj.fullpath,'/data/input.00/geometry','z_0'); obj.above1=h5readatt(obj.fullpath,'/data/input.00/geometry','above1'); obj.above2=h5readatt(obj.fullpath,'/data/input.00/geometry','above2'); obj.interior=h5readatt(obj.fullpath,'/data/input.00/geometry','interior'); obj.walltype=h5readatt(obj.fullpath,'/data/input.00/geometry','walltype'); try obj.L_r=h5readatt(obj.fullpath,'/data/input.00/geometry','L_r'); obj.L_z=h5readatt(obj.fullpath,'/data/input.00/geometry','L_z'); catch end obj.conformgeom=false; catch obj.conformgeom=true; obj.walltype=0; obj.r_a=obj.rgrid(1); obj.r_b=obj.rgrid(end); obj.above1=1; obj.above2=-1; obj.L_r=0; obj.L_z=0; end try obj.weight=h5readatt(obj.fullpath,'/data/part/','weight'); catch obj.weight=obj.msim/obj.me; end filesgrpinfo=h5info(obj.fullpath,'/files'); nbrst=h5readatt(obj.fullpath,'/files','jobnum'); obj.restartsteps(1)=0; obj.restarttimes(1)=0; grp=sprintf('/data/input.%02i/',0); obj.dt(1)=h5readatt(obj.fullpath,grp,'dt'); for i=1:nbrst grp=sprintf('/data/input.%02i/',i); obj.restartsteps(i+1)= h5readatt(obj.fullpath,grp,'startstep'); obj.restarttimes(i+1)= obj.restarttimes(i) + (obj.restartsteps(i+1)-obj.restartsteps(i))*obj.dt(i); obj.dt(i+1)=h5readatt(obj.fullpath,grp,'dt'); end obj.nplasma = h5readatt(obj.fullpath,'/data/input.00/','nplasma'); obj.potinn = h5readatt(obj.fullpath,'/data/input.00/','potinn'); obj.potout = h5readatt(obj.fullpath,'/data/input.00/','potout'); obj.B0 = h5readatt(obj.fullpath,'/data/input.00/','B0'); obj.Rcurv = h5readatt(obj.fullpath,'/data/input.00/','Rcurv'); obj.width = h5readatt(obj.fullpath,'/data/input.00/','width'); obj.n0 = h5readatt(obj.fullpath,'/data/input.00/','n0'); obj.temp = h5readatt(obj.fullpath,'/data/input.00/','temp'); try obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0d'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2d'); obj.it2 = h5readatt(obj.fullpath,'/data/input.00/','itparts'); catch obj.it0 = h5readatt(obj.fullpath,'/data/input.00/','it0'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it1'); obj.it1 = h5readatt(obj.fullpath,'/data/input.00/','it2'); end try try obj.nbspecies=h5readatt(obj.fullpath,'/data/part/','nbspecies'); catch obj.nbspecies=h5readatt(obj.fullpath,'/data/input.00/','nbspecies'); end obj.normalized=strcmp(h5readatt(obj.fullpath,'/data/input.00/','rawparts'),'y'); catch obj.nbspecies=1; obj.normalized=false; end try obj.nbcelldiag=h5readatt(obj.fullpath,'/data/celldiag/','nbcelldiag'); catch obj.nbcelldiag=0; end obj.omepe=sqrt(abs(obj.n0)*obj.qe^2/(obj.me*obj.eps_0)); obj.omece=obj.qe*obj.B0/obj.me; obj.npart= h5read(obj.fullpath, '/data/var0d/nbparts'); try obj.nudcol= h5read(obj.fullpath, '/data/var0d/nudcol'); catch end try obj.H0 = h5read(obj.fullpath,'/data/input.00/H0'); obj.P0 = h5read(obj.fullpath,'/data/input.00/P0'); catch obj.H0=3.2e-14; obj.P0=8.66e-25; end % Normalizations if old obj.tnorm=abs(1/obj.omepe); else obj.tnorm=min(abs(1/obj.omepe),abs(1/obj.omece)); end obj.rnorm=obj.vlight*obj.tnorm; obj.bnorm=obj.B0; obj.enorm=obj.vlight*obj.bnorm; obj.phinorm=obj.enorm*obj.rnorm; obj.vnorm=obj.vlight; % Grid data obj.rgrid= h5read(obj.fullpath, '/data/var1d/rgrid')*obj.rnorm; obj.zgrid= h5read(obj.fullpath, '/data/var1d/zgrid')*obj.rnorm; obj.dz=(obj.zgrid(end)-obj.zgrid(1))/double(obj.nz); rid=1; for i=1:length(obj.nnr) obj.dr(i)=(obj.rgrid(sum(obj.nnr(1:i))+1)-obj.rgrid(rid))/double(obj.nnr(i)); rid=rid+obj.nnr(i); end Br = h5read(obj.fullpath,'/data/fields/Br')*obj.bnorm; obj.Br= reshape(Br,length(obj.zgrid),length(obj.rgrid)); Bz = h5read(obj.fullpath,'/data/fields/Bz')*obj.bnorm; obj.Bz= reshape(Bz,length(obj.zgrid),length(obj.rgrid)); try Atheta = h5read(obj.fullpath,'/data/fields/Athet')*obj.bnorm; obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid)); [rmeshgrid,~]=meshgrid(obj.rgrid,obj.zgrid); obj.rAthet=(rmeshgrid.*obj.Athet)'; catch end obj.B=sqrt(obj.Bz.^2+obj.Br.^2); obj.costhet=(obj.Br./obj.B)'; obj.sinthet=(obj.Bz./obj.B)'; clear Br Bz try obj.t0d=h5read(obj.fullpath,'/data/var0d/time'); catch obj.t0d=obj.dt.*double(0:length(obj.epot)-1); end try for i=0:nbrst grp=sprintf('/data/input.%02i/',i); obj.Erxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'Erxt']),length(obj.zgrid),length(obj.rgrid))'*obj.enorm; obj.Ezxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'Ezxt']),length(obj.zgrid),length(obj.rgrid))'*obj.enorm; obj.potxt(:,:,i+1)=reshape(h5read(obj.fullpath,[grp,'potxt']),length(obj.zgrid),length(obj.rgrid))'*obj.phinorm; end catch end obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder'); obj.ngauss = h5read(obj.fullpath,'/data/input.00/ngauss'); obj.plasmadim = h5read(obj.fullpath,'/data/input.00/plasmadim'); obj.radii = h5read(obj.fullpath,'/data/input.00/radii'); obj.epot = h5read(obj.fullpath,'/data/var0d/epot'); obj.ekin = h5read(obj.fullpath,'/data/var0d/ekin'); obj.etot = h5read(obj.fullpath,'/data/var0d/etot'); try obj.etot0 = h5read(obj.fullpath,'/data/var0d/etot0'); obj.eerr = obj.etot-obj.etot0; catch obj.eerr = obj.etot-obj.etot(2); end if(obj.normalized) obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,1); obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,1); obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,1); else obj.pot=gridquantity(obj.fullpath,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,obj.phinorm); obj.Er=gridquantity(obj.fullpath,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,obj.enorm); obj.Ez=gridquantity(obj.fullpath,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,obj.enorm); end try obj.t2d = h5read(obj.fullpath,'/data/fields/time'); catch info=h5info(obj.fullpath,'/data/fields/partdensity'); obj.t2d=obj.dt*(0:info.objspace.Size(2)-1)*double(obj.it1); end try info=h5info(obj.fullpath,'/data/fields/moments'); obj.femorder = h5read(obj.fullpath,'/data/input.00/femorder'); kr=obj.femorder(2)+1; obj.knotsr=augknt(obj.rgrid,kr); kz=obj.femorder(1)+1; obj.knotsz=augknt(obj.zgrid,kz); try obj.CellVol= reshape(h5read(obj.fullpath,'/data/fields/volume'),length(obj.knotsz)-kz,length(obj.knotsr)-kr); obj.CellVol=permute(obj.CellVol,[2,1,3])*obj.rnorm^3; catch zvol=fnder(spmak(obj.knotsz,ones(1,length(obj.knotsz)-kz)), -1 ); rvol=fnder(spmak(obj.knotsr,2*pi*[obj.rgrid' 2*obj.rgrid(end)-obj.rgrid(end-1)]), -1 ); ZVol=diff(fnval(zvol,obj.knotsz)); RVol=diff(fnval(rvol,obj.knotsr)); obj.CellVol=RVol(3:end-1)*ZVol(3:end-1)'; obj.CellVol=padarray(obj.CellVol,[1,1],'replicate','post'); end try obj.geomweight = h5read(obj.fullpath,'/data/input.00/geometry/geomweight'); obj.geomweight= reshape(obj.geomweight,length(obj.zgrid),length(obj.rgrid),[]); obj.geomweight = permute(obj.geomweight,[2,1,3]); catch obj.geomweight=ones(length(obj.rgrid),length(obj.zgrid),3); end try obj.dirichletweight = h5read(obj.fullpath,'/data/input.00/geometry/dirichletweight'); obj.dirichletweight= reshape(obj.dirichletweight,length(obj.zgrid),length(obj.rgrid),[]); obj.dirichletweight = permute(obj.dirichletweight,[2,1,3]); catch obj.dirichletweight=obj.geomweight; end try obj.gtilde = h5read(obj.fullpath,'/data/input.00/geometry/gtilde'); obj.gtilde= reshape(obj.gtilde,length(obj.zgrid),length(obj.rgrid),[]); obj.gtilde = permute(obj.gtilde,[2,1,3]); catch obj.gtilde=zeros(length(obj.rgrid),length(obj.zgrid),3); end geomweight=ones(length(obj.rgrid),length(obj.zgrid)); if(obj.normalized) obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, 1, geomweight, 1); obj.phi=splinequantity(obj.fullpath,'/data/fields/phi', obj.knotsr, obj.knotsz, obj.femorder, 1, obj.geomweight(:,:,1), -1); else obj.N=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), geomweight, 1); end obj.fluidUR=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 2); obj.fluidUTHET=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 3); obj.fluidUZ=splinevelocity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, geomweight, 4); if(obj.normalized) obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me, geomweight); obj.fluidEkin=splineenergy(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.me*0.5, geomweight); else obj.Presstens=splinepressure(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim, geomweight); obj.fluidEkin=splineenergy(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim*0.5, geomweight); end try obj.celltype=h5read(obj.fullpath,'/data/input.00/geometry/ctype')'; obj.linked_s=h5read(obj.fullpath,'/data/input.00/geometry/linked_s'); obj.bsplinetype=h5read(obj.fullpath,'/data/input.00/geometry/bsplinetype'); obj.bsplinetype=reshape(obj.bsplinetype,length(obj.knotsz)-kz,length(obj.knotsr)-kr); catch obj.celltype=[]; obj.linked_s=[]; end catch obj.CellVol=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)'; obj.CellVol=obj.CellVol'; obj.N=griddensity(obj.fullpath, '/data/fields/partdensity', sum(obj.nnr)+1, obj.nz+1, obj.CellVol, abs(obj.qsim/obj.qe), true); obj.fluidUR=gridquantity(obj.fullpath, '/data/fields/fluidur', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); obj.fluidUTHET=gridquantity(obj.fullpath, '/data/fields/fluiduthet', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); obj.fluidUZ=gridquantity(obj.fullpath, '/data/fields/fluiduz', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true); end % If we have a maxwellian source, read its parameters try obj.maxwellsrce.rlim=h5read(obj.fullpath, '/data/input.00/maxwellsource/rlimits'); obj.maxwellsrce.zlim=h5read(obj.fullpath, '/data/input.00/maxwellsource/zlimits'); obj.maxwellsrce.frequency=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','frequency'); obj.maxwellsrce.radialtype=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','radialtype'); obj.maxwellsrce.temperature=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','temperature'); obj.maxwellsrce.time_end=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','time_end'); obj.maxwellsrce.time_start=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','time_start'); obj.maxwellsrce.vth=h5readatt(obj.fullpath, '/data/input.00/maxwellsource','vth'); obj.maxwellsrce.rate=obj.maxwellsrce.frequency*obj.weight/(pi*(diff(obj.maxwellsrce.rlim.^2))*diff(obj.maxwellsrce.zlim)); obj.maxwellsrce.current=obj.maxwellsrce.frequency*obj.weight*obj.qe; obj.maxwellsrce.present=true; catch obj.maxwellsrce.present=false; end %% load neutcol parameters try obj.neutcol.neutdens=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','neutdens')); obj.neutcol.neutpressure=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','neutpressure')); obj.neutcol.scatter_fac=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','scatter_fac')); obj.neutcol.Eion=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','Eion')); obj.neutcol.E0=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','E0')); obj.neutcol.Escale=double(h5readatt(obj.fullpath, '/data/input.00/neutcol','Escale')); try obj.neutcol.io_cross_sec=double(h5read(obj.fullpath, '/data/input.00/neutcol/io_cross_sec')); obj.neutcol.io_cross_sec(:,2)=obj.neutcol.io_cross_sec(:,2)*obj.rnorm^2; obj.neutcol.io_cross_sec(:,3)=[log(obj.neutcol.io_cross_sec(2:end,2)./obj.neutcol.io_cross_sec(1:end-1,2))... ./log(obj.neutcol.io_cross_sec(2:end,1)./obj.neutcol.io_cross_sec(1:end-1,1)); 0]; obj.neutcol.iom_cross_sec=zeros(500,3); obj.neutcol.iom_cross_sec(:,1)=logspace(log10(obj.neutcol.Eion+0.001),log10(5e4),size(obj.neutcol.iom_cross_sec,1)); obj.neutcol.iom_cross_sec(:,2)=obj.sigmiopre(obj.neutcol.iom_cross_sec(:,1),true); obj.neutcol.iom_cross_sec(:,3)=abs([log(obj.neutcol.iom_cross_sec(2:end,2)./obj.neutcol.iom_cross_sec(1:end-1,2))... ./log(obj.neutcol.iom_cross_sec(2:end,1)./obj.neutcol.iom_cross_sec(1:end-1,1)); 0]); catch obj.neutcol.io_cross_sec=[]; obj.neutcol.iom_cross_sec=[]; end try obj.neutcol.ela_cross_sec=double(h5read(obj.fullpath, '/data/input.00/neutcol/ela_cross_sec')); obj.neutcol.ela_cross_sec(:,2)=obj.neutcol.ela_cross_sec(:,2)*obj.rnorm^2; obj.neutcol.ela_cross_sec(:,3)=[log(obj.neutcol.ela_cross_sec(2:end,2)./obj.neutcol.ela_cross_sec(1:end-1,2))... ./log(obj.neutcol.ela_cross_sec(2:end,1)./obj.neutcol.ela_cross_sec(1:end-1,1)); 0]; catch obj.neutcol.ela_cross_sec=[]; end obj.neutcol.present=true; catch obj.neutcol.present=false; end %% load spline boundaries try obj.spl_bound.nbsplines=h5readatt(obj.fullpath, '/data/input.00/geometry_spl','nbsplines'); for i=1:obj.spl_bound.nbsplines splgroup=sprintf('/data/input.00/geometry_spl/%02d',i); obj.spl_bound.boundary(i).knots=h5read(obj.fullpath,sprintf('%s/knots',splgroup)); obj.spl_bound.boundary(i).Dval=h5readatt(obj.fullpath,splgroup,'Dirichlet_val'); obj.spl_bound.boundary(i).coefs=reshape(h5read(obj.fullpath,sprintf('%s/pos',splgroup)),2,[])'; obj.spl_bound.boundary(i).order=h5readatt(obj.fullpath,splgroup,'order'); obj.spl_bound.boundary(i).kind=h5readatt(obj.fullpath,splgroup,'kind'); + try + obj.spl_bound.boundary(i).type=h5readatt(obj.fullpath,splgroup,'type'); + catch + end obj.spl_bound.boundary(i).fun=spmak(obj.spl_bound.boundary(i).knots,obj.spl_bound.boundary(i).coefs'); end catch obj.spl_bound.nbsplines=0; end %% load non ideal power supply parameters try obj.psupply.targetbias=h5readatt(obj.fullpath, '/data/input.00/psupply','targetbias'); obj.psupply.expdens=h5readatt(obj.fullpath, '/data/input.00/psupply','expdens'); obj.psupply.PSresistor=h5readatt(obj.fullpath, '/data/input.00/psupply','PSresistor'); obj.psupply.geomcapacitor=h5readatt(obj.fullpath, '/data/input.00/psupply','geomcapacitor'); obj.psupply.nbhdt=h5readatt(obj.fullpath, '/data/input.00/psupply','nbhdt'); obj.psupply.biases=h5read(obj.fullpath, '/data/var0d/biases'); obj.psupply.current=h5read(obj.fullpath, '/data/var0d/current'); obj.psupply.tau=obj.psupply.PSresistor*obj.psupply.geomcapacitor*obj.psupply.expdens/obj.neutcol.neutdens; obj.psupply.active=true; obj.psupply.bdpos=h5read(obj.fullpath, '/data/input.00/psupply/bdpos'); catch obj.psupply.active=false; end % Read the main particles parameters if(readparts) if(obj.normalized) obj.R = h5partsquantity(obj.fullpath,'/data/part','R'); obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z'); else obj.R = h5partsquantity(obj.fullpath,'/data/part','R',obj.rnorm); obj.Z = h5partsquantity(obj.fullpath,'/data/part','Z',obj.rnorm); end try obj.THET = h5partsquantity(obj.fullpath,'/data/part','THET'); catch clear obj.THET end try obj.Rindex=h5partsquantity(obj.fullpath,'/data/part','Rindex'); obj.Zindex=h5partsquantity(obj.fullpath,'/data/part','Zindex'); catch clear obj.Rindex obj.Zindex end vscale=obj.vnorm; obj.VR = h5partsquantity(obj.fullpath,'/data/part','UR',vscale); obj.VZ = h5partsquantity(obj.fullpath,'/data/part','UZ',vscale); obj.VTHET= h5partsquantity(obj.fullpath,'/data/part','UTHET',vscale); if(obj.normalized) obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe); else obj.partepot = h5partsquantity(obj.fullpath,'/data/part','pot',sign(obj.qsim)*obj.qe*obj.phinorm); end try obj.partindex = h5partsquantity(obj.fullpath,'/data/part/','partindex'); catch end if(obj.nbspecies >1) obj.species=h5parts.empty(obj.nbspecies,0); for i=2:obj.nbspecies obj.species(i-1)=h5parts(obj.fullpath,sprintf('/data/part/%2d',i),obj); end end end try obj.tpart = h5read(obj.fullpath,'/data/part/time'); obj.nbparts = h5read(obj.fullpath,'/data/part/Nparts'); catch obj.nbparts=obj.npart; obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2); end if(obj.nbcelldiag > 0) obj.celldiag=h5parts.empty; j=0; for i=1:obj.nbcelldiag nbparts=h5read(obj.fullpath,sprintf('%s/Nparts',sprintf('/data/celldiag/%02d',i))); if (sum(nbparts)>0) j=j+1; obj.celldiag(j)=h5parts(obj.fullpath,sprintf('/data/celldiag/%02d',i),obj); obj.celldiag(j).rindex=double(h5readatt(obj.fullpath, sprintf('/data/celldiag/%02d',i),'rindex'))+(1:2); obj.celldiag(j).zindex=double(h5readatt(obj.fullpath, sprintf('/data/celldiag/%02d',i),'zindex'))+(1:2); end end end end %------------------------------------------ % Functions for accesing secondary simulation quantities function Atheta=Atheta(obj,R,Z) %% returns the magnetic vector potential at position R,Z interpolated from stored Athet in h5 file % halflz=(obj.zgrid(end)+obj.zgrid(1))/2; % Atheta=0.5*obj.B0*(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)... % .*besseli(1,2*pi*R/obj.width).*cos(2*pi*(Z-halflz)/obj.width)); Atheta=interp2(obj.rgrid,obj.zgrid,obj.Athet,R,Z); end function quantity=H(obj,indices) %% computes the total energy for the main specie % for the particle with index indices{1} at timepart step indices{2} % which is time obj.timepart(indices{2}) if strcmp(indices{1},':') p=1:obj.VR.nparts;% if nothing is defined we load all particles else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); %if nothing is defined all time steps are considered else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,1)>2 track=indices{3}; else track=false; end quantity=0.5*obj.me*(obj.VR(p,t,track).^2+obj.VTHET(p,t,track).^2+obj.VZ(p,t,track).^2)+obj.partepot(p,t,track); end function quantity=P(obj,indices) %P computes the canonical angular momentum for the main specie % for the particle with index indices{1} at timepart step indices{2} % which is time obj.timepart(indices{2}) if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,1)>2 track=indices{3}; else track=false; end quantity=obj.R(p,t,track).*(obj.VTHET(p,t,track)*obj.me+sign(obj.qsim)*obj.qe*obj.Atheta(obj.R(p,t,track),obj.Z(p,t,track))); end function quantity=Vpar(obj,varargin) %Vpar Computes the parallel velocity for the main specie % for the particle with index indices{1} at timepart step indices{2} % which is time obj.timepart(indices{2}) if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,1)>2 track=indices{3}; else track=false; end Zp=obj.Z(p,t,track);% get the particle axial positon Rp=obj.R(p,t,track);% get the particle radial position % interpolate the magnetic field at the particle position Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima'); Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima'); Bp=sqrt(Bzp.^2+Brp.^2); % calculate the projection angle of the radial and axial % directions on the magnetic field line Costhet=Bzp./Bp; Sinthet=Brp./Bp; % calculate the actuale parallel velocity quantity=obj.VR(p,t,track).*Sinthet+obj.VZ(p,t,track).*Costhet; end function quantity=Vperp(obj,varargin) %Vperp Computes the perpendicular velocity in the guidind center reference frame, % for the main specie particle indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,2)>2 track=indices{3}; else track=false; end % if gcs is true, gives the perpendicular velocity in the % guiding center system by substracting the EXB azimuthal % velocity % else gives the total perpendicular velocity if size(indices,2)>3 gcs=indices{4}; else gcs=false; end % get the particle position Zp=obj.Z(p,t,track); Rp=obj.R(p,t,track); % interpolate the magnetic field at the particle position Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima'); Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima'); Bp=sqrt(Bzp.^2+Brp.^2); % calculate the projecting angles Costhet=Bzp./Bp; Sinthet=Brp./Bp; Vdrift=zeros(size(Zp)); if gcs % for each particle and each timestep % calculate the azimuthal ExB drift velocity for j=1:length(t) [~, tfield]=min(abs(obj.t2d-obj.tpart(t(j)))); timeEr=obj.Er(:,:,tfield); timeEz=obj.Ez(:,:,tfield); %posindE=sub2ind(size(timeEr),Rind(:,j),Zind(:,j)); timeErp=interp2(obj.zgrid,obj.rgrid,timeEr,Zp(:,j),Rp(:,j)); timeEzp=interp2(obj.zgrid,obj.rgrid,timeEz,Zp(:,j),Rp(:,j)); Vdrift(:,j)=(timeEzp.*Brp(:,j)-timeErp.*Bzp(:,j))./Bp(:,j).^2; end end % calculate the perpendicular velocity quantity=sqrt((obj.VTHET(p,t,track)-Vdrift).^2+(obj.VR(p,t,track).*Costhet-obj.VZ(p,t,track).*Sinthet).^2); end function quantity=cyclphase(obj,varargin) %cyclphase Computes the cyclotronic phase for the main specie % for particles with indices{1} at time indices{2} if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,2)>2 track=indices{3}; else track=false; end Zp=obj.Z(p,t,track); Rp=obj.R(p,t,track); % [~, zind(1)]=min(abs(obj.zgrid-0.005262)); % [~, zind(2)]=min(abs(obj.zgrid-0.006637)); % [~, rind(1)]=min(abs(obj.rgrid-0.0784)); % [~, rind(2)]=min(abs(obj.rgrid-0.07861)); % indices=Zp=obj.zgrid(zind(1)) &... % Rp=obj.rgrid(rind(1)); %Zp=Zp(indices); %Rp=Rp(indices); %p=indices; Bzp=interp2(obj.zgrid,obj.rgrid,obj.Bz',Zp,Rp,'makima'); Brp=interp2(obj.zgrid,obj.rgrid,obj.Br',Zp,Rp,'makima'); Bp=sqrt(Bzp.^2+Brp.^2); Costhet=Bzp./Bp; Sinthet=Brp./Bp; % compute the projection of the perpendicular velocity in the % radial direction vr=(obj.VR(p,t,track).*Costhet-obj.VZ(p,t,track).*Sinthet); % Get the perpendicular velocity vperp=obj.Vperp(p,t,track,true); vr=vr(indices); vperp=vperp(indices); cospsi=vr./vperp; quantity=acos(cospsi); end function p=borderpoints(obj,subdiv) %borderpoints Return a cell array containing the curves %defining the boundary of the domain % for each boundary p(1,:) and p(2,:) give axial and radial position % for each boundary p(3,:) and p(4,:) give axial and radial normals %gw= contourc(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0]) p=cell(0,0); if nargin<2 subdiv=1; end ndiv=sum(subdiv); %outer cylinder if any(obj.geomweight(end,:,1)>=0) idp=ceil(length(obj.zgrid)/ndiv); imin=1; for j=1:length(subdiv) imax=min(imin+subdiv(j)*idp-1,length(obj.zgrid)); p{end+1}=[obj.zgrid(imin:imax)';obj.rgrid(end)*ones(imax-imin+1,1)'; zeros(imax-imin+1,1)';ones(imax-imin+1,1)']; imin=imax; end end %inner cylinder if any(obj.geomweight(1,:,1)>=0) idp=ceil(length(obj.zgrid)/ndiv); imin=1; for j=1:length(subdiv) imax=min(imin+subdiv(j)*idp-1,length(obj.zgrid)); p{end+1}=[obj.zgrid(imin:imax)';obj.rgrid(1)*ones(imax-imin+1,1)'; zeros(imax-imin+1,1)';-ones(imax-imin+1,1)']; imin=imax; end end if obj.walltype==2 % We have an elliptic insert that we want to isolate gw=obj.ellipseborder; zpos=obj.zgrid(obj.zgrid<(min(gw(1,:))) | obj.zgrid>(max(gw(1,:)))); p{2}=[zpos,obj.rgrid(end)*ones(size(zpos))]'; p{1}=[obj.zgrid';obj.rgrid(1)*ones(size(obj.zgrid))']; gw=obj.ellipseborder; p{3}=gw; elseif obj.walltype~=0 % extract all the walls gw=contourc(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0]); [x,y,~]=C2xyz(gw); for i=1:length(x) %subdiv=[4,1,2]; ndiv=sum(subdiv); idp=ceil(length(x{i})/ndiv); imin=1; for j=1:length(subdiv) imax=min(imin+subdiv(j)*idp-1,length(x{i})); p{end+1}=[x{i}(imin:imax);y{i}(imin:imax)]; imin=imax; end end end % figure % for i=1:length(p) % plot(p{i}(1,:),p{i}(2,:)) % hold on % end end function p=ellipseborder(obj) %ellipseborder returns the boundary points defining the %elliptic insert z=linspace(-0.998,0.998,1000)*obj.z_r; p=zeros(4,length(z)); for i=1:length(z) p(1,i)=z(i)+obj.z_0; p(2,i)=obj.r_0-obj.r_r*sqrt(1-(z(i)/obj.z_r)^2); p(3,i)=2/(obj.z_r^2)*(z(i)); p(4,i)=2/(obj.r_r^2)*(p(2,i)-obj.r_0); end norm=sqrt(p(3,:).^2+p(4,:).^2); p(3,:)=double(obj.interior)*p(3,:)./norm; p(4,:)=double(obj.interior)*p(4,:)./norm; end function charge=totcharge(obj,fieldstep) % Integrates the density profile over the full volume to obtain % the total number of electrons in the volume n=splinedensity(obj.fullpath, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder,ones(size(obj.CellVol)), 1, 1); charge=sum(sum(n(:,:,fieldstep))); end function Gamma=Axialflux(obj,timestep,zpos) % Computes the axial particle flux n*Uz at timestep timestep and axial position zpos Gamma=obj.fluidUZ(:,zpos,timestep).*obj.N(:,zpos,timestep); end function Gamma=Metallicflux(obj,timestep,subdiv) % Computes the particle flux at time obj.t2d(timestep) on the % metallic boundaries if nargin<3 subdiv=1; end % We find the borderpoints p=obj.borderpoints(subdiv); gamma=cell(size(p)); Nr=cell(size(p)); Nz=cell(size(p)); for i=1:length(p) bp=p{i}; if size(bp,1)==2 % We get the normals at these positions and normalise them Nr{i}=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,3),bp(1,:),bp(2,:)); Nz{i}=-interp2(obj.zgrid,obj.rgrid,obj.geomweight(:,:,2),bp(1,:),bp(2,:)); norm=sqrt(Nr{i}.^2+Nz{i}.^2); Nr{i}=Nr{i}./norm; Nz{i}=Nz{i}./norm; else Nr{i}=bp(4,:); Nz{i}=bp(3,:); end gamma{i}=zeros(size(bp,2),length(timestep)); end [z,r]=ndgrid(obj.zgrid,obj.rgrid); N=obj.N(:,:,timestep(1)); n=griddedInterpolant(z,r,N'); uz=griddedInterpolant(z,r,obj.fluidUZ(:,:,timestep(1))'); ur=griddedInterpolant(z,r,obj.fluidUR(:,:,timestep(1))'); % we get the density and fluid velocities at the desired time % steps and interpolate them at the boundary position for j=1:length(timestep) n.Values=obj.N(:,:,timestep(j))'; uz.Values=obj.fluidUZ(:,:,timestep(j))'; ur.Values=obj.fluidUR(:,:,timestep(j))'; for i=1:length(p) bp=p{i}; gamma{i}(:,j)=n(bp(1:2,:)').*(ur(bp(1:2,:)').*Nr{i}'+uz(bp(1:2,:)').*Nz{i}'); end end % return the boundary position p and the corresponding flux % gamma Gamma.p=p; Gamma.gamma=gamma; end - function [pot] = PotentialWell(obj,fieldstep) + function [pot] = PotentialWell(obj,fieldstep,fieldaligned) %PotentialWell Computes the potential well at the given timestep on the FEM grid points % interpolates the model data on rgrid and zgrid + if nargin<3 + fieldaligned=false; + end model=obj.potentialwellmodel(fieldstep); z=model.z; - r=model.r; modpot=model.pot; - [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); - pot=zeros(length(obj.zgrid),length(obj.rgrid),length(fieldstep)); + + if fieldaligned + r=model.rathet; + lvls=linspace(min(obj.rAthet(:)),max(obj.rAthet(:)),400); + [Zmesh,Rmesh]=meshgrid(obj.zgrid,lvls); + else + r=model.r; + [Zmesh,Rmesh]=meshgrid(obj.zgrid,obj.rgrid); + end + pot=zeros(size(Zmesh,2),size(Zmesh,1),length(fieldstep)); for i=1:length(fieldstep) pot(:,:,i)=griddata(z,r,modpot(:,i),Zmesh,Rmesh)'; end end function Epar = Epar(obj,fieldstep) % Computes the electric field component parallel to the magnetic field line Epar=obj.Er(:,:,fieldstep).*(obj.Br./obj.B)' + (obj.Bz./obj.B)'.*obj.Ez(:,:,fieldstep); end function Eperp = Eperp(obj,fieldstep) % Computes the electric field component perpendicular to the magnetic field line Eperp=obj.Er(:,:,fieldstep).*(obj.Bz./obj.B)' - (obj.Br./obj.B)'.*obj.Ez(:,:,fieldstep); end function Ekin = Ekin(obj,varargin) %Ekin Computes the classical kinetic energy of particles indices{1} at % time obj.tpart(indices{2}) in Joules if(~iscell(varargin)) indices=mat2cell(varargin); else indices=varargin; end if strcmp(indices{1},':') p=1:obj.R.nparts; else p=indices{1}; end if strcmp(indices{2},':') t=1:length(obj.tpart); else t=indices{2}; end % if track is true we look at specific particles with their % index and follow them in time % if it is false we just care about the distribution function % and specific particles can have different positions in the % resulting array for each timestep if size(indices,1)>2 track=indices{3}; else track=false; end Vr=obj.VR(p,t,track); Vthet= obj.VTHET(p,t,track); Vz=obj.VZ(p,t,track); Ekin=0.5*obj.msim/obj.weight*(Vr.^2+Vthet.^2+Vz.^2); end function sig=sigio(obj,E,init) %sigio returns the total ionisation cross-section in m^2 % at energy E[eV] % init is only used during the loading of the h5 file if nargin <3 init=false; end sig=zeros(size(E)); if(~init &&( ~obj.neutcol.present || isempty(obj.neutcol.io_cross_sec))) sig=zeros(size(E)); return end for i=1:length(E(:)) if(E(i)>obj.neutcol.Eion) sig(ind2sub(size(E),i))=obj.fit_cross_sec(E(ind2sub(size(E),i)),obj.neutcol.io_cross_sec); end end end function sig=sigmio(obj,E) %sigmio returns the total ionisation cross-section for momentum exchange for the incoming electron in m^2 % at energy E[eV] sig=zeros(size(E)); if(~obj.neutcol.present || isempty(obj.neutcol.iom_cross_sec)) return end for i=1:length(E(:)) if(E(i)>obj.neutcol.Eion) sig(ind2sub(size(E),i))=obj.fit_cross_sec(E(ind2sub(size(E),i)),obj.neutcol.iom_cross_sec); end end end function sigm=sigmela(obj,E) %sigmela returns the elastic collision cross-section for momentum exchange for the incoming electron in m^2 % at energy E[eV] sigm=zeros(size(E)); if(~obj.neutcol.present || isempty(obj.neutcol.ela_cross_sec)) return end for i=1:length(E(:)) sigm(ind2sub(size(E),i))=obj.fit_cross_sec(E(ind2sub(size(E),i)),obj.neutcol.ela_cross_sec); end end function sig=sigela(obj,E) %sigmela returns the elastic collision cross-section for the incoming electron in m^2 % at energy E[eV] % if used this will give the frequency of elastic collisions E0=obj.neutcol.E0; chi=E./(0.25*E0+E); sig=(2*chi.^2)./((1-chi).*((1+chi).*log((1+chi)./(1-chi))-2*chi)).*obj.sigmela(E); end function [Forces, Density]=Forcespline(obj,it,fdens,getmean) %Forcesplie calculates the fluid force terms in each direction %at time obj.t2d(it) % if fdens return the force density in N/m^3 othewise give % the force in N % if getmean return only the time averaged quanties over % time samples[it(1)...it(end] if strcmp(it,':') it=floor(0.95*size(obj.t2d)):size(obj.t2d)-1; end if nargin<3 fdens=true; end if nargin <4 getmean=false; end % To be able to calculate the centered finite difference in % time, we remove the first and last time indices it(it<2)=[]; it(it>length(obj.t2d)-1)=[]; m_e=obj.msim/obj.weight; q_e=obj.qsim/obj.weight; n=obj.N(:,:,it); [r,~]=meshgrid(obj.rgrid,obj.zgrid); Rinv=1./r'; Rinv(isinf(Rinv))=0; % get inverse of density to get the force in N Density.N=n; invn=1./n; invn(isnan(invn) | isinf(invn))=0; % Calculate electric forces Eforcer=q_e*obj.Er(:,:,it); Eforcez=q_e*obj.Ez(:,:,it); Dragforcer=zeros(size(n,1),size(n,2),size(n,3)); Dragforcethet=zeros(size(n,1),size(n,2),size(n,3)); Dragforcez=zeros(size(n,1),size(n,2),size(n,3)); time=obj.t2d(it); Forces.it=it; Forces.time=time; if getmean if ~fdens n=ones(size(n)); end % Electric forces Forces.Eforcer=mean(n.*q_e.*obj.Er(:,:,it),3); Forces.Eforcez=mean(n.*q_e.*obj.Ez(:,:,it),3); % Magnetic forces Forces.Bforcer=mean(q_e.*obj.fluidUTHET(:,:,it).*obj.Bz'.*n,3); Forces.Bforcethet=mean(q_e.*(obj.fluidUZ(:,:,it).*obj.Br'-obj.fluidUR(:,:,it).*obj.Bz').*n,3); Forces.Bforcez=mean(-q_e.*obj.fluidUTHET(:,:,it).*obj.Br'.*n,3); % Inertial forces Forces.inertforcer=mean(-m_e.*n.*(-obj.fluidUTHET(:,:,it).^2.*Rinv... +obj.fluidUR(:,:,it).*obj.fluidUR.der(:,:,it,[1 0])... +obj.fluidUZ(:,:,it).*obj.fluidUR.der(:,:,it,[0 1])),3); Forces.inertforcethet=mean(-m_e*n.*(obj.fluidUR(:,:,it).*obj.fluidUTHET(:,:,it).*Rinv... +obj.fluidUR(:,:,it).*obj.fluidUTHET.der(:,:,it,[1 0])... +obj.fluidUZ(:,:,it).*obj.fluidUTHET.der(:,:,it,[0 1])),3); Forces.inertforcez=mean(-m_e*n.*(obj.fluidUR(:,:,it).*obj.fluidUZ.der(:,:,it,[1 0])... +obj.fluidUZ(:,:,it).*obj.fluidUZ.der(:,:,it,[0 1])),3); % Pressure forces Forces.Pressforcer=mean(-n.*( squeeze(obj.Presstens.der(1,:,:,it,[1 0]))... + squeeze(obj.Presstens(1,:,:,it) - obj.Presstens(4,:,:,it)).*Rinv... + squeeze(obj.Presstens.der(3,:,:,it,[0 1])))... .*invn,3); Forces.Pressforcethet=mean(-n.*( squeeze(obj.Presstens.der(2,:,:,it,[1 0]))... + squeeze(obj.Presstens.der(5,:,:,it,[0 1])) ... + 2*squeeze(obj.Presstens(2,:,:,it)).*Rinv ... ).*invn,3); Forces.Pressforcez=mean(-n.*( squeeze(obj.Presstens.der(3,:,:,it,[1 0]))... + squeeze(obj.Presstens(3,:,:,it)).*Rinv... + squeeze(obj.Presstens.der(6,:,:,it,[0 1])) )... .*invn,3); % ellastic coll drag forces if( obj.neutcol.present) Ek=squeeze(obj.fluidEkin(1,:,:,it)+obj.fluidEkin(2,:,:,it)+obj.fluidEkin(3,:,:,it)); sigm=obj.sigmela(Ek/obj.qe)+obj.sigio(Ek/obj.qe)+obj.sigmio(Ek/obj.qe); dragfreq=obj.neutcol.neutdens.*sigm.*sqrt(2*obj.weight/obj.msim*Ek); Forces.Dragforcer=mean(-m_e*n.*dragfreq.*obj.fluidUR(:,:,it),3); Forces.Dragforcethet=mean(-m_e*n.*dragfreq.*obj.fluidUTHET(:,:,it),3); Forces.Dragforcez=mean(-m_e*n.*dragfreq.*obj.fluidUZ(:,:,it),3); else Forces.Dragforcer=0; Forces.Dragforcethet=0; Forces.Dragforcez=0; end % effective drag frequency due to the maxwellian source if( obj.maxwellsrce.present) dragfreqsrc=obj.maxwellsrce.frequency*obj.weight/(pi*diff(obj.maxwellsrce.zlim)*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)).*invn; dragfreqsrc(isinf(dragfreqsrc))=0; Forces.Dragforcer=Forces.Dragforcer+mean(-n.*m_e.*dragfreqsrc.*obj.fluidUR(:,:,it),3); Forces.Dragforcethet=Forces.Dragforcethet+mean(-n.*m_e.*dragfreqsrc.*obj.fluidUTHET(:,:,it),3); Forces.Dragforcez=Forces.Dragforcez+mean(-n.*m_e.*dragfreqsrc.*obj.fluidUZ(:,:,it),3); end % Time derivative for fluid accelleration cdt=(obj.t2d(it+1)-obj.t2d(it-1)); cdt=reshape(cdt,1,1,[]); Forces.durdt=mean(m_e*(obj.fluidUR(:,:,it+1)-obj.fluidUR(:,:,it-1))./cdt,3); Forces.duthetdt=mean(m_e*(obj.fluidUTHET(:,:,it+1)-obj.fluidUTHET(:,:,it-1))./cdt,3); Forces.duzdt=mean(m_e*(obj.fluidUZ(:,:,it+1)-obj.fluidUZ(:,:,it-1))./cdt,3); else % Allocate memory Bforcer=zeros(size(n,1),size(n,2),size(n,3)); Bforcez=zeros(size(n,1),size(n,2),size(n,3)); Bforcethet=zeros(size(n,1),size(n,2),size(n,3)); inertforcer=zeros(size(n,1),size(n,2),size(n,3)); inertforcez=zeros(size(n,1),size(n,2),size(n,3)); inertforcethet=zeros(size(n,1),size(n,2),size(n,3)); Pressforcer=zeros(size(n,1),size(n,2),size(n,3)); Pressforcethet=zeros(size(n,1),size(n,2),size(n,3)); Pressforcez=zeros(size(n,1),size(n,2),size(n,3)); durdt=zeros(size(n,1),size(n,2),size(n,3)); duthetdt=zeros(size(n,1),size(n,2),size(n,3)); duzdt=zeros(size(n,1),size(n,2),size(n,3)); fluiduThet=obj.fluidUTHET(:,:,it); Density.fluiduThet=fluiduThet; for j=1:size(n,3) % Magnetic forces Bforcer(:,:,j)=q_e.*fluiduThet(:,:,j).*obj.Bz'; Bforcethet(:,:,j)=q_e.*(obj.fluidUZ(:,:,it(j)).*obj.Br'-obj.fluidUR(:,:,it(j)).*obj.Bz'); Bforcez(:,:,j)=-q_e.*fluiduThet(:,:,j).*obj.Br'; % Inertial forces inertforcer(:,:,j)=-m_e.*(-fluiduThet(:,:,j).^2.*Rinv... +obj.fluidUR(:,:,it(j)).*obj.fluidUR.der(:,:,it(j),[1 0])... +obj.fluidUZ(:,:,it(j)).*obj.fluidUR.der(:,:,it(j),[0 1])); inert1=obj.fluidUR(:,:,it(j)).*fluiduThet(:,:,j).*Rinv; inert2=obj.fluidUR(:,:,it(j)).*obj.fluidUTHET.der(:,:,it(j),[1 0]); inert3=obj.fluidUZ(:,:,it(j)).*obj.fluidUTHET.der(:,:,it(j),[0 1]); inertforcethet(:,:,j)=-m_e.*(inert1... +inert2... +inert3); inertforcez(:,:,j)=-m_e.*(obj.fluidUR(:,:,it(j)).*obj.fluidUZ.der(:,:,it(j),[1 0])... +obj.fluidUZ(:,:,it(j)).*obj.fluidUZ.der(:,:,it(j),[0 1])); % Pressure forces Pr1=squeeze(obj.Presstens.der(1,:,:,it(j),[1 0])); Pr2=squeeze(obj.Presstens(1,:,:,it(j)) - obj.Presstens(4,:,:,it(j))).*Rinv; Pr3=squeeze(obj.Presstens.der(3,:,:,it(j),[0 1])); Pressforcer(:,:,j)=-( Pr1... + Pr2... + Pr3 )... .*invn(:,:,j); Pthet1=squeeze(obj.Presstens.der(2,:,:,it(j),[1 0])); Pthet2=squeeze(obj.Presstens.der(5,:,:,it(j),[0 1])); Pthet3=2*squeeze(obj.Presstens(2,:,:,it(j))).*Rinv; Pressforcethet(:,:,j)=-( Pthet1... + Pthet2 ... + Pthet3 ... ).*invn(:,:,j); Pz1=squeeze(obj.Presstens.der(3,:,:,it(j),[1 0])); Pz2=squeeze(obj.Presstens(3,:,:,it(j))).*Rinv; Pz3=squeeze(obj.Presstens.der(6,:,:,it(j),[0 1])); Pressforcez(:,:,j)=-( Pz1... + Pz2... + Pz3 )... .*invn(:,:,j); % ellastic coll drag forces if( obj.neutcol.present) Ek=squeeze(obj.fluidEkin(1,:,:,it(j))+obj.fluidEkin(2,:,:,it(j))+obj.fluidEkin(3,:,:,it(j))); sigm=obj.sigmela(Ek/obj.qe)+obj.sigio(Ek/obj.qe)+obj.sigmio(Ek/obj.qe); dragfreq=obj.neutcol.neutdens.*sigm.*sqrt(2*obj.weight/obj.msim*Ek); Dragforcer(:,:,j)=-m_e*dragfreq.*obj.fluidUR(:,:,it(j)); Dragforcethet(:,:,j)=-m_e*dragfreq.*obj.fluidUTHET(:,:,it(j)); Dragforcez(:,:,j)=-m_e*dragfreq.*obj.fluidUZ(:,:,it(j)); end % effective drag frequency due to the maxwellian source if( obj.maxwellsrce.present) dragfreqsrc=obj.maxwellsrce.frequency*obj.weight/(pi*diff(obj.maxwellsrce.zlim)*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2))*invn(:,:,j); dragfreqsrc(isinf(dragfreqsrc))=0; Dragforcer(:,:,j)=Dragforcer(:,:,j)+-m_e*dragfreqsrc.*obj.fluidUR(:,:,it(j)); Dragforcethet(:,:,j)=Dragforcethet(:,:,j)+-m_e*dragfreqsrc.*obj.fluidUTHET(:,:,it(j)); Dragforcez(:,:,j)=Dragforcez(:,:,j)+-m_e*dragfreqsrc.*obj.fluidUZ(:,:,it(j)); end % Time derivative cdt=(obj.t2d(it(j)+1)-obj.t2d(it(j)-1)); durdt(:,:,j)=m_e*(obj.fluidUR(:,:,it(j)+1)-obj.fluidUR(:,:,it(j)-1))/cdt; duthetdt(:,:,j)=m_e*(obj.fluidUTHET(:,:,it(j)+1)-obj.fluidUTHET(:,:,it(j)-1))/cdt; duzdt(:,:,j)=m_e*(obj.fluidUZ(:,:,it(j)+1)-obj.fluidUZ(:,:,it(j)-1))/cdt; end if(~fdens) Forces.Eforcer=Eforcer; Forces.Eforcez=Eforcez; Forces.Bforcer=Bforcer; Forces.Bforcethet=Bforcethet; Forces.Bforcez=Bforcez; Forces.inertforcer=inertforcer; Forces.inertforcethet=inertforcethet; Forces.inertforcez=inertforcez; Forces.Pressforcer=Pressforcer; Forces.Pressforcethet=Pressforcethet; Forces.Pressforcez=Pressforcez; Forces.durdt=durdt; Forces.duthetdt=duthetdt; Forces.duzdt=duzdt; Forces.Dragforcer=Dragforcer; Forces.Dragforcethet=Dragforcethet; Forces.Dragforcez=Dragforcez; else % multiply by density to have force density Forces.Eforcer=Eforcer.*n; Forces.Eforcez=Eforcez.*n; Forces.Bforcer=Bforcer.*n; Forces.Bforcethet=Bforcethet.*n; Forces.Bforcez=Bforcez.*n; Forces.inertforcer=inertforcer.*n; Forces.inertforcethet=inertforcethet.*n; Forces.inertforcez=inertforcez.*n; Forces.Pressforcer=Pressforcer.*n; Forces.Pressforcethet=Pressforcethet.*n; Forces.Pressforcez=Pressforcez.*n; Forces.durdt=durdt.*n; Forces.duthetdt=duthetdt.*n; Forces.duzdt=duzdt.*n; Forces.Dragforcer=Dragforcer.*n; Forces.Dragforcethet=Dragforcethet.*n; Forces.Dragforcez=Dragforcez.*n; end end end function [lr,rb,lz,zb]= clouddims(obj,it,zpos,fracn) % clouddims return the cloud axial and radial limit at time it % and axial position zpos % fracn defines the fraction of the maximum density below which % we consider to have a vacuum if nargin<4 fracn=0.1; end % get the density n=obj.N(:,:,it); lr=cell(1,length(it)); lz=lr; rb=lr; zb=rb; for i=1:size(n,3) nthresh=fracn*max(max(n(:,:,i))); % find the points outside of the cloud outside=find(n(:,zpos,i)2) rmpos=outside(j); rppos=outside(j+1); lr{i}(k)=obj.rgrid(rppos-1)-obj.rgrid(rmpos+1); rb{i}(:,k)=[max(rmpos+1,1) min(rppos-1,sum(obj.nnr))]; k=k+1; end end maxgap=2; k=1; for I=rmpos+1:rppos-1 outside=find(n(I,:,i)maxgap) maxgap=zgap(j); zmpos=outside(j); zppos=outside(j+1); lz{i}(k)=obj.zgrid(zppos-1)-obj.zgrid(zmpos+1); zb{i}(:,k)=[max(zmpos+1,1) min(zppos-1,obj.nz)]; k=k+1; end end end end end %------------------------------------------ % Functions for plotting evolving quantities function displaysplbound(obj,ax,rescale) %displaysplbound display on axis ax the boundary of the %simulation domain and the Dirichlet and Neumann walls defined %with spline curves if nargin<2 ax=gca; end if nargin<3 rescale=1; end hold on for i=1:obj.spl_bound.nbsplines knots=obj.spl_bound.boundary(i).knots(1:end); coeffs=obj.spl_bound.boundary(i).coefs'*rescale; pp=spmak(knots,coeffs); sizec=size(coeffs,2); order=length(knots)-sizec; s=linspace(knots(order),knots(sizec+1),1000); fittedpos=fnval(pp,s); plot(fittedpos(1,:),fittedpos(2,:),'-') plot(coeffs(1,:),coeffs(2,:),'rx','markersize',14) end end function displayraddim(obj,it,zpos,fracn) %displayraddim display the evolution of the radial dimension of the cloud in %time to find if the cloud size get below a critical radial %size at which the ionisation is not sufficient to compensate %the losses % also plot the well radial dimensions in time if nargin<3 zpos=floor(length(obj.zgrid)/2); end if nargin<4 fracn=0.1; end [lr,rb,lz,zb]=obj.clouddims(it,zpos,fracn); t=obj.t2d(it); Lr=zeros(size(lr)); er=obj.Er(:,:,it); r_min=Lr; r_minpred=r_min; well_r=Lr; nb=Lr; for i=1:length(lr) if ~isempty(lr{i}) && ~isempty(lz{i}) [Lr(i),id]=max(lr{i}); rm=rb{i}(1,id); rp=rb{i}(2,id); nb(i)=mean(obj.N(rm:rp,zpos,it(i))); Lp=min(lz{i}); Lm=mean(lz{i}); rpos=rm:rp; vperp=-er(rpos,zpos,i)./obj.Bz(zpos,rpos)'; Ek=0.5*obj.me*vperp.^2/obj.qe; sigio=obj.sigio(Ek); sigd=obj.sigmela(Ek)+obj.sigmio(Ek)+sigio; omegap2=obj.qe^2*obj.N(rpos,zpos,it(i))/obj.eps_0/obj.me; omegac2=(obj.qe*obj.Bz(zpos,rpos)'/obj.me).^2; ur=er(rpos,zpos,i)*obj.qe./((omegap2-omegac2)*obj.me).*sigd.*vperp*obj.neutcol.neutdens; r_minpred(i)=mean(obj.N(rp,zpos,it(i))*Lp*ur./(nb(i)*obj.neutcol.neutdens*sigio.*vperp*Lm));%mean(1./(-1/obj.rgrid(rm)+obj.neutcol.neutdens*sigio.*vperp./ur*(Lm/Lp)*nb(i)/obj.N(rp,zpos,it(i)))); rpos=rp; vperp=-er(rpos,zpos,i)./obj.Bz(zpos,rpos)'; Ek=0.5*obj.me*vperp.^2/obj.qe; sigio=obj.sigio(Ek); ur=obj.fluidUR(rpos,zpos,it(i)); r_min(i)=max(obj.N(rp,zpos,it(i))*Lp*ur/(nb(i)*obj.neutcol.neutdens*sigio.*vperp*Lm),0);%max(mean(1./(-1/obj.rgrid(rm)+obj.neutcol.neutdens*sigio.*vperp./ur*(Lm/Lp)*nb(i)/obj.N(rp,zpos,it(i)))),0); nb(i)=nb(i)*Lm*2*pi*obj.rgrid(rm)*Lr(i); else Lr(i)=NaN; r_min(i)=NaN; r_minpred(i)=NaN; end potwell=obj.PotentialWell(it(i))'; outside=find(isnan(potwell(:,zpos))); gap=diff(outside); for j=1:length(gap) if(gap(j)>2) rmpos=outside(j)+1; rppos=outside(j+1)-1; well_r(i)=obj.rgrid(rppos)-obj.rgrid(rmpos); end end end f=figure('Name', sprintf('%s rlims B=%f phi=%f',obj.name,obj.B0, obj.phinorm*(obj.potout-obj.potinn))); plot(t,Lr,'displayname','\Deltar_{cloud}','linewidth',1.3) hold on plot(t,r_min,'displayname','\Deltar_{min} (u_r simu)','linewidth',1.3) plot(t,r_minpred,'displayname','\Deltar_{min} (u_r pred)','linewidth',1.3) plot(t,well_r,'displayname','\Deltar_{well}','linewidth',1.3) ylabel('\Delta r [m]') yyaxis right plot(t,nb,'--','displayname','N') legend('location','eastoutside') xlabel('t [s]') ylabel('N') set(gca,'fontsize',12) yyaxis left ylimits=ylim; %ylim([ylimits(1) 1.1*max(Lr)]) title(sprintf('cloud radial limits at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%s/%s_%d_rlims',obj.folder,obj.name,zpos),[15 10]); end function displaypsi(obj,deltat) %% plot the initial and final radial profile at position z=0 and show the normalized enveloppe function Psi % relevant for Davidson annular distribution function f=figure('Name', sprintf('%s Psi',obj.name)); f.Name= sprintf('%s Psi',obj.name); zpos=floor(length(obj.zgrid)/2); tinit=1; tend=length(obj.t2d); if iscell(deltat) deltat=cell2mat(deltat); end if(obj.R.nt<2) h0=obj.H0; p0=obj.P0; else h0=mean(H(obj,{1:obj.VR.nparts,obj.VR.nt,false})); p0=mean(P(obj,{1:obj.VR.nparts,obj.VR.nt,false})); end lw=1.5; Mirrorratio=(obj.Rcurv-1)/(obj.Rcurv+1); locpot=mean(obj.pot(:,zpos,tend-deltat:tend),3); psi=1+obj.qe*locpot(:)/h0-1/(2*obj.me*h0)*(p0./obj.rgrid+obj.qe*0.5*obj.B0.*(obj.rgrid-obj.width/pi*Mirrorratio*cos(2*pi*obj.zgrid(zpos)/obj.width)*besseli(1,2*pi*obj.rgrid/obj.width))).^2; locdens=mean(obj.N(:,zpos,tend-deltat:tend),3); [maxn,In]=max(locdens);%M.N(:,zpos,tinit)); plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) hold on plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(tend-deltat)*1e9,obj.t2d(tend)*1e9),'linewidth',lw) plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end)*maxn*obj.rgrid(In),'DisplayName','N=c*1/r','linewidth',lw) plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end).^2*maxn*obj.rgrid(In)^2,'DisplayName','N=c*1/r^2','linewidth',lw) plot(obj.rgrid(In-2:end),1./obj.rgrid(In-2:end).^4*maxn*obj.rgrid(In)^4,'DisplayName','N=c*1/r^4','linewidth',lw) xlabel('r [m]') ylabel('n_e [m^{-3}]') I=find(psi>0); if (length(I)>1) I=[I(1)-2; I(1)-1; I; I(end)+1; I(end)+2]; else I=obj.nnr(1):length(psi); end rq=linspace(obj.rgrid(max(I(1),1)),obj.rgrid(I(end)),500); psiinterp=interp1(obj.rgrid(I),psi(I),rq,'pchip'); zeroindices=find(diff(psiinterp>=0),2); maxpsiinterp=max(psiinterp); plot(rq,maxn*psiinterp/abs(maxpsiinterp),'Displayname','normalized \Psi [a.u.]','linewidth',lw) ylim([0 inf]) for i=1:length(zeroindices) border=plot([rq(zeroindices(i)) rq(zeroindices(i))],[0 obj.rgrid(In)/obj.rgrid(In-2)*maxn],'k--','linewidth',lw); set(get(get(border,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); end legend xlim([0 0.02]) grid on title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%sPsi',obj.name),[15 10]); end function f=displayrprofile(obj,t,zpos,init) %% plot the initial and final radial profile at the axial center of the simulation space % also plot the azimuthal fluid rotation frequency profile % t: time index considered % zpos: axial position index % init: initial time considered for comparison f=figure('Name', sprintf('%s Prof',obj.name)); if nargin < 3 || length(zpos)<1 zpos=floor(length(obj.zgrid)/2); end if nargin<4 init=false; end if(iscell(t)) t=cell2mat(t); end lw=1.5; if init tinit=t(1); t=t(2:end); end locdens=mean(obj.N(:,zpos,t),3); %inverse of radius Rinv=1./obj.rgrid; Rinv(isnan(Rinv))=0; %azimuthal velocity and azimuthal rotation frequency in m/s and %1/s vthet=mean(obj.fluidUTHET(:,zpos,t),3); omegare=(vthet.*Rinv); % plot the initial density if(init) plot(obj.rgrid,obj.N(:,zpos,tinit),'bx-','DisplayName',sprintf('t=%1.2f[ns]',obj.t2d(tinit)*1e9),'linewidth',lw) end hold on %plot the time averaged current density plot(obj.rgrid,locdens,'rx-','DisplayName',sprintf('t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(t(1))*1e9,obj.t2d(t(end))*1e9),'linewidth',lw) xlabel('r [m]') ylabel('n_e [m^{-3}]') legend('location','Northwest') % limit the axis to the simulation domain if obj.conformgeom xlim([obj.rgrid(1) obj.rgrid(sum(obj.nnr(1:2)))]) else xlim([obj.rgrid(1) obj.rgrid(end)]) end grid on ylimits=ylim(); % plot the metallic walls for a constant radius coaxial % configuration if obj.conformgeom plot(obj.rgrid(1)*[1 1],ylimits,'k--') plot(obj.rgrid(end)*[1 1],ylimits,'k--') else plot(obj.r_a*[1 1],ylimits,'k--') if obj.walltype==0 plot(obj.r_b*[1 1],ylimits,'k--') elseif obj.walltype==1 rmax=obj.r_0-obj.r_r*sqrt(1-(obj.zgrid(zpos)-obj.z_0)^2/obj.z_r^2); plot(rmax*[1 1],ylimits,'k--') end end yyaxis right % plot the azimuthal fluid rotation frequency profile plot(obj.rgrid,omegare,'DisplayName',sprintf('<\\omega_{re}> t=[%1.2f-%1.2f] [ns] averaged',obj.t2d(t(1))*1e9,obj.t2d(t(end))*1e9),'linewidth',lw) ylabel('\omega_{re} [1/s]') title(sprintf('Radial density profile at z=%1.2e[m]',obj.zgrid(zpos))) obj.savegraph(f,sprintf('%srProf',obj.name),[15 10]); end function displayenergy(obj) %% Plot the time evolution of the system energy and number of simulated macro particles tmin=1; tmax=length(obj.ekin); f=figure('Name', sprintf('%s Energy',obj.name)); subplot(2,1,1) plot(obj.t0d(tmin:tmax),obj.ekin(tmin:tmax),'o-',... obj.t0d(tmin:tmax),obj.epot(tmin:tmax),'d-',... obj.t0d(tmin:tmax),obj.etot(tmin:tmax),'h-',... obj.t0d(tmin:tmax),obj.etot0(tmin:tmax),'h-',... obj.t0d(tmin:tmax),obj.eerr(tmin:tmax),'x--') %obj.t0d(tmin:tmax),obj.ekin(tmin:tmax)-obj.epot(tmin:tmax),'--', legend('E_{kin}', 'E_{pot}', 'E_{tot}','E_{ref}','E_{err}') xlabel('Time [s]') ylabel('Energies [J]') grid on xlimits=xlim(); subplot(2,1,2) try semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'h-') catch semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)/obj.etot(2)),'h-') end hold on xlabel('t [s]') ylabel('E_{err}/E_{tot}') xlim(xlimits) grid on try yyaxis right plot(obj.t0d(tmin:tmax),abs(obj.npart(tmin:tmax)./obj.npart(1)*100),'d--') ylabel('Nparts %') %ylim([0 110]) catch end ylimits=ylim; for i=1:length(obj.restarttimes) plot(obj.restarttimes(i)*[1 1],ylimits,'k--') end obj.savegraph(f,sprintf('%s/%sEnergy',obj.folder,obj.name)); end function f=displaycharge(obj,f,linelegend) %% Plot the time evolution of the system charge of electrons % f: figure handle if you want to stack several such curves % on the same figure % linelegend: legend for this charge evolution tmin=1; tmax=length(obj.ekin); if nargin<2 f=figure('Name', sprintf('%s Charge',obj.name)); end if nargin<3 linelegend=''; end ax=f.CurrentAxes; if isempty(ax) ax=axes(f); end try plot(ax,obj.t0d(tmin:tmax),abs(obj.npart(tmin:tmax)*obj.qsim),'linewidth',1.5,'displayname',linelegend) hold on ylabel(ax,'Total charge [C]') xlabel(ax,'t [s]') grid on if(nargin>2) legend end set(ax,'fontsize',12) catch end if nargin < 2 obj.savegraph(f,sprintf('%s/%scharge',obj.folder,obj.name)); end end function displaySimParticles(obj) %% Plot the time evolution of the number of simulated markers in the main specie f=figure('Name', sprintf('%s Trapped particles',obj.name)); plot(obj.t0d,obj.npart,'linewidth',1.5) xlabel('t [s]') ylabel('N particles') set(gca,'fontsize',12) obj.savegraph(f,sprintf('%s/%sntrapped',obj.folder,obj.name),[10 12]); end function displayLarmorRad(obj,time2d) if nargin<2 time2d=length(obj.t2d); end % Plot the larmor radius for created particles with low energy % the larmor radius is calculated by considering that the % initial perpendicular velocity \approx the ExB velocity if time2d>0 Er=obj.Er(:,:,time2d); Ez=obj.Ez(:,:,time2d); else Er=obj.Erxt(:,:,1); Ez=obj.Ezxt(:,:,1); end rl=abs(obj.me/obj.qe*(-Er.*obj.Bz'+Ez.*obj.Br')./(obj.B.^3)'); figure rl(obj.geomweight(:,:,1)<0)=0; contourf(obj.zgrid,obj.rgrid,rl) hold on contour(obj.zgrid,obj.rgrid,obj.geomweight(:,:,1),[0 0],'r-','linewidth',3) if time2d>0 n=obj.N(:,:,time2d); maxN=max(n (:)); n=n/maxN*mean(rl(:)); contour(obj.zgrid,obj.rgrid,n,linspace(0,1,6)*mean(rl(:)),'r:','linewidth',3) end c=colorbar; xlabel('z [m]') ylabel('r [m]') c.Label.String='r_L [m]'; end function displayHP(obj,tstart) % Plot the histogramm of the total energy and canonical angular momentum at time tstart and % end time of the simulation over the full simulation space for the main specie if(iscell(tstart)) tstart=cell2mat(tstart); end if(obj.R.nt>=2) tstart=obj.R.nt; f=figure('Name', sprintf('%s HP',obj.name)); legtext=sprintf("t=%2.1f - %2.1f [ns]",obj.tpart(tstart)*1e9,obj.tpart(end)*1e9); subplot(1,2,1) partsmax=min(obj.nbparts(end),obj.R.nparts); Hloc=H(obj,{1:obj.nbparts(1),1,false}); h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9)); hold on Hloc=H(obj,{1:partsmax,obj.R.nt,false}); %,'Binwidth',h1.BinWidth h1=histogram(Hloc,20,'BinLimits',[min(Hloc(:)) max(Hloc(:))],'DisplayName',legtext); ylabel('counts') xlabel('H [J]') legend subplot(1,2,2) Ploc=P(obj,{1:obj.nbparts(1),1,false}); h2=histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(1)*1e9)); hold on Ploc=P(obj,{1:partsmax,obj.R.nt,false}); histogram(Ploc,50,'BinLimits',[min(Ploc(:)) max(Ploc(:))],'DisplayName',legtext); ylabel('counts') xlabel('P [kg\cdotm^2\cdots^{-1}]') %clear P %clear H legend %xlim([0.95*h2.BinLimits(1) 1.05*h2.BinLimits(2)]) obj.savegraph(f,sprintf('%s/%sParts_HP',obj.folder,obj.name)); end end function displayaveragetemp(obj) % Computes and show the particles average temperature as a function of time f=figure('Name',sprintf('%s potinn=%f part temperature',obj.name,obj.potinn*obj.phinorm)); vr2=obj.VR(:,:,false); vr2=mean(vr2.^2,1)-mean(vr2,1).^2; vz2=obj.VZ(:,:,false); vz2=mean(vz2.^2,1)-mean(vz2,1).^2; vthet2=obj.VTHET(:,:,false); vthet2=mean(vthet2.^2,1)-mean(vthet2,1).^2; plot(obj.tpart,0.5*obj.me*vr2/obj.qe,'displayname','T_r') hold on plot(obj.tpart,0.5*obj.me*vz2/obj.qe,'displayname','T_z') plot(obj.tpart,0.5*obj.me*vthet2/obj.qe,'displayname','T_{thet}') xlabel('time [s]') ylabel('T [eV]') title(sprintf('\\phi_a=%.1f kV \\phi_b=%.1f kV R=%.1f',obj.potinn*obj.phinorm/1e3,obj.potout*obj.phinorm/1e3,obj.Rcurv)) legend grid obj.savegraph(f,sprintf('%s/%s_partstemp',obj.folder,obj.name)); end function displayCurrentsevol(obj,timesteps) % Computes and display the time evolution of the outgoing currents on each domain boundary % at timesteps timesteps if nargin<2 timesteps=1:length(obj.t2d); end currents=obj.OutCurrents(timesteps); f=figure('Name',sprintf('%s Currents',obj.name)); if(obj.B(1,1)>obj.B(end,1)) lname='HFS'; rname='LFS'; else lname='LFS'; rname='HFS'; end plot(obj.t2d(timesteps),currents(1,:),'Displayname',lname,'linewidth',1.8); hold on plot(obj.t2d(timesteps),currents(2,:),'Displayname',rname,'linewidth',1.8); plot(obj.t2d(timesteps),currents(3,:),'Displayname','outer cylinder','linewidth',1.8); plot(obj.t2d(timesteps),currents(4,:),'Displayname','inner cylinder','linewidth',1.8); if size(currents,1)>=5 plot(obj.t2d(timesteps),currents(5,:),'Displayname','ellipse','linewidth',1.8); end legend('location','Northeast') xlabel('time [s]') ylabel('I [A]') grid on set(gca,'fontsize',12) title(sprintf('\\phi_b-\\phi_a=%.2g kV, R=%.1f',(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv)) obj.savegraph(f,sprintf('%s/%s_outCurrents',obj.folder,obj.name),[16 12]); end function displayChargeLossevol(obj,timesteps,toptitle,scalet,dens) % Computes and display the time evolution of the outgoing currents on each domain boundary % at time obj.t2d(timesteps) %scalet=true scales the time by the ellastic collision %frequency %dens = true plot the time evolution of the maximum electron %density in the simulation domain otherwise plot the total %number of electrons in the domain if nargin<2 timesteps=1:length(obj.t2d); end if nargin<4 scalet=true; end if nargin <5 dens=true; end if scalet if obj.neutcol.present vexb0=(obj.Ez(:,:,1).*obj.Br'-obj.Er(:,:,1).*obj.Bz')./(obj.B'.^2); vexb0(obj.geomweight(:,:,1)<=0)=0; E=0.5*obj.msim/obj.weight*mean(abs(vexb0(:)))^2/obj.qe; taucol=1/(obj.neutcol.neutdens*mean(abs(vexb0(:)))*(obj.sigio(E)+obj.sigmela(E)+obj.sigmio(E))); try Sio_S=1e17*(obj.neutcol.neutdens*mean(abs(vexb0(:)))*obj.sigio(E))/(obj.maxwellsrce.frequency*obj.weight/(pi*(obj.maxwellsrce.rlim(2)^2-obj.maxwellsrce.rlim(1)^2)*diff(obj.maxwellsrce.zlim))) catch end tlabel='t/\tau_d [-]'; else taucol=2*pi/obj.omece; tlabel='t/\tau_ce [-]'; end else taucol=1e-9; tlabel='t [ns]'; end if dens N=obj.N(:,:,timesteps); geomw=obj.geomweight(:,:,1); geomw(geomw<0)=0; geomw(geomw>0)=1; N=N.*geomw; nmax=squeeze(max(max(N,[],1),[],2)); tn=(obj.t2d(timesteps)); nlabel='n_{e,max} [m^{-3}]'; ndlabel='n_{e,max}'; else t0dst=find(obj.t0d>=obj.t2d(timesteps(1)),1,'first'); t0dlst=find(obj.t0d<=obj.t2d(timesteps(end)),1,'last'); tn=obj.t0d(t0dst:t0dlst); nmax=obj.npart(t0dst:t0dlst)*obj.weight; nlabel='Nb e^-'; ndlabel='Nb e^-'; end [currents,pos]=obj.OutCurrents(timesteps); P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar currents=currents/P; f=figure('Name',sprintf('%s Charges',obj.name)); % Plot the evolution of nb of particles yyaxis right p=plot(tn/taucol,nmax,'b-.','linewidth',1.8,'Displayname',ndlabel); ylabel(nlabel) ax=gca; ax.YAxis(2).Color=p.Color; ylim([0 inf]) if(obj.B(1,1)>obj.B(end,1)) lname='HFS'; rname='LFS'; else lname='LFS'; rname='HFS'; end yyaxis left mincurr=max(currents(:))*5e-3; if (max(currents(1,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(1,:),'r:','Displayname',lname,'linewidth',1.8); end hold on if (max(currents(2,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(2,:),'r--','Displayname',rname,'linewidth',1.8); end if (max(currents(3,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(3,:),'r-','Displayname','outer cylinder','linewidth',1.8); end if (max(currents(4,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(4,:),'Displayname','inner cylinder','linewidth',1.8); end if (size(currents,1)>=5 && max(currents(5,:)>mincurr)) plot(obj.t2d(timesteps)/taucol,currents(5,:),'r-','Displayname','ellipse','linewidth',1.8); end xlabel(tlabel) ylabel('I/p_n [A/mbar]') grid on set(gca,'fontsize',12) ax.YAxis(1).Color='red'; legend('Orientation','horizontal','location','south','numcolumns',3) if nargin <3 title(sprintf('\\phi_b-\\phi_a=%.2g kV, B=%f T',(obj.potout-obj.potinn)*obj.phinorm/1e3,max(obj.B(:)))) elseif ~isempty(toptitle) title(toptitle) end obj.savegraph(f,sprintf('%s/%s_ChargeEvol%i%i',obj.folder,obj.name,scalet,dens),[16 14]); end function display1Dpotentialwell(obj,timestep,rpos) % Display the potential well along the magentic field line % passing by rgrid(rpos) at the center of the simulation space if iscell(timestep) timestep=cell2mat(timestep); end f=figure('Name',sprintf('%s 1D Potential well',obj.name)); model=obj.potentialwellmodel(timestep); z=model.z; r=model.r; Pot=model.pot; rathet=model.rathet; if (mod(rpos, 1) ~= 0) [~,rpos]=min(abs(M.rgrid-rpos)); end crpos=obj.rgrid(rpos); id=find(timestep==0); timestep(id)=[]; n=obj.N(:,:,timestep); if(~isempty(timestep==0)) N0=zeros(obj.N.nr+1,obj.N.nz+1); n=cat(3,n(:,:,1:id-1),N0,n(:,:,id:end)); end n=mean(n,3); linepot=zeros(length(obj.zgrid),length(timestep)); rathetpos=obj.rAthet(rpos,ceil(length(obj.zgrid)/2)); F=scatteredInterpolant(z',rathet',Pot(:,1)); for i=1:length(timestep) F=scatteredInterpolant(z',rathet',Pot(:,i)); linepot(:,i)=F(obj.zgrid,rathetpos*ones(size(obj.zgrid))); %linepot(:,i)=griddata(z,rathet,pot(:,i),obj.zgrid,rathetpos); end linepot=mean(linepot,2); [Zinit,~]=meshgrid(obj.zgrid,obj.rAthet(:,1)); n=griddata(Zinit,obj.rAthet,n,obj.zgrid,rathetpos); plot(obj.zgrid,linepot) ylabel('Potentiel [eV]') xlabel('z [m]') xlim([obj.zgrid(1) obj.zgrid(end)]) hold(gca, 'on') yyaxis right plot(obj.zgrid,n) ylabel('n [m^{-3}]') if length(timestep)==1 title(sprintf('Potential well t=%1.2f [ns] r=%1.2f [mm]',obj.t2d(timestep)*1e9,1e3*crpos)) else title(sprintf('Potential well t=[%1.2f-%1.2f] [ns] r=%1.2f [mm]',obj.t2d(timestep(1))*1e9,obj.t2d(timestep(end))*1e9,1e3*crpos)) end obj.savegraph(f,sprintf('%s/%s_well1Dr_%d',obj.folder,obj.name,rpos)); end function displayVdistribRThetZ(obj,timestep, rpos, zpos) %displayVdistribRThetZ plot the velocity distribution function % in m/s %extracted from the markers at position window from rpos(1) %rpos(end) and zpos(1) to zpos(end) % and at time obj.tpart(timestep) %rpos and zpos are given as grid indices if(obj.R.nt>=2) if nargin<2 timesteppart=length(obj.tpart); else timesteppart=timestep; end if nargin<3 || isempty(rpos) rpos=1:length(obj.rgrid); rspan=[obj.rgrid(1) obj.rgrid(end)]; else r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)]; rspan=[r(rpos) r(rpos+1)]; end if nargin<4 || isempty(zpos) zpos=1:length(obj.zgrid); zspan=[obj.zgrid(1) obj.zgrid(end)]; else z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)]; zspan=[z(zpos) z(zpos+1)]; end nbp=min(obj.nbparts(1),obj.R.nparts); R=obj.R(1:nbp,1,false); Z=obj.Z(1:nbp,1,false); Vr=obj.VR(1:nbp,1,false); Vz=obj.VZ(1:nbp,1,false); Vthet=obj.VTHET(1:nbp,1,false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Vr=Vr(ids); Vz=Vz(ids); Vthet=Vthet(ids); vTr=std(Vr,1); vTz=std(Vz,1); vTthet=std(Vthet,1); nbp=min(obj.nbparts(timesteppart),obj.R.nparts); Rend=obj.R(1:nbp,timesteppart,false); Zend=obj.Z(1:nbp,timesteppart,false); Vrend=obj.VR(1:nbp,timesteppart,false); Vzend=obj.VZ(1:nbp,timesteppart,false); Vthetend=obj.VTHET(1:nbp,timesteppart,false); ids=Rend>=rspan(1) & Rend<=rspan(2) & Zend>=zspan(1) & Zend<=zspan(2); nbtot=sum(ids) Vrend=Vrend(ids); Vzend=Vzend(ids); Vthetend=Vthetend(ids); vTrend=std(Vrend,1); vTzend=std(Vzend,1); vTthetend=std(Vthetend,1); binwidth=abs(max(Vrend)-min(Vrend))/sqrt(length(Vrend)); f=figure('Name',sprintf("%s vrz distrib",obj.file)); [~,time2did]=min(abs(obj.t2d-obj.tpart(timestep))); subplot(1,4,1); obj.dispV(Vr,Vrend,'V_r [m/s]',[1,timesteppart]) [~,time2did]=min(abs(obj.t2d-obj.tpart(timestep))); if length(rpos)==1 vexb=-obj.Er(rpos,zpos,time2did)/obj.Bz(zpos,rpos)'; vexb=mean(vexb(:)); if ~isempty(obj.neutcol.ela_cross_sec) % plot the radial drift velocity as nu_dE_r/(B\Omega_c) vdr=obj.neutcol.neutdens*obj.sigmela(vexb^2*obj.me*0.5/obj.qe)*vexb*-obj.Er(rpos,zpos,time2did)... ./(obj.B(zpos,rpos)'.*obj.B(zpos,rpos)'*obj.qe/obj.me); vdr=mean(vdr(:)); ylimits=ylim; plot(vdr*[1 1],ylimits,'k--','displayname',sprintf('V_{d,pred}=%1.2g [m/s]',vdr)) end end subplot(1,4,2); obj.dispV(Vthet,Vthetend,'V\theta [m/s]',[1,timesteppart]) hold on drawnow ylimits=ylim; if length(rpos)==1 if ~isempty(obj.Erxt) vexbext=-obj.Erxt(rpos,zpos)/obj.Bz(zpos,rpos)'; plot(vexbext*[1 1],ylimits,'k--','displayname',sprintf('V_{ExB,ext}=%1.2g [m/s]',vexbext)) end plot(vexb*[1 1],ylimits,'k-.','displayname',sprintf('V_{ExB,tot}=%1.2g [m/s]',vexb)) end subplot(1,4,3); obj.dispV(Vz,Vzend,'Vz [m/s]',[1,timesteppart]) subplot(1,4,4); obj.dispV(sqrt(Vr.^2+(Vthet).^2+Vz.^2),sqrt(Vrend.^2+(Vthetend).^2+Vzend.^2),'Vtot [m/s]',[1,timesteppart],'maxwell') sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.tpart(timestep)*1e9, rspan*1e3, zspan*1e3)) obj.savegraph(f,sprintf('%s/%sParts_V_RZ',obj.folder,obj.name),[25 14]); end end function displayEkin(obj,timestep, rpos, zpos) %displayEkin plot the kinetic energy distribution function in %eV %extracted from the markers at position window from rpos(1) %rpos(end) and zpos(1) to zpos(end) % and at time obj.tpart(timestep) %rpos and zpos are given as grid indices if(obj.R.nt>=2) if nargin<2 timesteppart=[1 length(obj.tpart)]; else if length(timestep)<2 timesteppart=[1 timestep]; else timesteppart=[timestep(1) timestep(end)]; end end if nargin<3 || isempty(rpos) rspan=[obj.rgrid(1) obj.rgrid(end)]; else r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)]; rspan=[r(rpos) r(rpos+1)]; end if nargin<4 || isempty(zpos) zspan=[obj.zgrid(1) obj.zgrid(end)]; else z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)]; zspan=[z(zpos) z(zpos+1)]; end nbp=min(obj.nbparts(timesteppart(1)),obj.R.nparts); R=obj.R(1:nbp,timesteppart(1),false); Z=obj.Z(1:nbp,timesteppart(1),false); Ekin=obj.Ekin(1:nbp,timesteppart(1),false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Ekin=Ekin(ids)/obj.qe; nbp=min(obj.nbparts(timesteppart(2)),obj.R.nparts); Rend=obj.R(1:nbp,timesteppart(2),false); Zend=obj.Z(1:nbp,timesteppart(2),false); Ekinend=obj.Ekin(1:nbp,timesteppart(2),false); ids=Rend>=rspan(1) & Rend<=rspan(2) & Zend>=zspan(1) & Zend<=zspan(2); Ekinend=Ekinend(ids)/obj.qe; f=figure('Name',sprintf("%s E_k distrib",obj.file)); obj.dispV(Ekin,Ekinend,'E_k',timesteppart,'maxwell') sgtitle(sprintf('dt=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.dt*1e9, rspan*1e3, zspan*1e3)) obj.savegraph(f,sprintf('%s/%sParts_E_kin',obj.folder,obj.name)); end end function displayVdistribParPer(obj,timestep, rpos, zpos, gcs) %displayVdistribParPer plot the velocity distribution function % in m/s for the parallel and perpendicular velocity %extracted from the markers at position window from rpos(1) %rpos(end) and zpos(1) to zpos(end) % and at time obj.tpart(timestep) %rpos and zpos are given as grid indices % gcs define if you give the perpendicular velocity in the % guiding center frame or in the laboratory frame if(obj.R.nt>=2) if nargin<2 timesteppart=length(obj.tpart); else timesteppart=timestep; end if nargin<3 || isempty(rpos) rspan=[obj.rgrid(1) obj.rgrid(end)]; else r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)]; rspan=[r(rpos) r(rpos+1)]; end if nargin<4 || isempty(zpos) zspan=[obj.zgrid(1) obj.zgrid(end)]; else z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)]; zspan=[z(zpos) z(zpos+1)]; end if nargin<5 gcs=false; % define if we look in the guiding center system end nbp=min(obj.nbparts(1),obj.R.nparts); R=obj.R(1:nbp,1,false); Z=obj.Z(1:nbp,1,false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Vperp=obj.Vperp(1:nbp,1,false,gcs); Vpar=obj.Vpar(1:nbp,1,false); Vperp=Vperp(ids); Vpar=Vpar(ids); nbp=min(obj.nbparts(timesteppart),obj.R.nparts); R=obj.R(1:nbp,timesteppart,false); Z=obj.Z(1:nbp,timesteppart,false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Vperpend=obj.Vperp(1:nbp,timesteppart,false,gcs); Vparend=obj.Vpar(1:nbp,timesteppart,false); Vperpend=Vperpend(ids); Vparend=Vparend(ids); %binwidth=abs(max(Vparend)-min(Vparend))/500; f=figure('Name',sprintf("%s v parper distrib",obj.file)); subplot(1,2,1) if gcs lgd='v_\perp gcs [m/s]'; else lgd='v_\perp [m/s]'; end obj.dispV(Vperp,Vperpend,lgd,[1,timesteppart], 'maxwell') subplot(1,2,2) obj.dispV(abs(Vpar),abs(Vparend),'v_{par} [m/s]',[1,timesteppart],'None') sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm]',obj.tpart(timestep)*1e9, rspan*1e3, zspan*1e3)) obj.savegraph(f,sprintf('%s/%sParts_V_parper',obj.folder,obj.name)); if gcs obj.savegraph(f,sprintf('%s/%sParts_V_parpergcs',obj.folder,obj.name)); else obj.savegraph(f,sprintf('%s/%sParts_V_parper',obj.folder,obj.name)); end end end function display2DVdistrib(obj,timestep, rpos, zpos, gcs) %display2DVdistrib plot the velocity distribution function % in m/s for the parallel and perpendicular velocity % and for the radial azimuthal velocity % as a 2D contour plot the show the velocity phase space distribution %extracted from the markers at position window from rpos(1) %rpos(end) and zpos(1) to zpos(end) % and at time obj.tpart(timestep) %rpos and zpos are given as grid indices % gcs define if you give the perpendicular velocity in the % guiding center frame or in the laboratory frame if(obj.R.nt>=2) if nargin<2 timesteppart=length(obj.tpart); else timesteppart=timestep; end if nargin<3 || isempty(rpos) rspan=[obj.rgrid(1) obj.rgrid(end)]; else r=[obj.rgrid(1);(obj.rgrid(1:end-1)+obj.rgrid(2:end))*0.5;obj.rgrid(end)]; rspan=[r(rpos(1)) r(rpos(end)+1)]; end if nargin<4 || isempty(zpos) zspan=[obj.zgrid(1) obj.zgrid(end)]; else z=[obj.zgrid(1);(obj.zgrid(1:end-1)+obj.zgrid(2:end))*0.5;obj.zgrid(end)]; zspan=[z(zpos(1)) z(zpos(end)+1)]; end if nargin<5 gcs=false; % define if we look in the guiding center system end nbp=min(obj.nbparts(timesteppart),obj.R.nparts); R=obj.R(1:nbp,timesteppart,false); Z=obj.Z(1:nbp,timesteppart,false); ids=R>=rspan(1) & R<=rspan(2) & Z>=zspan(1) & Z<=zspan(2); Vperp=obj.Vperp(1:nbp,timesteppart,false,gcs); Vpar=obj.Vpar(1:nbp,timesteppart,false); Vr=obj.VR(1:nbp,timesteppart,false); Vthet=obj.VTHET(1:nbp,timesteppart,false); Vper=Vperp(ids); Vpar=Vpar(ids); Vr=Vr(ids); Vthet=Vthet(ids); nbp=sum(ids(:)); f=figure('Name',sprintf("%s v parper distrib",obj.file)); subplot(2,1,1) [N,Xedges,Yedges] = histcounts2(Vpar,Vper,20); Xedges=(Xedges(1:end-1)+Xedges(2:end))/2; Yedges=(Yedges(1:end-1)+Yedges(2:end))/2; contourf(Xedges,Yedges,N') xlabel('v_{par} [m/s]') ylabel('v_{\perp} [m/s]') c=colorbar; c.Label.String='Counts'; subplot(2,1,2) [N,Xedges,Yedges] = histcounts2(Vthet,Vr,20); Xedges=(Xedges(1:end-1)+Xedges(2:end))/2; Yedges=(Yedges(1:end-1)+Yedges(2:end))/2; contourf(Xedges,Yedges,N') %histogram2(Vthet,Vr,'displaystyle','tile','binmethod','auto') %scatter(Vthet,Vr) xlabel('v_\theta [m/s]') ylabel('v_r [m/s]') c=colorbar; c.Label.String='Counts'; sgtitle(sprintf('t=%1.2e[ns] r=[%2.1f, %2.1f] [mm] z=[%2.1f, %2.1f] [mm] N=%3i',mean(obj.tpart(timestep))*1e9, rspan*1e3, zspan*1e3,nbp)) mkdir(sprintf('%s/vdist',obj.folder)) if gcs obj.savegraph(f,sprintf('%s/vdist/%sParts_V_2dparpergcs_r%iz%it%i',obj.folder,obj.name,floor(mean(rpos)),floor(mean(zpos)),floor(mean(timestep)))); else obj.savegraph(f,sprintf('%s/vdist/%sParts_V_2dparper_r%iz%it%i',obj.folder,obj.name,floor(mean(rpos)),floor(mean(zpos)),floor(mean(timestep)))); end end end function [p, maxnb, c]=displayPhaseSpace(obj,type,partsstep, Rindex, Zindex,legendtext, figtitle, f, maxnb, c, gcs) if nargin<8 f=figure; f=gca; end if nargin<7 figtitle=sprintf('r=%1.2f [mm] z=%1.2f [mm] \\Delta\\phi=%1.1f[kV] R=%1.1f',obj.rgrid(Rindex)*1e3,obj.zgrid(Zindex)*1e3,(obj.potout-obj.potinn)*obj.phinorm/1e3,obj.Rcurv); end if nargin <6 legendtext=sprintf('t=%1.3g [s]',obj.tpart(partsstep)); end fieldstep=find(obj.tpart(partsstep(end))==obj.t2d,1); if nargin>=10 ctemp=c; n=zeros(length(c{1}),length(c{2})); else nbins=15; n=zeros(nbins); end if nargin <11 gcs=true; end for i=1:length(partsstep) odstep=find(obj.tpart(partsstep(i))==obj.t0d); nbp=min(obj.R.nparts,obj.nbparts(partsstep(i))); Rp=obj.R(1:nbp,partsstep(i),false); Zp=obj.Z(1:nbp,partsstep(i),false); deltar=obj.dr(2)/2; deltarm=obj.rgrid(Rindex)-sqrt(obj.rgrid(Rindex)^2-deltar^2-2*obj.rgrid(Rindex)*deltar); deltaz=obj.dz/2; Indices=Rp>=obj.rgrid(Rindex)-deltarm & Rp=obj.zgrid(Zindex)-deltaz & Zp0)),max(Blines(obj.geomweight(:,:,1)>0)),20); Blines(obj.geomweight(:,:,1)<0)=NaN; [~,h1]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'-.','color','k','linewidth',1.5,'Displayname','Magnetic field lines'); % Draw the metallic boundaries and the geometry itself [c1,hContour]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,-geomw,[0,0],'linewidth',1.5); drawnow; xlim(ax1,[obj.zgrid(1)*1000 obj.zgrid(end)*1000]) % Change the color of the metallic boundaries to grey hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' try hFills(end).ColorData = uint8([150;150;150;255]); for idx = 1 : numel(hFills)-1 hFills(idx).ColorData(4) = 0; % default=255 end catch end grid on; hold on; f.PaperOrientation='landscape'; f.PaperUnits='centimeters'; papsize=[16 14]; f.PaperSize=papsize; set(ax1,'fontsize',14) %axis equal obj.savegraph(f,sprintf('%sfluid_dens',obj.name)) end function displaymagfield(obj) %displaymagfield display the magnetic field lines and the %amplitude of the magnetic field using a contour % also show the domain boundaries B=obj.B'; f=figure('Name', sprintf('%s B field',obj.name)); B(obj.geomweight(:,:,1)<0)=NaN; ax1=gca; title(sprintf('Configuration')) h=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,B,linspace(min(B(:)),max(B(:)),50),'Displayname','B [T]', 'linestyle','none'); hold on; %% Magnetic field lines Blines=obj.rAthet; levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),50); Blines(obj.geomweight(:,:,1)<0)=NaN; [~,h1]=contour(ax1,obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'r-.','linewidth',1.5,'Displayname','Magnetic field lines'); colormap(ax1,'parula') % Grey outline geomw=obj.geomweight(:,:,1); geomw(geomw>0)=NaN; geomw(geomw<0)=min(B(:)); [c1,hContour]=contourf(ax1,obj.zgrid*1000,obj.rgrid*1000,geomw, [0 0]); drawnow; xlim(ax1,[obj.zgrid(1)*1000 obj.zgrid(end)*1000]) if(obj.conformgeom) ylim([ax1 ],[obj.rgrid(1)*1000 obj.rgrid(rgridend)*1000]) else ylim([ax1],[obj.rgrid(1)*1000 obj.rgrid(end)*1000]) end legend([h1],{'Magnetic field lines'},'location','northwest') xlabel(ax1,'z [mm]') ylabel(ax1,'r [mm]') %title(ax1,sprintf('Density t=[%1.2g-%1.2g]s n_e=%1.2gm^{-3}',M.t2d(fieldstart),M.t2d(fieldend),double(maxdens))) c = colorbar(ax1); c.Label.String= 'B [T]'; view(ax1,2) %set(h,'edgecolor','none'); grid on; hFills=hContour.FacePrims; [hFills.ColorType] = deal('truecoloralpha'); % default = 'truecolor' %caxis([min(B(:)) max(B(:))]) try hFills(1).ColorData = uint8([150;150;150;255]); for idx = 2 : numel(hFills) hFills(idx).ColorData(4) = 0; % default=255 end catch end [~, name, ~] = fileparts(obj.file); if( obj.walltype >=2 && obj.walltype<=4) rectangle('Position',[obj.zgrid(1) obj.r_b obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none') ylimits=ylim; ylim([ylimits(1),ylimits(2)+1]) end if(isempty(obj.spl_bound)) rectangle('Position',[obj.zgrid(1) obj.r_a-0.001 obj.zgrid(end)-obj.zgrid(1) 0.001]*1e3,'FaceColor',[150 150 150]/255,'Edgecolor','none') ylimits=ylim; ylim([ylimits(1)-1,ylimits(2)]) end f.PaperUnits='centimeters'; %axis equal papsize=[14 5 ]; pos=f.Position; pos(3)=floor(1.3*pos(3)); f.Position=pos; obj.savegraph(f,sprintf('%s/%s_Bfield',obj.folder,obj.name),papsize); end function displaySurfFlux(obj,timestep, subdiv) %displaySurfFlux plot the current densities %on the domain boundaries for time t2d(timestep) %directly on the boundaries themselves %make it easier to see where the currents are collected if nargin<3 subdiv=1; end mflux= obj.Metallicflux(timestep,subdiv); lflux= -squeeze(obj.Axialflux(timestep,1))'; rflux= squeeze(obj.Axialflux(timestep,length(obj.zgrid)))'; time=obj.t2d(timestep); if nargin<3 ids=1:length(mflux); end %% P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar f=figure('name','fluxevol'); linew=5; %obj.displaysplbound(gca,1e3); contour(obj.zgrid*1e3,obj.rgrid*1e3,obj.geomweight(:,:,1),[0 0],'b-','linewidth',1.5); hold on for i=1:length(mflux.p) x=mflux.p{i}(1,:)*1000; y=mflux.p{i}(2,:)*1000; y(end)=NaN; c=mean(mflux.gamma{i},2)'*obj.qe/(100^2)/P; c(c<=0)=NaN; patch(x,y,c,'EdgeColor','interp','LineWidth',linew); hold on end x=obj.zgrid(1)*ones(size(obj.rgrid))*1000; y=obj.rgrid*1000; y(end)=NaN; c=mean(lflux,1)*obj.qe/(100^2)/P; c(c<=0)=NaN; patch(x,y,c,'EdgeColor','interp','LineWidth',linew); x=obj.zgrid(end)*ones(size(obj.rgrid))*1e3; y=obj.rgrid*1000; y(end)=NaN; c=mean(rflux,1)*obj.qe/(100^2)/P; c(c<=0)=NaN; patch(x,y,c,'EdgeColor','interp','LineWidth',linew); title(sprintf('t=%4.2f [ns]',mean(time)*1e9)) c=colorbar; c.Label.String= 'j\cdotn [A/(cm^2 mbar)]'; xlabel('z [mm]') ylabel('r [mm]') colormap(jet) set(gca,'colorscale','log') %% Magnetic field lines Blines=obj.rAthet; levels=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),20); Blines(obj.geomweight(:,:,1)<0)=NaN; [~,h1]=contour(obj.zgrid*1000,obj.rgrid*1000,Blines,real(levels),'m-.','linewidth',1.5,'Displayname','Magnetic field lines'); %axis equal obj.savegraph(f,sprintf('%s/%s_surfFlux_it2d_%i',obj.folder,obj.name,floor(mean(timestep))),[16 14]); end % interactive window to display the terms of the pressure tensor dispespicPressure(obj,logdensity,showgrid,fixed,temperature) % interactive window to display the electron density, magnetic % field lines, electric potential and field at given time steps dispespicFields(obj,logdensity,showgrid,fixed,parper) function displaycollfreq(obj) %displaycollfreq plot the collision frequencies in Hz/mbar for a range of %electron kinetic energies in eV for the different collision %processes considered: ionisation and elastic collisions E=logspace(1,4,1000); v=sqrt(2/obj.msim*obj.weight*E*obj.qe); P=obj.neutcol.neutdens*obj.kb*300/100;% pressure at room temperature in mbar tauio=P./(obj.neutcol.neutdens*obj.sigio(E).*v); tauiom=P./(obj.neutcol.neutdens*obj.sigmio(E).*v); tauelam=P./(obj.neutcol.neutdens*obj.sigmela(E).*v); f=figure('name','t scales coll'); loglog(E,1./tauio,'displayname','ionisation','linewidth',1.5) hold on loglog(E,1./tauiom,'displayname','ionisation momentum','linewidth',1.5) loglog(E,1./tauelam,'displayname','elastic','linewidth',1.5) loglog(E,1./tauio+1./tauiom+1./tauelam,'displayname','total drag','linewidth',1.5) loglog([E(1) E(end)],1./(2*pi/obj.omece)* [1 1],'--','displayname','cyclotronic','linewidth',1.5) xlabel('Electron kinetic energy [eV]') ylabel('\nu [Hz/mbar]') legend('location','southeast') grid on obj.savegraph(f,sprintf('%s/collfreqscales',obj.folder),[14 12]); end function displaycrosssec(obj) %displaycrosssec plot the collision crosssections in m^2 for a range of %electron kinetic energies in eV for the different collision %processes considered: ionisation and elastic collisions E=logspace(1,4,1000); sig_io=obj.sigio(E); sig_iom=obj.sigmio(E); sig_elam=obj.sigmela(E); f=figure('name','t scales coll'); loglog(E,sig_io,'displayname','ionisation','linewidth',1.5) hold on loglog(E,sig_iom,'displayname','ionisation momentum','linewidth',1.5) loglog(E,sig_elam,'displayname','elastic','linewidth',1.5) loglog(E,sig_io+sig_elam+sig_iom,'displayname','total drag','linewidth',1.5) xlabel('Energy [eV]') ylabel('\sigma [m^{2}]') legend('location','southeast') grid on obj.savegraph(f,sprintf('%s/coll_cross_sec_scales',obj.folder),[14 12]); end %------------------------------------------ % Helper functions needed for other functions function [zpos,rpos]=getpos(obj,tstep) % interactive window to return an specific axial and radial % position picked from the cloud density if nargin<2 tstep=length(obj.t2d); end n=obj.N(:,:,tstep); n(obj.geomweight(:,:,1)<0)=NaN; figure contourf(obj.zgrid,obj.rgrid,n); xlabel('z [m]') ylabel('r [m]') [x,y]=ginput(1); zpos=find(x>obj.zgrid,1,'last'); rpos=find(y>obj.rgrid,1,'last'); hold on plot(obj.zgrid(zpos),obj.rgrid(rpos),'rx') fprintf('zpos=%i rpos=%i z=%1.4f r=%1.4f\n',zpos,rpos,obj.zgrid(zpos),obj.rgrid(rpos)) end function changed=ischanged(obj) %ischanged Check if the file has been changed since the initial loading of the file %and if some data must be reloaded try filedata=dir(obj.fullpath); checkedtimestamp=filedata.date; if (max(checkedtimestamp > obj.timestamp) ) changed=true; return end changed=false; return catch changed=true; return end end function dispV(obj,V,Vend,label,t, dist, vd) %dispV generic functio to plot the velocity distribution and %comapare two timesteps V and Vend at time t(1) and t(2) if nargin<6 dist='gaussian'; end if nargin<7 vd=0; end vmean=mean(V(~isnan(V))); vtherm=std(V(~isnan(V)),1); vmeanend=mean(Vend(~isnan(Vend))); vthermend=std(Vend(~isnan(Vend)),1); if(length(V)>1) [Counts,edges]=histcounts(V,'binmethod','sqrt'); binwidth=mean(diff(edges)); plot([edges(1) 0.5*(edges(2:end)+edges(1:end-1)) edges(end)],[0 Counts 0],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(t(1))*1e9)); hold on end hold on [Counts,edges]=histcounts(Vend,'binmethod','sqrt'); plot([edges(1) 0.5*(edges(2:end)+edges(1:end-1)) edges(end)],[0 Counts 0],'DisplayName',sprintf("t=%2.3d [ns]",obj.tpart(t(2))*1e9)); if strcmp(dist,'maxwell') vfit=linspace(0,edges(end),300); a=vmeanend/sqrt(2); dist=sqrt(2/pi)*vfit.^2.*exp(-((vfit).^2-vd^2)/2/a^2)/a^3; dist=dist/max(dist); plot(vfit,max(Counts)*dist,'displayname',sprintf('Maxw mu=%2.2g sigma=%2.2g',vmeanend,vthermend)) elseif strcmp(dist,'gaussian') vfit=linspace(edges(1),edges(end),300); dist=exp(-(vfit-vmeanend).^2/2/vthermend^2); plot(vfit,max(Counts)*dist,'displayname',sprintf('gauss mu=%2.2g sigma=%2.2g',vmeanend,vthermend)) end ylabel('counts') xlabel(label) grid on legend('location','southoutside','orientation','vertical') end function cross_sec=fit_cross_sec(obj,energy,crosssec_table) %Interpolate the cross-section at the given energy using the %crosssec_table and an exponential fitting cross_sec=0; if (energy<=0 || isnan(energy) || isinf(energy)) return end id=find(energy>crosssec_table(:,1),1,'last'); if(isempty(id)) id=1; end id=min(size(crosssec_table,1)-1,id); id=max(1,id); cross_sec=crosssec_table(id,2)*(energy/crosssec_table(id,1))^crosssec_table(id,3); end function fighandle=savegraph(obj, fighandle, name, papsize) %% Saves the given figure as a pdf a .fig and an eps using export_fig fighandle.PaperUnits='centimeters'; if (nargin < 4) papsize=[14 16]; end set(fighandle, 'Color', 'w'); fighandle.PaperSize=papsize; %export_fig(fighandle,name,'-png','-r300') exportgraphics(fighandle,sprintf('%s.png',name),'Resolution',300) print(fighandle,name,'-dpdf','-fillpage') savefig(fighandle,name) set(fighandle, 'Color', 'w'); exportgraphics(fighandle,sprintf('%s.eps',name)) %export_fig(fighandle,name,'-eps','-painters') end function sig=dsigmaio(obj,Ekin, Ebar, Ei, E0, chi, gamma) % calculates the integrand used for the ionisation collision % cross section for momentum exchange for the incoming electron % it is only used by obj.sigmiopre gamma=reshape(gamma,1,[],1); chi=reshape(chi,1,1,[]); siggamma=sin(gamma).*(E0^2+8*(1-chi)*(Ekin-Ei)*E0)./(E0+4*(1-chi)*(Ekin-Ei)-4*(1-chi)*(Ekin-Ei).*cos(gamma)).^2/2; sigchi=(Ekin-Ei)./(Ebar*atan((Ekin-Ei)/(2*Ebar)).*(1+(chi*(Ekin-Ei)/Ebar).^2)); dp=1- trapz(gamma,sqrt((1-chi).*(1-Ei/Ekin)).*cos(gamma).*siggamma,2);%- trapz(gamma,sqrt((1-chi).*(1-Ei/Ekin)).*cos(gamma).*siggamma,2); sig=sigchi.*dp; end function sigm=sigmiopre(obj,E, init) % returns the precalculated values used for the interpolation % of the ionisation collision cross-section for momentum % exchange for the incoming electron if nargin <3 init=false; end if(~init &&( ~obj.neutcol.present || isempty(obj.neutcol.io_cross_sec))) sigm=zeros(size(E)); return end Ebar=obj.neutcol.scatter_fac; Ei=obj.neutcol.Eion; E0=obj.neutcol.E0; nE=numel(E); nchi=300; ngamma=300; gamma=linspace(0,pi,ngamma); chi=linspace(0,0.5,nchi); %sigm2=zeros(nE,nchi); sigm=zeros(size(E)); for i=1:nE if(E(i)>=Ei) sigm2=zeros(nchi,1); for j=1:nchi %sigm2(j)=trapz(alpha,trapz(gamma,obj.dsigmaio(E(i),Ebar,Ei,E0,chi(j),alpha,gamma),2),1); sigm2(j)=obj.dsigmaio(E(i),Ebar,Ei,E0,chi(j),gamma); end sigm(i)=trapz(chi,sigm2)*obj.sigio(E(i),init); %sigm(i)=trapz(chi,trapz(alpha,trapz(gamma,dsigmaio(obj,E(i),Ebar,Ei,E0,chi,alpha,gamma),2),1),3)*obj.sigio(E(i),init); end end end end end diff --git a/matlab/@espic2dhdf5/potentialwellmodel.m b/matlab/@espic2dhdf5/potentialwellmodel.m index 6e798b7..d115fa5 100644 --- a/matlab/@espic2dhdf5/potentialwellmodel.m +++ b/matlab/@espic2dhdf5/potentialwellmodel.m @@ -1,133 +1,138 @@ function model=potentialwellmodel(obj,timestep,calcfull,mag,nblv) % Computes the potential well at the given timestep and return the model to be able to % interpolate either using grid coordinates or magnetic field line coordinates % the potential well is calculated along each magnetic field % line by taking the difference between a local potential well maximum % and the maximum local minimum along the line on % each side of the maximum % model.z and model.r are the axial and radial positions where % the potential depth has been calculated and is used with a % scatteredinterpolant to calculate the well depth at the % desired position % model.pot is the potential well depth % model.rathet is the magnetic field vector value at the % corresponding position multiplied by the radial position if iscell(timestep) timestep=cell2mat(timestep); end % Defines if only the well region is kept or if the hills are also kept if nargin<3 calcfull=false; end if nargin<4 mag=[]; end if nargin<5 nblv=400; end +outsidethreshold=-5e-3; % if one of the timesteps is 0 we take the external potential id=find(timestep==0); if(~isempty(id)) timestep(id)=[]; end Phi=-obj.pot(:,:,timestep); if(~isempty(id)) phiext=-obj.potxt(:,:,1); if isempty(obj.potxt) phiext=-obj.pot(:,:,1); end Phi=cat(3,Phi(:,:,1:id-1),phiext,Phi(:,:,id:end)); end % We get the magnetic field lines rA_theta values %lvls=sort(unique([obj.rAthet(:,1)',obj.rAthet(:,end)', obj.rAthet(1,:),obj.rAthet(end,:)])); if isempty(mag) Blines=obj.rAthet; -lvls=logspace(log10(min(Blines(obj.geomweight(:,:,1)>0))),log10(max(Blines(obj.geomweight(:,:,1)>0))),nblv); +vals=obj.rAthet(obj.geomweight(:,:,1)>=outsidethreshold&obj.rAthet>0); +lvls=logspace(log10(min(vals)),log10(max(vals)),nblv); %lvls=linspace(min(Blines(obj.geomweight(:,:,1)>0)),max(Blines(obj.geomweight(:,:,1)>0)),200); -Blines(obj.geomweight(:,:,1)<=0)=NaN; +Blines(obj.geomweight(:,:,1)<=outsidethreshold)=NaN; lvls(isnan(lvls))=[]; % We obtain the magnetic field lines coordinates for each % level in r and z contpoints=contourc(obj.zgrid,obj.rgrid,Blines,lvls(1:end)); else [z,r]=ndgrid(obj.zgrid,obj.rgrid); locw=griddedInterpolant(z,r,obj.geomweight(:,:,1)','linear','none'); idrmin=find(obj.rgrid(1)mag.r,1,'last'); idzmin=find(obj.zgrid(1)mag.z,1,'last'); [z,r]=ndgrid(mag.z(idzmin:idzmax),mag.r(idrmin:idrmax)); geomwB=locw(z,r)'; geomwB(isnan(geomwB))=-1; Blines=mag.Athet.*mag.r; Blines=Blines((idrmin:idrmax),(idzmin:idzmax)); lvls=logspace(log10(min(Blines(geomwB>0))),log10(max(Blines(geomwB>0))),nblv); Blines(geomwB<0)=NaN; lvls(isnan(lvls))=[]; % We obtain the magnetic field lines coordinates for each % level in r and z contpoints=contourc(mag.z(idzmin:idzmax),mag.r(idrmin:idrmax),Blines,lvls(1:end)); end [x,y,zcont]=C2xyz(contpoints); % memory allocation for spped potfinal=zeros(numel(cell2mat(x)),length(timestep)); Pot=cell(1,size(x,2)-1); rathet=cell(1,size(x,2)-1); [z,r]=ndgrid(obj.zgrid,obj.rgrid); % for each timestep we calculate the well for i=1:length(timestep)+~isempty(id) locPhi=Phi(:,:,i); % We delete points outside of the domain - locPhi(obj.geomweight(:,:,1)<0)=NaN; + %locPhi(obj.geomweight(:,:,1)=3)% We neglect if the field length is too small Pot{j}=locPhi(xloc,yloc); Pot{j}=movmean(Pot{j},1); locpot=Pot{j}; + if any(isnan(locpot),'all') + disp( "bad line") + end %locpot=locpot-min(locpot); % along the given field line j we calculate the % minimum on the left and right side of the % position k and calculate the local well depth for k=2:length(locpot)-1 left=max(locpot(1:k-1)); right=max(locpot(k+1:end)); Pot{j}(k)=locpot(k)-min(left,right); end Pot{j}(1)=0; Pot{j}(end)=0; else Pot{j}=zeros(size(xloc)); end end potfinal(:,i)=cell2mat(Pot); end %potfinal=cell2mat(potfinal); if ~calcfull - potfinal(potfinal>=0)=NaN; + %potfinal(potfinal>=0)=NaN; end model.z=cell2mat(x); model.r=cell2mat(y); model.pot=-potfinal; model.rathet=cell2mat(rathet); end \ No newline at end of file diff --git a/matlab/Magnetic_Field_GLS_2020/refurb_modif.data b/matlab/Magnetic_Field_GLS_2020/refurb_modif.data index f09b79c..8077c12 100644 --- a/matlab/Magnetic_Field_GLS_2020/refurb_modif.data +++ b/matlab/Magnetic_Field_GLS_2020/refurb_modif.data @@ -1,242 +1,244 @@ # z [mm] # r [mm] # Cathode -cthP99 -75 110 -cthP98 -75 100 -cthP97 -75 95 -cthP96 -75 90.020801 +cthP99 -82.45 110 +cthP98 -82.45 100 +cthP97 -82.45 95 +cthP96 -82.45 90.020801 cthP95 -80 90.020801 cthP94 -82.846074 90.020801 +cthP94 -85 90.020801 cthP93 -93.945889 90.020801 +cthP93 -93.945889 85 cthP92 -93.945889 80.010989 cthP911 -90 80.010989 cthP911 -85 80.010989 cthP91 -77.522162 80.010989 cthEP0 -19.9797 63.1981 cthEP1 -17.4845 63.0446 cthEP2 -14.9889 62.8954 cthEP3 -12.3933 62.7449 cthEP4 -9.8972 62.6048 cthEP5 -7.40087 62.4694 cthEP6 -4.90429 62.3388 cthEP7 -2.40746 62.213 cthEP8 0.189511 62.0874 cthEP9 2.68684 61.9718 cthEP10 5.1844 61.8615 cthEP11 7.6822 61.7565 cthEP12 10.2801 61.6531 cthEP13 12.7784 61.5594 cthEP14 15.2768 61.4714 cthEP15 17.7755 61.3892 cthEP16 20.2743 61.313 cthEP17 22.8733 61.2403 cthEP18 25.3725 61.1768 cthEP19 27.8718 61.1197 cthEP20 30.3713 61.0691 cthEP21 32.8709 61.0253 cthEP22 35.4707 60.987 cthEP23 37.9705 60.9572 cthEP24 40.4704 60.9346 cthEP25 42.9703 60.9193 cthEP26 45.4703 60.9114 cthEP27 48.0703 60.9112 cthEP28 50.5703 60.9189 cthEP29 53.0703 60.9343 cthEP30 55.5701 60.9576 cthEP31 58.1699 60.9903 cthEP32 60.6696 61.03 cthEP33 63.1692 61.0778 cthEP34 65.6685 61.1338 cthEP35 68.1677 61.1979 cthEP36 70.7666 61.2733 cthEP37 73.2653 61.354 cthEP38 75.7637 61.4427 cthEP39 78.2619 61.5392 cthEP40 80.7597 61.6434 cthEP41 83.3571 61.7595 cthEP42 85.8543 61.8784 cthEP43 88.3511 62.0038 cthEP44 90.8477 62.1353 cthEP45 93.4438 62.2778 cthEP46 95.9397 62.4195 cthEP47 98.4355 62.565 cthEP48 100.931 62.7131 cthEP49 103.526 62.8685 cthP2 105.72253 63 cthP4 106.72253 63 cthP5 109.57202 62.864366 cthP6 112.39575 62.458692 cthP7 115.16818 61.786644 cthP8 117.86425 60.854301 emtrP2 120.25 59.9 emtrP1 124.75 58.1 cthP13 131.1696 55.53216 cthP14 132.70522 54.781069 cthP15 134.11253 53.81064 cthP16 135.36045 52.642314 cthP17 136.42138 51.301909 cthP18 137.27189 49.819043 cthP19 137.89319 48.226481 cthP20 138.27155 46.559414 cthP21 138.39861 44.854678 cthP24 138.2641 43.469046 cthP25 137.86576 42.155894 cthP26 137.21889 40.945686 cthP27 136.34835 39.884931 cthP28 135.2876 39.014391 cthP29 134.07739 38.367521 cthP30 132.76424 37.969181 cthP31 131.39861 37.834678 cthP11 122.5 37.834678 cthP32 -100 57.5 cthP40 -154.21892 64.40626 cthP42 -158.26726 65.357779 cthP44 -162.03129 67.125946 cthP46 -165.34829 69.634314 cthP48 -168.07484 72.77443 cthP52 -168.90702 74.602138 cthP54 -169.02759 76.606762 cthP56 -168.42046 78.521037 cthP58 -167.16659 80.089745 cthP60 -165.43313 81.103738 cthP62 -163.45121 81.427828 cthP64 -161.48505 81.018804 cthP78 -159.79681 79.9312 cthP66 -149.8737 79.9312 cthP67 -149.8737 89.9172 cthP68 -161.0308 89.9172 cthP69 -161.0308 95 cthP70 -161.0308 100 cthP71 -161.0308 110 # Body andP13 300 36.3 andP11 203 37.3 andP1 196 40.4 andP4 190 43.3 andP5 187.60079 45.8374 andP5 185.20158 48.3748 andP5 182.80237 50.9122 andP5 180.40316 53.4496 andP5 178.00395 55.9870 andP5 175.60474 58.5244 andP5 173.20553 61.0618 andP5 170.80632 63.5992 andP5 168.40711 66.1366 andP5 166.00790 68.6740 andP5 163.60869 71.2114 andP5 161.20948 73.7488 andP5 158.81027 76.2862 andP12 156.41106 78.8236 andP13 151.912526 80.623013 andP14 147.413992 82.422426 andP15 142.915458 84.221839 andP16 138.416924 86.021252 andP2 133.91839 87.820665 crnEP50 131.995 87.5419 crnEP49 129.846 87.0707 crnEP48 127.792 86.6332 crnEP47 125.637 86.1882 crnEP46 123.48 85.7572 crnEP45 121.418 85.3589 crnEP44 119.255 84.9557 crnEP43 117.188 84.5842 crnEP42 115.021 84.2092 crnEP41 112.949 83.8648 crnEP40 110.777 83.5181 crnEP39 108.701 83.2007 crnEP38 106.524 82.8822 crnEP37 104.444 82.5915 crnEP36 102.363 82.3139 crnEP35 100.18 82.0368 crnEP34 98.0952 81.7854 crnEP33 96.0088 81.5466 crnEP32 93.8216 81.3098 crnEP31 91.7324 81.0964 crnEP30 89.6421 80.8952 crnEP29 87.451 80.6973 crnEP28 85.3585 80.5204 crnEP27 83.265 80.3552 crnEP26 81.0709 80.1943 crnEP25 78.9757 80.0522 crnEP24 76.8798 79.921 crnEP23 74.7833 79.8007 crnEP22 72.6861 79.6908 crnEP21 70.4886 79.5869 crnEP20 68.3905 79.4979 crnEP19 66.292 79.4189 crnEP18 64.1931 79.3494 crnEP17 61.994 79.2868 crnEP16 59.8946 79.2364 crnEP15 57.795 79.195 crnEP14 55.6953 79.1623 crnEP13 53.5954 79.1382 crnEP12 51.3955 79.1218 crnEP11 49.2955 79.1145 crnEP10 47.1955 79.1151 crnEP9 45.0955 79.1235 crnEP8 42.9956 79.1394 crnEP7 40.7957 79.164 crnEP6 38.696 79.1947 crnEP5 36.5963 79.2325 crnEP4 34.4968 79.2771 crnEP3 32.3974 79.3284 crnEP2 30.1982 79.3891 crnEP1 28.0992 79.4536 crnP0 26.7 79.5 crnP4 24.365737 79.996789 crnP8 22.336075 81.2522 crnP12 22.011506 82.151259 crnP16 22.590076 82.91212 crnP20 23.543174 82.83962 crnP24 24 82 crnP24 28 82 crnP24 32 82 crnP24 36 82 crnP31 39.053729 82 crnP311 39.053729 84 crnP312 39.053729 86 crnP313 39.053729 88 crnP32 39.053729 90.020801 crnP321 36 90.020801 crnP322 32 90.020801 crnP33 28 90.020801 crnP32 27.953913 90.020801 crnP38 27.953913 100 crnP39 27.953913 110 # Coaxial insert cthP71 -271.6998 110 cthP70 -271.6998 100 cthP69 -271.6998 95 insP14 -271.6998 89.9172 insP13 -282.884 89.9172 insP12 -282.884 69.9453 insP11 -296.9384 69.9453 insP10 -296.9384 49.9742 insP9 -259.5 49.9735 insP89 -250 49.1295 insP88 -208.75 45.4648 insP87 -167.5 41.8002 insP86 -126.25 38.1355 insP85 -85 34.4708 insP84 -43.75 30.8062 insP83 -2.5 27.1415 insP82 38.75 23.4768 insP81 80 19.8121 insP8 92.5184 18.7 insP72 110 18.7 insP71 120 18.7 insP71 130 18.7 insP71 140 18.7 insP71 150 18.7 insP71 160 18.7 insP71 170 18.7 insP71 180 18.7 insP71 200 18.7 insP3 217.3 18.7 insP4 250 12 insP5 250 8.4 insP6 300 8.4 \ No newline at end of file diff --git a/src/psupply_mod.f90 b/src/psupply_mod.f90 index bbec101..798387e 100644 --- a/src/psupply_mod.f90 +++ b/src/psupply_mod.f90 @@ -1,325 +1,325 @@ module psupply use constants implicit none type power_supply logical :: active = .false. ! is the power supply active real(kind=db):: geomcapacitor = 1 ! capacitance of the metalic vessel normalised to the neutral density in vessel used in the neutcol module real(kind=db):: PSresistor = 1 ! internal resistance of the power supply normalised to the actual neutral density in vessel real(kind=db):: targetbias = 0 ! Set voltage on the power supply real(kind=db):: bias = 0 ! current voltage on the power supply integer :: nbbounds = 2 ! number of boundaries defined in the geometry integer :: lststp = 0 ! previous step on which the bias was updated real(kind=db):: current(3) = 0 ! current collected on the boundaries normalised to the simulated collision neutral density set in neutcol module ! 1 is at time i-2nbhdt, 2 is at time i-nbhdt and 3 is at time i integer, allocatable:: bdpos(:) ! sign of each boundary for collected charge to determine direction of current real(kind=db),allocatable:: charge(:) ! Charge collected on each boundary and per nbdt real(kind=db),allocatable:: biases(:) ! Actual potentials at each boundary integer :: nbhdt = 10 ! half of the number of time steps between each calls to RK4 real(kind=db):: expdens ! [m-3] experimental neutral density real(kind=db):: neutcoldens ! [m-3] neutral density used in neutcol module real(kind=db):: frequency ! [Hz] frequency of an imposed oscillation in bias real(kind=db):: deltabias ! [V] Amplitude of the oscillations around targetbias end type type(power_supply):: the_ps contains ! read the input parameters from the input file and setup the necesary variables for the module to work subroutine psupply_init(fileid,cstep,nbbounds,neutcoldens,rstbias) use splinebound use basic, only: phinorm, tnorm, mpirank, qnorm, potinn, potout, dt use constants use mpi use geometry use weighttypes use fields integer:: fileid, cstep, nbbounds, istat, nbhdt, ierr,i real(kind=db),OPTIONAL, INTENT(IN):: rstbias real(kind=db):: neutcoldens ! [m-3] real(kind=db):: expneutdens = 1 ! [m-3] real(kind=db):: PsResistor = 1 ! [Ohm] real(kind=db):: geomcapacitor = 1 ! [F] real(kind=db):: targetbias = 0 ! [V] real(kind=db):: frequency = 0 ! [Hz] real(kind=db):: deltabias = 0 ! [V] integer, allocatable:: bdpos(:) character(len=1000) :: line logical :: active = .false. NAMELIST /psupplyparams/ expneutdens, PsResistor, geomcapacitor, targetbias, nbhdt, active, bdpos, frequency, deltabias the_ps%lststp=cstep the_ps%nbbounds=nbbounds allocate(the_ps%bdpos(nbbounds),bdpos(nbbounds)) allocate(the_ps%charge(nbbounds),the_ps%biases(nbbounds)) the_ps%bdpos=0 bdpos=0 the_ps%charge=0 the_ps%biases=0 the_ps%current=0 ! read the input parameters from file Rewind(fileid) READ(fileid, psupplyparams, iostat=istat) if (istat.gt.0) then backspace(fileid) read(fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in pssupplyparams: '//trim(line) call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if ! save the parameters on output IF(mpirank .eq. 0) THEN WRITE(*, psupplyparams) END IF ! rescale the targetbias set on the power supply the_ps%targetbias=abs(targetbias)/phinorm the_ps%bdpos=bdpos the_ps%frequency=frequency*tnorm*2*pi the_ps%deltabias=deltabias/phinorm ! save the experimental neutral density the_ps%expdens=expneutdens ! save the neutral collision density the_ps%neutcoldens=neutcoldens if(present(rstbias))then ! Initialize the current bias from the restart value the_ps%bias=rstbias else ! initialize with the file input parameters if (the_domain%nbsplines.gt.0) then do i=1,the_ps%nbbounds if(the_ps%bdpos(i) .lt. 0)then the_ps%bias=-the_domain%boundaries(i)%Dirichlet_val exit end if end do else the_ps%bias=(potout-potinn) end if end if ! set the initial bias where(the_ps%bdpos .lt. 0) - the_ps%biases=the_ps%bias*the_ps%bdpos + the_ps%biases=-the_ps%bias end where ! Normalise resistor and capacitor to adapt to experimental pressure the_ps%PSresistor = PSresistor*the_ps%expdens/the_ps%neutcoldens*qnorm/(tnorm*phinorm) the_ps%geomcapacitor = geomcapacitor*phinorm/qnorm the_ps%nbhdt = nbhdt the_ps%active = active if( .not. the_ps%active) return ! Initialize the biases if (the_domain%nbsplines.gt.0) then do i=1,the_ps%nbbounds the_domain%boundaries(i)%Dirichlet_val=the_ps%biases(i)+the_ps%deltabias*sin(the_ps%frequency*cstep*dt) end do else potinn=the_ps%biases(1)+the_ps%deltabias*sin(the_ps%frequency*cstep*dt) Potout=0 Phidown=Potinn Phiup=Potout end if ! recalculate gtilde to adapt for the new biases CALL total_gtilde(vec1, vec2, gtilde, gridwdir) call comp_gradgtilde ! Recompute the vacuum field call vacuum_field end subroutine ! save to the result file the parameters of this module read from the input Subroutine psupply_diag(File_handle, str) use mpi Use futils use basic, only: tnorm, phinorm, qnorm implicit none Integer:: File_handle Character(len=*):: str CHARACTER(len=256):: grpname Integer:: ierr, mpirank CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0 .and. the_ps%active) THEN Write(grpname,'(a,a)') trim(str),"/psupply" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "expdens", the_ps%expdens) Call attach(File_handle, trim(grpname), "targetbias", the_ps%targetbias*phinorm) Call attach(File_handle, trim(grpname), "PSresistor", the_ps%PSresistor/the_ps%expdens*the_ps%neutcoldens/qnorm*(tnorm*phinorm)) Call attach(File_handle, trim(grpname), "geomcapacitor", the_ps%geomcapacitor/phinorm*qnorm) Call attach(File_handle, trim(grpname), "nbhdt", the_ps%nbhdt) Call putarr(File_handle,trim(grpname)//"/bdpos", the_ps%bdpos) END IF End subroutine psupply_diag ! gneral routine called from stepon to update the psupply bias subroutine psupply_step(ps,p,cstep) use particletypes use geometry use weighttypes use fields use basic, only: Potinn, potout type(power_supply):: ps type(particles):: p(:) integer:: cstep, i if (.not. ps%active ) return ! calculate the charge collected on each boundary due to the contribution of each specie call add_charge(ps,p) ! calculate the current flowing between the electrodes due to the cloud call calc_current(ps,cstep) ! calculate the bias at the new time step call updt_bias(ps,cstep) if(mod(cstep-ps%lststp,2*ps%nbhdt) .ne. 0) return ! update the bias on the geometry for the Dirichlet b.c. if (the_domain%nbsplines.gt.0) then do i=1,ps%nbbounds the_domain%boundaries(i)%Dirichlet_val=ps%biases(i) end do else potinn=ps%biases(1) Potout=0 Phidown=Potinn Phiup=Potout end if ! recalculate gtilde to adapt for the new biases CALL total_gtilde(vec1, vec2, gtilde, gridwdir) call comp_gradgtilde end subroutine ! calculates the current flowing between the electrodes due to the cloud subroutine calc_current(ps,cstep) use geometry use basic, only: phinorm, dt use fields type(power_supply):: ps integer:: cstep if(mod(cstep-ps%lststp,ps%nbhdt).eq.0) then ! communicate the charge accumulation in this timestep call reduce_charge(ps) if(mod(cstep-ps%lststp,ps%nbhdt*2).eq.0)then ! calculates the current by adding the contribution of each boundary if (mpirank .eq. 0)then ps%current(3)=sum(-ps%charge*ps%bdpos)/(ps%nbhdt*dt) end if ps%lststp=cstep else ! calculates the current by adding the contribution of each boundary if (mpirank .eq. 0)then ps%current(2)=sum(-ps%charge*ps%bdpos)/(ps%nbhdt*dt) end if end if ps%charge=0 end if end subroutine ! calculate the charge deposited by each specie on the electrodes (used to calculate the resulting current) subroutine add_charge(ps,p) use particletypes use basic, only: qnorm type(power_supply):: ps type(particles):: p(:) integer:: i do i=1,size(p,1) if(.not. p(i)%is_field) cycle !Add the normalised contribution of each specie ps%charge=ps%charge+p(i)%nblost(5:)*p(i)%weight*p(i)%q/qnorm end do end subroutine ! Time integrate the ODE of the actual bias between the accelerating electrodes ! and broadcast it to all the workers subroutine updt_bias(ps,cstep) use basic, only: dt implicit none type(power_supply):: ps integer:: cstep real(kind=db):: bias,k1,k2,k3,k4, hdeltat if(mod(cstep-ps%lststp,2*ps%nbhdt) .ne. 0) return ! half delta t if (ps%PSresistor.gt.0)then hdeltat=dt*ps%nbhdt/(ps%PSresistor*ps%geomcapacitor) bias=ps%bias ! we update the bias using RK4 k1=-(bias+ps%current(1)*ps%PSresistor-ps%targetbias) k2=-(bias+hdeltat*k1+ps%current(2)*ps%PSresistor-ps%targetbias) k3=-(bias+hdeltat*k2+ps%current(2)*ps%PSresistor-ps%targetbias) k4=-(bias+2*hdeltat*k3+ps%current(3)*ps%PSresistor-ps%targetbias) ps%bias=bias+(k1+2*k2+2*k3+k4)*2*hdeltat/6 end if !Write(*,*) " new bias ", ps%bias*phinorm where (ps%bdpos .lt. 0) - ps%biases=-ps%bias + ps%biases=-ps%bias+ps%deltabias*sin(ps%frequency*cstep*dt) end where where (ps%bdpos .eq. 2) ps%biases=ps%deltabias*sin(ps%frequency*cstep*dt) end where ! broadcast the bias to all the mpi processes call bcast_bias(ps) ps%current(1)=ps%current(3) end subroutine updt_bias ! gather on node 0 the collected charge on each metallic boundary subroutine reduce_charge(ps) use mpi use mpihelper use basic, ONLY: mpirank type(power_supply):: ps integer:: ierr if(mpirank .eq. 0) then call MPI_REDUCE(MPI_IN_PLACE,ps%charge,ps%nbbounds,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) !Write(*,*) "curr charge ", ps%charge else call MPI_REDUCE(ps%charge,ps%charge,ps%nbbounds,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) end if end subroutine ! broadcast to all the nodes the new bias imposed by the power supply on the electrodes subroutine bcast_bias(ps) use mpi use mpihelper type(power_supply):: ps integer:: ierr call MPI_BCAST(ps%biases,ps%nbbounds,db_type,0,MPI_COMM_WORLD,ierr) end subroutine end module psupply diff --git a/src/splinebound_mod.f90 b/src/splinebound_mod.f90 index 051ae73..59ccdd8 100644 --- a/src/splinebound_mod.f90 +++ b/src/splinebound_mod.f90 @@ -1,930 +1,931 @@ MODULE splinebound USE constants USE bsplines USE forSISL, only: newCurve, freeCurve, freeIntCurve, writeSISLcurve, writeSISLpoints Use forSISLdata IMPLICIT NONE - INTEGER(SELECTED_INT_KIND(4)), PARAMETER :: bd=-1, bd_Dirichletconst=0, bd_Dirichletvar=1, bd_Neumann=2 + INTEGER, PARAMETER :: bd=-1, bd_Dirichletconst=0, bd_Dirichletvar=1, bd_Neumann=2 type cellkind integer:: spldirkind=0 !< -1 outside (return -1) no dist to calculate; 0 boundary calculate dist with linked boundaries; 1 inside (return 1) no dist to calculate integer:: spltotkind=0 !< -1 outside (return -1) no dist to calculate; 0 boundary calculate dist with linked boundaries; 1 inside (return 1) no dist to calculate integer:: linkedboundaries(2)=0 !< stores the spline curve indices in the spline_domain of the spline boundaries that are the closest and at a distance lower than dist_extent (1) !< (1) is for dirichlet boundaries !< (2) is for domain boundaries integer:: leftknot(4)=0 !< knots pointer for s1424 in wtot then wdir real(kind=db):: lguess(2)=-1 !< Spline parameter left limit as start guess real(kind=db):: rguess(2)=-1 !< Spline parameter right limit as start guess end type cellkind TYPE spline_boundary ! all curves assume right handedness to set which side of the curve is inside or outside type(SISLCurve):: curve Real(kind=db):: Dirichlet_val !< Value for the dirichlet boundary condition created by this boundary Real(kind=db):: epsge=1.0e-5 !< geometric resolution used for calculating distances Real(kind=db):: epsce=1.0e-9 !< value of weight below which it is 0 INTEGER(kind(bd)):: type=bd_Dirichletconst !< type of boundary conditions END TYPE spline_boundary type spline_domain integer:: nbsplines = 0 !< number of spline boundaries in the domain type(spline_boundary), allocatable:: boundaries(:) !< List of boundaries in the domain Real(kind=db):: dist_extent=0.1 !< distance used for the merging with the plateau function for the weight type(cellkind), ALLOCATABLE:: cellk(:,:) !< Precomputed parameters at each cell for faster weight computation type(spline2d), pointer:: splrz => null() !< Pointer to the main spline grid used for the FEM solver Integer:: nb1 !< Number of grid points in the 1st dimension Integer:: nb2 !< Number of grid points in the 2nd dimension real(kind=db), ALLOCATABLE:: x1(:) !< Grid points in first direction for weight interpolation real(kind=db), ALLOCATABLE:: x2(:) !< Grid points in 2nd direction for weight interpolation real(kind=db), ALLOCATABLE:: dx1(:) !< inverse cell width in first direction for weight interpolation real(kind=db), ALLOCATABLE:: dx2(:) !< inverse cell width in 2nd direction for weight interpolation !type(SISLsurf):: Dirdomweight !< structure storing precalculated geometric weight for faster evaluation !type(SISLsurf):: totdomweight !< structure storing precalculated total weight for faster evaluation type(spline2d):: Dirdomweightspl !< structure storing precalculated geometric weight for faster evaluation type(spline2d):: totdomweightspl !< structure storing precalculated total weight for faster evaluation end type spline_domain CONTAINS !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Reads a spline domain from the namelist or from a h5 file. Needs to be called for the initialization of the module !> @param[in] Fileid Text file id of the input file containing namelists !> @param[out] spldom spline domain !> @param[in] splrz bspline structure used by the FEM comming form bspline library !> @param[in] rnorm distance normalization constant !> @param[in] phinorm electric potential normalization constant !--------------------------------------------------------------------------- subroutine read_splinebound(Fileid, spldom, splrz, rnorm, Phinorm) use mpi Integer:: Fileid type(spline_domain):: spldom type(spline2d):: splrz real(kind=db):: rnorm, phinorm Integer:: nbsplines, istat, mpirank, ierr real(kind=db):: dist_extent Character(len=128):: h5fname="", line real(kind=db) :: Dvals(30)=0 integer:: i namelist /spldomain/ nbsplines, dist_extent, h5fname, Dvals CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) REWIND(fileid) READ(fileid,spldomain, iostat=istat) if (istat.gt.0) then if(mpirank .eq. 0) then backspace(fileid) read(fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in geomparams: '//trim(line) end if call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if if(mpirank .eq. 0) WRITE(*, spldomain) Dvals=Dvals/phinorm dist_extent=dist_extent/rnorm if (.not. trim(h5fname)=='' ) then call setspline_domain(spldom, splrz, dist_extent, 0) call splinebound_readh5domain(h5fname,spldom, rnorm, phinorm) call classifycells(spldom) do i=1,spldom%nbsplines spldom%boundaries(i)%Dirichlet_val=Dvals(i) end do return else WRITE(*,*) "Error the filename h5fname is not defined. No boundary has been set!" call mpi_Abort(MPI_COMM_WORLD, -1, ierr) end if end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Saves the spline boundaries to the result file !> @param[in] File_handle futils h5 file id !> @param[in] curr_grp groupname under which the boundaries must be saved !> @param[in] spldom spline domain !--------------------------------------------------------------------------- Subroutine splinebound_diag(File_handle, curr_grp, spldom) use mpi Use futils Use basic, ONLY: rnorm, phinorm Integer:: File_handle type(spline_domain):: spldom Character(len=*):: curr_grp CHARACTER(len=128):: grpname Integer:: ierr, mpirank, i CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0) THEN Write(grpname,'(a,a)') trim(curr_grp),"/geometry_spl" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "dist_extent",spldom%dist_extent) Call attach(File_handle, trim(grpname), "nbsplines", spldom%nbsplines) do i=1,spldom%nbsplines Write(grpname,'(a,a,i2.2)') trim(curr_grp),"/geometry_spl/",i If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "Dirichlet_val", spldom%boundaries(i)%Dirichlet_val*phinorm) Call attach(File_handle, trim(grpname), "order", spldom%boundaries(i)%curve%ik) Call attach(File_handle, trim(grpname), "kind", spldom%boundaries(i)%curve%ikind) + Call attach(File_handle, trim(grpname), "type", spldom%boundaries(i)%type) Call attach(File_handle, trim(grpname), "dim", spldom%boundaries(i)%curve%idim) CALL putarr(File_handle, TRIM(grpname)//"/pos", spldom%boundaries(i)%curve%ecoef*rnorm) CALL putarr(File_handle, TRIM(grpname)//"/knots", spldom%boundaries(i)%curve%et) end do END IF End subroutine splinebound_diag !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Read a spline boundary domain from an h5 file structure !> @param[out] spldom new spline domain !> @param[in] filename filename of the h5 file !> @param[in] rnorm distance normalization constant !> @param[in] phinorm electric potential normalization constant !--------------------------------------------------------------------------- subroutine splinebound_readh5domain(filename, spldom, rnorm, phinorm) use futils use forSISL implicit none Character(len=*),intent(in) :: filename type(spline_domain),intent(inout) :: spldom integer:: h5id, i real(kind=db):: rnorm, phinorm CHARACTER(len=128):: grpname integer:: periodic integer:: order, dim, bdtype INTEGER:: posrank, posdim(2), err real(kind=db):: Dval, epsge, epsce real(kind=db),allocatable:: points(:,:) call openf(filename, h5id,'r','d') call getatt(h5id, '/geometry_spl/','nbsplines', spldom%nbsplines) ! prepare memory if (allocated(spldom%boundaries)) then do i=1,size(spldom%boundaries,1) call free_bsplinecurve(spldom%boundaries(i)) end do DEALLOCATE(spldom%boundaries) end if allocate(spldom%boundaries(spldom%nbsplines)) ! Read each boundary curve individually do i=1,spldom%nbsplines Write(grpname,'(a,i2.2)') "/geometry_spl/",i If(.not. isgroup(h5id, trim(grpname))) THEN Write(*,*) "Error the geometry definition file is invalid" END IF periodic=0 Call getatt(h5id, trim(grpname), "Dirichlet_val", Dval) Call getatt(h5id, trim(grpname), "epsge", epsge) Call getatt(h5id, trim(grpname), "epsce", epsce) Call getatt(h5id, trim(grpname), "order", order) Call getatt(h5id, trim(grpname), "dim", dim) err=0 Call getatt(h5id, trim(grpname), "periodic", periodic,err) if(err .lt.0) periodic=0 CALL getdims(h5id, TRIM(grpname)//"/pos", posrank, posdim) allocate(points(posdim(1),posdim(2))) CALL getarr(h5id, TRIM(grpname)//"/pos", points) points=points/rnorm Call setspline_boundary(spldom%boundaries(i),transpose(points), order-1, Dval/phinorm, epsge,epsce, periodic) bdtype=bd err=0 Call getatt(h5id, trim(grpname), "type", bdtype,err) if(err.ge.0) spldom%boundaries(i)%type=bdtype deallocate(points) end do call closef(h5id) end subroutine splinebound_readh5domain !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> initialize a spline domain and allocate the necessary memory !> @param[out] spldom new spline domain !> @param[in] splrz bspline structure used by the FEM comming form bspline library !> @param[in] dist_extent normalized characteristic fall lenght of the weight !> @param[in] nb_splines number of boundary splines to allocate !--------------------------------------------------------------------------- subroutine setspline_domain(spldom,splrz,dist_extent, nb_splines) type(spline_domain):: spldom type(spline2d), TARGET:: splrz real(kind=db):: dist_extent integer:: nb_splines, nb1, nb2 ! Store the grid parameters to speed-up calculations nb1=splrz%sp1%nints nb2=splrz%sp2%nints spldom%nb1=nb1 spldom%nb2=nb2 spldom%splrz=>splrz allocate(spldom%cellk(0:nb1-1,0:nb2-1)) allocate(spldom%x1(0:nb1)) allocate(spldom%x2(0:nb2)) allocate(spldom%dx1(0:nb1-1)) allocate(spldom%dx2(0:nb2-1)) spldom%x1(0:)=splrz%sp1%knots(0:nb1) spldom%x2(0:)=splrz%sp2%knots(0:nb2) spldom%dx1(0:)=1/(spldom%x1(1:nb1)-spldom%x1(0:nb1-1)) spldom%dx2(0:)=1/(spldom%x2(1:nb2)-spldom%x2(0:nb2-1)) !Prepare structures to host singular spline boundaries spldom%nbsplines=nb_splines if(spldom%nbsplines.gt. 0) allocate(spldom%boundaries(nb_splines)) spldom%dist_extent=dist_extent end subroutine setspline_domain !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief initialize a spline boundary and allocate the necessary memory !> @param[out] b_curve new spline boundary !> @param[in] cpoints control points at the node positions !> @param[in] degree degree of the spline polynomia defining the boundary curve !> @param[in] D_val Normalized value of the Dirichlet boundary condition for this curve !> @param[in] epsge geometric precision used by SISL !> @param[in] epsce arithmetic precision used by SISL !> @param[in] periodic set if the spline curve is periodic !--------------------------------------------------------------------------- subroutine setspline_boundary(b_curve, cpoints, degree, D_val, epsge, epsce, periodic) Use bsplines use forSISL,ONLY: newcurve, s1630 use mpi type(spline_boundary):: b_curve Real(kind=db):: cpoints(:,:) Real(REAL64),ALLOCATABLE:: points(:) Real(REAL64):: astpar integer:: degree integer, optional:: periodic Integer:: order, ierr Real(kind=db):: D_val Real(kind=db),OPTIONAL :: epsge, epsce Integer:: nbpoints, dim, jstat, bsptype integer:: period period=0 if(present(periodic))period=periodic nbpoints= size(cpoints,2) dim=size(cpoints,1) order=degree+1 if(nbpoints .lt. order) then WRITE(*,'(a,i3,a,i5)') "Error: the number of points", nbpoints, " is insuficient for the required order ", order CALL mpi_finalize(ierr) call EXIT(-1) end if allocate(points(dim*nbpoints)) points=reshape(cpoints,(/dim*nbpoints/)) bsptype=1 ! open boundaries b-spline if(period.gt.0) bsptype=-1 ! closed periodic curve astpar=0.0 ! starting parameter for the knots vector ! initialize a new curve using SISL CALL s1630(points, nbpoints, astpar, bsptype, dim, order, b_curve%curve, jstat) if (jstat > 0 ) WRITE(*,*) "Warning ", jstat," in curve initialisation s1630 for splineweight" if (jstat < 0 ) WRITE(*,*) "Error ", jstat," in curve initialisation s1630 for splineweight" b_curve%Dirichlet_val=D_val if(present(epsge)) b_curve%epsge=epsge if(present(epsce)) b_curve%epsce=epsce end subroutine setspline_boundary !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the Dirichlet boundary weight from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] x1(:) array of axial positions where the weights are evaluated !> @param[in] x2(:) array of radial positions where the weights are evaluated !> @param[out] w(:,0:) matrix of weights with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_w(spldom,x1,x2,w) use forSISL,ONLY: s1424 use bsplines type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: w(:,0:) Integer,allocatable::i(:),j(:) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) if (size(w,2).gt.1) then CALL speval(spldom%Dirdomweightspl, x1, x2, i, j, w(:,0), w(:,1), w(:,2)) else CALL speval(spldom%Dirdomweightspl, x1, x2, i, j, w(:,0)) end if End SUBROUTINE spline_w !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the total geometric weight from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] x1(:) array of axial positions where the weights are evaluated !> @param[in] x2(:) array of radial positions where the weights are evaluated !> @param[out] w(:,0:) matrix of weights with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_wtot(spldom,x1,x2,w,idwall) use forSISL,ONLY: s1424 use bsplines type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: w(:,0:) INTEGER, optional, INTENT(OUT):: idwall(:) Integer:: k Integer,allocatable::i(:),j(:) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) if(present(idwall)) then Do k=1,size(x2,1) idwall(k)=spldom%cellk(i(k),j(k))%linkedboundaries(2) END DO end if if (size(w,2).gt.1) then CALL speval(spldom%totdomweightspl, x1, x2, i, j, w(:,0), w(:,1), w(:,2)) else CALL speval(spldom%totdomweightspl, x1, x2, i, j, w(:,0)) end if End SUBROUTINE spline_wtot !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the interpolation in the domain of the Dirichlet boundary conditions from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] z(:) array of axial positions where the weights are evaluated !> @param[in] r(:) array of radial positions where the weights are evaluated !> @param[out] g(:,0:) matrix of boundary interpolations g with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_g(spldom,x1,x2,g,w) use forSISL,ONLY: s1424 use bsplines type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: g(:,0:) Real(kind=db), INTENT(IN),OPTIONAL::w(:,0:) REAL(real64),allocatable:: gtmp(:,:) Integer:: k Integer,allocatable::i(:),j(:) !type(cellkind):: cellk allocate(gtmp(size(x2,1),0:size(g,2)-1)) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) if(present(w)) then gtmp=w else CALL speval(spldom%Dirdomweightspl, x1, x2,i,j, gtmp(1,:), gtmp(2,:), gtmp(3,:)) end if Do k=1,size(x2,1) if(spldom%cellk(i(k),j(k))%spldirkind.eq.0)then if(gtmp(k,0) .ge. 0) then if(size(g,2) .gt. 1) then g(k,1:2)=-gtmp(k,1:2)*spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val end if g(k,0)=(1-gtmp(k,0))*spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val else g(k,0)=spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val if(size(g,2).gt. 1) then g(k,1:2)=0 end if end if else g(k,:)=0 end if end DO End SUBROUTINE spline_g !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Evaluates the geometric weight induced by the spline curve defined by b_curve at position (z,r) !> @param[in] b_curve spline_boundary containing the spline curve parameters !> @param[in] z axial position where the weight is evaluated !> @param[in] r radial position where the weight is evaluated !> @param[out] weight(:) weight index defines the order of derivation by r or z !> @param[in] h distance from the spline at which the weight is 1 !> @param[out] distance unscaled distance between evaluation point and spline b_curve !> @param[inout] leftknot initial guess for the closest spline knot of the points (r,z) !--------------------------------------------------------------------------- subroutine splineweight(b_curve, z, r, weight, h, distance, guess, lguess, rguess) Use forSISL, ONLY: s1227,s1221, s1774 type(spline_boundary):: b_curve Real(kind=db)::r,z Real(kind=db):: weight(0:) Real(kind=db),OPTIONAL:: distance real(kind=db),OPTIONAL:: guess real(kind=db),OPTIONAL:: lguess real(kind=db),OPTIONAL:: rguess integer:: sstatus, der, left,siz real(kind=db):: h, d, tpos, proj, norm real(kind=real64):: curvepos(2*b_curve%curve%idim) real(kind=db):: leftpar, rightpar,guesspar weight=0 der=1 sstatus=-1 guesspar=-1.0_db if(present(lguess) .and. present(rguess)) then leftpar=lguess rightpar=rguess guesspar=(lguess+rguess)/2 call s1774(b_curve%curve,(/z,r/),b_curve%curve%idim,b_curve%epsge,leftpar,rightpar,guesspar,tpos,sstatus) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1774 for splineweight at ", z, r else call dist(b_curve,(/z,r/),d,tpos) end if ! position and derivative wrt r,z call s1221(b_curve%curve,der,tpos,left,curvepos,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1227 for splineweight at ", z, r if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1227 for splineweight at ", z, r d=sqrt((curvepos(1)-z)**2+(curvepos(2)-r)**2) weight(0)=1-max((h-d)/h,0.0_db)**3 norm=sqrt(curvepos(3)**2+curvepos(4)**2) if(norm.gt.0) curvepos(3:4)=curvepos(3:4)/norm ! if the projection of the distance vector on the normal is negative, the weight is negative proj=(-(z-curvepos(1))*curvepos(4)+(r-curvepos(2))*curvepos(3)) if (proj .lt. 0 .or. abs(abs(proj) -sqrt((z-curvepos(1))**2+(r-curvepos(2))**2)).gt.1e-10) weight(0)=-weight(0) !if (proj .lt. 0 ) weight(0)=-weight(0) siz=size(weight,1) if (size(weight,1).gt.1 .and. abs(weight(0)) .lt. 1) then weight(1)=-3*curvepos(4)*(h-d)**2/h**3 weight(2)=+3*curvepos(3)*(h-d)**2/h**3 end if if(present(distance)) distance=d if(present(guess)) guess=tpos end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the closest distance between the point and the selected spline b_curve !> @param[in] b_curve spline_boundary containing the spline curve parameters !> @param[in] point(:) array containing the position from which to calculate the distance !> @param[out] distance distance from the point to the spline !> @param[in] pos parameter value of the closest point on the spline !--------------------------------------------------------------------------- subroutine dist(b_curve, point, distance, pos) Use forSISL, ONLY: s1957,s1953, s1221,s1227 type(spline_boundary):: b_curve Real(kind=db):: point(:) real(kind=db):: distance Real(kind=db),optional::pos REAL(real64):: posres, epsco, epsge,curvepos(2),d,distmin REAL(real64),allocatable::intpar(:) integer:: numintpt, numintcu,i,left,sstatus type(SISLIntCurve),ALLOCATABLE:: intcurve(:) epsco=1.0e-15 epsge=1.0e-15 !epsco=0 !epsge=b_curve%epsge numintpt=0 sstatus=0 distmin=HUGE(d) call s1953(b_curve%curve,point,b_curve%curve%idim,epsco,epsge,numintpt,intpar,numintcu,intcurve,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1953 for splineweight at ", point(1), point(2) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1953 for splineweight at ",point(1), point(2) if(numintpt .gt. 1) then Do i=1,numintpt call s1227(b_curve%curve,0,intpar(i),left,curvepos,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1221 for splineweight at ", point(1), point(2) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1221 for splineweight at ",point(1), point(2) d=(curvepos(1)-point(1))**2+(curvepos(2)-point(2))**2 if(d .lt. distmin) then distmin=d posres=intpar(i) end if end do else if(numintpt .gt. 0) then posres=intpar(1) end if distance=distmin if(numintcu.ge.1) then posres=0.5*(intcurve(1)%epar1(1)+intcurve(1)%epar1(2)) end if call s1221(b_curve%curve,0,posres,left,curvepos,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1227 for splineweight at ", point(1), point(2) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1227 for splineweight at ", point(1), point(2) distance=sqrt((curvepos(1)-point(1))**2+(curvepos(2)-point(2))**2) if (present(pos)) pos=posres END subroutine SUBROUTINE classify(x1, x2, cellk, spldom, wpredir, wpretot) real(kind=db), INTENT(IN):: x2(2), x1(2) type(cellkind), intent(INOUT):: cellk type(spline_domain)::spldom Real(kind=db):: zeval(4),reval(4), wpretot, wpredir real(kind=db), allocatable:: guess(:,:), w(:,:,:) Real(kind=db):: dmin, insidedir, insidetot, distance integer:: i,k allocate(guess(spldom%nbsplines,4)) allocate(w(0:2,spldom%nbsplines,4)) w=0 cellk%spldirkind=0 guess=-1.0_db dmin=HUGE(spldom%dist_extent) cellk%linkedboundaries=0 zeval=(/ x1(1),x1(2),x1(1),x1(2) /) reval=(/ x2(1),x2(1),x2(2),x2(2) /) insidedir=1 insidetot=1 do i=1,spldom%nbsplines do k=1,4 ! calculate the weight for each spline boundaries at each cell corner call splineweight(spldom%boundaries(i),zeval(k),reval(k),w(:,i,k),spldom%dist_extent,distance,guess(i,k)) ! We find the closest boundary to this point if(distance .lt. dmin) then ! If we are close enough we check if we are below dist_extent and need to calculate the distance each time if(distance .lt. spldom%dist_extent) then if(spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) then cellk%linkedboundaries(1)=i cellk%spldirkind=0 end if cellk%linkedboundaries(2)=i cellk%spltotkind=0 end if dmin=distance ! Otherwise we define the interior by the closest spline if(spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) then insidedir=w(0,i,k) end if insidetot=w(0,i,k) end if ! The neumann boundaries take precedence over the dirichlet boundaries ! this is important when they define what is outside of the simulation domain. if(spldom%boundaries(i)%type.eq. bd_Neumann.and. w(0,i,k).lt.0)then insidetot=w(0,i,k) if(distance.lt.spldom%dist_extent)then cellk%linkedboundaries(2) =i else cellk%linkedboundaries(2) =0 end if end if end do end do if(cellk%linkedboundaries(1) .gt. 0) then i=cellk%linkedboundaries(1) cellk%lguess(1)=minval(guess(i,:),1,guess(i,:).ge.0) cellk%rguess(1)=maxval(guess(i,:),1) wpredir=w(0,i,1) else cellk%spldirkind=sign(1,int(insidedir)) wpredir=insidedir end if if(cellk%linkedboundaries(2) .gt. 0) then i=cellk%linkedboundaries(2) wpretot=w(0,i,1) else cellk%spltotkind=sign(1,int(insidetot)) wpretot=insidetot end if end subroutine subroutine classifycells(spldom) use forSISL, ONLY: s1537, s1424 use bsplines type(spline_domain):: spldom integer:: i,j, dims(2), nbeval1, nbeval2,k,l real(kind=db)::val type(cellkind):: cellk real(kind=db), allocatable:: wpretot(:,:,:), wpredir(:,:,:), c(:,:), x1(:), x2(:) allocate(wpretot(1:1,0:spldom%nb1,0:spldom%nb2)) allocate(wpredir(1:1,0:spldom%nb1,0:spldom%nb2)) nbeval1=spldom%nb1+3 nbeval2=spldom%nb2+3 ! We set the interpolation points such that the spline interpolation of the weight uses the same knots as the spline interpolation of the electric potential allocate(x1(0:nbeval1-1),x2(0:nbeval2-1)) x1(0)=spldom%x1(0) x1(1)=(spldom%x1(0)+spldom%x1(1))/2.0_db j=0 do i=2,spldom%nb1 j=j+1 x1(i)=spldom%x1(j) end do x1(nbeval1-2)=(spldom%x1(spldom%nb1-1)+3*spldom%x1(spldom%nb1))/2.0_db x1(nbeval1-1)=spldom%x1(spldom%nb1) ! We do the same for x2 x2(0)=spldom%x2(0) x2(1)=(spldom%x2(0)+spldom%x2(1))/2.0_db j=0 do i=2,spldom%nb2 j=j+1 x2(i)=spldom%x2(j) end do x2(nbeval2-2)=(spldom%x2(spldom%nb2-1)+spldom%x2(spldom%nb2))/2.0_db x2(nbeval2-1)=spldom%x2(spldom%nb2) wpretot=0 wpredir=0 !$OMP PARALLEL DO private(i,j) do i=0,spldom%nb1-1 !DIR$ UNROLL do j=0,spldom%nb2-1 call classify(spldom%x1(i:i+1),spldom%x2(j:j+1),spldom%cellk(i,j),spldom, wpredir(1,i,j),wpretot(1,i,j)) end do end do !$OMP END PARALLEL DO deallocate(wpretot) deallocate(wpredir) allocate(wpretot(1:1,0:nbeval1-1,0:nbeval2-1)) allocate(wpredir(1:1,0:nbeval1-1,0:nbeval2-1)) !$OMP PARALLEL DO private(i,j,cellk,k,l) do i=0,nbeval1-1 call locintv(spldom%splrz%sp1,x1(i),k) do j=0,nbeval2-1 call locintv(spldom%splrz%sp2,x2(j),l) cellk=spldom%cellk(k,l) If(abs(cellk%spldirkind) .eq. 1) Then wpredir(1,i,j)=cellk%spldirkind else call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(i),x2(j), wpredir(:,i,j),spldom%dist_extent) end IF If(abs(cellk%spltotkind) .eq. 1) Then wpretot(1,i,j)=cellk%spltotkind else call splineweight(spldom%boundaries(cellk%linkedboundaries(2)), x1(i),x2(j), wpretot(:,i,j),spldom%dist_extent) end IF end do end do !$OMP END PARALLEL DO ! Set the approximated spline weight for the Dirichlet boundary conditions CALL set_splcoef((/3,3/),x1,x2,spldom%Dirdomweightspl) call get_dim(spldom%Dirdomweightspl,dims) allocate(c(dims(1),dims(2))) call get_splcoef(spldom%Dirdomweightspl, wpredir(1,:,:), c) CALL gridval(spldom%Dirdomweightspl,spldom%x1(1),spldom%x2(1), val ,(/0,0/),c) ! Set the approximated spline weight for the Neumann boundary conditions CALL set_splcoef((/3,3/),x1,x2,spldom%totdomweightspl) call get_splcoef(spldom%totdomweightspl, wpretot(1,:,:), c) CALL gridval(spldom%totdomweightspl,spldom%x1(1),spldom%x2(1), val ,(/0,0/),c) deallocate(c) end subroutine subroutine getindex(x1,x2,spldom, i, j) use distrib, ONLY: closest type(spline_domain):: spldom real(kind=db):: x1(:), x2(:) integer:: i(:),j(:) call locintv(spldom%splrz%sp1,x1, i) call locintv(spldom%splrz%sp2,x2, j) end subroutine SUBROUTINE speval(sp, xp, yp, leftx, lefty, f00, f10, f01) ! ! Compute the function f00 and its derivatives ! f10 = d/dx f ! f01 = d/dy f ! assuming that its PPFORM/BCOEFSC was already computed! ! TYPE(spline2d), INTENT(inout) :: sp DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp INTEGER, DIMENSION(:), INTENT(in) :: leftx, lefty DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: f00 DOUBLE PRECISION, DIMENSION(:), INTENT(out), OPTIONAL :: f10, f01 ! INTEGER :: np DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp)) INTEGER :: i, nidbas(2) DOUBLE PRECISION :: temp0(SIZE(xp),sp%sp2%order), temp1(SIZE(xp),sp%sp2%order) LOGICAL :: nlppform ! ! Apply periodicity if required ! np = SIZE(xp) nidbas(1) = sp%sp1%order-1 nidbas(2) = sp%sp2%order-1 nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform ! ! Locate the interval containing x, y ! x(:) = xp(:) - sp%sp1%knots(leftx(:)) y(:) = yp(:) - sp%sp2%knots(lefty(:)) ! ! Compute function/derivatives ! ! Using PPFORM !---------- DO i=1,np CALL my_ppval1(nidbas(1), x(i), sp%ppform(:,leftx(i)+1,:,lefty(i)+1), & & temp0(i,:), temp1(i,:)) END DO ! CALL my_ppval0(nidbas(2), y, temp0, 0, f00) if(present(f01))then CALL my_ppval0(nidbas(2), y, temp0, 1, f01) end if if(present(f10))then CALL my_ppval0(nidbas(2), y, temp1, 0, f10) end if !----------- CONTAINS !+++ SUBROUTINE my_ppval0(p, x, ppform, jder, f) ! ! Compute function and derivatives from the PP representation ! for many points x(:) INTEGER, INTENT(in) :: p DOUBLE PRECISION, INTENT(in) :: x(:) DOUBLE PRECISION, INTENT(in) :: ppform(:,:) INTEGER, INTENT(in) :: jder DOUBLE PRECISION, INTENT(out) :: f(:) DOUBLE PRECISION :: fact INTEGER :: j SELECT CASE (jder) CASE(0) ! function value SELECT CASE(p) CASE(1) f(:) = ppform(:,1) + x(:)*ppform(:,2) CASE(2) f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*ppform(:,3)) !!$ CASE(3) !!$ f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*(ppform(:,3)+x(:)*ppform(:,4))) CASE(3:) f(:) = ppform(:,p+1) DO j=p,1,-1 f(:) = f(:)*x(:) + ppform(:,j) END DO END SELECT CASE(1) ! 1st derivative SELECT CASE(p) CASE(1) f(:) = ppform(:,2) CASE(2) f(:) = ppform(:,2) + x(:)*2.d0*ppform(:,3) !!$ CASE(3) !!$ f(:) = ppform(:,2) + x(:)*(2.d0*ppform(:,3)+x(:)*3.0d0*ppform(:,4)) CASE(3:) f(:) = p*ppform(:,p+1) DO j=p-1,1,-1 f(:) = f(:)*x(:) + j*ppform(:,j+1) END DO END SELECT CASE default ! 2nd and higher derivatives f(:) = ppform(:,p+1) fact = p-jder DO j=p,jder+1,-1 f(:) = f(:)/fact*j*x(:) + ppform(:,j) fact = fact-1.0d0 END DO DO j=2,jder f(:) = f(:)*j END DO END SELECT END SUBROUTINE my_ppval0 !+++ SUBROUTINE my_ppval1(p, x, ppform, f0, f1) ! ! Compute function and first derivative from the PP representation INTEGER, INTENT(in) :: p DOUBLE PRECISION, INTENT(in) :: x DOUBLE PRECISION, INTENT(in) :: ppform(:,:) DOUBLE PRECISION, INTENT(out) :: f0(:) DOUBLE PRECISION, INTENT(out) :: f1(:) DOUBLE PRECISION :: fact INTEGER :: j SELECT CASE(p) CASE(1) f0(:) = ppform(1,:) + x*ppform(2,:) f1(:) = ppform(2,:) CASE(2) f0(:) = ppform(1,:) + x*(ppform(2,:)+x*ppform(3,:)) f1(:) = ppform(2,:) + x*2.d0*ppform(3,:) CASE(3) f0(:) = ppform(1,:) + x*(ppform(2,:)+x*(ppform(3,:)+x*ppform(4,:))) f1(:) = ppform(2,:) + x*(2.d0*ppform(3,:)+x*3.0d0*ppform(4,:)) CASE(4:) f0 = ppform(p+1,:) f1 = f0 DO j=p,2,-1 f0(:) = ppform(j,:) + x*f0(:) f1(:) = f0(:) + x*f1(:) END DO f0(:) = ppform(1,:) + x*f0(:) END SELECT END SUBROUTINE my_ppval1 !+++ END SUBROUTINE speval subroutine free_bsplinecurve(b_curve) type(spline_boundary):: b_curve call freeCurve(b_curve%curve) !call freeIntCurve(b_curve%intcurve) end subroutine END MODULE splinebound