momentumgrid_dke_jd

 
 Generates non uniform pitch angle and momentum grids 
 In the case of a non-uniform grid, the code ensures that a minimum number of points are present 
 in the trapped region and near mhu = 1. Also, generates fine grids near the trapped/passing boundary
 for each radial position. 

   INPUTS: 

     - uniform: (0) Non-uniform grid, (1) uniform grid
       - grid_power: Power for step variation. (grid_power > 1) NOTE: grid_power = 0 for old method with linear step variation 
    - dke_mode: Resolution of the (1) DKE equation ,(0) FP equation
    - r_dke: Indices of points for DKE calculation 
    - xmhu0T2: Trapped/passing boundaries in mhu space for all radial positions
    - nmhu_S_in: Number of angular values for the flux grid
    - np_S_in: Number of momentum values for the flux grid
    - pnmax_S: Maximum value of the flux momentum grid
    - sfac: reduction factors
    - snmhu: numbers of fine grid points in mhu
    - snp: reduction factor for the fine grid points in p
    - np_S_ref: position in the mpomentum grid where dpn_S is twice the minimum value

   OUTPUTS:

       - pn_f: momentum half-grid (distribution function)
       - dpn_f: 
       - pnm_S: 
       - dpnm_S: 
       - pnp_S: 
       - dpnp_S: 
       - pn_S: momentum full grid (flux)
    - mhu_f: pitch-angle half-grid (distribution function)
       - dmhu_f: 
       - mhum_S: 
       - dmhum_S: 
       - mhup_S: 
       - dmhup_S: 
    - mhu_S: pitch-angle full-grid (flux)

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

