matrix3Dgridbuilder_dke_yp

   Calculate matrix 3D grids for the 3D electron relativistic drift kinetic solver.

   by Y. Peysson (CEA-DRFC) (yves.peysson@cea.fr) and J. Decker (CEA-DRFC) (joan.decker@cea.fr)

0001 function [gridDKE] = matrix3Dgridbuilder_dke_yp(dkeparam,display,equilDKE,radialDKE,momentumDKE)
0002 %
0003 %   Calculate matrix 3D grids for the 3D electron relativistic drift kinetic solver.
0004 %
0005 %   by Y. Peysson (CEA-DRFC) (yves.peysson@cea.fr) and J. Decker (CEA-DRFC) (joan.decker@cea.fr)
0006 %
0007 xmhubounce2 = equilDKE.xmhu0T2;
0008 xrho = equilDKE.xrho;
0009 %
0010 mhu = momentumDKE.mhu;
0011 dmhu = momentumDKE.dmhu;
0012 mhum = momentumDKE.mhum;
0013 dmhum = momentumDKE.dmhum;
0014 mhup = momentumDKE.mhup;
0015 dmhup = momentumDKE.dmhup;
0016 mhu_S = momentumDKE.mhu_S;
0017 pn = momentumDKE.pn;
0018 dpn = momentumDKE.dpn;
0019 pnm = momentumDKE.pnm;
0020 dpnm = momentumDKE.dpnm;
0021 pnp = momentumDKE.pnp;
0022 dpnp = momentumDKE.dpnp;
0023 pn_S = momentumDKE.pn_S;
0024 spfac = dkeparam.spfac;
0025 %
0026 xpsin = radialDKE.xpsin_f;
0027 xpsin_dke = radialDKE.xpsin_f_dke;
0028 xdpsin_dke = radialDKE.xdpsin_f_dke;
0029 xpsinm_dke = radialDKE.xpsinm_S_dke;
0030 xdpsinm_dke = radialDKE.xdpsinm_S_dke;
0031 xpsinp_dke = radialDKE.xpsinp_S_dke;
0032 xdpsinp_dke = radialDKE.xdpsinp_S_dke;
0033 xpsin_S_dke = radialDKE.xpsin_S_dke;
0034 rdke = radialDKE.r_dke;
0035 ir_display_out = radialDKE.ir_display_out;
0036 %
0037 npn = length(pn);%Size of the momentum grid for the distribution function (half grid)
0038 nmhu = length(mhu);%Size of the pitch-angle grid for the distribution function (half grid)
0039 nr = length(xpsin);%Size of the radial grid for the distribution function (half grid)
0040 nr_dke = length(rdke);%Size of the effective radial grid for the drift kinetic equation (nr_dke = nr for FPE equation only)
0041 %
0042 masku = triu(ones(npn,npn)) - 0.5*diag(ones(1,npn));
0043 maskl = tril(ones(npn,npn)) - 0.5*diag(ones(1,npn));
0044 %
0045 %Angular grid
0046 %
0047 mhu2 = mhu.*mhu;
0048 Xmhu = ones(npn,1)*mhu;
0049 XXmhu = repmat(Xmhu,[1,1,nr]);
0050 Xmhu2 = ones(npn,1)*mhu2;
0051 X1mmhu2 = 1 - Xmhu2;
0052 XX1mmhu2 = repmat(X1mmhu2,[1,1,nr]);
0053 Xdmhu = ones(npn,1)*dmhu;
0054 XXdmhu = repmat(Xdmhu,[1,1,nr]);
0055 Xmhum = ones(npn,1)*mhum;
0056 Xmhum(find(Xmhum==0)) = eps;
0057 XXmhum = repmat(Xmhum,[1,1,nr]);
0058 Xmhumm = ones(npn,1)*[mhu(1),mhu(1:nmhu-1)];
0059 Xmhu2m = Xmhum.*Xmhum;
0060 X1mmhu2m = 1 - Xmhu2m;
0061 X1mmhu2m(find(X1mmhu2m==0)) = eps;
0062 XX1mmhu2m = repmat(X1mmhu2m,[1,1,nr]);
0063 Xdmhum = ones(npn,1)*dmhum;
0064 XXdmhum = repmat(Xdmhum,[1,1,nr]);
0065 Xdmhumm = ones(npn,1)*[dmhu(1),dmhu(1:nmhu-1)];
0066 Xmhup = ones(npn,1)*mhup;
0067 Xmhup(find(Xmhup==0)) = eps;
0068 XXmhup = repmat(Xmhup,[1,1,nr]);
0069 Xmhupp= ones(npn,1)*[mhu(2:nmhu),mhu(nmhu)];
0070 Xmhu2p = Xmhup.*Xmhup;
0071 X1mmhu2p = 1 - Xmhu2p;
0072 X1mmhu2p(find(X1mmhu2p==0)) = eps;
0073 XX1mmhu2p = repmat(X1mmhu2p,[1,1,nr]);
0074 Xdmhup = ones(npn,1)*dmhup;
0075 XXdmhup = repmat(Xdmhup,[1,1,nr]);
0076 Xdmhupp= ones(npn,1)*[dmhu(2:nmhu),dmhu(nmhu)];
0077 %
0078 %Momentum grid
0079 %
0080 ip_th_ref = min(find(pn >= 1));
0081 Xdpn = dpn'*ones(1,nmhu);
0082 pn2 = pn.*pn;
0083 pn3 = pn2.*pn;
0084 xpn = pn'*ones(1,npn);
0085 xpn2 = xpn.*xpn;
0086 xpn3 = xpn2.*xpn;
0087 Xpn = pn'*ones(1,nmhu);
0088 Xdpn = dpn'*ones(1,nmhu);
0089 Xpn2 = Xpn.*Xpn;
0090 Xpn3 = Xpn2.*Xpn;
0091 Xdpnm = dpnm'*ones(1,nmhu);
0092 XXdpnm = repmat(Xdpnm,[1,1,nr]);
0093 pn2m = pnm.*pnm;
0094 Xpnm = pnm'*ones(1,nmhu);
0095 XXpnm = repmat(Xpnm,[1,1,nr]);
0096 Xpn2m = Xpnm.*Xpnm;
0097 XXpn2m = XXpnm.*XXpnm;
0098 Xdpnp   = dpnp'*ones(1,nmhu);
0099 XXdpnp = repmat(Xdpnp,[1,1,nr]);
0100 pn2p = pnp.*pnp;  
0101 Xpnp = pnp'*ones(1,nmhu);
0102 XXpnp = repmat(Xpnp,[1,1,nr]);
0103 Xpn2p = Xpnp.*Xpnp;
0104 XXpn2p = XXpnp.*XXpnp;
0105 pnmm = [eps,pn(1:npn-1)];
0106 Xpnmm = pnmm'*ones(1,nmhu);
0107 XXpnmm = repmat(Xpnmm,[1,1,nr]);
0108 XXpn2mm = XXpnmm.*XXpnmm;
0109 Xdpnmm  = [dpn(1),dpn(1:npn-1)]'*ones(1,nmhu);
0110 pnpp = [pn(2:npn),pn(npn)];
0111 Xpnpp = pnpp'*ones(1,nmhu);
0112 XXpnpp = repmat(Xpnpp,[1,1,nr]);
0113 XXpn2pp = XXpnpp.*XXpnpp;
0114 Xdpnpp  = [dpn(2:npn),dpn(npn)]'*ones(1,nmhu);
0115 XXpn = repmat(Xpn,[1,1,nr]);
0116 XXpn2 = repmat(Xpn2,[1,1,nr]);
0117 %
0118 % Fill structure gridDKE for output
0119 %
0120 gridDKE.masku = masku;
0121 gridDKE.maskl = maskl;
0122 %
0123 gridDKE.mhu2 = mhu2;
0124 gridDKE.Xmhu = Xmhu;
0125 gridDKE.XXmhu = XXmhu;
0126 gridDKE.Xmhu2 = Xmhu2;
0127 gridDKE.X1mmhu2 = X1mmhu2;
0128 gridDKE.XX1mmhu2 = XX1mmhu2;
0129 gridDKE.Xdmhu = Xdmhu;
0130 gridDKE.XXdmhu = XXdmhu;
0131 gridDKE.Xmhum = Xmhum;
0132 gridDKE.XXmhum = XXmhum;
0133 gridDKE.Xmhumm = Xmhumm;
0134 gridDKE.Xmhu2m = Xmhu2m;
0135 gridDKE.X1mmhu2m = X1mmhu2m;
0136 gridDKE.XX1mmhu2m = XX1mmhu2m;
0137 gridDKE.Xdmhum = Xdmhum;
0138 gridDKE.XXdmhum = XXdmhum;
0139 gridDKE.Xdmhumm = Xdmhumm;
0140 gridDKE.Xmhup = Xmhup;
0141 gridDKE.XXmhup = XXmhup;
0142 gridDKE.Xmhupp = Xmhupp;
0143 gridDKE.Xmhu2p = Xmhu2p;
0144 gridDKE.X1mmhu2p = X1mmhu2p;
0145 gridDKE.XX1mmhu2p = XX1mmhu2p;
0146 gridDKE.Xdmhup = Xdmhup;
0147 gridDKE.XXdmhup = XXdmhup;
0148 gridDKE.Xdmhupp = Xdmhupp;
0149 gridDKE.mhu = mhu;
0150 gridDKE.dmhu = dmhu;
0151 gridDKE.mhum = mhum;
0152 gridDKE.dmhum = dmhum;
0153 gridDKE.mhup = mhup;
0154 gridDKE.dmhup = dmhup;
0155 gridDKE.mhu_S = mhu_S;
0156 %
0157 %Momentum grid
0158 %
0159 gridDKE.ip_th_ref = ip_th_ref;
0160 gridDKE.Xdpn = Xdpn;
0161 gridDKE.pn2 = pn2;
0162 gridDKE.pn3 = pn3;
0163 gridDKE.xpn = xpn;
0164 gridDKE.xpn2 = xpn2;
0165 gridDKE.xpn3 = xpn3;
0166 gridDKE.Xpn = Xpn;
0167 gridDKE.Xdpn = Xdpn;
0168 gridDKE.Xpn2 = Xpn2;
0169 gridDKE.Xpn3 = Xpn3;
0170 gridDKE.Xdpnm = Xdpnm;
0171 gridDKE.XXdpnm = XXdpnm;
0172 gridDKE.pn2m = pn2m;
0173 gridDKE.Xpnm = Xpnm;
0174 gridDKE.XXpnm = XXpnm;
0175 gridDKE.Xpn2m = Xpn2m;
0176 gridDKE.XXpn2m = XXpn2m;
0177 gridDKE.Xdpnp = Xdpnp;
0178 gridDKE.XXdpnp = XXdpnp;
0179 gridDKE.pn2p = pn2p;  
0180 gridDKE.Xpnp = Xpnp;
0181 gridDKE.XXpnp = XXpnp;
0182 gridDKE.Xpn2p = Xpn2p;
0183 gridDKE.XXpn2p = XXpn2p;
0184 gridDKE.pnmm = pnmm;
0185 gridDKE.Xpnmm = Xpnmm;
0186 gridDKE.XXpnmm = XXpnmm;
0187 gridDKE.XXpn2mm = XXpn2mm;
0188 gridDKE.Xdpnmm = Xdpnmm;
0189 gridDKE.pnpp = pnpp;
0190 gridDKE.Xpnpp = Xpnpp;
0191 gridDKE.XXpnpp = XXpnpp;
0192 gridDKE.XXpn2pp = XXpn2pp;
0193 gridDKE.Xdpnpp = Xdpnpp;
0194 gridDKE.XXpn = XXpn;
0195 gridDKE.XXpn2 = XXpn2;
0196 gridDKE.pn = pn;
0197 gridDKE.dpn = dpn;
0198 gridDKE.pnm = pnm;
0199 gridDKE.dpnm = dpnm;
0200 gridDKE.pnp = pnp;
0201 gridDKE.dpnp = dpnp;
0202 gridDKE.pn_S = pn_S;
0203 gridDKE.spfac = spfac;
0204 %
0205 gridDKE.xmhubounce2 = xmhubounce2;
0206 gridDKE.xrho = xrho;
0207 gridDKE.xpsin = xpsin;
0208 gridDKE.xpsin_dke = xpsin_dke;
0209 gridDKE.xdpsin_dke = xdpsin_dke;
0210 gridDKE.xpsinm_dke = xpsinm_dke;
0211 gridDKE.xdpsinm_dke = xdpsinm_dke;
0212 gridDKE.xpsinp_dke = xpsinp_dke;
0213 gridDKE.xdpsinp_dke = xdpsinp_dke;
0214 gridDKE.rdke = rdke;
0215 gridDKE.xpsin_S_dke = xpsin_S_dke;
0216 gridDKE.ir_display_out = ir_display_out;
0217 %
0218 gridDKE.npn = npn;
0219 gridDKE.nmhu = nmhu;
0220 gridDKE.nr = nr;
0221 gridDKE.nr_dke = nr_dke;
0222 %
0223 if display.display_mode == 2, 
0224     Fig_PitchGrid = figure('Name','Pitch-angle grid');
0225     figure(Fig_PitchGrid),plot(1:nmhu,mhu,'bo',1:nmhu,mhum,'g*',1:nmhu,mhup,'r+');grid on;zoom on;
0226     axis([1,nmhu+1,-1,1])
0227     title('Pitch-angle grid');
0228     xlabel('Index');
0229     ylabel('mhu');
0230     drawnow
0231 %
0232     Fig_MomGrid = figure('Name','Momentum grid');
0233     figure(Fig_MomGrid),plot(1:npn,pn,'bo',1:npn,pnm,'g*',1:npn,pnp,'r+');grid on;zoom on; 
0234     axis([1,npn+1,0,max(pnp)])
0235     title('Momentum grid');
0236     xlabel('Index');
0237     ylabel('pn');
0238     drawnow
0239 %
0240     Fig_MomStepGrid = figure('Name','Momentum step grid');
0241     figure(Fig_MomStepGrid),semilogx(pn,dpn,'bo',pnm,dpnm,'g*',pnp,dpnp,'r+');grid on;zoom on; 
0242     title('Momentum mesh (grid step)');
0243     xlabel('pn');
0244     ylabel('dpn');
0245     drawnow
0246 %
0247     Fig_RadGrid = figure('Name','Radial grid');
0248     figure(Fig_RadGrid),plot(1:nr,xpsin,'bo',rdke,xpsin(rdke),'r+',rdke,xrho(rdke),'bx');grid on;zoom on; 
0249     axis([0,nr+1,0,1])
0250     title('Radial mesh');
0251     xlabel('Index');
0252     ylabel('(psin: +,rho: x)');
0253     drawnow
0254 %
0255     Fig_MetRel = figure('Name','Metric relation');
0256     figure(Fig_MetRel),subplot(121),plot(xpsin,xrho,'r-.',xpsin,xrho,'b+');grid on;zoom on
0257     xlabel('psi/psia');ylabel('rho');title('Metric relation');
0258     axis('square');axis([0,1,0,1]);
0259     figure(Fig_MetRel),subplot(122),plot(xrho,xpsin,'r-.',xrho,xpsin,'b+');grid on;zoom on
0260     xlabel('rho');ylabel('psi/psia');title('Metric relation');
0261     axis('square');axis([0,1,0,1]);
0262     drawnow
0263 end    
0264 
0265 
0266 
0267 
0268 
0269 
0270 
0271 
0272 
0273