fpengine_dke_yp

PURPOSE ^

SYNOPSIS ^

function [Xa,Xb,Xc,Xd,Xe,Xf,Xg,Xh,Xi] =fpengine_dke_yp(dkeparam,ftp_mode,Xdeltap_imj,Xdeltap_ipj,Xdeltam_ijm,Xdeltam_ijp,Xdeltap_imjmm,Xdeltap_imjpp,Xdeltap_ipjmm,Xdeltap_ipjpp,Xdeltam_immjm,Xdeltam_ippjm,Xdeltam_immjp,Xdeltam_ippjp,XRlambdam,XRlambdap,Xpnmm,Xpnm,Xpn,Xpnp,Xpnpp,Xpn2m,Xpn2p,Xdpnmm,Xdpnm,Xdpn,Xdpnp,Xdpnpp,Xmhumm,Xmhum,Xmhu,Xmhup,Xmhupp,X1mmhu2,X1mmhu2m,X1mmhu2p,Xdmhumm,Xdmhum,Xdmhu,Xdmhup,Xdmhupp,XDpp_ipj,XDpp_imj,XDpm_ipj,XDpm_imj,XFp_ipj,XFp_imj,XDmm_ijp,XDmm_ijm,XDmp_ijp,XDmp_ijm,XFm_ijp,XFm_ijm,XDpp_ij,XDpm_ij,XDmp_ij,XFp_ij);

DESCRIPTION ^

