Generating ASCOT5 input entirely from TRANSP

[Ian-Dolby]

TRANSP outputs the plasma and field data that it uses internally during a given time window, in a "plasma state file", provided that the namelist setting FI_OUTTIM is set. The markers generated by TRANSP's NBI module, NUBEAM, are output in the marker birth file which is part of standard output.

It should be possible to use this to generate completely self-consistent input for ASCOT5 from TRANSP.
(currently trying to fork and add a feature branch to include existing scripts which do this)

[miekkasarki]

Hi @Ian-Dolby and thanks for providing these! We'll need to transition those from scripts to templates. I can help with that once I have the stable release 5.5 ready.

Did you have anything on the output available? We have a couple of Matlab scripts in ASCOT4 format and I could include those into this issue. (I can do the work of converting them to ASCOT5 python functions). Meanwhile I'll paste the Matlab scripts here:

transp2ascot4distribution.m

  function a=transp2ascot4distribution(transpCDF,ascotH5)
%
% We cheat by putting all the variuous RZ-points into a single column in 
%
%
   distName='rzPitchEdist';

   disp('Reading transp...')
   tra=read_transp_netcdf_aug(transpCDF);
   
   if(mod(numel(tra.grid.r),2) ~= 0)
      error('Our scheme to have two R columns and many lines failed!') 
   end
   
   disp('Converting...')
   d.version    = 3;
   d.ndims      = 6;
   d.nslots     = round([2 numel(tra.grid.r)/2 numel(tra.pitch) numel(tra.energy) 1 1]);
   d.vectlength = 1;
   
   ordinate = permute(tra.dist,[3 2 1]);
   d.ordinate(1,1,:,:,:) = ordinate(1:(numel(tra.grid.r)/2),    :,:);
   d.ordinate(1,2,:,:,:) = ordinate((numel(tra.grid.r)/2+1):end,:,:);
    
   d.abscissa_edges={               ...
       [1.0 2.0 3.0],                   ...
       linspace(1,2,d.nslots(2)+1), ...
       tra.pitch_boundary,          ...
       tra.energy_boundary,         ...
       [0.0 1.0],                   ...
       [0.5 1.5]                    };
   d.abscissa_units =  {'m',  'm',  'vpar/vtot',  'eV',     's',    ''};
   d.abscissa_names =  {'R',  'z',  'pitch',      'energy', 'time', 'species'};
   d.ordinate_name  =  {'density'};
   d.ordinate_unit  =  {'m^{-3}eV^{-1}'};

   
   a.species.testParticle.anum   = 2;
   a.species.testParticle.znum   = 1;
   a.species.testParticle.mass   = 1.6605e-27;
   a.species.testParticle.charge = 1.6022e-19;

   
   

   if(nargout > 0)
       a.dists.(distName) = d;
   end
   
   if(nargin > 1)
        disp('Writing...')
        write_ascot4_distribution(d,distName,ascotH5);
        write_ascot4_species(a.species.testParticle, ascotH5);
   end

end

read_transp_netcdf_aug.m

function tra=read_transp_netcdf_aug(filename)
% This file reads a nubeam file, such as: 29226A02_fi_1.cdf


tra.filemame=filename;
tra.grid.z    =ncread(filename,'Z2D')/100;
tra.grid.r    =ncread(filename,'R2D')/100;
tra.grid.theta=ncread(filename,'TH2D');
tra.grid.X    =ncread(filename,'X2D');
tra.grid.volumes =ncread(filename,'BMVOL');

tra.pitch     =ncread(filename,'A_D_NBI');
tra.pitch_boundary=ncread(filename,'AB_D_NBI'); 
tra.energy    =ncread(filename,'E_D_NBI'); %eV
tra.energy_boundary=ncread(filename,'EB_D_NBI'); 
tra.energy_units='eV';

