thetastar_dke_en

PURPOSE ^

SYNOPSIS ^

function [XXthetastar1,XXthetastar2,XXalpha] = thetastar_dke_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:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [XXthetastar1,XXthetastar2,XXalpha] = thetastar_dke_en(dkepath,dkeparam,dkedisplay,equilDKE,gridDKE,Zmomcoef,ohm,radialDKE)
0002 %
0003 % Calculation of the critical poloidal angle value theta_star for runaway avalanches
0004 %
0005 % INPUT:
0006 %
0007 %   - equil: magnetic equilibrium structure
0008 %   - gridDKE: grid structure
0009 %   - Zmomcoef: momentum dynamics structure
0010 %   - radialDKE: radial grid structure
0011 %
0012 % OUTPUT:
0013 %
0014 %   - XXthetastar1: poloidal angle value theta_star between 0 and pi
0015 %       [np,nmhu,nr] on each flux surface ie [np,nmhu,length(B)]
0016 %   - XXthetastar2: poloidal angle value theta_star between pi and 2*pi
0017 %   -the sum sig [npn,nmhu,l]
0018 %
0019 % By E. Nilsson (CEA-IRFM, emelie.nilsson@cea.fr) Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker
0020 % (CEA-DRFC, joan.decker@cea.fr)
0021 %
0022 ntheta = length(equilDKE.mtheta);
0023 Bmax=max(sqrt(sqrt(equilDKE.XBx'.^2+equilDKE.XBy'.^2)+equilDKE.XBphi'.^2));
0024 %
0025 jj = 1;
0026 %
0027 for ir = 1:gridDKE.nr,
0028     for ipn = 1:gridDKE.npn,
0029         for imhu = 1:gridDKE.nmhu,
0030             grid_index(jj).ir = ir;
0031             grid_index(jj).ipn = ipn;
0032             grid_index(jj).imhu = imhu;
0033             jj = jj + 1;
0034         end
0035     end
0036 end
0037 %
0038 mdce_mode = 3;%3 for parfor, i or 3+i for GPU
0039 mdce_remote = 0;%0 local, x remotely (see dkepath structure)
0040 mdce_enforce = 1;
0041 mdce_display = 1;
0042 %
0043 %
0044 clustermode.loop_thetastar_dke_en.scheduler.mode = real(mdce_mode);
0045 if ~isreal(mdce_mode),
0046     clustermode.loop_thetastar_dke_en.scheduler.gpu = imag(mdce_mode);
0047 else
0048     clustermode.loop_thetastar_dke_en.scheduler.gpu = 0;
0049 end
0050 %
0051 clustermode.loop_thetastar_dke_en.scheduler.enforce = mdce_enforce;
0052 %
0053 clustermode.ftest_mdce.scheduler.display = mdce_display;
0054 %
0055 [dkecluster] = clustermode_luke(clustermode,'loop_thetastar_dke_en',dkepath,mdce_remote);
0056 %
0057 tic
0058 if isfield(dkeparam,'options_fzero'),
0059     options_fzero = dkeparam.options_fzero;
0060 else
0061     options_fzero = [];
0062 end
0063 %
0064 %for igrid = 1:length(grid_index),
0065 %    [YYalpha(igrid),YYthetastar1(igrid),YYthetastar2(igrid)] = loop_thetastar_dke_en(igrid,grid_index,equilDKE.xrho,equilDKE.equil_fit,equilDKE.xB0,Bmax,gridDKE.mhu,Zmomcoef.gamma,options_fzero);
0066 %    igrid
0067 %end
0068 %
0069 [flag,YYalpha,YYthetastar1,YYthetastar2] = mdce_run(@loop_thetastar_dke_en,{NaN,grid_index,equilDKE.xrho,equilDKE.equil_fit,equilDKE.xB0,Bmax,gridDKE.mhu,Zmomcoef.gamma,options_fzero},1,1:length(grid_index),dkecluster);
0070 %
0071 toc
0072 %
0073 for ir = 1:gridDKE.nr,
0074     XXalpha(ir,:,:) = reshape(YYalpha(gridDKE.npn*gridDKE.nmhu*(ir-1)+1:gridDKE.npn*gridDKE.nmhu*ir),gridDKE.nmhu,gridDKE.npn)';
0075     XXthetastar1(ir,:,:) = reshape(YYthetastar1(gridDKE.npn*gridDKE.nmhu*(ir-1)+1:gridDKE.npn*gridDKE.nmhu*ir),gridDKE.nmhu,gridDKE.npn)';
0076     XXthetastar2(ir,:,:) = reshape(YYthetastar2(gridDKE.npn*gridDKE.nmhu*(ir-1)+1:gridDKE.npn*gridDKE.nmhu*ir),gridDKE.nmhu,gridDKE.npn)';
0077 end
0078 %
0079 if iscell(XXalpha) && iscell(XXthetastar1) && iscell(XXthetastar2)
0080     XXalpha = cell2mat(XXalpha);
0081     XXthetastar1 = cell2mat(XXthetastar1);
0082     XXthetastar2 = cell2mat(XXthetastar2);
0083 end
0084 %
0085 %keyboard
0086 %check discrepancies between thetastar1 and 2
0087 for ir=1:gridDKE.nr
0088 a=isnan(XXthetastar1(ir,:,:));
0089 b=isnan(XXthetastar2(ir,:,:));
0090 a=squeeze(a);
0091 b=squeeze(b);
0092 [I,J]=find(abs(a-b)>0)
0093 end
0094 % for k=1:length(I)
0095 % k1(k)=XXthetastar1(1,I(k),J(k));
0096 % k2(k)=XXthetastar2(1,I(k),J(k));
0097 % f(k)=2.0/(1 - gridDKE.XXmhu(1,I(k)).^2)./(Zmomcoef.gamma(J(k))+ 1.0);
0098 % end
0099 % %
0100 % for k=1:601
0101 % theta(k)=(k-1)*2*pi/600;
0102 % equilval1 = equilval_yp(equilDKE.equil_fit,0.4,theta(k));
0103 % Bx1 = equilval1.Bx;
0104 % By1 = equilval1.By;
0105 % B = equilval1.B;
0106 % dB1dtheta = equilval1.dBdtheta;
0107 % Bstar(k)=B/equilDKE.xB0;
0108 % dBdth(k)=dB1dtheta;
0109 % Bfield(k)=B;
0110 %
0111 % end
0112 
0113 %
0114 %figure
0115 %hold on
0116 %plot(theta,Bstar,'r')
0117 % if(~isempty(k2))
0118 % plot(k2,f,'ko')
0119 % end
0120 % %
0121 %
0122  
0123 
0124 
0125 
0126 
0127 if 0,
0128 
0129 %
0130 for ipsi = 1:npsi,
0131     for ipn = 1:npn,
0132         for imhu = 1:nmhu,
0133             ymin = fthetastar_dke_en(0,equilDKE.xrho(ipsi),equilDKE.equil_fit,equilDKE.xB0(ipsi),gridDKE.mhu(imhu),Zmomcoef.gamma(ipn));
0134             ymax = fthetastar_dke_en(pi,equilDKE.xrho(ipsi),equilDKE.equil_fit,equilDKE.xB0(ipsi),gridDKE.mhu(imhu),Zmomcoef.gamma(ipn));
0135             if ymin*ymax <= 0,
0136                 thetastar1 = fzero(@(theta) fthetastar_dke_en(theta,equilDKE.xrho(ipsi),equilDKE.equil_fit,equilDKE.xB0(ipsi),gridDKE.mhu(imhu),Zmomcoef.gamma(ipn)),pi);
0137                 thetastar2 = fzero(@(theta) fthetastar_dke_en(theta,equilDKE.xrho(ipsi),equilDKE.equil_fit,equilDKE.xB0(ipsi),gridDKE.mhu(imhu),Zmomcoef.gamma(ipn)),pi/2 + pi);
0138                 %
0139                 if thetastar1 >= 0 && thetastar1 <= 2*pi,
0140                     XXthetastar1(ipn,imhu,ipsi) = thetastar1;
0141                     %
0142                     equilval1 = equilval_yp(equilDKE.equil_fit,equilDKE.xrho(ipsi),thetastar1);
0143                     %
0144                     x1 = equilval1.x;
0145                     y1 = equilval1.y;
0146                     Bx1 = equilval1.Bx;
0147                     By1 = equilval1.By;
0148                     dB1dtheta = equilval1.dBdtheta;
0149                     %
0150                 end
0151                 %
0152                 if thetastar2 >= 0 && thetastar2 <= 2*pi,
0153                     XXthetastar2(ipn,imhu,ipsi) = thetastar2;
0154                     %
0155                     equilval2 = equilval_yp(equilDKE.equil_fit,equilDKE.xrho(ipsi),thetastar2);
0156                     %
0157                     x2 = equilval2.x;
0158                     y2 = equilval2.y;
0159                     Bx2 = equilval2.Bx;
0160                     By2 = equilval2.By;
0161                     dB2dtheta = equilval2.dBdtheta;                    
0162                     %
0163                 end
0164                 %
0165                 if ~isnan(XXthetastar1(ipn,imhu,ipsi)) && ~isnan(XXthetastar2(ipn,imhu,ipsi)),
0166                     XXalpha(ipn,imhu,ipsi) = equilDKE.xB0(ipsi)*((x1*x1+y1*y1)/abs(-By1*x1+Bx1*y1)/abs(dB1dtheta) + (x2*x2+y2*y2)/abs(-By2*x2+Bx2*y2)/abs(dB2dtheta));
0167                 end
0168                 %
0169                 disp([' - imhu: ',int2str(imhu),'/',int2str(nmhu),' - ipn: ',int2str(ipn),'/',int2str(npn),' - ipsi: ',int2str(ipsi),'/',int2str(npsi)])
0170             end
0171         end
0172     end
0173 end
0174 %
0175 end
0176 %
0177 toc
0178 %
0179 % keyboard
0180 % for ir=1:gridDKE.nr
0181 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXthetastar1(ir,:,:)),'p','mhu','theta^*_1',NaN,NaN,0,NaN,30); colorbar;
0182 % figure;graph2D_jd(gridDKE.pn,gridDKE.mhu,squeeze(XXthetastar2(ir,:,:)),'p','mhu','theta^*_2',NaN,NaN,0,NaN,30); colorbar;
0183 % end
0184 %
0185 end
0186 
0187 
0188 function [alpha,thetastar1,thetastar2] = loop_thetastar_dke_en(igrid,grid_index,xrho,equil_fit,xB0,Bmax,mhu,gamma,options_fzero)
0189     %
0190     if nargin < 9,options_fzero = [];end
0191     %
0192     ymin = fthetastar_dke_en(0,xrho(grid_index(igrid).ir),equil_fit,xB0(grid_index(igrid).ir),mhu(grid_index(igrid).imhu),gamma(grid_index(igrid).ipn));
0193     ymax = fthetastar_dke_en(pi,xrho(grid_index(igrid).ir),equil_fit,xB0(grid_index(igrid).ir),mhu(grid_index(igrid).imhu),gamma(grid_index(igrid).ipn));
0194     if ymin*ymax <= 0,
0195         %keyboard
0196         %theta_guess = 2*asin((2/(1-mhu(grid_index(igrid).imhu)^2)/(gamma(grid_index(igrid).ipn)+1)-1)./(Bmax(grid_index(igrid).ir)/xB0(grid_index(igrid).ir)-1));
0197         %del = (pi-theta_guess)/10;
0198         thetastar1 = fzero(@(theta) fthetastar_dke_en(theta,xrho(grid_index(igrid).ir),equil_fit,xB0(grid_index(igrid).ir),mhu(grid_index(igrid).imhu),gamma(grid_index(igrid).ipn)),[0 pi],options_fzero);
0199         thetastar2 = fzero(@(theta) fthetastar_dke_en(theta,xrho(grid_index(igrid).ir),equil_fit,xB0(grid_index(igrid).ir),mhu(grid_index(igrid).imhu),gamma(grid_index(igrid).ipn)),[pi 2*pi],options_fzero);       
0200         %
0201         if thetastar1 >= 0 && thetastar1 <= 2*pi,
0202             %
0203             equilval1 = equilval_yp(equil_fit,xrho(grid_index(igrid).ir),thetastar1);
0204             %
0205             x1 = equilval1.x;
0206             y1 = equilval1.y;
0207             Bx1 = equilval1.Bx;
0208             By1 = equilval1.By;
0209             dB1dtheta = equilval1.dBdtheta;
0210             %
0211         else
0212             thetastar1 = NaN;
0213         end
0214         %
0215         if thetastar2 >= 0 && thetastar2 <= 2*pi,
0216             %
0217             equilval2 = equilval_yp(equil_fit,xrho(grid_index(igrid).ir),thetastar2);
0218             %
0219             x2 = equilval2.x;
0220             y2 = equilval2.y;
0221             Bx2 = equilval2.Bx;
0222             By2 = equilval2.By;
0223             dB2dtheta = equilval2.dBdtheta;
0224             %
0225         else
0226             thetastar2 = NaN;
0227         end
0228         %
0229         if ~isnan(thetastar1) && ~isnan(thetastar2),
0230             alpha = xB0(grid_index(igrid).ir)*((x1*x1+y1*y1)/abs(-By1*x1+Bx1*y1)/abs(dB1dtheta) + (x2*x2+y2*y2)/abs(-By2*x2+Bx2*y2)/abs(dB2dtheta));
0231         else
0232             alpha = NaN;
0233         end
0234     else
0235         thetastar1 = NaN;
0236         thetastar2 = NaN;
0237         alpha = NaN;
0238     end
0239 end
0240 %
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.