Fokker-Planck engine for zero (usual FP code) or first order (neoclassical corrections for bootstrap current).
It calculates the matrix coefficients starting from diffusion and convection in momentum space. Full conservative
form. All processes are therefore considered in the same way

    Input:

    Output:

    - Xa: Matrix coefficient for f0(i+1/2,j-1/2)
    - Xb: Matrix coefficient for f0(i-1/2,j+1/2)
    - Xc: Matrix coefficient for f0(i+1/2,j+1/2)
    - Xd: Matrix coefficient for f0(i+3/2,j+1/2)
    - Xe: Matrix coefficient for f0(i+1/2,j+3/2)
    - Xf: Matrix coefficient for f0(i-1/2,j-1/2)
    - Xg: Matrix coefficient for f0(i+3/2,j+3/2)
    - Xh: Matrix coefficient for f0(i-1/2,j+3/2)
    - Xi: Matrix coefficient for f0(i+3/2,j-1/2)
      (Notations: _imjp -> (i,j+1);_imjm -> (i,j);_ijp -> (i+1/2,j+1);_ijm -> (i+1/2,j);_immjp -> (i-1/2,j+1);_immjm -> (i-1/2,j)

By Joan Decker (MIT-RLE, jodecker@mit.edu) and Yves Peysson (CEA-DRFC, yves.peysson@cea.fr)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function      [Xa,Xb,Xc,Xd,Xe,Xf,Xg,Xh,Xi] = ...
0002                 fpengine_dke_yp(dkeparam,ftp_mode,...
0003                 Xdeltap_imj,Xdeltap_ipj,Xdeltam_ijm,Xdeltam_ijp,...
0004                 Xdeltap_imjmm,Xdeltap_imjpp,Xdeltap_ipjmm,Xdeltap_ipjpp,...
0005                 Xdeltam_immjm,Xdeltam_ippjm,Xdeltam_immjp,Xdeltam_ippjp,...
0006                 XRlambdam,XRlambdap,...
0007                 Xpnmm,Xpnm,Xpn,Xpnp,Xpnpp,Xpn2m,Xpn2p,Xdpnmm,Xdpnm,Xdpn,Xdpnp,Xdpnpp,...
0008                 Xmhumm,Xmhum,Xmhu,Xmhup,Xmhupp,X1mmhu2,X1mmhu2m,X1mmhu2p,Xdmhumm,Xdmhum,Xdmhu,Xdmhup,Xdmhupp,...
0009                 XDpp_ipj,XDpp_imj,XDpm_ipj,XDpm_imj,XFp_ipj,XFp_imj,...
0010                 XDmm_ijp,XDmm_ijm,XDmp_ijp,XDmp_ijm,XFm_ijp,XFm_ijm,...
0011                 XDpp_ij,XDpm_ij,XDmp_ij,XFp_ij);
0012 %
0013 %Fokker-Planck engine for zero (usual FP code) or first order (neoclassical corrections for bootstrap current).
0014 %It calculates the matrix coefficients starting from diffusion and convection in momentum space. Full conservative
0015 %form. All processes are therefore considered in the same way
0016 %
0017 %    Input:
0018 %
0019 %    Output:
0020 %
0021 %    - Xa: Matrix coefficient for f0(i+1/2,j-1/2)
0022 %    - Xb: Matrix coefficient for f0(i-1/2,j+1/2)
0023 %    - Xc: Matrix coefficient for f0(i+1/2,j+1/2)
0024 %    - Xd: Matrix coefficient for f0(i+3/2,j+1/2)
0025 %    - Xe: Matrix coefficient for f0(i+1/2,j+3/2)
0026 %    - Xf: Matrix coefficient for f0(i-1/2,j-1/2)
0027 %    - Xg: Matrix coefficient for f0(i+3/2,j+3/2)
0028 %    - Xh: Matrix coefficient for f0(i-1/2,j+3/2)
0029 %    - Xi: Matrix coefficient for f0(i+3/2,j-1/2)
0030 %      (Notations: _imjp -> (i,j+1);_imjm -> (i,j);_ijp -> (i+1/2,j+1);_ijm -> (i+1/2,j);_immjp -> (i-1/2,j+1);_immjm -> (i-1/2,j)
0031 %
0032 %By Joan Decker (MIT-RLE, jodecker@mit.edu) and Yves Peysson (CEA-DRFC, yves.peysson@cea.fr)
0033 %
0034 if nargin < 47,
0035     error('Wrong number of input arguments in fpengine_dke_yp')
0036 end
0037 %
0038 if ftp_mode >= 1, ftp_mode = 1;end%Force ftp_mode to be 0 (FP) or 1 (DKE)
0039 %
0040 if ~isfield(dkeparam,'oldcross_mode'),dkeparam.oldcross_mode = 0;end
0041 %
0042 %Center cross-derivatives
0043 %
0044 %(as done by M. Shoucri and I. Shkarofsky,Comp. Phys. Comm., 82 (1994) 287
0045 %Same results for LHCD but not very accurate for particle conservation, and
0046 %unstable solution when D// > 1). It is therefore only given for a link
0047 %with the past ! WARNING: IT MUST BE NEVER USED EXCEPT TESTS AND THE
0048 %REFERENCE PUBLICATION OF LUKE
0049 %
0050 if dkeparam.oldcross_mode == 1,%************ Must never be used **************
0051    Xa = - Xpnp.*XDpm_ipj.*sqrt(X1mmhu2)./(Xdmhum + Xdmhup)./Xdpn;%grid point (i+1/2,j-1/2)
0052    Xa = Xa + Xpnm.*XDpm_imj.*sqrt(X1mmhu2)./(Xdmhum + Xdmhup)./Xdpn;
0053    Xa = Xa - XDmm_ijm.*X1mmhu2m.*XRlambdam./Xdmhu./Xdmhum;
0054    Xa = Xa + Xpn.*sqrt(X1mmhu2m).*XFm_ijm.*Xdeltam_ijm.*(Xmhum./Xmhumm).^ftp_mode.*XRlambdam./Xdmhu;
0055 %
0056    Xb = -Xpn2m.*XDpp_imj./Xdpn./Xdpnm;%grid point (i-1/2,j+1/2)
0057    Xb = Xb - Xpn.*sqrt(X1mmhu2p).*XDmp_ijp.*XRlambdap./Xdmhu./(Xdpnm + Xdpnp);
0058    Xb = Xb + Xpn.*sqrt(X1mmhu2m).*XDmp_ijm.*XRlambdam./Xdmhu./(Xdpnm + Xdpnp);
0059    Xb = Xb - Xpn2m.*XFp_imj.*Xdeltap_imj.*(Xpnm./Xpnmm).^ftp_mode./Xdpn;    
0060 %
0061    Xc = Xpn2p.*XDpp_ipj./Xdpn./Xdpnp;%grid point (i+1/2,j+1/2)
0062    Xc = Xc + Xpn2p.*XFp_ipj.*Xdeltap_ipj.*(Xpnp./Xpn).^ftp_mode./Xdpn;
0063    Xc = Xc + Xpn2m.*XDpp_imj./Xdpn./Xdpnm;
0064    Xc = Xc - Xpn2m.*XFp_imj.*(1 - Xdeltap_imj).*(Xpnm./Xpn).^ftp_mode./Xdpn;
0065    Xc = Xc + XDmm_ijp.*X1mmhu2p.*XRlambdap./Xdmhu./Xdmhup;    
0066    Xc = Xc - Xpn.*sqrt(X1mmhu2p).*XFm_ijp.*XRlambdap.*Xdeltam_ijp.*(Xmhup./Xmhu).^ftp_mode./Xdmhu;
0067    Xc = Xc + XDmm_ijm.*X1mmhu2m.*XRlambdam./Xdmhu./Xdmhum;
0068    Xc = Xc + Xpn.*sqrt(X1mmhu2m).*XFm_ijm.*(1 - Xdeltam_ijm).*(Xmhum./Xmhu).^ftp_mode.*XRlambdam./Xdmhu;
0069 %
0070    Xd = - Xpn2p.*XDpp_ipj./Xdpn./Xdpnp;%grid point (i+3/2,j+1/2)
0071    Xd = Xd + Xpn.*sqrt(X1mmhu2p).*XDmp_ijp.*XRlambdap./Xdmhu./(Xdpnm + Xdpnp);
0072    Xd = Xd - Xpn.*sqrt(X1mmhu2m).*XDmp_ijm.*XRlambdam./Xdmhu./(Xdpnm + Xdpnp);
0073    Xd = Xd + Xpn2p.*XFp_ipj.*(1 - Xdeltap_ipj).*(Xpnp./Xpnpp).^ftp_mode./Xdpn;    
0074 %
0075    Xe = Xpnp.*XDpm_ipj.*sqrt(X1mmhu2)./(Xdmhum + Xdmhup)./Xdpn;%grid point (i+1/2,j+3/2)
0076    Xe = Xe - Xpnm.*XDpm_imj.*sqrt(X1mmhu2)./(Xdmhum + Xdmhup)./Xdpn;
0077    Xe = Xe - XDmm_ijp.*X1mmhu2p.*XRlambdap./Xdmhu./Xdmhup;
0078    Xe = Xe - Xpn.*sqrt(X1mmhu2p).*XFm_ijp.*(1 - Xdeltam_ijp).*(Xmhup./Xmhupp).^ftp_mode.*XRlambdap./Xdmhu;    
0079 %
0080    Xf = XDpm_ij.*Xpn.*sqrt(X1mmhu2)./(Xdmhum + Xdmhup)./(Xdpnm + Xdpnp);%grid point (i-1/2,j-1/2)
0081    Xf = Xf + XDmp_ij.*Xpn.*sqrt(X1mmhu2)./(Xdmhum + Xdmhup)./(Xdpnm + Xdpnp);
0082 %
0083    Xg = Xf;%grid point (i+3/2,j+3/2)
0084 %
0085    Xh = -XDpm_ij.*Xpn.*sqrt(X1mmhu2)./(Xdmhum + Xdmhup)./(Xdpnm + Xdpnp);%grid point (i-1/2,j+3/2)
0086    Xh = Xh - XDmp_ij.*Xpn.*sqrt(X1mmhu2)./(Xdmhum + Xdmhup)./(Xdpnm + Xdpnp);    
0087 %
0088    Xi = Xh;%grid point (i+3/2,j-1/2)
0089 elseif dkeparam.oldcross_mode == 0,%************ Must be always used **************
0090 %
0091 %Cross derivatives with exact boundary conditions (WARNING: MUST BE ALWAYS USED)
0092 %
0093    Xa = - Xpnp.*XDpm_ipj.*sqrt(X1mmhu2).*Xdeltap_ipjmm./(Xdmhum + Xdmhup)./Xdpn;%grid point (i+1/2,j-1/2)
0094    Xa = Xa + Xpnm.*XDpm_imj.*sqrt(X1mmhu2).*(1 - Xdeltap_imjmm)./(Xdmhum + Xdmhup)./Xdpn;
0095    Xa = Xa - XDmm_ijm.*X1mmhu2m.*XRlambdam./Xdmhu./Xdmhum;
0096    Xa = Xa + Xpn.*sqrt(X1mmhu2m).*XFm_ijm.*Xdeltam_ijm.*(Xmhum./Xmhumm).^ftp_mode.*XRlambdam./Xdmhu;
0097 %
0098    Xb = -Xpn2m.*XDpp_imj./Xdpn./Xdpnm;%grid point (i-1/2,j+1/2) OK
0099    Xb = Xb - Xpn.*sqrt(X1mmhu2p).*XDmp_ijp.*XRlambdap.*Xdeltam_immjp./Xdmhu./(Xdpnm + Xdpnp);
0100    Xb = Xb + Xpn.*sqrt(X1mmhu2m).*XDmp_ijm.*XRlambdam.*(1 - Xdeltam_immjm)./Xdmhu./(Xdpnm + Xdpnp);
0101    Xb = Xb - Xpn2m.*XFp_imj.*Xdeltap_imj.*(Xpnm./Xpnmm).^ftp_mode./Xdpn;    
0102 %
0103    Xc = Xpn2p.*XDpp_ipj./Xdpn./Xdpnp;%grid point (i+1/2,j+1/2)
0104    Xc = Xc + Xpn2p.*XFp_ipj.*Xdeltap_ipj.*(Xpnp./Xpn).^ftp_mode./Xdpn;
0105    Xc = Xc + Xpn2m.*XDpp_imj./Xdpn./Xdpnm;
0106    Xc = Xc - Xpn2m.*XFp_imj.*(1 - Xdeltap_imj).*(Xpnm./Xpn).^ftp_mode./Xdpn;
0107    Xc = Xc + XDmm_ijp.*X1mmhu2p.*XRlambdap./Xdmhu./Xdmhup;    
0108    Xc = Xc - Xpn.*sqrt(X1mmhu2p).*XFm_ijp.*XRlambdap.*Xdeltam_ijp.*(Xmhup./Xmhu).^ftp_mode./Xdmhu;
0109    Xc = Xc + XDmm_ijm.*X1mmhu2m.*XRlambdam./Xdmhu./Xdmhum;
0110    Xc = Xc + Xpn.*sqrt(X1mmhu2m).*XFm_ijm.*(1 - Xdeltam_ijm).*(Xmhum./Xmhu).^ftp_mode.*XRlambdam./Xdmhu;
0111 %
0112    Xd = - Xpn2p.*XDpp_ipj./Xdpn./Xdpnp;%grid point (i+3/2,j+1/2) OK
0113    Xd = Xd + Xpn.*sqrt(X1mmhu2p).*XDmp_ijp.*XRlambdap.*Xdeltam_ippjp./Xdmhu./(Xdpnm + Xdpnp);
0114    Xd = Xd - Xpn.*sqrt(X1mmhu2m).*XDmp_ijm.*XRlambdam.*(1 - Xdeltam_ippjm)./Xdmhu./(Xdpnm + Xdpnp);
0115    Xd = Xd + Xpn2p.*XFp_ipj.*(1 - Xdeltap_ipj).*(Xpnp./Xpnpp).^ftp_mode./Xdpn;    
0116 %
0117    Xe = Xpnp.*XDpm_ipj.*sqrt(X1mmhu2).*Xdeltap_ipjpp./(Xdmhum + Xdmhup)./Xdpn;%grid point (i+1/2,j+3/2) OK
0118    Xe = Xe - Xpnm.*XDpm_imj.*sqrt(X1mmhu2).*(1 - Xdeltap_imjpp)./(Xdmhum + Xdmhup)./Xdpn;
0119    Xe = Xe - XDmm_ijp.*X1mmhu2p.*XRlambdap./Xdmhu./Xdmhup;
0120    Xe = Xe - Xpn.*sqrt(X1mmhu2p).*XFm_ijp.*(1 - Xdeltam_ijp).*(Xmhup./Xmhupp).^ftp_mode.*XRlambdap./Xdmhu;    
0121 %
0122    Xf = XDpm_imj.*Xpnm.*sqrt(X1mmhu2).*Xdeltap_imjmm./(Xdmhum + Xdmhup)./Xdpn;%grid point (i-1/2,j-1/2) OK
0123    Xf = Xf + XDmp_ijm.*Xpn.*sqrt(X1mmhu2m).*XRlambdam.*Xdeltam_immjm./Xdmhu./(Xdpnm + Xdpnp);
0124 %
0125    Xg = XDpm_ipj.*Xpnp.*sqrt(X1mmhu2).*(1 - Xdeltap_ipjpp)./(Xdmhum + Xdmhup)./Xdpn;%grid point (i+3/2,j+3/2) OK
0126    Xg = Xg + XDmp_ijp.*Xpn.*sqrt(X1mmhu2p).*XRlambdap.*(1 - Xdeltam_ippjp)./Xdmhu./(Xdpnm + Xdpnp);
0127 %
0128    Xh = -XDpm_imj.*Xpnm.*sqrt(X1mmhu2).*Xdeltap_imjpp./(Xdmhum + Xdmhup)./Xdpn;%grid point (i-1/2,j+3/2) OK
0129    Xh = Xh - XDmp_ijp.*Xpn.*sqrt(X1mmhu2p).*XRlambdap.*(1 - Xdeltam_immjp)./Xdmhu./(Xdpnm + Xdpnp);    
0130 %
0131    Xi = -XDpm_ipj.*Xpnp.*sqrt(X1mmhu2).*(1 - Xdeltap_ipjmm)./(Xdmhum + Xdmhup)./Xdpn;%grid point (i+3/2,j-1/2) OK
0132    Xi = Xi - XDmp_ijm.*Xpn.*sqrt(X1mmhu2m).*XRlambdam.*Xdeltam_ippjm./Xdmhu./(Xdpnm + Xdpnp);
0133 end
Community support and wiki are available on Redmine. Last update: 18-Apr-2019.