0001 function [momentumDKE] = momentumgrid_dke_jd(dkeparam,equilDKE,radialDKE);
0002 %
0003 % Generates non uniform pitch angle and momentum grids
0004 % In the case of a non-uniform grid, the code ensures that a minimum number of points are present
0005 % in the trapped region and near mhu = 1. Also, generates fine grids near the trapped/passing boundary
0006 % for each radial position.
0007 %
0008 %   INPUTS:
0009 %
0010 %     - uniform: (0) Non-uniform grid, (1) uniform grid
0011 %       - grid_power: Power for step variation. (grid_power > 1) NOTE: grid_power = 0 for old method with linear step variation
0012 %    - dke_mode: Resolution of the (1) DKE equation ,(0) FP equation
0013 %    - r_dke: Indices of points for DKE calculation
0014 %    - xmhu0T2: Trapped/passing boundaries in mhu space for all radial positions
0015 %    - nmhu_S_in: Number of angular values for the flux grid
0016 %    - np_S_in: Number of momentum values for the flux grid
0017 %    - pnmax_S: Maximum value of the flux momentum grid
0018 %    - sfac: reduction factors
0019 %    - snmhu: numbers of fine grid points in mhu
0020 %    - snp: reduction factor for the fine grid points in p
0021 %    - np_S_ref: position in the mpomentum grid where dpn_S is twice the minimum value
0022 %
0023 %   OUTPUTS:
0024 %
0025 %       - pn_f: momentum half-grid (distribution function)
0026 %       - dpn_f:
0027 %       - pnm_S:
0028 %       - dpnm_S:
0029 %       - pnp_S:
0030 %       - dpnp_S:
0031 %       - pn_S: momentum full grid (flux)
0032 %    - mhu_f: pitch-angle half-grid (distribution function)
0033 %       - dmhu_f:
0034 %       - mhum_S:
0035 %       - dmhum_S:
0036 %       - mhup_S:
0037 %       - dmhup_S:
0038 %    - mhu_S: pitch-angle full-grid (flux)
0039 %
0040 %  by Joan Decker <jodecker@mit.edu> (MIT/RLE) and Yves Peysson <yves.peysson@cea.fr> (CEA/DRFC)
0041 %
0042 if nargin < 3,
0043     error('Wrong number of input arguments in momentumgrid_dke_jd');
0044 end
0045 %
0046 dke_mode = dkeparam.dke_mode;
0047 bounce_mode = dkeparam.bounce_mode;
0048 uniform = dkeparam.grid_uniformity;
0049 grid_power = dkeparam.grid_power;
0050 grid_method = dkeparam.grid_method;
0051 nmhu_S_in = dkeparam.nmhu_S;
0052 np_S_in = dkeparam.np_S;
0053 pnmax_S = dkeparam.pnmax_S;
0054 snp = dkeparam.snp; 
0055 np_S_ref = dkeparam.np_S_ref;
0056 sfac0 = dkeparam.sfac0;% reduction factor for fine grid near mhu = 0
0057 sfac1 = dkeparam.sfac1; % reduction factor for fine grid near mhu = 1
0058 sfact = dkeparam.sfact; % reduction factor for fine grid near mhu = mhut
0059 snmhu0 = dkeparam.snmhu0;% number of fine grid points near mhu = 0
0060 snmhu1 = dkeparam.snmhu1;% number of fine grid points near mhu = 1
0061 snmhut = dkeparam.snmhut;% number of fine grid points near mhu = mhut
0062 %
0063 r_dke = radialDKE.r_dke;
0064 xmhu0T2 =  equilDKE.xmhu0T2;
0065 %
0066 if uniform == 1,%if uniform grid called, the old method is used for the mhu-grid
0067     %grid_power = 0;
0068 end
0069 %
0070 %Pitch-angle grid
0071 %
0072 if ~isfield(dkeparam,'mhu_S') || isempty(dkeparam.mhu_S),
0073     if grid_power == 0,
0074         nmhumin = nmhu_S_in;
0075         dmhu = 2/nmhumin; % grid step for uniform grid
0076     %
0077         sdmhu0 = dmhu/sfac0; % fine grid step near mhu = 0
0078         sdmhu1 = dmhu/sfac1; % fine grid step near mhu = 1
0079         sdmhut = dmhu/sfact; % fine grid step near mhu = mhut
0080     %
0081     % ---------- fine grid near mhu = 0 and  mhu = 1 ----------------
0082     %
0083         xmhut = sqrt(xmhu0T2);
0084         xmhut = xmhut(xmhut > 0 & xmhut < 1);
0085         mhutmin = min([xmhut,1]);
0086         mhutmax = max([0,xmhut]);
0087     %
0088     % ------ ensures that there are always enough points near mhu = 0 and  mhu = 1 --------
0089     %
0090         if ((1/4 + snmhu0)*sdmhu0) > mhutmin 
0091             sdmhu0 = mhutmin/(1/4 + snmhu0);
0092         else
0093             sdmhu0 = sdmhu0;
0094         end
0095     %
0096         if ((3/4 + snmhu1)*sdmhu1) > (1 - mhutmax)
0097             sdmhu1 = (1 - mhutmax)/(3/4 + snmhu1);
0098         else
0099             sdmhu1 = sdmhu1;
0100         end
0101     %
0102         smhu0min = 0;
0103         smhu0max = snmhu0*sdmhu0;
0104         smhu0 = [smhu0min:sdmhu0:smhu0max]; % fine grid at the mhu=0 vicinity
0105     %
0106         smhu1min = 1 - snmhu1*sdmhu1;
0107         smhu1max = 1;
0108         smhu1 = [smhu1min:sdmhu1:smhu1max]; % fine grid at the mhu=1 vicinity
0109     %
0110     % --------------- generates half mhu grid ---------------------------
0111     %
0112         if xrho_f == 0
0113             lmhu = trapeze_dke_jd(smhu0max,smhu1min,sdmhut,sfact);
0114             mhu = [smhu0,lmhu,smhu1];
0115         else
0116             mhu = smhu0;
0117             if dke_mode == 0
0118                 sdmhut = min([sdmhut,diff(xmhut)/3]);
0119             else
0120                 sdmhut = min([sdmhut,diff(xmhut)/6]);
0121             end        
0122             N = length(xmhut);
0123             i = 1;
0124             while ~isempty(xmhut(i:N))
0125                 if smhu0max > (xmhut(i) - (1/4 + snmhut)*sdmhut)
0126                     sdmhuti = sdmhu0;
0127                 elseif smhu1min < (xmhut(i) + (3/4 + snmhut)*sdmhut)
0128                     sdmhuti = sdmhu1;
0129                 else
0130                     sdmhuti = sdmhut;
0131                 end
0132     %
0133                 rmhumin = xmhut(i) - (1/4 + snmhut)*sdmhuti;
0134                 rmhumax =  xmhut(i) + (3/4 + snmhut)*sdmhuti;
0135     %
0136                 if i > 1 & max(mhu) > rmhumin % when two consecutive fine radial grid overlap
0137                     d = max(mhu) - rmhumin;
0138                     while d > 0
0139                         mhu = mhu(1:length(mhu)-1);
0140                         rmhumin = rmhumin + sdmhuti;
0141                         d = max(mhu) - rmhumin;
0142                     end
0143                     if (-d) > sdmhuti
0144                         lmhu = max(mhu) + (-d)/2;
0145                     else
0146                         lmhu = [];
0147                     end
0148                 else
0149                     mhu = mhu(mhu <= rmhumin);
0150                 lmhumin = max(mhu);
0151                 lmhumax = rmhumin;       
0152                 lmhu = trapeze_dke_jd(lmhumin,lmhumax,sdmhuti,dmhu);    
0153              end 
0154              rmhu = [rmhumin:sdmhuti:rmhumax];
0155              mhu = [mhu,lmhu,rmhu];
0156              i = i+1;
0157           end
0158     %
0159     % fine grid near mhu = 1
0160     %
0161           smhu1 = smhu1(smhu1 >= max(mhu));
0162           lmhumin = max(mhu);
0163           lmhumax = min(smhu1);  
0164           lmhu = trapeze_dke_jd(lmhumin,lmhumax,sdmhuti/2,dmhu);   
0165           mhu = [mhu,lmhu,smhu1];
0166        end
0167     %
0168        mhum = mhu(2:length(mhu));
0169        mhu = [mhu(1),mhum(diff(mhu) > eps)];
0170        mhu_S = [-fliplr(mhu(2:length(mhu))),mhu];
0171     %
0172     %   Case of an uniform mhu grid
0173     %
0174        if uniform == 1,
0175         mhu_S = [-1:2/(nmhumin+1):1];
0176        end
0177     %
0178     % calculates the radial positions for dke solutions
0179     %
0180        if sum(xrho_f == 0) > 0
0181         dke_mode = 0
0182        end
0183        if dke_mode == 1
0184         i_dke = [];
0185         i_rho = [];
0186         r_dke =[];
0187         for i = 1:length(xrho_f),
0188             i_mhu = max(find(mhu < xmhut(i)));
0189             i_dke = [i_dke, i_mhu];
0190             if i > 1 & i_mhu == max(i_rho)
0191                 r_dke = [r_dke, length(i_rho)];
0192                 i_rho = [i_rho, i_mhu + 1];
0193              elseif  i > 1 & i_mhu == max(i_rho) + 1
0194                 r_dke = [r_dke, length(i_rho) + 1];
0195                 i_rho = [i_rho, i_mhu, i_mhu + 1];
0196              else
0197                 r_dke = [r_dke, length(i_rho) + 2];
0198                 i_rho = [i_rho, i_mhu - 1, i_mhu, i_mhu + 1];
0199              end            
0200           end
0201           xrho_dke = 1./(2./(3/4*mhu(i_rho)+1/4*mhu(i_rho+1)).^2 - 1)*Rp_norm;
0202           xrho_S = NaN*xrho_S;
0203        else 
0204         r_dke = [1:1:length(xrho_f)];
0205         xrho_dke = xrho_f;
0206        end
0207     %
0208        xrho_f = xrho_dke;%Final grid with extra points for solving the drift kinetic equation
0209     %
0210        if nr_S > 1,
0211         ir_display_out = find(min(abs(xrho_f(r_dke) - rho_display_in)) == abs(xrho_f(r_dke) - rho_display_in));
0212         if length(ir_display_out) == 2, ir_display_out = min(ir_display_out); end
0213        else
0214         ir_display_out = 1;
0215        end
0216     else
0217     %
0218     %New optimized method for reducing the number of pitch-angle points.
0219     %
0220         nmhumin = nmhu_S_in;
0221         dmhu = 2/(nmhumin-1); % grid step for uniform grid
0222     %
0223         sdmhu0 = dmhu/sfac0; % fine grid step near mhu = 0
0224         sdmhu1 = dmhu/sfac1; % fine grid step near mhu = 1
0225         sdmhut = dmhu/sfact; % fine grid step near mhut
0226     %
0227         xmhut = sqrt(xmhu0T2);
0228         N = length(xmhut);
0229     %
0230         if bounce_mode > 0,
0231             mhu = makegrid_dke_jd(0,snmhu0,sdmhu0,xmhut(1),snmhut,sdmhut,snmhu0,dmhu,grid_power);
0232             for ix=1:N-1,
0233                 %sdmhutmin = diff(mhu(length(mhu) - 1:length(mhu)));
0234                 %snmhut_new = max([snmhut - fix(log(sdmhut/sdmhutmin)/log(grid_power)),1]);
0235                 mhu = [mhu,makegrid_dke_jd(xmhut(ix),snmhut,sdmhut,xmhut(ix+1),snmhut,sdmhut,1,dmhu,grid_power)];
0236             end
0237             %sdmhutmin = diff(mhu(length(mhu) - 1:length(mhu)));
0238             mhu = [mhu,makegrid_dke_jd(xmhut(N),snmhut,sdmhut,1,snmhu1,sdmhu1,snmhu1,dmhu,grid_power)];
0239         else
0240             mhu = makegrid_dke_jd(0,snmhu0,sdmhu0,1,snmhu1,sdmhu1,snmhu0+snmhu1,dmhu,grid_power);
0241         end
0242     %
0243         mhu = sort(mhu(find(diff([-1,mhu]) > 10*eps)));%Sorting for RFP since the trapped domain is not monotonic with the minor radius
0244         mhu_S = [-fliplr(mhu(2:length(mhu))),mhu];
0245     %
0246     end
0247     %
0248 else% pitch angle grid specified in dkeparam
0249     %
0250     mhu_S = dkeparam.mhu_S;
0251     %
0252 end
0253 %
0254 nmhu_S = length(mhu_S);
0255 nmhu_f = nmhu_S - 1;
0256 %
0257 % added to avoid imaginary numbers later
0258 %
0259 mhu_S(1) = -1;
0260 mhu_S(nmhu_S) = 1;
0261 %
0262 mhum_S = mhu_S(1:nmhu_S-1);%point j (flux grid)
0263 mhup_S = mhu_S(2:nmhu_S);%point j+1 (flux grid)
0264 mhum_S(mhum_S == 0) = eps;%Avoid singularities (NaN)
0265 mhup_S(mhup_S == 0) = eps;%Avoid singularities (NaN)
0266 mhu_f = (mhum_S + mhup_S)/2;%point j+1/2 (f grid)
0267 %
0268 dmhu_f = diff(mhu_S);%point j+1/2 (f grid)
0269 dmhu_S = mhu_f(2:nmhu_f) - mhu_f(1:nmhu_f-1);
0270 dmhum_S = [2*(mhu_f(1)+1),dmhu_S];%point j (flux grid)
0271 dmhup_S = [dmhu_S,2*(1-mhu_f(nmhu_f))];%point j+1 (flux grid)
0272 %
0273 dmhum_S(dmhum_S == 0) = eps;%Avoid singularities (NaN)
0274 dmhup_S(dmhup_S == 0) = eps;%Avoid singularities (NaN)
0275 %
0276 %Momentum grid
0277 %
0278 dpn_threshold = 0.1;
0279 %
0280 if ~isfield(dkeparam,'pn_S') || isempty(dkeparam.pn_S),
0281     if uniform == 1,
0282         pgridfactor = 1;%Uniform grid: pgridfactor = 1
0283         np_S_ref = 20;
0284     else
0285         pgridfactor = snp;%Non-uniform grid pgridfactor > 1
0286     end
0287     %
0288     %if np_S_in < nmhu_S,
0289     %    np_S = nmhu_S;
0290     %else
0291         np_S = np_S_in;
0292     %end
0293     %
0294     pn_S(1) = 0;
0295     dpn_S_ref = pnmax_S(1)/(np_S(1) - 1);
0296     dpn_S_0 = dpn_S_ref/pgridfactor;%Momentum step close to p = 0
0297     ip = 1;
0298     while pn_S(ip) <= pnmax_S(1) + dpn_S_ref*1e-3,%Small steps close to p = 0
0299     %   dpn_S(ip) = (dpn_S_ref-dpn_S_0)*tanh((np_S(1) - ip - np_S_ref)/25)/2 + (dpn_S_ref+dpn_S_0)/2;
0300        dpn_S(ip) = (dpn_S_ref-dpn_S_0)*tanh((ip-np_S_ref))/2 + (dpn_S_ref+dpn_S_0)/2;
0301         pn_S(ip+1) = pn_S(ip) + dpn_S(ip);
0302         ip = ip+1;
0303     end
0304     pn_S = pn_S(1:length(pn_S) - 1);
0305     %
0306     if length(np_S) > 1 && length(pnmax_S) > 1,% high energy sparse grid (runaway studies)
0307         if imag(pnmax_S(2)) == 0,
0308             pn_S_tail = linspace(pnmax_S(1),pnmax_S(2),np_S(2));
0309         else
0310             pn_S_tail = linspace(pnmax_S(1),imag(pnmax_S(2))*pnmax_S(1),np_S(2));
0311         end
0312         %
0313         pn_S = unique([pn_S,pn_S_tail]);
0314     end
0315     %
0316 else% momentum grid specified in dkeparam
0317     %
0318     pn_S = dkeparam.pn_S;
0319     %
0320 end
0321 %
0322 npn_S = length(pn_S);
0323 %
0324 npn_f = npn_S - 1;
0325 pn_f = (pn_S(1:npn_S-1) + pn_S(2:npn_S))/2;%point i+1/2 (f grid)
0326 dpn_f = diff(pn_S);%point i+1/2 (f grid)
0327 %
0328 dpn_S = pn_f(2:npn_f) - pn_f(1:npn_f-1);
0329 pnm_S = pn_S(1:npn_S-1);%point i (flux grid)
0330 pnm_S(pnm_S == 0) = eps;%Avoid singularities (NaN)
0331 pnp_S = pn_S(2:npn_S);%point i+1 (flux grid)
0332 pnp_S(pnp_S == 0) = eps;%Avoid singularities (NaN)
0333 dpnm_S = [2*pn_f(1),dpn_S];
0334 dpnp_S = [dpn_S,2*(pn_S(npn_S) - pn_f(npn_f))];
0335 %
0336 testdpn_f = find(dpn_f > dpn_threshold);
0337 %if isempty(testdpn_f) == 0, info_dke_yp(2,'WARNING: momentum step too large, results may be inaccurate !');end
0338 %
0339 if abs(uniform) == 1,%for test only with uniform grid
0340     mhu_S_pos = linspace(0,1,floor((nmhu_S_in + 1)/2));
0341     mhu_S = [-fliplr(mhu_S_pos(2:length(mhu_S_pos))),mhu_S_pos];
0342     nmhu_S = length(mhu_S);
0343      mhu_f = (mhu_S(1:nmhu_S-1) + mhu_S(2:nmhu_S))/2;
0344     nmhu_f = length(mhu_f);
0345     dmhu_f = diff(mhu_S);
0346     dmhu_S = [2*(mhu_f(1)+1),diff(mhu_f),2*(1-mhu_f(nmhu_f))];
0347     mhum_S = mhu_S(1:nmhu_f);
0348     mhup_S = mhu_S(2:nmhu_f+1);
0349     dmhum_S = dmhu_S(1:nmhu_f);
0350     dmhup_S = dmhu_S(2:nmhu_f+1);
0351 end
0352 %
0353 momentumDKE.pn = pn_f;
0354 momentumDKE.dpn = dpn_f;
0355 momentumDKE.pnm = pnm_S;
0356 momentumDKE.dpnm = dpnm_S;
0357 momentumDKE.pnp = pnp_S;
0358 momentumDKE.dpnp = dpnp_S;
0359 momentumDKE.pn_S = pn_S;
0360 momentumDKE.mhu = mhu_f;
0361 momentumDKE.dmhu = dmhu_f;
0362 momentumDKE.mhum = mhum_S;
0363 momentumDKE.dmhum = dmhum_S;
0364 momentumDKE.mhup = mhup_S;
0365 momentumDKE.dmhup = dmhup_S;
0366 momentumDKE.mhu_S = mhu_S;
0367 
0368 
0369