radialgrid_dke_jd

 
 Generates a non uniform radial grid. 

   INPUTS: 

        - dke_mode : (1) resolution of the dke (0) of the FP equation
        - xpsin_S_dke : flux radial grid
        - dpsin_dke : maximum space between two consecutive points for DKE
        radial derivative calculations
        - psin_display_in: radial location at which data are displayed
        (graphic mode)

   OUTPUTS:

       - xpsin_f: positions where FP equation is to be solved [1,nr]
       - xpsin_f_dke: positions where DKE equation is to be solved
       [1,nr_dke]
       - xdpsin_f_dke: 
       - xpsinm_S_dke: 
       - xdpsinm_S_dke: 
       - xpsinp_S_dke: 
       - xdpsinp_S_dke: 
       - xpsin_S_dke: 
          - r_dke: index of the radial positions for DKE calculations
        - rho_display_out: radial location at which data are displayed
        (graphic mode)
        - ir_display_out: index of the radial location at which data are
        displayed (graphic mode)

 by Joan Decker <jodecker@alum.mit.edu> (MIT/RLE) and Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)

0001 function [radialDKE] = radialgrid_dke_jd(equil,dkeparam,display)
0002 %
0003 % Generates a non uniform radial grid.
0004 %
0005 %   INPUTS:
0006 %
0007 %        - dke_mode : (1) resolution of the dke (0) of the FP equation
0008 %        - xpsin_S_dke : flux radial grid
0009 %        - dpsin_dke : maximum space between two consecutive points for DKE
0010 %        radial derivative calculations
0011 %        - psin_display_in: radial location at which data are displayed
0012 %        (graphic mode)
0013 %
0014 %   OUTPUTS:
0015 %
0016 %       - xpsin_f: positions where FP equation is to be solved [1,nr]
0017 %       - xpsin_f_dke: positions where DKE equation is to be solved
0018 %       [1,nr_dke]
0019 %       - xdpsin_f_dke:
0020 %       - xpsinm_S_dke:
0021 %       - xdpsinm_S_dke:
0022 %       - xpsinp_S_dke:
0023 %       - xdpsinp_S_dke:
0024 %       - xpsin_S_dke:
0025 %          - r_dke: index of the radial positions for DKE calculations
0026 %        - rho_display_out: radial location at which data are displayed
0027 %        (graphic mode)
0028 %        - ir_display_out: index of the radial location at which data are
0029 %        displayed (graphic mode)
0030 %
0031 % by Joan Decker <jodecker@alum.mit.edu> (MIT/RLE) and Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)
0032 %
0033 %
0034 dke_mode = dkeparam.dke_mode;
0035 dpsin_dke = dkeparam.dpsin_dke;%for drift kinetic calculations
0036 %
0037 xpsin_S_dke = unique([dkeparam.psin_S,rho2psi_jd(equil,dkeparam.rho_S)]);
0038 %
0039 if isnan(display.psin_display), 
0040    psin_display_in=rho2psi_jd(equil,display.rho_display);
0041 else
0042    psin_display_in = display.psin_display;
0043 end
0044 %
0045 xpsin_S_dke = sort(xpsin_S_dke);
0046 nr_S = length(xpsin_S_dke);
0047 psinmin = 1e-5;%May be reduced but often need larger number of pitch-angle grid points for DKE (increasing nmhu may help)
0048 psinmax = 1-psinmin;
0049 %
0050 if nr_S == 1,%Case of local studies, no fast electron transport
0051     if xpsin_S_dke < psinmin, 
0052         xpsin_S_dke = psinmin;
0053         info_dke_yp(2,['WARNING: psin = 0 not allowed, psin = ',num2str(psinmin),' enforced !']);
0054     end
0055     if xpsin_S_dke > psinmax, 
0056         xpsin_S_dke = psinmax;
0057         info_dke_yp(2,['WARNING: psin = 0 not allowed, psin = ',num2str(psinmax),' enforced !']);
0058     end
0059 %
0060     xpsin_f_dke = xpsin_S_dke;
0061     xpsin_S_dke = [0,1];
0062     xdpsin_f_dke = NaN;
0063     xpsinm_S_dke = NaN;
0064     xdpsinm_S_dke = NaN;
0065     xpsinp_S_dke = NaN;
0066     xdpsinp_S_dke = NaN;
0067 %
0068     nr_f = 1;
0069 else%Case of global studies, fast electron transport allowed
0070     if xpsin_S_dke(1) ~= 0,%add extra point for the flux grid if not entered
0071         xpsin_S_dke = [0,xpsin_S_dke];
0072         nr_S = nr_S + 1;
0073     end
0074     if xpsin_S_dke(nr_S) ~= 1,%add extra point for the flux grid if not entered
0075         xpsin_S_dke = [xpsin_S_dke,1];
0076         nr_S = nr_S + 1;
0077     end
0078     maskmin = max(find((xpsin_S_dke(1) + xpsin_S_dke(2:nr_S))/2 <= psinmin)) + 1;
0079     if ~isempty(maskmin), %remove points too close from 0: ensures xpsin_f_dke >= psimin;
0080         xpsin_S_dke = [xpsin_S_dke(1), 2*psinmin - xpsin_S_dke(1), xpsin_S_dke(maskmin+1:nr_S)];
0081     end
0082     nr_S = length(xpsin_S_dke);
0083     maskmax = min(find((xpsin_S_dke(1:nr_S-1) + xpsin_S_dke(nr_S))/2 >= psinmax));
0084     if ~isempty(maskmax), %remove points too close from 1: ensures xpsin_f_dke <= psimax;
0085         xpsin_S_dke = [xpsin_S_dke(1:maskmax-1), 2*psinmax - xpsin_S_dke(nr_S), xpsin_S_dke(nr_S)];
0086     end  
0087     nr_S = length(xpsin_S_dke);
0088     nr_f = nr_S - 1;
0089 %
0090     xpsin_f_dke = (xpsin_S_dke(1:nr_S-1) + xpsin_S_dke(2:nr_S))/2;%point l+1/2 (f grid)
0091     xdpsin_f_dke = diff(xpsin_S_dke);
0092     xdpsin_S_dke = xpsin_f_dke(2:nr_f) - xpsin_f_dke(1:nr_f-1);
0093     xpsinm_S_dke = xpsin_S_dke(1:nr_S-1);%point i (flux grid)
0094     xdpsinm_S_dke = [2*xpsin_f_dke(1),xdpsin_S_dke];
0095     xpsinp_S_dke = xpsin_S_dke(2:nr_S);%point i+1 (flux grid)
0096     xdpsinp_S_dke = [xdpsin_S_dke,2*(1 - xpsin_f_dke(nr_f))];
0097     xdpsinm_S_dke(xdpsinm_S_dke == 0) = eps;
0098     xdpsinp_S_dke(xdpsinp_S_dke == 0) = eps;
0099 end
0100 %
0101 ir_dke_display_out = find(abs(xpsin_f_dke - psin_display_in) == min(abs(xpsin_f_dke - psin_display_in)));
0102 %
0103 % for DKE calculations, extra radial points are needed close to each
0104 %   position for the purpose of radial derivative
0105 %
0106 if dke_mode == 1%DKE calculations
0107     if (xpsin_f_dke(1) - dpsin_dke < psinmin) | (xpsin_f_dke(nr_f) + dpsin_dke > psinmax)
0108         dpsin_dke = psinmin/2;
0109     end
0110    % first left dke point
0111     xpsin_f = [xpsin_f_dke(1) - dpsin_dke,xpsin_f_dke(1)];
0112     ir_dke = 2;
0113     r_dke = ir_dke;
0114     if (nr_f > 1) & (xpsin_f_dke(2) > xpsin_f_dke(1) + dpsin_dke + eps)
0115         xpsin_f = [xpsin_f,xpsin_f_dke(1) + dpsin_dke];
0116         ir_dke = ir_dke + 1;
0117     end
0118     for ir_f = 2:nr_f
0119        % left dke point
0120         if max(xpsin_f) < xpsin_f_dke(ir_f) - dpsin_dke - eps,
0121             xpsin_f = [xpsin_f,xpsin_f_dke(ir_f) - dpsin_dke];
0122             ir_dke = ir_dke + 1;
0123         elseif max(r_dke) < ir_dke %if one and only one fp point between two dke point, it is centered
0124             xpsin_f(ir_dke) = (xpsin_f_dke(ir_f - 1) + xpsin_f_dke(ir_f))/2;
0125         end
0126        % central dke point
0127         xpsin_f = [xpsin_f,xpsin_f_dke(ir_f)];
0128         ir_dke = ir_dke + 1;
0129         r_dke = [r_dke,ir_dke];       
0130        % right dke point
0131         if (ir_f < nr_f) & (xpsin_f_dke(ir_f + 1) > xpsin_f_dke(ir_f) + dpsin_dke + eps) 
0132             xpsin_f = [xpsin_f,xpsin_f_dke(ir_f) + dpsin_dke];
0133             ir_dke = ir_dke + 1;
0134         end
0135     end
0136    % last right dke point
0137     xpsin_f = [xpsin_f,xpsin_f_dke(nr_f) + dpsin_dke];
0138     ir_dke = ir_dke + 1;
0139 else%FP calculations only
0140     xpsin_f = xpsin_f_dke;
0141     r_dke = cumsum(ones(1,nr_f)); 
0142 end
0143 %
0144 ir_display_out = r_dke(ir_dke_display_out);
0145 %
0146 radialDKE.xrho_S_dke = psi2rho_jd(equil,xpsin_S_dke);
0147 radialDKE.xpsin_f = xpsin_f;
0148 radialDKE.xpsin_f_dke = xpsin_f_dke;
0149 radialDKE.xdpsin_f_dke = xdpsin_f_dke;
0150 radialDKE.xpsinm_S_dke = xpsinm_S_dke;
0151 radialDKE.xdpsinm_S_dke = xdpsinm_S_dke;
0152 radialDKE.xpsinp_S_dke = xpsinp_S_dke;
0153 radialDKE.xdpsinp_S_dke = xdpsinp_S_dke;
0154 radialDKE.xpsin_S_dke = xpsin_S_dke;
0155 radialDKE.r_dke = r_dke;
0156 radialDKE.ir_display_out = ir_display_out;
0157 
0158 
0159 
0160