tra.dist      =ncread(filename,'F_D_NBI') * 0.5 ;
% Caveat about the normalisation of FBM w.r.t. the pitch angle
%The TRANSP coordinate for the distribution function FBM is dA = dOmega/4pi as reported in the netCDF file label.
%Since there is no FBM dependence on theta, the integral over all pitch angles is
%2pi/4pi \int_0^2pi FBM*sin(theta) d(theta) 
%which is equivalent to consider FBM as the probability density with respect to the coordinate g = cos(theta)/2, i.e.
%df/dg = FBM
%Other codes use p = cos theta = 2*g as direct pitch angle coordinate. In that case their pitch angle probability density is:
%df/dp = df/d(2g) = 0.5 * df/dg = 0.5 * FBM 
% https://www.aug.ipp.mpg.de/aug/manuals/transp/fbm/pitch.html

% d(cos(theta))=sin(theta)d(theta)


tra.dist_units=ncreadatt(filename,'F_D_NBI','units');

tra.time=ncread(filename,'TIME');

tra.runid=ncread(filename,'TRANSP_RUNID');
disp(['Transp runid: "' tra.runid' '"'])


tra.n_fi=squeeze(sum(sum(tra.dist,2),1));
tra.n_fi=tra.n_fi*100^3; % 1/cm^3 => 1/m^3
tra.n_fi=tra.n_fi*(tra.pitch_boundary(2)-tra.pitch_boundary(1));
tra.n_fi=tra.n_fi*(tra.energy_boundary(2)-tra.energy_boundary(1));



r=linspace(min(tra.grid.r),max(tra.grid.r),30);
z=linspace(min(tra.grid.z),max(tra.grid.z),50);

d = griddata(tra.grid.r,tra.grid.z,tra.n_fi,r,z');
set(pcolor(r,z,d),'linestyle','none');
hold on;
scatter(tra.grid.r,tra.grid.z,50,tra.n_fi,'filled');
axis image
box on
xlabel('R (m)')
ylabel('z (m)')
cbh=colorbar;
%ylabel(cbh,'1/(m^3 d(omega/4pi))')
ylabel(cbh,'1/(m^3)')
title('Fast ion density')

colormap(jetW)

tra.volumeAverage = sum(tra.n_fi .* tra.grid.volumes) ./ sum(tra.grid.volumes);
fprintf('Volume averaged density is %e 1/m^3\n',tra.volumeAverage)

tra.vdist = sum(tra.dist.*permute(repmat(tra.grid.volumes,[1,75,50]),[2 3 1]),3);
end

[miekkasarki]

Actually it might be this script that converts the distribution properly:

transp_4dPE_dist_2_ascot_4dPE_dist.m

function [a,tra]=transp_4dPE_dist_2_ascot_4dPE_dist(transpCDF,ascotH5)

    %% PARAMETERS  
    
    constants
  
    
    nR       = 30;
    nZ       = 42;
  
    
    mass   = m_d;
    charge = e;
    anum   = 2;
    znum   = 1;
    origin = 1;
    
    %% DISTRIBUTION

   distName='rzPitchEdist';

   disp('Reading transp...')
   tra=read_transp_netcdf_aug(transpCDF);
   drawnow   
   
   disp('Initializing...')
   
   
    nE     = numel(tra.energy);
    nP     = numel(tra.pitch);
    nT       = 1;
    nSpecies = 1;
    nOrdDim  = 1;

     
   tra.grid.rmin=min(tra.grid.r);
   tra.grid.rmax=max(tra.grid.r);
   rmin =  tra.grid.rmin - 0.1*(tra.grid.rmax - tra.grid.rmin);
   rmax =  tra.grid.rmax + 0.1*(tra.grid.rmax - tra.grid.rmin);

   tra.grid.zmin=min(tra.grid.z);
   tra.grid.zmax=max(tra.grid.z);
   zmin =  tra.grid.zmin - 0.1*(tra.grid.zmax - tra.grid.zmin);
   zmax =  tra.grid.zmax + 0.1*(tra.grid.zmax - tra.grid.zmin);
   
   
   d.version    = 3;
   d.ndims      = 6;
   d.nslots     = round([nR nZ nP nE  nT nSpecies]');
   d.vectlength = nOrdDim;
   
    
   d.abscissa_edges={               ...
       linspace( rmin,rmax,d.nslots(1)+1),               ...
       linspace( zmin,zmax,d.nslots(2)+1),               ...
       tra.pitch_boundary,                               ...
       tra.energy_boundary * e,                          ...
       [0.0 1.0],                   ...
       [0.5 1.5]                    };
   d.abscissa_units =  {'m',  'm',  '',  'J',     's',        ''};
   d.abscissa_names =  {'R',  'z',  'pitch',  'energy',  'time', 'species'};
   d.ordinate_name  =  {'density'};
   d.ordinate_unit  =  {'s^3/m^6'};

   for iabs=1:d.ndims
    d.abscissa{iabs}=mean(...
        [ d.abscissa_edges{iabs}(1:(end-1)) ... 
          d.abscissa_edges{iabs}(2:end) ],2);
    d.abscissa_size{iabs}=...   
        diff(d.abscissa_edges{iabs});            
   end
   
   
   disp('Precalculating r,z weights for each transp slot for each ascot rz-gridpoint')
   
   r  = 0.5* ( d.abscissa_edges{1}(1:(end-1)) + d.abscissa_edges{1}(2:(end)));
   z  = 0.5* ( d.abscissa_edges{2}(1:(end-1)) + d.abscissa_edges{2}(2:(end)));
   
   nRZ=numel(tra.grid.r);
   w=zeros(nRZ,nR,nZ);
   for iRZ=1:nRZ;
       values=zeros(size(tra.grid.r));
       values(iRZ)=1.0;
       w(iRZ,:,:) = griddata(tra.grid.r,tra.grid.z,values,r,z')';
       if(false)
           clf
           pcolor(r,z,squeeze(w(iRZ,:,:))')
           hold on
           plot(tra.grid.r(iRZ),tra.grid.z(iRZ),'rx')
           caxis([0 1])
           colormap(jetW)
           drawnow
           pause(1)
       end
   end
   
   w(isnan(w)) = 0.0;
   
   for iR=nR:-1:1
       for iZ=nZ:-1:1
           rz_ind{iR,iZ}=find(w(:,iR,iZ) ~= 0);
           rz_wei{iR,iZ}=w( rz_ind{iR,iZ} ,iR,iZ);
       end
   end
   
   if(false)
       for iR=1:2:nR
           for iZ=1:7:nZ
               clf
               scatter(tra.grid.r,tra.grid.z,20,w(:,iR,iZ),'filled');
               hold on
               plot(r(iR),z(iZ),'rx');
               caxis([0 1])
               %           colormap(jetW)
               title(sum(w(:,iR,iZ)))
               drawnow          
               pause(1)
           end
       end
   end

   disp('Converting transp to ASCOT4 Pitch-Energy...')

   
   
   
   
   d.ordinate = zeros(d.nslots');
   
   
   
    for iZ=1:nZ
       fprintf('.');
        for iR=1:nR
%            for iE=1:nE
%                for iP=1:nP

            if(~isempty(rz_wei{iR,iZ}))
                rzw=repmat(reshape(rz_wei{iR,iZ},1,1,[]),[nE nP 1]);
                d.ordinate(    iR,iZ,:,:,1,1) = ...
                        sum( (tra.dist(:,:,rz_ind{iR,iZ})  ...
                        .* rzw),3 )' ;
            end
%                end
%            end
        end
    end
   fprintf('\n');

   
   
     
   
   
   %% Units and Jacobian
   
   J=1/(100^-3 .* e) ; % 1/(cm^3 eV) => 1/(m^3 eV)
   %keyboard
   

   
   
   d.ordinate = reshape(d.ordinate.* J,[1 size(d.ordinate) ]);
   
  
   
   
  
   
   
   %% TEST PARTICLE SPECIES   
   
   a.species.testParticle.anum   = anum;
   a.species.testParticle.znum   = znum;
   a.species.testParticle.mass   = mass;
   a.species.testParticle.charge = charge;
   a.species.testParticle.origin = origin;

   %% Also create vpa,vpe distribution
   
   
   a.dists.(distName) = d;
   disp('Converting ASCOT4 Pitch-Energy to ASCOT4 Vpa-Vpe...')

   so=convert_ascot4_4dDist_PE_2_VpaVpe(a);
   %so.dists
   
   %% OUTPUT VARIABLES

   if(nargout > 0)
       a.dists.rzVDist = so.dists.rzVDist;
   end
   
   if(nargin > 1)
        disp(['Writing... to file "' ascotH5  '".'])
        write_ascot4_distribution(d,distName,ascotH5);
        write_ascot4_distribution(so.dists.rzVDist,'rzVDist',ascotH5);
        write_ascot4_species(a.species.testParticle, ascotH5);
   end

   %   keyboard


   dv2=diff(d.abscissa_edges{3}(1:2)) *... 
       diff(d.abscissa_edges{4}(1:2));
   
   figure()
   plot2dgrid(d.abscissa_edges{1},d.abscissa_edges{2}, ...
              dv2*squeeze(sum(sum(d.ordinate,5),4)));
   colormap(jetW)
   colorbar
   axis image


end


function [w,ind]=peWeightsInd(x0,xvec)
  
    if(x0 < xvec(1))
        ind= 1;
        w  = 1;
        return
    end
    if(x0 > xvec(end))
        ind= numel(xvec);
        w  = 1;
        return
    end
        
    
    [~,in]=min(abs(xvec-x0));
    
    xn=xvec(in);
    
    if( xn == x0)
      ind = in;
      w   = 1.0;
      return
    end

    if(x0 > xn)
      ind = [ in in+1 ] ;
    else
      ind = [ in-1 in ] ;        
    end
    
    if(min(ind) == 0)
      keyboard  
    end
    
    w(2)  = (x0 - xvec(ind(1))) / ( xvec(ind(2)) - xvec(ind(1)))  ;
    w(1)  = 1.0 - w(2);
    
    
end

[Ian-Dolby]

Sure, I'll start changing them once I have fixed some known issues with the scripts (e.g some of the bulk neutrals are missed out).
No, I have not yet set up the output scripts in a way that is ready to be uploaded.

I notice that those Matlab scripts do not seem to automatically adjust for TRANSP's subtly different pitch convention, in which the pitch values are multiplied by -1 if the plasma current is in the opposite direction to the toroidal B field, which presumably does not matter for AUG. Maybe I missed where it does that, though.

[asnicker]

Hi Ian, I guess the main reason is indeed that in AUG (where these scripts have been mainly used) B and Ip are in the same direction.

[miekkasarki]

I pushed a tool to read TRANSP distributions to ASCOT5 format. The tool is in no way complete, mainly because it is missing:

• The sign of pitch issue Ian brought up
• I've no idea what the units in TRANSP are so it is likely missing some scaling factor
• It's now hard-coded to read Deuterium beams only - or at least that's what I think based on the names of the variables.

I'll address these once I have the data to make 1:1 comparisons between ASCOT5 and TRANSP. Meanwhile here's a plot with some test data showing the output of this function (top row is the converted ASCOT5 distribution and the bottom row is raw data from TRANSP).

By the way @Ian-Dolby, I noticed that in the matlab script the TRANSP distribution was multiplied with 0.5 with no comments explaining why. I recall that at the end of the training camp you got a match between ASCOT5 and TRANSP except for a factor of 2...

[Ian-Dolby]

The pitch sign issue can be sorted using the variable 'nsjdotb' in the FI output cdf file, which equals -1 if the plasma current and toroidal B field are anti-parallel. I did not know this variable existed in the FI output file.

There is a handy document containing more information about the FI distribution output from TRANSP: Hyun-Tae Kim | 2015 TRANSP users meeting | Culham | 22 Jan 2015 . It shows that the factor of 2 comes from a conversion of the distribution from solid angle to pitch, although I haven't gone through the calculation myself.

Regarding the discrepancy I mentioned, it turned out that TRANSP and ASCOT5 use different conditions for a marker colliding with the wall: in TRANSP, a marker is counted as lost if it comes within one Larmor radius of the wall, even when FLR effects are turned off. ASCOT in GC mode will have very few collisions, if any, with the wall, and this is similar for hybrid mode (from what I have seen).
This did not affect Andrea's similar work on MAST because MAST did not have a physical wall very close to the plasma, unlike NSTX. My results from TRANSP and ASCOT5 are now within 10% of each other.

[miekkasarki]

Thanks! I didn't know the presentation you mentioned existed but it seems really helpful.

10 % difference already sounds quite good to me. My guess would be that the codes use different values for the Coulomb logarithm, and that would largely explain the remaining difference.