thetastar_dke_interp_en

PURPOSE ^

Calculation of the critical poloidal angle value theta_star for runaway avalanches

SYNOPSIS ^

function [XXthetastar1,XXthetastar2,XXalpha] = thetastar_dke_interp_en(dkepath,dkeparam,dkedisplay,equilDKE,gridDKE,Zmomcoef,ohm,radialDKE)

DESCRIPTION ^

 Calculation of the critical poloidal angle value theta_star for runaway avalanches

 INPUT:

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

 OUTPUT:

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

 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_dke_interp_en(dkepath,dkeparam,dkedisplay,equilDKE,gridDKE,Zmomcoef,ohm,radialDKE)
0002 % Calculation of the critical poloidal angle value theta_star for runaway avalanches
0003 %
0004 % INPUT:
0005 %
0006 %   - equil: magnetic equilibrium structure
0007 %   - gridDKE: grid structure
0008 %   - Zmomcoef: momentum dynamics structure
0009 %   - radialDKE: radial grid structure
0010 %
0011 % OUTPUT:
0012 %
0013 %   - XXthetastar1: poloidal angle value theta_star between 0 and pi
0014 %       [np,nmhu,nr] on each flux surface ie [np,nmhu,length(B)]
0015 %   - XXthetastar2: poloidal angle value theta_star between pi and 2*pi
0016 %   -the sum sig [npn,nmhu,l]
0017 %
0018 % By E. Nilsson (CEA-IRFM, emelie.nilsson@cea.fr) Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker
0019 % (CEA-DRFC, joan.decker@cea.fr)
0020 %
0021 method='spline';%'linear';
0022 npn = size(gridDKE.XXmhu,1);
0023 nmhu = size(gridDKE.XXmhu,2);
0024 npsi=size(equilDKE.XBx,1); %nr of radial positions
0025 theta=equilDKE.mtheta;
0026 ntheta = length(theta);
0027 xpsin = equilDKE.xpsin;
0028 ap=equilDKE.ap;
0029 Rp=equilDKE.Rp;Zp=equilDKE.Zp;
0030 B0=equilDKE.xB0;
0031 Bmax = max((sqrt(equilDKE.XBx.^2 + equilDKE.XBy.^2 + equilDKE.XBphi.^2)),[],2)';
0032 XXB0 = permute(repmat(B0',[1,npn,nmhu]), [2,3,1]);
0033 XXBmax = permute(repmat(Bmax',[1,npn,nmhu]), [2,3,1]);
0034 %
0035 mhu0min=sqrt(1-2./(Zmomcoef.gamma+1));
0036 mhu0max=sqrt(1-2*B0./(Zmomcoef.gamma+1)./Bmax);
0037 %
0038 for ip=1:npn
0039     [dummy,I]=min(abs(gridDKE.mhu-mhu0min(ip)))
0040     mhu0minLUKE(ip)=gridDKE.mhu(I);
0041     min_index(ip)=I;
0042     [dummy,J]=min(abs(gridDKE.mhu-mhu0max(ip)))
0043     max_index(ip)=J;
0044     mhu0maxLUKE(ip)=gridDKE.mhu(J);
0045     gridmhu_S(ip,:)=linspace(mhu0min(ip),mhu0max(ip),nmhu/2+1); %cells
0046     %
0047     for im=1:(length(gridmhu_S)-1) %half grid
0048         XXmhuip(ip,im)=(gridmhu_S(ip,im)+gridmhu_S(ip,im+1))./2;
0049         dmhuip(ip,im)=gridmhu_S(ip,im+1)-gridmhu_S(ip,im);
0050     end
0051     %
0052 end
0053 %
0054 XXmhuim = -fliplr(XXmhuip);
0055 dmhuim = fliplr(dmhuip);
0056 %
0057 XXmhui = [XXmhuim,XXmhuip];
0058 dmhui = [dmhuim, dmhuip];
0059 XXBstarnew = 2*XXB0./(1 + Zmomcoef.XXgamma)./(1 - XXmhui.^2);
0060 %
0061 XXBstarnew=permute(XXBstarnew,[3,1,2]);
0062 %
0063 keyboard
0064 Br=equilDKE.XBx;
0065 Bz=equilDKE.XBy;
0066 Bphi=equilDKE.XBphi;
0067 for ipsi=1:npsi
0068     Bp=sqrt(Br(ipsi,:).^2+Bz(ipsi,:).^2);        %poloidal Bfield
0069     Bt=abs(Bphi(ipsi,:));                %toroidal Bfield
0070     B=sqrt(Bt.^2+Bp.^2);
0071     XXthetastar1new(ipsi,:,:)=interp1(B(1:(ntheta+1)/2),theta(1:(ntheta+1)/2),XXBstarnew(ipsi,:,:),method);     %calc theta_k* <pi
0072     XXthetastar2new(ipsi,:,:)=interp1(B((ntheta+1)/2:end),theta((ntheta+1)/2:end),XXBstarnew(ipsi,:,:),method);%calc theta_k* >pi
0073     [x1,y1,Bx1,By1,Bp1,dB1dtheta] = equilval_B_dB_en(equilDKE.equil_fit,equilDKE.xrho(ipsi),XXthetastar1new(ipsi,:,:));
0074     [x2,y2,Bx2,By2,Bp2,dB2dtheta] = equilval_B_dB_en(equilDKE.equil_fit,equilDKE.xrho(ipsi),XXthetastar2new(ipsi,:,:));
0075     XXalphaNew(ipsi,:,:)=B0(ipsi)^2.*equilDKE.xrho(ipsi)./abs(dB1dtheta)./Bp1 + B0(ipsi)^2.*equilDKE.xrho(ipsi)./abs(dB2dtheta)./Bp2;
0076 end
0077 %XalphaNewINT = integral_dke_jd(dmhui(25,:),squeeze(XXalphaNew(1,25,:))',2)./sum(dmhui(25,:)); %integrate over mhu, pos n neg : /2
0078 figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXalphaNew(1,:,:)),'p','mhu','alphaNew',NaN,NaN,0,NaN,100); colorbar;
0079 for ip=1:npn
0080     XalphaNewINT(ip) = integral_dke_jd(dmhui(ip,:),squeeze(XXalphaNew(1,ip,:))',2); %integrate over mhu, pos n neg : /2
0081 end
0082 figure
0083 graph1D_jd(gridDKE.pn,XalphaNewINT,0,0,'p','XalphaNew','',NaN,NaN,NaN,'-','.','b')
0084 %XalphaNewINT = integral_dke_jd(dmhui,squeeze(XXalphaNew(1,:,:)),2)./sum(dmhui(25,:)); %integrate over mhu, pos n neg : /2
0085 %
0086 b=1;
0087 
0088 ind=0;
0089 
0090 
0091 %%
0092 %put points back onto the LUKE grid, no 'cell on edge' effect taken into
0093 %account
0094 alphajj=NaN(1,npn,nmhu);
0095 
0096 for ip = 1:npn
0097     for ik=1:(max_index(ip)+1-min_index(ip))
0098         dxij = -gridDKE.dmhu(ip)/2 +gridDKE.dmhu(ip+ik)*(ik-1);
0099         dxij2 = gridDKE.dmhu(ip)/2 + gridDKE.dmhu(ip+ik)*(ik-1);
0100         dmhuif = ones(1,nmhu);
0101         dmhuif = dmhui(ip,:).*dmhuif;
0102         alphajj(1,ip,min_index(ip)+ik-1) = sum(XXalphaNew(1,ip,(gridDKE.mhu(min_index(ip))+dxij)<=XXmhui(ip,:) ...
0103             & XXmhui(ip,:) <=(gridDKE.mhu(min_index(ip))+dxij2))*dmhuif(min_index(ip)+ik-1))/gridDKE.dmhu(min_index(ip)+ik-1);
0104         alphajj(1,ip,nmhu-min_index(ip)+1-ik+1) = sum(XXalphaNew(1,ip,(gridDKE.mhu(nmhu-min_index(ip)+1)-dxij2)<=XXmhui(ip,:) ...
0105             & XXmhui(ip,:) <=(gridDKE.mhu(nmhu-min_index(ip)+1)-dxij))*dmhuif(nmhu-min_index(ip)+1))/gridDKE.dmhu(nmhu-min_index(ip)+1);
0106     end
0107 end
0108 
0109 %%
0110 %try to divide the cells on the edge correctly
0111 alphajj=NaN(1,npn,nmhu);
0112 b=1;
0113 ind=1;
0114 for ip = 1:npn
0115     for ik=1:(max_index(ip)+1-min_index(ip))
0116         dxij = -gridDKE.dmhu(ip)/2 +gridDKE.dmhu(ip+ik)*(ik-1);
0117         dxij2 = gridDKE.dmhu(ip)/2 + gridDKE.dmhu(ip+ik)*(ik-1);
0118         edge_test = (abs(mhu0minLUKE(ip)+dxij2)>(XXmhui(ip,:)-dmhui(ip,ik)/2)).*(abs(mhu0minLUKE(ip)+dxij2)<(XXmhui(ip,:)+dmhui(ip,ik)/2)); %check if a cell is in 2 LUKE cells
0119         dmhuif = ones(1,nmhu);
0120         dmhuif(1,ind) = abs(b);
0121         dmhuif(1,nmhu-ind+1) = abs(b);
0122         ind = find(edge_test>0);
0123         a = abs(mhu0minLUKE(ip)+dxij2-(XXmhui(ip,ind)+dmhui(ip,ind)/2))/dmhui(ip,ind); %fraction belonging to
0124         dmhuif(1,ind) = 1-a;
0125         dmhuif(1,nmhu-ind+1) = 1-a;
0126         dmhuif = dmhui(ip,:).*dmhuif;
0127         xip=(gridDKE.mhu(min_index(ip))+dxij)<=XXmhui(ip,:) & XXmhui(ip,:) <=(gridDKE.mhu(min_index(ip))+dxij2);
0128         xim=(gridDKE.mhu(nmhu-min_index(ip)+1)-dxij2)<=XXmhui(ip,:) & XXmhui(ip,:) <=(gridDKE.mhu(nmhu-min_index(ip)+1)-dxij);
0129        %
0130        alphajj(1,ip,min_index(ip)+ik-1) = sum(XXalphaNew(1,ip,(gridDKE.mhu(min_index(ip))+dxij)<=XXmhui(ip,:) & XXmhui(ip,:) <=(gridDKE.mhu(min_index(ip))+dxij2))*dmhuif(min_index(ip)+ik-1))/gridDKE.dmhu(min_index(ip)+ik-1);
0131        alphajj(1,ip,nmhu-min_index(ip)+1-ik+1) = sum(XXalphaNew(1,ip,(gridDKE.mhu(nmhu-min_index(ip)+1)-dxij2)<=XXmhui(ip,:) & XXmhui(ip,:) <=(gridDKE.mhu(nmhu-min_index(ip)+1)-dxij))*dmhuif(nmhu-min_index(ip)+1))/gridDKE.dmhu(nmhu-min_index(ip)+1);
0132        %
0133        % alphajj(1,ip,min_index(ip)+ik-1) = sum(squeeze(XXalphaNew(1,ip,xip)).*dmhuif(ik:ik+sum(xip>0)-1)')/gridDKE.dmhu(min_index(ip)+ik-1);
0134        % alphajj(1,ip,nmhu-min_index(ip)+1-ik+1) = sum(squeeze(XXalphaNew(1,ip,xim)).*dmhuif(ik:nmhu-min_index(ip)+1-ik+1+sum(xim>0)-1)')/gridDKE.dmhu(nmhu-min_index(ip)+1-ik+1);
0135         b =a;
0136     end
0137 end
0138 %
0139 figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(alphajj(1,:,:)),'p','mhu','alphajj',NaN,NaN,0,NaN,100); colorbar;
0140 %
0141 figure;
0142 Xalphajj = squeeze(alphajj(1,:,1:100));
0143 Xalphajj(isnan(Xalphajj)) = 0;
0144 Xalphajjint = integral_dke_jd(gridDKE.dmhu,Xalphajj,2);
0145 graph1D_jd(gridDKE.pn,Xalphajjint,0,0,'p','Xalphajj','',NaN,NaN,NaN,'-','.','b')
0146 legend(sprintf( '%s: %d', 'rho'));
0147 %nmhu/2+2+(nmhu/2-min_index(ip)-1)+ik-1
0148 %nmhu-min_index(ip)+1)
0149 %
0150 XXBstar = 2*XXB0./(1 + Zmomcoef.XXgamma)./(1 - gridDKE.XXmhu.^2);
0151 XXmaskBstar = XXBstar>=(XXB0) & XXBstar<(XXBmax);
0152 % XXmaskBstar2=XXBstar>=(XXB0) & XXBstar<(XXB0);
0153 % %
0154 % [I,J,K]=ind2sub(size(XXmaskBstar2), find(XXmaskBstar2>0));
0155 % XXBstar2=zeros(npn,nmhu,npsi);
0156 % for ik=1:length(I)
0157 %     XXBstar2(I(ik),J(ik),K(ik))=2*B0(K(ik))./(1 + (Zmomcoef.gamma(I(ik))-1*(gridDKE.dpn(I(ik))/10)*Zmomcoef.xz(I(ik),1)/(2*Zmomcoef.gamma(I(ik)))))./(1 - (gridDKE.mhu(J(ik))-0*sign(gridDKE.mhu(J(ik)).*gridDKE.dmhu(J(ik)+1)/2)).^2);
0158 % % g=Zmomcoef.gamma(I(ik))
0159 % % gcorr=Zmomcoef.gamma(I(ik))-1*(gridDKE.dpn(I(ik))/10)*Zmomcoef.xz(I(ik),1)/(2*Zmomcoef.gamma(I(ik)))
0160 % end
0161 %XXmaskBstar2 = double(XXmaskBstar2);
0162 %XXBstar2=permute(XXBstar2,[3,1,2]);
0163 XXmaskBstar = double(XXmaskBstar);
0164 %XXmaskBstar(XXmaskBstar==0) = NaN;
0165 XXBstar=(XXBstar.*XXmaskBstar);%+XXBstar2.*XXmaskBstar2);
0166 XXBstar(XXBstar==0) = NaN;
0167 XXBstar=permute(XXBstar,[3,1,2]);
0168 %
0169 XXthetastar1=NaN(npsi,npn,nmhu);
0170 XXthetastar2=NaN(npsi,npn,nmhu);
0171 XXalpha=NaN(npsi,npn,nmhu);
0172 tic
0173     Br=equilDKE.XBx;
0174     Bz=equilDKE.XBy;
0175     Bphi=equilDKE.XBphi;
0176     for ipsi=1:npsi
0177         Bp=sqrt(Br(ipsi,:).^2+Bz(ipsi,:).^2);        %poloidal Bfield
0178         Bt=abs(Bphi(ipsi,:));                %toroidal Bfield
0179         B=sqrt(Bt.^2+Bp.^2);
0180         XXthetastar1(ipsi,:,:)=interp1(B(1:(ntheta+1)/2),theta(1:(ntheta+1)/2),XXBstar(ipsi,:,:),method);     %calc theta_k* <pi
0181         XXthetastar2(ipsi,:,:)=interp1(B((ntheta+1)/2:end),theta((ntheta+1)/2:end),XXBstar(ipsi,:,:),method);%calc theta_k* >pi
0182         [x1,y1,Bx1,By1,Bp1,dB1dtheta] = equilval_B_dB_en(equilDKE.equil_fit,equilDKE.xrho(ipsi),XXthetastar1(ipsi,:,:));
0183         [x2,y2,Bx2,By2,Bp2,dB2dtheta] = equilval_B_dB_en(equilDKE.equil_fit,equilDKE.xrho(ipsi),XXthetastar2(ipsi,:,:));
0184         XXalpha(ipsi,:,:)=B0(ipsi)^2.*equilDKE.xrho(ipsi)./abs(dB1dtheta)./Bp1 + B0(ipsi)^2.*equilDKE.xrho(ipsi)./abs(dB2dtheta)./Bp2;
0185     end
0186 toc
0187 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXthetastar1(1,:,:)),'p','mhu','theta^*_1',NaN,NaN,0,NaN,100); colorbar;
0188 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXthetastar2(1,:,:)),'p','mhu','theta^*_2',NaN,NaN,0,NaN,100); colorbar;
0189 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXalpha(1,:,:)),'p','mhu','XXalpha',NaN,NaN,0,NaN,100); colorbar;
0190 %
0191 % figure; hold on
0192 % for k=1:100
0193 %     plot(XXthetastar2(1,:,k),'b*')
0194 %  plot(XXthetastar2new(1,:,k),'r*')
0195 % end
0196 figure; hold on
0197 for k=1:50
0198 graph1D_jd(gridDKE.mhu,squeeze(XXalpha(1,k,:)),0,0,'mhu','alpha','',NaN,NaN,NaN,'none','.','b')
0199 graph1D_jd(XXmhui(k,:),squeeze(XXalphaNew(1,k,:)),0,0,'mhu','alpha','',NaN,NaN,NaN,'none','.','r')
0200 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.