build_helmex_equil_yp

PURPOSE ^

SYNOPSIS ^

function [equil_magnetic,prho,pspsin,pspsinT,pp,li,jout,qpsi,ppsi] = build_helmex_equil_yp(nit,ppsi_in,ptot_in,pj_in,Rp_in,Zp_in,Bt_in,Ip_in,ap_in,R_in,Z_in,ntheta,display_opt)

DESCRIPTION ^

   Create a magnetic equilibrium structure using HELENA code (mex file version)

   INPUTS

       - psipin : flux poloidal (V*s, doit etre nul au centre) [Xpts,1]
       - ptot : ne Te + ni Ti (keV * m-3) [Xpts,1] 
       - jin : densite de courant (A/m2) [Xpts,1]
       - R0 : grand rayon (m) [1,1]
       - B0 : champ magnetique (T) [1,1]
       - Ip : courant plasma (A) [1,1], Ip < 0 -> utilisation de psip pour le calcul du courant
       - lsc : type d'input pour la derniere surface magnetique fermee, 0 = plasma sym,1 = plasma asymetrique, donnee de a, elong, triangh,triangb, 2 = plasma asymetrique, donnee de la LSC en (R,Z) [1,1]
       - a  : petit rayon (m) [1,1]
       - elong : ellipticite de la derniere surface magnetique fermee (1 pour TS) (-) [1,1]
       - triangh : triangularite haute de la derniere surface magnetique fermee (0 pour TS) (-) [1,1]
       - triangb : triangularite basse de la derniere surface magnetique fermee (0 pour TS) (-) [1,1]
       - init : si =1, alors on part de l'equilibre initial specifie par df2,b,kkbig [1,1] si =0, alors on part d'un equilibre initial quelconque
       - b : parametre b de l'equilibre initial (-) [1,1]
       - df2 : vecteur df2 de l'equilibre initial (-) [101,1]
       - kkbig : matrice kkbig deja calculee (-) [KKLDA,4*MAXNODE]
       - ecrit : matrice kkbig deja calculee (-) [KKLDA,4*MAXNODE]
       - Rext  : grand rayon de la derniere surface de flux [250,1]
       - Zext  : altitude de la derniere surface de flux [250,1]
       - xrsig = [XR1,SIG1] -> defini la position et la largeur de la gaussienne pour un maillage plus fin autour de XR1
       - dpr : vecteur dpr de l'equilibre initial (-) [101,1] (P')
       - FR     : Fourier coefficient of LCMS transform

   OUTPUTS

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [equil_magnetic,prho,pspsin,pspsinT,pp,li,jout,qpsi,ppsi] = build_helmex_equil_yp(nit,ppsi_in,ptot_in,pj_in,Rp_in,Zp_in,Bt_in,Ip_in,ap_in,R_in,Z_in,ntheta,display_opt)
0002 %
0003 %   Create a magnetic equilibrium structure using HELENA code (mex file version)
0004 %
0005 %   INPUTS
0006 %
0007 %       - psipin : flux poloidal (V*s, doit etre nul au centre) [Xpts,1]
0008 %       - ptot : ne Te + ni Ti (keV * m-3) [Xpts,1]
0009 %       - jin : densite de courant (A/m2) [Xpts,1]
0010 %       - R0 : grand rayon (m) [1,1]
0011 %       - B0 : champ magnetique (T) [1,1]
0012 %       - Ip : courant plasma (A) [1,1], Ip < 0 -> utilisation de psip pour le calcul du courant
0013 %       - lsc : type d'input pour la derniere surface magnetique fermee, 0 = plasma sym,1 = plasma asymetrique, donnee de a, elong, triangh,triangb, 2 = plasma asymetrique, donnee de la LSC en (R,Z) [1,1]
0014 %       - a  : petit rayon (m) [1,1]
0015 %       - elong : ellipticite de la derniere surface magnetique fermee (1 pour TS) (-) [1,1]
0016 %       - triangh : triangularite haute de la derniere surface magnetique fermee (0 pour TS) (-) [1,1]
0017 %       - triangb : triangularite basse de la derniere surface magnetique fermee (0 pour TS) (-) [1,1]
0018 %       - init : si =1, alors on part de l'equilibre initial specifie par df2,b,kkbig [1,1] si =0, alors on part d'un equilibre initial quelconque
0019 %       - b : parametre b de l'equilibre initial (-) [1,1]
0020 %       - df2 : vecteur df2 de l'equilibre initial (-) [101,1]
0021 %       - kkbig : matrice kkbig deja calculee (-) [KKLDA,4*MAXNODE]
0022 %       - ecrit : matrice kkbig deja calculee (-) [KKLDA,4*MAXNODE]
0023 %       - Rext  : grand rayon de la derniere surface de flux [250,1]
0024 %       - Zext  : altitude de la derniere surface de flux [250,1]
0025 %       - xrsig = [XR1,SIG1] -> defini la position et la largeur de la gaussienne pour un maillage plus fin autour de XR1
0026 %       - dpr : vecteur dpr de l'equilibre initial (-) [101,1] (P')
0027 %       - FR     : Fourier coefficient of LCMS transform
0028 %
0029 %   OUTPUTS
0030 %
0031 
0032 
0033 %   by Y. Peysson (CEA/DSM/IRFM,yves.peysson@cea.fr) and J. Decker (EPFL/CRPP,joan.decker@epfl.ch)
0034 %
0035 %Note: the toroidal magnetic field must be always positive for helmex77
0036 %and the plasma current must be positive for helmex77 so that the current
0037 %output is correctly normalized
0038 %
0039 init = [0,1e-3,1e-7];%Second argument is the accuracy on the current, the third on the global equilibrium
0040 %
0041 ppsi_init = abs(ppsi_in)/max(abs(ppsi_in));
0042 pj_init = pj_in;
0043 ptot_init = ptot_in;
0044 %
0045 if exist(['helmex77.',mexext]) == 3,
0046      [C2,C3,rho2m,Rout,Zout,rho,Vp,drhoor,qpsi,ppsi,ftra,Fdia,psiT,B2,invB2,Sp, ...
0047             rav,oor,drhoav,xshift,xrad,xell,xtriapos,xtrianeg,jout,ipout,pout, ...
0048             r2,oor2,r2tau2,drho2ob2,r3otau3,r3otau,BR,BZ,b,df2,dpr,kkbig,li, ...
0049             dimerc,drmerc,balcrit,frmode,ifail,nchi,cpsurf,radius,gem11,gem12,gem33,rspot,zspot,brspot,bzspot]= ...
0050             helmex77(ppsi_init,ptot_init,pj_init,Rp_in,abs(Bt_in),abs(Ip_in),2,ap_in,0,0,0,R_in,Z_in-Zp_in,init,0,0*ppsi_in,[],2);
0051     %
0052     if nit > 1,
0053         %
0054         disp('Start iteration for numerical consistency.')
0055         %
0056         pj_in = jout;
0057         %ppsi_in = ppsi/max(ppsi);
0058         %ptot_in = pout;
0059         %
0060         figure('name','HELENA convergence');
0061         %
0062         for it = 1:nit,
0063             %
0064             [C2,C3,rho2m,Rout,Zout,rho,Vp,drhoor,qpsi,ppsi,ftra,Fdia,psiT,B2,invB2,Sp, ...
0065                 rav,oor,drhoav,xshift,xrad,xell,xtriapos,xtrianeg,jout,ipout,pout, ...
0066                 r2,oor2,r2tau2,drho2ob2,r3otau3,r3otau,BR,BZ,b,df2,dpr,kkbig,li, ...
0067                 dimerc,drmerc,balcrit,frmode,ifail,nchi,cpsurf,radius,gem11,gem12,gem33,rspot,zspot,brspot,bzspot]= ...
0068                 helmex77(ppsi_init,ptot_init,pj_in,Rp_in,abs(Bt_in),abs(Ip_in),2,ap_in,0,0,0,R_in,Z_in-Zp_in,init,0,0*ppsi_in,[],2);
0069             %
0070             delta(it) = sum(abs(jout - pj_in))/max(abs(jout));%current convergence criterion
0071             %
0072             graph1D_jd([1:it],delta,0,1,'Iteration','Current convergence','HELENA',NaN,'','','-','none','b',2);hold on       
0073             %
0074             pj_in = 0.7*jout + 0.3*pj_in;%weighting for better stability of the convergence
0075             ppsi_in = ppsi;
0076         end
0077     else
0078         disp('WARNING: no numerical consistency performed.')
0079     end
0080     %
0081     [qe,me,mp,mn,e0,mu0] = pc_dke_yp;
0082     Ip_helmex = ipout.* max(xrad).*abs(Bt_in)/mu0*sign(Ip_in);%For test only (MA)
0083     %
0084     BPHI = (Fdia*ones(1,ntheta))./Rout*sign(Bt_in);%Recover also the correct sign
0085     BR(1,:) = 0;%avoid problems at the plasma center. BR must be zero exactly (sometimes spurious values)
0086     BZ(1,:) = 0;%avoid problems at the plasma center. BZ must be zero exactly (sometimes spurious values)
0087     BR = BR*sign(Ip_in);%Recover also the correct sign
0088     BZ = BZ*sign(Ip_in);%Recover also the correct sign
0089     ppsi = ppsi*sign(Ip_in);%Recover also the correct sign
0090     %
0091     BP = sqrt(BR.*BR+BZ.*BZ);%Poloidal magnetic field
0092     %
0093     if display_opt >= 1,
0094         disp('------------------------------------------------------');
0095         disp(['Magnetic field on axis guess (T) : ',num2str(Bt_in),', from HELMEX77 : ',num2str(BPHI(1,1))])
0096         disp(['Plasma current guess (MA) : ',num2str(Ip_in/1e6),', from HELMEX77 : ',num2str(Ip_helmex/1e6)])
0097         disp('------------------------------------------------------');
0098     end
0099     %
0100     ppsi = ppsi/2/pi;%In order to match the definition of ppsi in LUKE
0101     ppsin = ppsi'/ppsi(end);%Normalized flux grid
0102     theta = linspace(0,2*pi,ntheta);% uniform theta grid
0103     %
0104     %        Transfer to the equil_magnetic structure
0105     %
0106     if display_opt  >= 2,
0107         %
0108         figure('name','Separatrices'),
0109         [ax] = graph1D_jd(R_in - Rp_in,Z_in - Zp_in,0,0,'R - R_{p} (m)','Z - Z_{p} (m)','HELENA separatrix',NaN,'','','-','none','b',2);        
0110         [ax] = graph1D_jd(Rout(end,:) - Rp_in,Zout(end,:) - Zp_in,0,0,'R - R_{p} (m)','Z - Z_{p} (m)','HELENA separatrix',NaN,'','','--','none','r',2);
0111         legend('Guess','HELENA');
0112         [ax] = graph1D_jd(R_in(1,1) - Rp_in,Z_in(1,1) - Zp_in,0,0,'R - R_{p} (m)','Z - Z_{p} (m)','HELENA separatrix',NaN,'','','none','+','b',2);
0113         [ax] = graph1D_jd(Rout(end,1) - Rp_in,Zout(end,1) - Zp_in,0,0,'R - R_{p} (m)','Z - Z_{p} (m)','HELENA separatrix',NaN,'','','none','*','r',2);
0114         [ax] = graph1D_jd(Rout(1,1) - Rp_in,Zout(1,1) - Zp_in,0,0,'R - R_{p} (m)','Z - Z_{p} (m)','HELENA separatrix',NaN,'','','none','+','r',2);
0115         axis('equal')
0116         %
0117         figure('name','HELENA separatrix')
0118         [ax] = graph1D_jd(Rout(end,:),Zout(end,:) + Zp_in,0,0,'R (m)','Z (m)','HELENA separatrix',NaN,'','','--','none','r',2);
0119         [ax] = graph1D_jd(Rout(end,1),Zout(end,1) + Zp_in,0,0,'R (m)','Z (m)','HELENA separatrix',NaN,'','','none','*','r',2);
0120         [ax] = graph1D_jd(Rout(1,1),Zout(end,1) + Zp_in,0,0,'R (m)','Z (m)','HELENA separatrix',NaN,'','','none','+','r',2);
0121         axis('equal')
0122     end
0123     %
0124     equil_magnetic.psi_apRp = ppsi*(Rout(end,1) - Rp_in)/Rout(1,1);%psi_apRp = psi(poloidal usual)*ap/Rp for the cylindrical limit
0125     equil_magnetic.theta = theta;
0126     equil_magnetic.ptx = Rout - Rout(1,1);
0127     equil_magnetic.pty = Zout - Zp_in;
0128     equil_magnetic.ptBx = BR;
0129     equil_magnetic.ptBy = BZ;
0130     equil_magnetic.ptBPHI = BPHI;
0131     equil_magnetic.Rp = Rout(1,1);
0132     equil_magnetic.Zp = Zp_in;
0133     equil_magnetic.ap = Rout(end,1) - Rout(1,1);
0134     %
0135     equil_magnetic.helmex.qpsi = qpsi;
0136     equil_magnetic.helmex.ppsi = ppsi;
0137     equil_magnetic.helmex.jout = jout;
0138     equil_magnetic.helmex.li = li;
0139     %
0140     % Other output
0141     %
0142     prho = equil_magnetic.ptx(:,1).'/equil_magnetic.ptx(end,1);
0143     %
0144     pq_h = (qpsi(1:end-1) + qpsi(2:end))/2;
0145     pdpsin = diff(ppsin);
0146     ppsiT = [0,cumsum(pq_h'.*pdpsin)];
0147     pspsinT = sqrt(ppsiT/ppsiT(end));
0148     %
0149     pspsin = sqrt(ppsin);
0150     %
0151     pp = pout;%pressure given by HELENA
0152 else
0153     equil_magnetic = [];
0154     prho = [];
0155     pspsin = [];
0156     pspsinT = [];
0157     pp = [];
0158     li = [];
0159     disp([' ---> WARNING: The mexfile helmex77.',mexext,' does not exists ! Compile it for the computer you are using...']);
0160 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.