efield_dke_jd

PURPOSE ^

SYNOPSIS ^

function [ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,xepsi,xEfield_validity] = efield_dke_jd(dkepath,dkeparam,display,equilDKE,mksa,gridDKE,Zmomcoef,Zbouncecoef,ohm,ZXXF_c)

DESCRIPTION ^

 Calculate matrix coefficient for the Ohmic electric field in LUKE

 by Y. PEYSSON, CEA-DRFC (yves.peysson@cea.fr) and J. DECKER (joan.decke@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ZXXD_e,ZXXF_e,ZXXD_e_tp,ZXXF_e_tp,xepsi,xEfield_validity] = efield_dke_jd(dkepath,dkeparam,display,equilDKE,mksa,gridDKE,Zmomcoef,Zbouncecoef,ohm,ZXXF_c)
0002 %
0003 % Calculate matrix coefficient for the Ohmic electric field in LUKE
0004 %
0005 % by Y. PEYSSON, CEA-DRFC (yves.peysson@cea.fr) and J. DECKER (joan.decke@cea.fr)
0006 %
0007 display_mode = display.display_mode;
0008 %
0009 dke_mode = dkeparam.dke_mode;
0010 bounce_mode = dkeparam.bounce_mode;
0011 dke_mode = dkeparam.dke_mode;
0012 temppath = dkepath.temppath;
0013 %
0014 etime_efield = 0;
0015 time0 = clock;
0016 %
0017 %load([fullfile(temppath,'temp_DKE.mat'],'xepsi');
0018 %
0019 if exist('xepsi') == 1, %The bounce averaging coefficients were calculated before.
0020     %
0021       %load([fullfile(temppath,'temp_DKE.mat'],'ZXXD_e','ZXXF_e','ZXXD_e_tp','ZXXF_e_tp','xepsi','xEfield_validity')
0022     %
0023     if display_mode >= 1,info_dke_yp(2,['Electric field flux coefficients retrieved.']);end    
0024     %
0025 else
0026 %
0027   xrho = equilDKE.xrho;
0028   xrhoT = equilDKE.xrhoT;
0029   xrhoP = equilDKE.xrhoP;
0030   xrhoV = equilDKE.xrhoV;
0031 %
0032   if dkeparam.nit_ohm > 1,% frozen current mode
0033       %
0034       sigma = - sign(equilDKE.XBphi(1,1));%epsi is defined versus B, electron charge included : epsi = -E.B
0035       %
0036       if ~isnan(dkeparam.Itot_ohm),% uniform electric field enforced to reach scalar current target
0037           ohm.id = 'Iconst';
0038           ohm.rho = equilDKE.xrho;
0039           if ~isfield(ohm,'epsi') || mean(ohm.epsi) == 0,
0040               ohm.epsi = sigma*dkeparam.efield_ohm*ones(size(ohm.rho));
0041           else
0042               ohm.epsi = mean(ohm.epsi)*ones(size(ohm.rho));
0043           end
0044       elseif isempty(ohm) || any(ohm.epsi == 0),% a non-zero initial field is needed to calculate consistent field
0045           ohm.id = 'jconst';
0046           ohm.rho = equilDKE.xrho;
0047           ohm.epsi = sigma*dkeparam.efield_ohm*ones(size(ohm.rho));
0048       end
0049   end    
0050 %
0051   xepsi_fsav = zeros(size(equilDKE.xrho));
0052   if ~isempty(ohm) && isstruct(ohm) && ~isempty(fieldnames(ohm)),
0053       if length(ohm.rho) == 1,%local calculation or uniform field
0054           xepsi_fsav = ohm.epsi*ones(size(xrho))/mksa.Edreicer_ref;%Normalized to the Dreicer value
0055       elseif ~isfield(ohm,'rhotype') || strcmp(ohm.rhotype,'g'),%geometric rho
0056           xepsi_fsav = interp1(ohm.rho,ohm.epsi,xrho,'spline')/mksa.Edreicer_ref;%Normalized to the Dreicer value
0057       elseif strcmp(ohm.rhotype,'p'),
0058           xepsi_fsav = interp1(ohm.rho,ohm.epsi,xrhoP,'spline')/mksa.Edreicer_ref;%Normalized to the Dreicer value
0059       elseif strcmp(ohm.rhotype,'t'),
0060           xepsi_fsav = interp1(ohm.rho,ohm.epsi,xrhoT,'spline')/mksa.Edreicer_ref;%Normalized to the Dreicer value
0061       elseif strcmp(ohm.rhotype,'v'),
0062           xepsi_fsav = interp1(ohm.rho,ohm.epsi,xrhoV,'spline')/mksa.Edreicer_ref;%Normalized to the Dreicer value
0063       else
0064           error('unknown rho type in ohm structure')
0065       end
0066   end
0067   %
0068   xB0 = equilDKE.xB0;
0069   xBT0 = equilDKE.xBT0;
0070   xx0 = equilDKE.xx0;
0071   Rp = equilDKE.Rp;
0072   Xmhu = gridDKE.Xmhu;
0073   X1mmhu2p = gridDKE.X1mmhu2p;
0074   X1mmhu2m = gridDKE.X1mmhu2m;
0075   ip_th_ref = gridDKE.ip_th_ref; 
0076   rdke =gridDKE.rdke;
0077   npn = gridDKE.npn;%Size of the momentum grids for the distribution function (half grid)
0078   nmhu = gridDKE.nmhu;
0079   nr = gridDKE.nr;%Size of the radial grid for the distribution function (half grid)
0080   nr_dke =gridDKE.nr_dke;%Size of the effective radial grid for the drift kinetic equation (nr_dke = nr for FPE equation only)
0081 %
0082   XXlambda_b_p1m3p4  = Zbouncecoef.XXlambda_b_p1m3p4;
0083   XXlambda_b_p1m2p2  = Zbouncecoef.XXlambda_b_p1m2p2;
0084   XXRlambda_b_p1m1p2 = Zbouncecoef.XXRlambda_b_p1m1p2;
0085   XXRlambda_b_p1m1p2_p = Zbouncecoef.XXRlambda_b_p1m1p2_p;
0086   XXRlambda_b_p1m1p2_m = Zbouncecoef.XXRlambda_b_p1m1p2_m;
0087   XXRlambda_b_p2m2p2  = Zbouncecoef.XXRlambda_b_p2m2p2;
0088   XXRlambda_b_p2m2p2_p = Zbouncecoef.XXRlambda_b_p2m2p2_p;
0089   XXRlambda_b_p2m2p2_m = Zbouncecoef.XXRlambda_b_p2m2p2_m;
0090   xqbar = Zbouncecoef.xqbar;
0091   xqtilde = Zbouncecoef.xqtilde;
0092   xqhat = Zbouncecoef.xqhat;
0093     %
0094     % Output files initialization
0095     %
0096     %
0097     % Terms from Sp at l+1/2 (Notations: _imj -> (i,j+1+2);_ipj -> (i+1,j+1/2)
0098     %
0099     ZXXD_e.pp_ipj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0100     ZXXD_e.pp_imj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0101     ZXXD_e.pm_ipj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0102     ZXXD_e.pm_imj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0103     ZXXD_e.pr_ipj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0104     ZXXD_e.pr_imj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0105     %
0106     ZXXF_e.p_ipj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0107     ZXXF_e.p_imj = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0108     %
0109     % Terms from Sksi at l+1/2 (Notations: _ijp -> (i+1/2,j+1);_ijm -> (i+1/2,j)
0110     %
0111     ZXXD_e.mm_ijp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0112     ZXXD_e.mm_ijm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0113     ZXXD_e.mp_ijp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0114     ZXXD_e.mp_ijm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0115     ZXXD_e.mr_ijp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0116     ZXXD_e.mr_ijm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0117     %
0118     ZXXF_e.m_ijp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0119     ZXXF_e.m_ijm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%momentum dynamics only
0120     %
0121     % Terms from Spsi at (i+1/2,j+1/2) (Notations: _lm -> (l);_lp -> (l+1);
0122     %
0123     ZXXD_e.rr_lm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%radial dynamics only
0124     ZXXD_e.rr_lp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%radial dynamics only
0125     ZXXD_e.rp_lm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0126     ZXXD_e.rp_lp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0127     ZXXD_e.rm_lm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0128     ZXXD_e.rm_lp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0129     %
0130     ZXXF_e.r_lm = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%radial dynamics only
0131     ZXXF_e.r_lp = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);%radial dynamics only
0132     %
0133     % Terms used by the old scheme only for momentum cross-derivatives (M.
0134     % Shoucri and I. Shkarofsky,Comp. Phys. Comm., 82 (1994) 287)
0135     %
0136     ZXXD_e.pp_ij = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0137     ZXXD_e.pm_ij = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0138     ZXXD_e.mp_ij = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0139     ZXXF_e.p_ij = zeros(gridDKE.npn,gridDKE.nmhu,gridDKE.nr);
0140 %
0141   if dke_mode == 1
0142 %
0143     ZXXD_e_tp.pp_ippj = zeros(npn,nmhu,nr_dke);
0144     ZXXD_e_tp.pp_ipj = zeros(npn,nmhu,nr_dke);
0145     ZXXD_e_tp.pp_ij = zeros(npn,nmhu,nr_dke);
0146     ZXXD_e_tp.pp_imj = zeros(npn,nmhu,nr_dke);
0147     ZXXD_e_tp.pp_immj = zeros(npn,nmhu,nr_dke);
0148     ZXXD_e_tp.pm_ipj = zeros(npn,nmhu,nr_dke);
0149     ZXXD_e_tp.pm_ij = zeros(npn,nmhu,nr_dke);
0150     ZXXD_e_tp.pm_imj = zeros(npn,nmhu,nr_dke);
0151     ZXXD_e_tp.mp_ijp = zeros(npn,nmhu,nr_dke);
0152     ZXXD_e_tp.mp_ij = zeros(npn,nmhu,nr_dke);
0153     ZXXD_e_tp.mp_ijm = zeros(npn,nmhu,nr_dke);
0154     ZXXD_e_tp.mm_ijp = zeros(npn,nmhu,nr_dke);
0155     ZXXD_e_tp.mm_ijm = zeros(npn,nmhu,nr_dke);
0156 %
0157     ZXXF_e_tp.p_ippj = zeros(npn,nmhu,nr_dke);
0158     ZXXF_e_tp.p_ipj = zeros(npn,nmhu,nr_dke);
0159     ZXXF_e_tp.p_ij = zeros(npn,nmhu,nr_dke);
0160     ZXXF_e_tp.p_imj = zeros(npn,nmhu,nr_dke);
0161     ZXXF_e_tp.p_immj = zeros(npn,nmhu,nr_dke);
0162     ZXXF_e_tp.m_ijp = zeros(npn,nmhu,nr_dke);
0163     ZXXF_e_tp.m_ijm = zeros(npn,nmhu,nr_dke);
0164 %
0165   else
0166 %
0167     ZXXD_e_tp = NaN;
0168     ZXXF_e_tp = NaN;
0169 %
0170   end
0171 
0172   for ir = 1:nr,
0173     if bounce_mode == 0 || xepsi_fsav(ir) == 0,
0174         xepsi(ir) = xepsi_fsav(ir);
0175     elseif bounce_mode == 1,
0176         if isreal(xepsi_fsav),%flux surface-averaging definition of LUKE using surface densities
0177             %xepsi(ir) = xepsi_fsav(ir)*(xqbar(ir)/xqtilde(ir))*(xB0(ir)/xBT0(ir))*(1+xx0(ir)/Rp)/XXlambda_b_p1m3p4(1,1,ir);%Value at Bmin for FP or DKE calculations
0178             xepsi(ir) = xepsi_fsav(ir)*(xqhat(ir)/xqtilde(ir))/XXlambda_b_p1m2p2(1,1,ir);%Value at Bmin for FP or DKE calculations
0179         else,%flux surface-averaging definition of CRONOS using volume densities
0180             xepsi(ir) = imag(xepsi_fsav(ir))*(xqhat(ir)/xqtilde(ir))/XXlambda_b_p1m2p2(1,1,ir);%Value at Bmin for FP or DKE calculations
0181         end
0182     end
0183     if xepsi(ir) > 0.05, 
0184 %        info_dke_yp(2,'WARNING: Ohmic electric field is likely too large. See the documentation');%To be consistent with the assumption of the linearized Fokker-Planck equation
0185     end 
0186 %
0187 % Electric field flux coefficients calculations for f0 and g
0188 %
0189     ZXXF_e.p_ipj(:,:,ir) = +Xmhu.*XXRlambda_b_p1m1p2(:,:,ir)*xepsi(ir);
0190     ZXXF_e.p_ij(:,:,ir) = +Xmhu.*XXRlambda_b_p1m1p2(:,:,ir)*xepsi(ir); 
0191     ZXXF_e.p_imj(:,:,ir) = +Xmhu.*XXRlambda_b_p1m1p2(:,:,ir)*xepsi(ir);
0192     ZXXF_e.m_ijp(:,:,ir) = -sqrt(X1mmhu2p).*XXRlambda_b_p1m1p2_p(:,:,ir)*xepsi(ir);
0193     ZXXF_e.m_ijm(:,:,ir) = -sqrt(X1mmhu2m).*XXRlambda_b_p1m1p2_m(:,:,ir)*xepsi(ir);    
0194 %
0195 % Electric field flux coefficients calculations for ftp
0196 %
0197     if dke_mode == 1
0198 %
0199         ir_dke = find(ir == rdke);
0200         if isempty(ir_dke)
0201             continue
0202         end
0203 %
0204         ZXXF_e_tp.p_ippj(:,:,ir_dke) = +Xmhu.*XXRlambda_b_p2m2p2(:,:,ir)*xepsi(ir);
0205         ZXXF_e_tp.p_ipj(:,:,ir_dke) = +Xmhu.*XXRlambda_b_p2m2p2(:,:,ir)*xepsi(ir);
0206         ZXXF_e_tp.p_ij(:,:,ir_dke) = +Xmhu.*XXRlambda_b_p2m2p2(:,:,ir)*xepsi(ir);    
0207         ZXXF_e_tp.p_imj(:,:,ir_dke) = +Xmhu.*XXRlambda_b_p2m2p2(:,:,ir)*xepsi(ir);
0208         ZXXF_e_tp.p_immj(:,:,ir_dke) = +Xmhu.*XXRlambda_b_p2m2p2(:,:,ir)*xepsi(ir);
0209 %
0210         ZXXF_e_tp.m_ijp(:,:,ir_dke) = -sqrt(X1mmhu2p).*XXRlambda_b_p2m2p2_p(:,:,ir)*xepsi(ir);
0211         ZXXF_e_tp.m_ijm(:,:,ir_dke) = -sqrt(X1mmhu2m).*XXRlambda_b_p2m2p2_m(:,:,ir)*xepsi(ir);
0212     end
0213 %
0214   end
0215 %
0216 % Test for E < ED
0217 %
0218 xEfield_validity = max(squeeze(ZXXF_e.p_ij(ip_th_ref,[1,nmhu],:)./ZXXF_c.p_ij(ip_th_ref,[1,nmhu],:)));%Comparison between electric field and collisions pinches at pth_ref (if above 1, the code validity breaks close to p = 0 leading to numerical particle leak)
0219 if display_mode >= 2,
0220     figure('Name','Code validity with Ohmic electric field'),semilogy(xrho,abs(xEfield_validity),'r+')
0221     xlabel('r/a')
0222     ylabel('Fp_e/Fp_coll at pth_ref')
0223     zoom on, grid on
0224     title('Code validity with Ohmic electric field')
0225 end
0226 %
0227 %  save([fullfile(temppath,'temp_DKE.mat'],'ZXXD_e','ZXXF_e','ZXXD_e_tp','ZXXF_e_tp','xepsi','xEfield_validity','-append');
0228   %
0229   etime_efield = etime_efield + etime(clock,time0);
0230   if display_mode >= 1,info_dke_yp(2,['Electric field flux coefficients calculations done in ',num2str(etime_efield),' (s)']);end    
0231   %
0232 %
0233 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.