diff --git a/maxima/g0/gk-sheath/ms-gk-sheath.mac b/maxima/g0/gk-sheath/ms-gk-sheath.mac index 29ead044..fb3054db 100644 --- a/maxima/g0/gk-sheath/ms-gk-sheath.mac +++ b/maxima/g0/gk-sheath/ms-gk-sheath.mac @@ -11,6 +11,8 @@ load("gk-sheath/sheath-3x2v"); /* ...... USER INPUTS........ */ +outDir : "~/max-out/"$ /* Output directory */ + /* Serendipity basis. */ minPolyOrder_Ser : 1$ maxPolyOrder_Ser : 1$ @@ -43,7 +45,7 @@ for bInd : 1 thru length(bName) do ( minPolyOrderB : minPolyOrder[bInd], maxPolyOrderB : maxPolyOrder[bInd], for polyOrder : minPolyOrderB thru maxPolyOrderB do ( - fname : sconcat("~/max-out/bc_sheath_gyrokinetic_",bName[bInd],"_p",polyOrder,".c"), + fname : sconcat(outDir,"/bc_sheath_gyrokinetic_",bName[bInd],"_p",polyOrder,".c"), fh : openw(fname), printf(fh, "#include ~%"), @@ -69,7 +71,7 @@ for bInd : 1 thru length(bName) do ( /* Create a header file for sheath kernels. */ -fh : openw("~/max-out/gkyl_bc_sheath_gyrokinetic_kernels.h")$ +fh : openw(sconcat(outDir,"/gkyl_bc_sheath_gyrokinetic_kernels.h"))$ printf(fh, "#pragma once~%")$ printf(fh, "~%")$ printf(fh, "#include ~%")$ @@ -94,7 +96,7 @@ for bInd : 1 thru length(bName) do ( for vI : 1 thru length(vDims[cdim]) do ( for sI : 1 thru length(edge) do ( - printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl); ~%", edge[sI], cdim, vDims[cdim][vI], bName[bInd], polyOrder) + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl); ~%", edge[sI], cdim, vDims[cdim][vI], bName[bInd], polyOrder) ) ) diff --git a/maxima/g0/gk-sheath/sheath-1x1v.mac b/maxima/g0/gk-sheath/sheath-1x1v.mac index 92d54efd..a046fd57 100644 --- a/maxima/g0/gk-sheath/sheath-1x1v.mac +++ b/maxima/g0/gk-sheath/sheath-1x1v.mac @@ -54,7 +54,7 @@ genGkSheathKer1x1v(fh, cdim, vdim, basisFun, polyOrder) := block( for sI : 1 thru length(edge) do ( printf(fh, "~%"), - printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), /* Calculate expansion for deltaPhi = phiSheath - phiWall, evaluated at z=zVal, where zVal=-1 or 1 depending on if we're on the left or right edge of the global domain, respectively. */ diff --git a/maxima/g0/gk-sheath/sheath-1x2v.mac b/maxima/g0/gk-sheath/sheath-1x2v.mac index 54f41043..7fc512f5 100644 --- a/maxima/g0/gk-sheath/sheath-1x2v.mac +++ b/maxima/g0/gk-sheath/sheath-1x2v.mac @@ -13,7 +13,8 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( surfVarsC,bSurf,bV,bVpar,bMu,numBsurf,zvSurfVars,nodesMu,numNodesMu,bNMu,f_e,edge, surfPerpVal,sI,phiSheath_e,phiWall_e,deltaPhi_e,deltaPhi_q,vcutSq_q,fSurf_c,fSurf_e, fSurfMu_q,vparLo_e,vparUp_e,fReflSurfMu_e,j,nodeIdx,fSurfMu_c,fSurfMu_e,intLims,fReflSurfMu_c, - xBar_v,xSqBar_v,fReflSurf_c,fReflSurf_e,fRefl_c], + xBar_v,xSqBar_v,fReflSurf_c,fReflSurf_e,fRefl_c, + vcutFact_e,vcutFact_vpar0_e,vcutFact_q,vcutFactVpar_e,nodesVparMu,vcutFact_c,tempVars1,tempVars2], pdim : cdim+vdim, @@ -35,12 +36,12 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( bSurf : basisFromVars("gkhyb",surfVars,polyOrder), bV : basisFromVars("gkhyb",[vpar,mu],polyOrder), bVpar : basisFromVars("gkhyb",[vpar],polyOrder), - bMu : basisFromVars("gkhyb",[mu],polyOrder) + bMu : basisFromVars("ser",[mu],polyOrder) ) else ( bSurf : basisFromVars(basisFun,surfVars,polyOrder), bV : basisFromVars(basisFun,[vpar,mu],polyOrder), bVpar : basisFromVars(basisFun,[vpar],polyOrder), - bMu : basisFromVars(basisFun,[mu],polyOrder) + bMu : basisFromVars("ser",[mu],polyOrder) ), numBsurf : length(bSurf), @@ -62,6 +63,12 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( /* Expand input distribution function f on phase-space basis using input coeffs. */ f_e : doExpand1(f, bP), + + /* Expand vcut_fact on mu basis (bMu). */ + vcutFactVpar_e : doExpand1(vcut_fact, bMu), + + /* Evaluate vcut_fact at (mu_j) */ + vcutFact_q : evAtNodes(vcutFactVpar_e, nodesMu, zvSurfVars), /* Generate separate kernels for lower/upper boundaries. */ edge : ["lower", "upper"], @@ -70,7 +77,7 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( for sI : 1 thru length(edge) do ( printf(fh, "~%"), - printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), /* Calculate expansion for deltaPhi = phiSheath - phiWall, evaluated at z=zVal, where zVal=-1 or 1 depending on if we're on the left or right edge of the global domain, respectively. */ @@ -78,10 +85,9 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phiWall, bC)), deltaPhi_e : phiSheath_e - phiWall_e, - /* Evaluate vcutSqQ = vcut^2. q2Dm = q*2/m. */ + /* Evaluate base potential difference: deltaPhi = phiSheath - phiWall. */ deltaPhi_q : deltaPhi_e, - vcutSq_q : gcfac(float(fullratsimp(-q2Dm*deltaPhi_q))), - + /* Evaluate f at the z surface. */ fSurf_c : calcInnerProdList(surfVars, 1, bSurf, subst(surfPerpVar=surfPerpVal[sI], f_e)), fSurf_e : doExpand(fSurf_c, bSurf), @@ -94,7 +100,7 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " double vcutSq;~%"), printf(fh, " double fReflSurf[~a] = {0.}; ~%", length(bV)), printf(fh, "~%"), - + vparLo_e : subst(vpar=-1, vmap_e[1]), vparUp_e : subst(vpar=1, vmap_e[1]), printf(fh, " double vparLo = ~a;~%", float(expand(vparLo_e))), @@ -105,47 +111,20 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " double vparAbsSqUp = vmap[0]>0.? vparUp*vparUp : vparLo*vparLo;~%"), printf(fh, "~%"), - /* Write vcut^2 at boundary. */ - printf(fh, " vcutSq = ~a; ~%", vcutSq_q), - printf(fh, "~%"), - - /* If vcut^2 at this node is below all vpar^2 in this cell, BC at this node should be absorbing - so set coefficients of fRefl to 0 (no reflection from this node). */ - printf(fh, " if (vcutSq <= vparAbsSqLo) { // absorb (no reflection) ~%"), - printf(fh, "~%"), - - writeCExprsWithZeros1(fReflSurf, makelist(0.,j,1,length(bV))), - printf(fh, "~%"), - - /* If vcut^2 at this node is above all vpar^2 in this cell, BC at this node should be full reflection. */ - /* So set coefficients of fRefl to coefficients of f. */ - printf(fh, " } else if (vcutSq > vparAbsSqUp) { // full reflection ~%"), - printf(fh, "~%"), - - /* Project f onto vpar,mu basis (bV). */ - fSurf_c : gcfac(fullratsimp(calcInnerProdList(varsV, 1, bV, fSurf_e))), - /* Full reflection: set fRefl coefficients to f coefficients. */ - writeCExprsNoExpand1(fReflSurf, fSurf_c), - printf(fh, "~%"), - - /* If vcut^2 at this node is in this cell, BC at this node is partial reflection. */ - printf(fh, " } else { // partial reflection ~%"), - printf(fh, "~%"), - /* Reflected f at (mu)_j nodes. */ fReflSurfMu_e : makelist(0,j,1,numNodesMu), - printf(fh, " double wv = ~a;~%", float(fullratsimp((vparUp + vparLo)/2))), - printf(fh, " double dv = ~a;~%", float(fullratsimp(vparUp - vparLo))), + printf(fh, " double wv = ~a;~%", float(expand((vparUp + vparLo)/2))), + printf(fh, " double dv = ~a;~%", float(expand(vparUp - vparLo))), printf(fh, "~%"), - printf(fh, " double xBar, xSqBar;~%"), - printf(fh, " double fReflSurfMu[~a][~a] = {0.}; ~%", numNodesMu, length(bVpar)), + printf(fh, " double xBar, xSqBar;~%"), + printf(fh, " double fReflSurfMu[~a][~a] = {0.}; ~%", numNodesMu, length(bVpar)), printf(fh, "~%"), /* Loop over (mu)_j nodes. */ for j : 1 thru numNodesMu do ( - printf(fh, " // node (mu)_~a ~%", j-1), + printf(fh, " // node (mu)_~a ~%", j-1), nodeIdx : sublist_indices(nodesMu, lambda([x], x[1]=nodesMu[j][1]))[1], @@ -153,6 +132,25 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfMu_c : gcfac(fullratsimp(calcInnerProdList([vpar], 1, bVpar, fSurfMu_q[nodeIdx]))), fSurfMu_e : doExpand(fSurfMu_c, bVpar), + /* Compute vcut^2 at this mu node: vcutSq = alpha(mu) * q2Dm * deltaPhi. */ + printf(fh, " vcutSq = ~a; ~%", gcfac(float(expand(-vcutFact_q[nodeIdx]*q2Dm*deltaPhi_q)))), + printf(fh, "~%"), + + /* If vcut^2 at this node is below all vpar^2 in this cell, BC at this node should be absorbing + so set coefficients of fRefl at this (x,y)_i node to 0 (no reflection from this node). */ + printf(fh, " if (vcutSq <= vparAbsSqLo) { // absorb (no reflection) ~%"), + printf(fh, "~%"), + writeCExprsWithZeros1(fReflSurfMu[j-1], makelist(0,i,1,length(bVpar))), + printf(fh, "~%"), + + printf(fh, " } else if (vcutSq > vparAbsSqUp) { // full reflection ~%"), + printf(fh, "~%"), + writeCExprsNoExpand1(fReflSurfMu[j-1], fSurfMu_c), + printf(fh, "~%"), + + printf(fh, " } else { // partial reflection ~%"), + printf(fh, "~%"), + printf(fh, " if (wv > 0.) {~%"), printf(fh, " // vcut in logical space.~%"), printf(fh, " double vcut_l = 2.*(sqrt(vcutSq)-wv)/dv;~%"), @@ -193,6 +191,9 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " }~%"), printf(fh, "~%"), + printf(fh, " }~%"), /* End of partial reflection. */ + printf(fh, "~%"), + /* Expand vpar coefficients in vpar basis at this (mu)_j node. */ fReflSurfMu_e[j] : doExpand1(fReflSurfMu[j-1], bVpar) @@ -208,9 +209,6 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( writeCExprsNoExpand1(fReflSurf, fReflSurf_c), printf(fh, "~%"), - printf(fh, " }~%"), /* End of partial reflection else. */ - printf(fh, "~%"), - /* Expansion in vpar,mu of fReflSurf. */ fReflSurf_e : doExpand1(fReflSurf, bV), diff --git a/maxima/g0/gk-sheath/sheath-2x2v.mac b/maxima/g0/gk-sheath/sheath-2x2v.mac index 6f2f0597..1d5b87df 100644 --- a/maxima/g0/gk-sheath/sheath-2x2v.mac +++ b/maxima/g0/gk-sheath/sheath-2x2v.mac @@ -13,7 +13,8 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( bSurf,bV,bVpar,bMu,bX,numBsurf,zvSurfVars,nodesX,nodesXMu,nodesMu,numNodesX,numNodesXMu,numNodesMu, bNMu,bNX,f_e,edge,surfPerpVal,sI,phiSheath_e,phiWall_e,deltaPhi_e,deltaPhi_q,vcutSq_q,fSurf_c, fSurf_e,fSurfX_q,fSurfXMu_q,fReflSurfX_e,vparLo_e,vparUp_e,i,fSurfX_c,fReflSurfXMu_e,j,nodeIdx, - fSurfXMu_c,fSurfXMu_e,intLims,fReflSurfXMu_c,xBar_v,xSqBar_v,fReflSurfX_c,fReflSurf_c,fReflSurf_e,fRefl_c], + fSurfXMu_c,fSurfXMu_e,intLims,fReflSurfXMu_c,xBar_v,xSqBar_v,fReflSurfX_c,fReflSurf_c,fReflSurf_e,fRefl_c, + vcutFact_q,vcutFactVpar_e,nodesVparMu,vcutFact_c,tempVars1,tempVars2], pdim : cdim+vdim, @@ -36,12 +37,14 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( bV : basisFromVars("gkhyb",[vpar,mu],polyOrder), bVpar : basisFromVars("gkhyb",[vpar],polyOrder), bMu : basisFromVars("gkhyb",[mu],polyOrder), + bXMu : basisFromVars("ser",[x,mu],polyOrder), bX : basisFromVars("gkhyb",[x],polyOrder) ) else ( bSurf : basisFromVars(basisFun,surfVars,polyOrder), bV : basisFromVars(basisFun,[vpar,mu],polyOrder), bVpar : basisFromVars(basisFun,[vpar],polyOrder), bMu : basisFromVars(basisFun,[mu],polyOrder), + bXMu : basisFromVars("ser",[x,mu],polyOrder), bX : basisFromVars(basisFun,[x],polyOrder) ), numBsurf : length(bSurf), @@ -73,6 +76,12 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( /* Expand input distribution function f on phase-space basis using input coeffs. */ f_e : doExpand1(f, bP), + /* Expand vcut_fact on mu basis (bXMu). */ + vcutFactVpar_e : doExpand1(vcut_fact, bXMu), + + /* Evaluate vcut_fact at (x,mu)_ij nodes. */ + vcutFact_q : evAtNodes(vcutFactVpar_e, nodesXMu, zvSurfVars), + /* Generate separate kernels for lower/upper boundaries. */ edge : ["lower", "upper"], surfPerpVal : [-1, 1], @@ -80,7 +89,7 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( for sI : 1 thru length(edge) do ( printf(fh, "~%"), - printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), /* Calculate expansion for deltaPhi = phiSheath - phiWall, evaluated at z=zVal, where zVal=-1 or 1 depending on if we're on the left or right edge of the global domain, respectively. */ @@ -88,9 +97,8 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phiWall, bC)), deltaPhi_e : phiSheath_e - phiWall_e, - /* Evaluate vcutSq = vcut^2 at x nodes. */ + /* Evaluate deltaPhi at x nodes. */ deltaPhi_q : evAtNodes(deltaPhi_e, nodesX, surfVarsC), - vcutSq_q : gcfac(float(fullratsimp(-q2Dm*deltaPhi_q))), /* Evaluate f at the z surface. */ fSurf_c : calcInnerProdList(surfVars, 1, bSurf, subst(surfPerpVar=surfPerpVal[sI], f_e)), @@ -117,51 +125,28 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " double vparAbsSqUp = vmap[0]>0.? vparUp*vparUp : vparLo*vparLo;~%"), printf(fh, "~%"), - for i : 1 thru numNodesX do ( - printf(fh, " // node (x)_~a ~%", i-1), - - /* Write vcut^2 at this node. */ - printf(fh, " vcutSq = ~a;~%", vcutSq_q[i]), - printf(fh, "~%"), - - /* If vcut^2 at this node is below all vpar^2 in this cell, BC at this node should be absorbing - so set coefficients of fRefl at this (x) node to 0 (no reflection from this node). */ - printf(fh, " if (vcutSq <= vparAbsSqLo) { // absorb (no reflection)~%"), - printf(fh, "~%"), - - writeCExprsWithZeros1(fReflSurfX[i-1], makelist(0.,j,1,length(bV))), - printf(fh, "~%"), - - /* If vcut^2 at this node is above the max vpar^2 in this cell, BC at this node should be full reflection */ - /* so set coefficients of fRefl at this (x) node to coefficients of f. */ - printf(fh, " } else if (vcutSq > vparAbsSqUp) { // full reflection~%"), - printf(fh, "~%"), + printf(fh, " double wv;~%"), + printf(fh, " double dv;~%"), + printf(fh, "~%"), - /* Project f at this (x) node onto vpar,mu basis (bV). */ - fSurfX_c : gcfac(fullratsimp(calcInnerProdList(varsV, 1, bV, fSurfX_q[i]))), - /* Full reflection: set fRefl bZVpMu coefficients to f bZVpMu coefficients at this (x)_i node. */ - writeCExprsNoExpand1(fReflSurfX[i-1], fSurfX_c), - printf(fh, "~%"), + printf(fh, " double xBar, xSqBar;~%"), + printf(fh, " double fReflSurfXMu[~a][~a] = {0.}; ~%", numNodesMu, length(bVpar)), + printf(fh, "~%"), - /* If vcut^2 at this node is in this cell, BC at this node is partial reflection. */ - printf(fh, " } else { // partial reflection~%"), + for i : 1 thru numNodesX do ( + printf(fh, " // node (x)_~a ~%", i-1), printf(fh, "~%"), - /* Reflected f at (x)_i,(mu)_j nodes. Only need to store it for one - (x)_i node at a time. */ + /* Reflected f at (x)_i,(mu)_j nodes. */ fReflSurfXMu_e : makelist(0,j,1,numNodesMu), - printf(fh, " double wv = ~a;~%", float(fullratsimp((vparUp + vparLo)/2))), - printf(fh, " double dv = ~a;~%", float(fullratsimp(vparUp - vparLo))), - printf(fh, "~%"), - - printf(fh, " double xBar, xSqBar;~%"), - printf(fh, " double fReflSurfXMu[~a][~a] = {0.}; ~%", numNodesMu, length(bVpar)), + printf(fh, " wv = ~a;~%", float(expand((vparUp + vparLo)/2))), + printf(fh, " dv = ~a;~%", float(expand(vparUp - vparLo))), printf(fh, "~%"), /* Loop over (mu)_j nodes. */ for j : 1 thru numNodesMu do ( - printf(fh, " // node (mu)_~a ~%", j-1), + printf(fh, " // node (mu)_~a ~%", j-1), nodeIdx : sublist_indices(nodesXMu, lambda([x], x[1]=nodesX[i][1] and x[2]=nodesMu[j][1]))[1], @@ -169,6 +154,25 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfXMu_c : gcfac(fullratsimp(calcInnerProdList([vpar], 1, bVpar, fSurfXMu_q[nodeIdx]))), fSurfXMu_e : doExpand(fSurfXMu_c, bVpar), + /* Compute vcut^2 at this mu node: vcutSq = alpha(mu) * q2Dm * deltaPhi. */ + printf(fh, " vcutSq = ~a; ~%", gcfac(float(expand(-vcutFact_q[nodeIdx]*q2Dm*deltaPhi_q[i])))), + printf(fh, "~%"), + + /* If vcut^2 at this node is below all vpar^2 in this cell, BC at this node should be absorbing + so set coefficients of fRefl at this (x,y)_i node to 0 (no reflection from this node). */ + printf(fh, " if (vcutSq <= vparAbsSqLo) { // absorb (no reflection) ~%"), + printf(fh, "~%"), + writeCExprsWithZeros1(fReflSurfXMu[j-1], makelist(0,k,1,length(bVpar))), + printf(fh, "~%"), + + printf(fh, " } else if (vcutSq > vparAbsSqUp) { // full reflection ~%"), + printf(fh, "~%"), + writeCExprsNoExpand1(fReflSurfXMu[j-1], fSurfXMu_c), + printf(fh, "~%"), + + printf(fh, " } else { // partial reflection ~%"), + printf(fh, "~%"), + printf(fh, " if (wv > 0.) {~%"), printf(fh, " // vcut in logical space.~%"), printf(fh, " double vcut_l = 2.*(sqrt(vcutSq)-wv)/dv;~%"), @@ -179,7 +183,7 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( tempVars1 : [], tempVars1 : writeCExprs1noPowers(fReflSurfXMu[j-1], fReflSurfXMu_c, [vcut_l], tempVars1), printf(fh, "~%"), - + printf(fh, " } else {~%"), printf(fh, " // vcut in logical space.~%"), printf(fh, " double vcut_l = 2.*(-sqrt(vcutSq)-wv)/dv;~%"), @@ -198,15 +202,18 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( |bar{x}|<=1 and |bar{x^2}|<=1, where bar{g}=int g*f dx/int f dx. */ xBar_v : gcfac(float(fReflSurfXMu[j-1][1]/(sqrt(3)*fReflSurfXMu[j-1][0]))), xSqBar_v : gcfac(float((2*sqrt(5)*fReflSurfXMu[j-1][2]+5*fReflSurfXMu[j-1][0])/(15*fReflSurfXMu[j-1][0]))), - printf(fh, " // If the cut distribution f(vpar), where vpar \in [-1,vcut_l] or [vcut_l/1],~%"), - printf(fh, " // has a cell average < 0, set to 0. If it's >0 but not realizable, set to p=0.~%"), - printf(fh, " xBar = ~a;~%", xBar_v), - printf(fh, " xSqBar = ~a;~%", xSqBar_v), - printf(fh, " if (fReflSurfXMu[~a][0]<0.) {~%",j-1), - writeCExprsWithZeros1(fReflSurfXMu[j-1], makelist(0,i,1,length(bVpar))), - printf(fh, " } else if (fabs(xBar)>=1. || fabs(xSqBar)>=1.) {~%"), - writeCExprsWithZeros1(fReflSurfXMu[j-1], makelist(0,i,1,length(bVpar))), - printf(fh, " }~%"), + printf(fh, " // If the cut distribution f(vpar), where vpar \in [-1,vcut_l] or [vcut_l/1],~%"), + printf(fh, " // has a cell average < 0, set to 0. If it's >0 but not realizable, set to p=0.~%"), + printf(fh, " xBar = ~a;~%", xBar_v), + printf(fh, " xSqBar = ~a;~%", xSqBar_v), + printf(fh, " if (fReflSurfXMu[~a][0]<0.) {~%",j-1), + writeCExprsWithZeros1(fReflSurfXMu[j-1], makelist(0,k,1,length(bVpar))), + printf(fh, " } else if (fabs(xBar)>=1. || fabs(xSqBar)>=1.) {~%"), + writeCExprsWithZeros1(fReflSurfXMu[j-1], makelist(0,k,1,length(bVpar))), + printf(fh, " }~%"), + printf(fh, "~%"), + + printf(fh, " }~%"), /* End of partial reflection. */ printf(fh, "~%"), /* Expand vpar coefficients in vpar basis at this (x)_i (mu)_j node. */ @@ -224,9 +231,6 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( writeCExprsNoExpand1(fReflSurfX[i-1], fReflSurfX_c), printf(fh, "~%"), - printf(fh, " }~%"), /* End of partial reflection else. */ - printf(fh, "~%"), - /* Expansion in vpar,mu of fReflSurf at each (x)_i node. */ fReflSurfX_e[i] : doExpand1(fReflSurfX[i-1], bV) diff --git a/maxima/g0/gk-sheath/sheath-3x2v.mac b/maxima/g0/gk-sheath/sheath-3x2v.mac index 826a9a7a..26247580 100644 --- a/maxima/g0/gk-sheath/sheath-3x2v.mac +++ b/maxima/g0/gk-sheath/sheath-3x2v.mac @@ -14,7 +14,8 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( numNodesMu,bNMu,bNXY,f_e,edge,surfPerpVal,sI,phiSheath_e,phiWall_e,deltaPhi_e,deltaPhi_q, vcutSq_q,fSurf_c,fSurf_e,fSurfXY_q,fSurfXYMu_q,fReflSurfXY_e,vparLo_e,vparUp_e,i,fSurfXY_c, fReflSurfXYMu_e,j,nodeIdx,fSurfXYMu_c,fSurfXYMu_e,intLims,fReflSurfXYMu_c,xBar_v,xSqBar_v, - fReflSurfXY_c,fReflSurf_c,fReflSurf_e,fRefl_c], + fReflSurfXY_c,fReflSurf_c,fReflSurf_e,fRefl_c, + vcutFact_e,vcutFact_vpar0_e,vcutFact_q,vcutFactVpar_e,nodesVparMu,vcutFact_c,tempVars1,tempVars2], pdim : cdim+vdim, @@ -37,12 +38,14 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( bV : basisFromVars("gkhyb",[vpar,mu],polyOrder), bVpar : basisFromVars("gkhyb",[vpar],polyOrder), bMu : basisFromVars("gkhyb",[mu],polyOrder), + bXYMu : basisFromVars("ser",[x,y,mu],polyOrder), bXY : basisFromVars("gkhyb",[x,y],polyOrder) ) else ( bSurf : basisFromVars(basisFun,surfVars,polyOrder), bV : basisFromVars(basisFun,[vpar,mu],polyOrder), bVpar : basisFromVars(basisFun,[vpar],polyOrder), bMu : basisFromVars(basisFun,[mu],polyOrder), + bXYMu : basisFromVars("ser",[x,y,mu],polyOrder), bXY : basisFromVars(basisFun,[x,y],polyOrder) ), numBsurf : length(bSurf), @@ -74,6 +77,12 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( /* Expand input distribution function f on phase-space basis using input coeffs. */ f_e : doExpand1(f, bP), + /* Expand vcut_fact on vpar,mu basis (bV). */ + vcutFactVpar_e : doExpand1(vcut_fact, bXYMu), + + /* Evaluate vcut_fact at (x,y,mu)_ij nodes. */ + vcutFact_q : evAtNodes(vcutFactVpar_e, nodesXYMu, zvSurfVars), + /* Generate separate kernels for lower/upper boundaries. */ edge : ["lower", "upper"], surfPerpVal : [-1, 1], @@ -81,7 +90,7 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( for sI : 1 thru length(edge) do ( printf(fh, "~%"), - printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, double q2Dm, const double *phi, const double *phiWall, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, double q2Dm, const double *phi, const double *phiWall, const double *vcut_fact, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), /* Calculate expansion for deltaPhi = phiSheath - phiWall, evaluated at z=zVal, where zVal=-1 or 1 depending on if we're on the left or right edge of the global domain, respectively. */ @@ -89,9 +98,8 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phiWall, bC)), deltaPhi_e : phiSheath_e - phiWall_e, - /* Evaluate vcutSq = vcut^2 at x,y nodes. */ + /* Evaluate deltaPhi at x,y nodes. */ deltaPhi_q : evAtNodes(deltaPhi_e, nodesXY, surfVarsC), - vcutSq_q : gcfac(float(fullratsimp(-q2Dm*deltaPhi_q))), /* Evaluate f at the z surface. */ fSurf_c : calcInnerProdList(surfVars, 1, bSurf, subst(surfPerpVar=surfPerpVal[sI], f_e)), @@ -118,51 +126,26 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " double vparAbsSqUp = vmap[0]>0.? vparUp*vparUp : vparLo*vparLo;~%"), printf(fh, "~%"), + printf(fh, " double wv;~%"), + printf(fh, " double dv;~%"), + printf(fh, " double xBar, xSqBar;~%"), + printf(fh, " double fReflSurfXYMu[~a][~a] = {0.}; ~%", numNodesMu, length(bVpar)), + printf(fh, "~%"), + for i : 1 thru numNodesXY do ( printf(fh, " // node (x,y)_~a ~%", i-1), - - /* Write vcut^2 at this node. */ - printf(fh, " vcutSq = ~a;~%", vcutSq_q[i]), - printf(fh, "~%"), - - /* If vcut^2 at this node is below all vpar^2 in this cell, BC at this node should be absorbing - so set coefficients of fRefl at this (x,y) node to 0 (no reflection from this node). */ - printf(fh, " if (vcutSq <= vparAbsSqLo) { // absorb (no reflection)~%"), - printf(fh, "~%"), - - writeCExprsWithZeros1(fReflSurfXY[i-1], makelist(0.,j,1,length(bV))), printf(fh, "~%"), - /* If vcut^2 at this node is above the max vpar^2 in this cell, BC at this node should be full reflection */ - /* so set coefficients of fRefl at this (x,y) node to coefficients of f. */ - printf(fh, " } else if (vcutSq > vparAbsSqUp) { // full reflection~%"), - printf(fh, "~%"), - - /* Project f at this (x,y) node onto vpar,mu basis (bV). */ - fSurfXY_c : gcfac(fullratsimp(calcInnerProdList(varsV, 1, bV, fSurfXY_q[i]))), - /* Full reflection: set fRefl bZVpMu coefficients to f bZVpMu coefficients at this (x,y)_i node. */ - writeCExprsNoExpand1(fReflSurfXY[i-1], fSurfXY_c), - printf(fh, "~%"), - - /* If vcut^2 at this node is in this cell, BC at this node is partial reflection. */ - printf(fh, " } else { // partial reflection~%"), - printf(fh, "~%"), - - /* Reflected f at (x,y)_i,(mu)_j nodes. Only need to store it for one - (x,y)_i node at a time. */ + /* Reflected f at (x,y)_i,(mu)_j nodes. */ fReflSurfXYMu_e : makelist(0,j,1,numNodesMu), - printf(fh, " double wv = ~a;~%", float(fullratsimp((vparUp + vparLo)/2))), - printf(fh, " double dv = ~a;~%", float(fullratsimp(vparUp - vparLo))), - printf(fh, "~%"), - - printf(fh, " double xBar, xSqBar;~%"), - printf(fh, " double fReflSurfXYMu[~a][~a] = {0.}; ~%", numNodesMu, length(bVpar)), + printf(fh, " wv = ~a;~%", float(expand((vparUp + vparLo)/2))), + printf(fh, " dv = ~a;~%", float(expand(vparUp - vparLo))), printf(fh, "~%"), /* Loop over (mu)_j nodes. */ for j : 1 thru numNodesMu do ( - printf(fh, " // node (mu)_~a ~%", j-1), + printf(fh, " // node (mu)_~a ~%", j-1), nodeIdx : sublist_indices(nodesXYMu, lambda([x], x[1]=nodesXY[i][1] and x[2]=nodesXY[i][2] and x[3]=nodesMu[j][1]))[1], @@ -170,6 +153,25 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfXYMu_c : gcfac(fullratsimp(calcInnerProdList([vpar], 1, bVpar, fSurfXYMu_q[nodeIdx]))), fSurfXYMu_e : doExpand(fSurfXYMu_c, bVpar), + /* Compute vcut^2 at this mu node: vcutSq = alpha(mu) * q2Dm * deltaPhi. */ + printf(fh, " vcutSq = ~a; ~%", gcfac(float(expand(-vcutFact_q[nodeIdx]*q2Dm*deltaPhi_q[i])))), + printf(fh, "~%"), + + /* If vcut^2 at this node is below all vpar^2 in this cell, BC at this node should be absorbing + so set coefficients of fRefl at this (x,y)_i node to 0 (no reflection from this node). */ + printf(fh, " if (vcutSq <= vparAbsSqLo) { // absorb (no reflection) ~%"), + printf(fh, "~%"), + writeCExprsWithZeros1(fReflSurfXYMu[j-1], makelist(0,k,1,length(bVpar))), + printf(fh, "~%"), + + printf(fh, " } else if (vcutSq > vparAbsSqUp) { // full reflection ~%"), + printf(fh, "~%"), + writeCExprsNoExpand1(fReflSurfXYMu[j-1], fSurfXYMu_c), + printf(fh, "~%"), + + printf(fh, " } else { // partial reflection ~%"), + printf(fh, "~%"), + printf(fh, " if (wv > 0.) {~%"), printf(fh, " // vcut in logical space.~%"), printf(fh, " double vcut_l = 2.*(sqrt(vcutSq)-wv)/dv;~%"), @@ -180,7 +182,7 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( tempVars1 : [], tempVars1 : writeCExprs1noPowers(fReflSurfXYMu[j-1], fReflSurfXYMu_c, [vcut_l], tempVars1), printf(fh, "~%"), - + printf(fh, " } else {~%"), printf(fh, " // vcut in logical space.~%"), printf(fh, " double vcut_l = 2.*(-sqrt(vcutSq)-wv)/dv;~%"), @@ -194,20 +196,23 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " }~%"), printf(fh, "~%"), - + /* Evaluate realizability. In 1D a function is realization if |bar{x}|<=1 and |bar{x^2}|<=1, where bar{g}=int g*f dx/int f dx. */ xBar_v : gcfac(float(fReflSurfXYMu[j-1][1]/(sqrt(3)*fReflSurfXYMu[j-1][0]))), xSqBar_v : gcfac(float((2*sqrt(5)*fReflSurfXYMu[j-1][2]+5*fReflSurfXYMu[j-1][0])/(15*fReflSurfXYMu[j-1][0]))), - printf(fh, " // If the cut distribution f(vpar), where vpar \in [-1,vcut_l] or [vcut_l/1],~%"), - printf(fh, " // has a cell average < 0, set to 0. If it's >0 but not realizable, set to p=0.~%"), - printf(fh, " xBar = ~a;~%", xBar_v), - printf(fh, " xSqBar = ~a;~%", xSqBar_v), - printf(fh, " if (fReflSurfXYMu[~a][0]<0.) {~%",j-1), - writeCExprsWithZeros1(fReflSurfXYMu[j-1], makelist(0,i,1,length(bVpar))), - printf(fh, " } else if (fabs(xBar)>=1. || fabs(xSqBar)>=1.) {~%"), - writeCExprsWithZeros1(fReflSurfXYMu[j-1], makelist(0,i,1,length(bVpar))), - printf(fh, " }~%"), + printf(fh, " // If the cut distribution f(vpar), where vpar \in [-1,vcut_l] or [vcut_l/1],~%"), + printf(fh, " // has a cell average < 0, set to 0. If it's >0 but not realizable, set to p=0.~%"), + printf(fh, " xBar = ~a;~%", xBar_v), + printf(fh, " xSqBar = ~a;~%", xSqBar_v), + printf(fh, " if (fReflSurfXYMu[~a][0]<0.) {~%",j-1), + writeCExprsWithZeros1(fReflSurfXYMu[j-1], makelist(0,k,1,length(bVpar))), + printf(fh, " } else if (fabs(xBar)>=1. || fabs(xSqBar)>=1.) {~%"), + writeCExprsWithZeros1(fReflSurfXYMu[j-1], makelist(0,k,1,length(bVpar))), + printf(fh, " }~%"), + printf(fh, "~%"), + + printf(fh, " }~%"), /* End of partial reflection. */ printf(fh, "~%"), /* Expand vpar coefficients in vpar basis at this (x,y)_i (mu)_j node. */ @@ -225,9 +230,6 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( writeCExprsNoExpand1(fReflSurfXY[i-1], fReflSurfXY_c), printf(fh, "~%"), - printf(fh, " }~%"), /* End of partial reflection else. */ - printf(fh, "~%"), - /* Expansion in vpar,mu of fReflSurf at each (x,y)_i node. */ fReflSurfXY_e[i] : doExpand1(fReflSurfXY[i-1], bV)