thetastar_calc_dke_en

PURPOSE ^

Calculation of the equilibrium dependent factor alpha of the bouncefactor C*alpha for runaway

SYNOPSIS ^

function [XXthetastar1,XXthetastar2,XXalpha] = thetastar_calc_dke_en(equilDKE,gridDKE,Zmomcoef,pnRA,dpnRA,npnRA, gammaRA)

DESCRIPTION ^

 Calculation of the equilibrium dependent factor alpha of the bouncefactor C*alpha for runaway
 avalanches, by applying a refined mhu grid, then redistributing it on the
 LUKE mhu grid.

 INPUT:

   - equil: magnetic equilibrium structure
   - gridDKE: grid structure
   - Zmomcoef: momentum dynamics structure

 OUTPUT:

   - XXalpha: poloidal angle value theta_star between 0 and pi
       [nr,np,nmhu] on each flux surface
   - XXthetastar1: poloidal angle value theta_star between 0 and pi
       [nr,np,nmhu] on each flux surface
   - XXthetastar2: poloidal angle value theta_star between pi and 2*pi
   -the sum sig [nr,npn,nmhu]

 By E. Nilsson (CEA-IRFM, emelie.nilsson@cea.fr) Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker
 (CEA-DRFC, joan.decker@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [XXthetastar1,XXthetastar2,XXalpha] = thetastar_calc_dke_en(equilDKE,gridDKE,Zmomcoef,pnRA,dpnRA,npnRA, gammaRA)
0002 % Calculation of the equilibrium dependent factor alpha of the bouncefactor C*alpha for runaway
0003 % avalanches, by applying a refined mhu grid, then redistributing it on the
0004 % LUKE mhu grid.
0005 %
0006 % INPUT:
0007 %
0008 %   - equil: magnetic equilibrium structure
0009 %   - gridDKE: grid structure
0010 %   - Zmomcoef: momentum dynamics structure
0011 %
0012 % OUTPUT:
0013 %
0014 %   - XXalpha: poloidal angle value theta_star between 0 and pi
0015 %       [nr,np,nmhu] on each flux surface
0016 %   - XXthetastar1: poloidal angle value theta_star between 0 and pi
0017 %       [nr,np,nmhu] on each flux surface
0018 %   - XXthetastar2: poloidal angle value theta_star between pi and 2*pi
0019 %   -the sum sig [nr,npn,nmhu]
0020 %
0021 % By E. Nilsson (CEA-IRFM, emelie.nilsson@cea.fr) Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker
0022 % (CEA-DRFC, joan.decker@cea.fr)
0023 %
0024 method='spline';%'linear';
0025 %
0026 npsi = gridDKE.nr;
0027 npn = gridDKE.npn;
0028 nmhu = gridDKE.nmhu;
0029 %
0030 theta=equilDKE.mtheta;
0031 ntheta = length(theta);
0032 xpsin = equilDKE.xpsin;
0033 ap=equilDKE.ap;
0034 Rp=equilDKE.Rp;Zp=equilDKE.Zp;
0035 B0=equilDKE.xB0;
0036 Bmax = max((sqrt(equilDKE.XBx.^2 + equilDKE.XBy.^2 + equilDKE.XBphi.^2)),[],2)';
0037 XXB0 = permute(repmat(B0',[1,npnRA,nmhu]), [2,3,1]);
0038 XXBmax = permute(repmat(Bmax',[1,npnRA,nmhu]), [2,3,1]);
0039 %
0040 mhu0min=repmat(sqrt(1-2./(gammaRA+1)),[npsi,1]);
0041 for ipsi=1:npsi
0042     mhu0max(ipsi,:)=sqrt(1-2*B0(ipsi)./(gammaRA+1)./Bmax(ipsi));
0043 end
0044 %
0045 for ipo = 2;
0046 mhuF = ipo; %factor of increase in mhu grid points >=1
0047 noP = mhuF*nmhu/2+1; %nr of cells on HALF the refined grid
0048 %
0049 gridmhu_S = NaN(npsi,npnRA,noP);
0050 XXthetastar1new = NaN(npsi,npnRA,mhuF*nmhu);
0051 XXthetastar2new = NaN(npsi,npnRA,mhuF*nmhu);
0052 XXalphaNew = NaN(npsi,npnRA,mhuF*nmhu);
0053 %
0054 for ipsi=1:npsi
0055     for ip=1:npnRA
0056         gridmhu_S(ipsi,ip,:)=linspace(mhu0min(ipsi,ip),mhu0max(ipsi,ip),noP); % refined grid cells
0057     end
0058 end
0059 %
0060 XXmhuip=(gridmhu_S(:,:,1:(end-1))+gridmhu_S(:,:,(2:end)))./2;
0061 dmhuip=gridmhu_S(:,:,(2:end))-gridmhu_S(:,:,1:(end-1));
0062 %
0063 XXmhuim = -flipdim(XXmhuip,3);
0064 dmhuim =  flipdim(dmhuip,3);
0065 %
0066 XXmhui = cat(3, XXmhuim, XXmhuip); %concatenate [-xi,0] and [0,xi]
0067 dmhui = cat(3, dmhuim, dmhuip);
0068 %
0069 XXmhuiper = permute(XXmhui,[2,3,1]); %permute for Bstar calc
0070 XXgammai=permute(repmat(gammaRA,[noP*2-2,1,npsi]),[2,1,3]);
0071 XXB0i=permute(repmat(B0',[1,npnRA,noP*2-2]),[2,3,1]);
0072 XXBstarnew = 2*XXB0i./(1 + XXgammai)./(1 - XXmhuiper.^2);
0073 %
0074 XXBstarnew=permute(XXBstarnew,[3,1,2]);
0075 %
0076 tic
0077 Br=equilDKE.XBx;
0078 Bz=equilDKE.XBy;
0079 Bphi=equilDKE.XBphi;
0080 for ipsi=1:npsi
0081     Bp=sqrt(Br(ipsi,:).^2+Bz(ipsi,:).^2);        %poloidal Bfield
0082     Bt=abs(Bphi(ipsi,:));                %toroidal Bfield
0083     B=sqrt(Bt.^2+Bp.^2);
0084     XXthetastar1new(ipsi,:,:)=interp1(B(1:(ntheta+1)/2),theta(1:(ntheta+1)/2),XXBstarnew(ipsi,:,:),method);     %calc theta_k* <pi
0085     XXthetastar2new(ipsi,:,:)=interp1(B((ntheta+1)/2:end),theta((ntheta+1)/2:end),XXBstarnew(ipsi,:,:),method);%calc theta_k* >pi
0086     %
0087     equilval1 = equilval_yp(equilDKE.equil_fit,equilDKE.xrho(ipsi),squeeze(XXthetastar1new(ipsi,:,:)));
0088     %
0089     x1 = equilval1.x;
0090     y1 = equilval1.y;
0091     Bx1 = equilval1.Bx;
0092     dBx1dtheta = equilval1.dBxdtheta;
0093     By1 = equilval1.By;
0094     dBy1dtheta = equilval1.dBydtheta;
0095     Bz1 = equilval1.Bz;
0096     dB1zdtheta = equilval1.dBzdtheta;
0097     Bp1 = equilval1.BP;
0098     %
0099     B1 = sqrt(Bx1.^2 + By1.^2 + Bz1.^2);
0100     dB1dtheta = (Bx1.*dBx1dtheta + By1.*dBy1dtheta + Bz1.*dB1zdtheta)./B1;
0101     %
0102     equilval2 = equilval_yp(equilDKE.equil_fit,equilDKE.xrho(ipsi),squeeze(XXthetastar2new(ipsi,:,:)));
0103     %
0104     x2 = equilval2.x;
0105     y2 = equilval2.y;
0106     Bx2 = equilval2.Bx;
0107     dBx2dtheta = equilval2.dBxdtheta;
0108     By2 = equilval2.By;
0109     dBy2dtheta = equilval2.dBydtheta;
0110     Bz2 = equilval2.Bz;
0111     dB2zdtheta = equilval2.dBzdtheta;
0112     Bp2 = equilval2.BP;
0113     %
0114     B2 = sqrt(Bx2.^2 + By2.^2 + Bz2.^2);
0115     dB2dtheta = (Bx2.*dBx2dtheta + By2.*dBy2dtheta + Bz2.*dB2zdtheta)./B2;
0116     r1Rp=sqrt(x1.^2+y1.^2)./equilDKE.ap; %LUKE definition with ap
0117     r2Rp=sqrt(x2.^2+y2.^2)./equilDKE.ap;
0118     XXalphaNew(ipsi,:,:) = B0(ipsi)^2.*r1Rp./abs(dB1dtheta)./Bp1 + B0(ipsi)^2.*r2Rp./abs(dB2dtheta)./Bp2;
0119 end
0120 %
0121 XXalpha=NaN(npsi,npnRA,nmhu);
0122 XXthetastar1=NaN(npsi,npnRA,nmhu);
0123 XXthetastar2=NaN(npsi,npnRA,nmhu);
0124 %
0125 yp=gridDKE.mhu(nmhu/2:end)+gridDKE.dmhup(nmhu/2:end)/2;
0126 ym=gridDKE.mhu(1:nmhu/2+1)-gridDKE.dmhup(1:nmhu/2+1)/2;
0127 %
0128 for ipsi = 1:npsi;
0129     for ip = 1:npnRA;
0130         xp = squeeze(gridmhu_S(ipsi,ip,:))';
0131         xm = -fliplr(xp);
0132         fp = XXalphaNew(ipsi,ip,noP:end);
0133         fm = XXalphaNew(ipsi,ip,1:noP-1);
0134         XXalpha(ipsi,ip,1:nmhu/2) = redistrib_jd(xm,fm,ym,1);
0135         XXalpha(ipsi,ip,nmhu/2+1:end) = redistrib_jd(xp,fp,yp,1);
0136         fpt1 = XXthetastar1new(ipsi,ip,noP:end);
0137         fmt1 = XXthetastar1new(ipsi,ip,1:noP-1);
0138         XXthetastar1(ipsi,ip,1:nmhu/2) = redistrib_jd(xm,fmt1,ym,0);
0139         XXthetastar1(ipsi,ip,nmhu/2+1:end) = redistrib_jd(xp,fpt1,yp,0);
0140         fpt2 = XXthetastar2new(ipsi,ip,noP:end);
0141         fmt2 = XXthetastar2new(ipsi,ip,1:noP-1);
0142         XXthetastar2(ipsi,ip,1:nmhu/2) = redistrib_jd(xm,fmt2,ym,0);
0143         XXthetastar2(ipsi,ip,nmhu/2+1:end) = redistrib_jd(xp,fpt2,yp,0);
0144     end
0145 end
0146 %
0147 dmhu = permute(repmat(gridDKE.dmhu',[1,1,npnRA,]), [2,3,1]);
0148 XXalpha = XXalpha;
0149 %
0150 t(ipo) = toc;
0151 %
0152 for ip = 1:npnRA
0153     ind = find(XXalpha(1,ip,1:nmhu/2));
0154     alphmin(ipo,ip) = XXalpha(1,ip,ind(1));
0155     alphmax(ipo,ip) = XXalpha(1,ip,ind(end));
0156 end
0157 end
0158 %plot the convergence test
0159 % figure; graph1D_jd((1:ipo)*nmhu,alphmin,0,0,'nr of points in refined grid','alpha(ximin,p=const)','',NaN,NaN,NaN,'-','.','b')
0160 % figure; graph1D_jd((1:ipo)*nmhu,alphmax,0,0,'nr of points in refined grid','alpha(ximax,p=const)','',NaN,NaN,NaN,'-','.','r')
0161 % figure; graph1D_jd((1:ipo)*nmhu,t,0,0,'nr of points in refined grid','time (s) for calc.','',NaN,NaN,NaN,'-','.','b')
0162 %
0163 %figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXalphaNew(1,:,:)),'p','mhu','alphaNew',NaN,NaN,0,NaN,100); colorbar;
0164 %
0165 figure
0166 hold on
0167 for ipsi=1:npsi
0168     for ip=1:npnRA
0169         XalphaNewINT(ip) = integral_dke_jd(dmhui(ipsi,ip,:),squeeze(XXalphaNew(ipsi,ip,:))',2); %integrate over mhu
0170     end
0171     graph1D_jd(pnRA,XalphaNewINT,0,0,'p','XalphaNew','',NaN,NaN,NaN,'-','.','b')
0172 end
0173 %
0174 figure;graph2D_jd(pnRA,gridDKE.mhu,squeeze(XXalpha(1,:,:)),'p','mhu','alpha LUKE',NaN,NaN,0,NaN,100); colorbar;
0175 figure;graph2D_jd(pnRA,gridDKE.mhu,squeeze(XXthetastar1(1,:,:)),'p','mhu','theta*1 LUKE',NaN,NaN,0,NaN,100); colorbar;
0176 figure;graph2D_jd(pnRA,gridDKE.mhu,squeeze(XXthetastar2(1,:,:)),'p','mhu','theta*2 LUKE',NaN,NaN,0,NaN,100); colorbar;
0177 %
0178 figure;
0179 hold on
0180 for ipsi=1:npsi
0181 Xalpha = squeeze(XXalpha(ipsi,:,:));
0182 Xalpha(isnan(Xalpha)) = 0;
0183 Xalpha_int = integral_dke_jd(gridDKE.dmhu,Xalpha,2);
0184 graph1D_jd(pnRA,Xalpha_int,0,0,'p','Xalpha LUKE','',NaN,NaN,NaN,'-','.','b')
0185 end
0186 hold off
0187 %
0188 % figure
0189 % hold on
0190 % for ip=1:npn
0191 % plot(plot(XXmhui(1,ip,:),squeeze(XXalphaNew(1,ip,:)),'o'))
0192 % plot(gridDKE.mhu,squeeze(XXalpha(1,ip,:)),'r*')
0193 % end
0194 
0195 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXthetastar1(1,:,:)),'p','mhu','theta^*_1',NaN,NaN,0,NaN,100); colorbar;
0196 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXthetastar2(1,:,:)),'p','mhu','theta^*_2',NaN,NaN,0,NaN,100); colorbar;
0197 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXalpha(1,:,:)),'p','mhu','XXalpha',NaN,NaN,0,NaN,100); colorbar;
0198 %
0199 % figure; hold on
0200 % for k=1:100
0201 %     plot(XXthetastar2(1,:,k),'b*')
0202 %  plot(XXthetastar2new(1,:,k),'r*')
0203 % end
0204 % figure; hold on
0205 % for k=1:50
0206 % graph1D_jd(XXmhui(k,:),squeeze(XXalphaNew(1,k,:)),0,0,'mhu','alpha','',NaN,NaN,NaN,'none','.','r')
0207 % end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.