diff --git a/maxima/g0/gk-sheath/ms-gk-sheath.mac b/maxima/g0/gk-sheath/ms-gk-sheath.mac index 29ead044..7b5db1f7 100644 --- a/maxima/g0/gk-sheath/ms-gk-sheath.mac +++ b/maxima/g0/gk-sheath/ms-gk-sheath.mac @@ -8,9 +8,14 @@ load("gk-sheath/sheath-1x1v"); load("gk-sheath/sheath-1x2v"); load("gk-sheath/sheath-2x2v"); load("gk-sheath/sheath-3x2v"); +load("gk-sheath/sheath_surrogate-1x2v"); +load("gk-sheath/sheath_surrogate-2x2v"); +load("gk-sheath/sheath_surrogate-3x2v"); /* ...... USER INPUTS........ */ +outDir : "~/max-out"$ /* Output directory */ + /* Serendipity basis. */ minPolyOrder_Ser : 1$ maxPolyOrder_Ser : 1$ @@ -43,7 +48,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 ~%"), @@ -63,15 +68,45 @@ for bInd : 1 thru length(bName) do ( ) ) ), + + for useSurrogate : 0 thru 1 do ( + for cdim : minDim[bInd] thru maxDim[bInd] do ( + print(sconcat("Creating vcut kernels for serendipity basis, surrogate = ", useSurrogate, " for ", cdim, "x2v, p", polyOrder)), + if (polyOrder = 1) then ( + if cdim=1 then ( + genGkSheathVcutCalcKer1x2v(fh, cdim, 2, bName[bInd], polyOrder, useSurrogate) + ) else if cdim=2 then ( + genGkSheathVcutCalcKer2x2v(fh, cdim, 2, bName[bInd], polyOrder, useSurrogate) + ) else if cdim=3 then ( + genGkSheathVcutCalcKer3x2v(fh, cdim, 2, bName[bInd], polyOrder, useSurrogate) + ) + ) + ) + ), + + for cdim : minDim[bInd] thru maxDim[bInd] do ( + print(sconcat("Creating input builder kernels for serendipity basis for ", cdim, "x2v, p", polyOrder)), + if (polyOrder = 1) then ( + if cdim=1 then ( + genGkSheathInputKer1x2v(fh, cdim, 2, bName[bInd], polyOrder) + ) else if cdim=2 then ( + genGkSheathInputKer2x2v(fh, cdim, 2, bName[bInd], polyOrder) + ) else if cdim=3 then ( + genGkSheathInputKer3x2v(fh, cdim, 2, bName[bInd], polyOrder) + ) + ) + ), close(fh) ) )$ /* 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 ~%")$ +printf(fh, "#include ~%")$ printf(fh, "#include ~%")$ printf(fh, "#include ~%")$ printf(fh, "~%")$ @@ -83,7 +118,31 @@ printf(fh, " // from Kroger ~%")$ printf(fh, " return (3.*x-x*x*x*(6. + x*x - 2.*x*x*x*x)/5.)/(1.-x*x); ~%")$ printf(fh, "}~%")$ printf(fh, "~%")$ - + +printf(fh, "// Linear interpolation of the surrogate output vcut[ng] (defined on the mu-grid~%")$ +printf(fh, "// mu_grid[ng], which the surrogate model supplies as the second half of each~%")$ +printf(fh, "// node's output block) onto mu_new[n], normalising mu by mu_ref (= T/B). Results~%")$ +printf(fh, "// are written into out[n]; clamps at the grid boundaries.~%")$ +printf(fh, "GKYL_CU_DH~%")$ +printf(fh, "static inline void~%")$ +printf(fh, "bc_sheath_gyrokinetic_surr_interpf(const float *vcut, const float *mu_grid, int ng, const double *mu_new, int n, double mu_ref, double *out)~%")$ +printf(fh, "{~%")$ +printf(fh, " for (int i = 0; i < n; i++) {~%")$ +printf(fh, " double mu = mu_new[i]/mu_ref;~%")$ +printf(fh, " if (mu <= mu_grid[0]) { out[i] = vcut[0]; continue; }~%")$ +printf(fh, " if (mu >= mu_grid[ng - 1]) { out[i] = vcut[ng - 1]; continue; }~%")$ +printf(fh, " // binary search for the bracketing interval~%")$ +printf(fh, " int lo = 0, hi = ng - 1;~%")$ +printf(fh, " while (hi - lo > 1) {~%")$ +printf(fh, " int mid = (lo + hi) >> 1;~%")$ +printf(fh, " if (mu_grid[mid] <= mu) lo = mid; else hi = mid;~%")$ +printf(fh, " }~%")$ +printf(fh, " double t = (mu - mu_grid[lo]) / (mu_grid[hi] - mu_grid[lo]);~%")$ +printf(fh, " out[i] = vcut[lo] + t * (vcut[hi] - vcut[lo]);~%")$ +printf(fh, " }~%")$ +printf(fh, "}~%")$ +printf(fh, "~%")$ + printf(fh, "EXTERN_C_BEG~%")$ printf(fh, "~%")$ @@ -94,7 +153,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 *vcutsq, const double *f, double *fRefl); ~%", edge[sI], cdim, vDims[cdim][vI], bName[bInd], polyOrder) ) ) @@ -102,6 +161,27 @@ for bInd : 1 thru length(bName) do ( ) )$ +/* Eval vcutconst kernels. */ +for cdim : minDim_Ser thru maxDim_Ser do ( + for sI : 1 thru length(edge) do ( + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_vcutsq_const_~a_~ax2v_ser_p1(const double *phi, const double *phi_wall, double q2Dm, double *vcutSq_n) ; ~%", edge[sI], cdim) + ) +)$ + +/* Eval vcutsq surrogate kernels. */ +for cdim : minDim_Ser thru maxDim_Ser do ( + for sI : 1 thru length(edge) do ( + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_vcutsq_surr_~a_~ax2v_ser_p1(const double *vmap, const float *nn_out, int n_out, const double *temperature, const double *bmag, double *vcutsq_out) ; ~%", edge[sI], cdim) + ) +)$ + +/* input builder kernels. */ +for cdim : minDim_Ser thru maxDim_Ser do ( + for sI : 1 thru length(edge) do ( + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_build_input_~a_~ax2v_ser_p1(const double *phi, const double *phi_wall, const double *density, const double *temperature, const double *bmag, const double *bimpact_angle, int n_inp, float *nn_inp_out) ; ~%", edge[sI], cdim) + ) +)$ + printf(fh, "~%")$ printf(fh, "EXTERN_C_END~%")$ printf(fh, "~%")$ diff --git a/maxima/g0/gk-sheath/sheath-1x1v.mac b/maxima/g0/gk-sheath/sheath-1x1v.mac index 92d54efd..953f8067 100644 --- a/maxima/g0/gk-sheath/sheath-1x1v.mac +++ b/maxima/g0/gk-sheath/sheath-1x1v.mac @@ -11,7 +11,7 @@ fpprec : 24$ genGkSheathKer1x1v(fh, cdim, vdim, basisFun, polyOrder) := block( [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,vmap_e,vmapSq_e,vmap_prime_e,vparLo_e,vparUp_e,surfPerpVar, surfVars,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, + surfPerpVal,sI,phiSheath_e,phiWall_e,deltaPhi_e,deltaPhi_q,vcutsq_q,fSurf_c,fSurf_e,fSurfMu_q, fReflSurfMu_e,j,nodeIdx,fSurfMu_c,fSurfMu_e,intLims,fReflSurfMu_c, xBar_v,xSqBar_v,fReflSurf_c,fReflSurf_e,fRefl_c], @@ -54,24 +54,13 @@ 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), - - /* 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. */ - phiSheath_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), - phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phiWall, bC)), - deltaPhi_e : phiSheath_e - phiWall_e, - - /* Evaluate vcutSqQ = vcut^2. q2Dm = q*2/m. */ - deltaPhi_q : deltaPhi_e, - vcutSq_q : gcfac(float(fullratsimp(-q2Dm*deltaPhi_q))), - + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double *vcutsq, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + /* Evaluate f at the z surface. */ fSurf_c : calcInnerProdList(surfVars, 1, bSurf, subst(surfPerpVar=surfPerpVal[sI], f_e)), fSurf_e : doExpand(fSurf_c, bSurf), /* Variable declarations/allocations. */ - printf(fh, " double vcutSq;~%"), printf(fh, " double fReflSurf[~a] = {0.}; ~%", length(bV)), printf(fh, "~%"), @@ -85,13 +74,13 @@ genGkSheathKer1x1v(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), + /* vcutsq is passed as a pre-computed DG field; read the single surface value. */ + printf(fh, " double vcutsq_n = vcutsq[0]; ~%"), 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, " if (vcutsq_n <= vparAbsSqLo) { // absorb (no reflection) ~%"), printf(fh, "~%"), writeCExprsWithZeros1(fReflSurf, makelist(0.,j,1,length(bV))), @@ -99,7 +88,7 @@ genGkSheathKer1x1v(fh, cdim, vdim, basisFun, polyOrder) := block( /* 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, " } else if (vcutsq_n > vparAbsSqUp) { // full reflection ~%"), printf(fh, "~%"), /* Project f onto vpar basis (bV). */ @@ -124,7 +113,7 @@ genGkSheathKer1x1v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " if (wv > 0.) {~%"), printf(fh, " // vcut in logical space.~%"), - printf(fh, " double vcut_l = 2.*(sqrt(vcutSq)-wv)/dv;~%"), + printf(fh, " double vcut_l = 2.*(sqrt(vcutsq_n)-wv)/dv;~%"), printf(fh, "~%"), intLims : [[-1,vcut_l]], @@ -135,7 +124,7 @@ genGkSheathKer1x1v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " } else {~%"), printf(fh, " // vcut in logical space.~%"), - printf(fh, " double vcut_l = 2.*(-sqrt(vcutSq)-wv)/dv;~%"), + printf(fh, " double vcut_l = 2.*(-sqrt(vcutsq_n)-wv)/dv;~%"), printf(fh, "~%"), intLims : [[vcut_l,1]], diff --git a/maxima/g0/gk-sheath/sheath-1x2v.mac b/maxima/g0/gk-sheath/sheath-1x2v.mac index 54f41043..ab43f903 100644 --- a/maxima/g0/gk-sheath/sheath-1x2v.mac +++ b/maxima/g0/gk-sheath/sheath-1x2v.mac @@ -11,9 +11,10 @@ fpprec : 24$ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,surfVars, 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, + 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, + vcutsq_e,nodesVparMu,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 vcutsq on mu basis (bMu). */ + vcutsq_e : doExpand1(vcutsq, bMu), + + /* Evaluate vcutsq at (mu_j) */ + vcutsq_q : evAtNodes(vcutsq_e, nodesMu, zvSurfVars), /* Generate separate kernels for lower/upper boundaries. */ edge : ["lower", "upper"], @@ -70,18 +77,8 @@ 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), - - /* 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. */ - phiSheath_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), - phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phiWall, bC)), - deltaPhi_e : phiSheath_e - phiWall_e, - - /* Evaluate vcutSqQ = vcut^2. q2Dm = q*2/m. */ - deltaPhi_q : deltaPhi_e, - vcutSq_q : gcfac(float(fullratsimp(-q2Dm*deltaPhi_q))), - + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double *vcutsq, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + /* Evaluate f at the z surface. */ fSurf_c : calcInnerProdList(surfVars, 1, bSurf, subst(surfPerpVar=surfPerpVal[sI], f_e)), fSurf_e : doExpand(fSurf_c, bSurf), @@ -91,10 +88,10 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfMu_q : evAtNodes(fSurf_e, nodesMu, zvSurfVars), /* Variable declarations/allocations. */ - printf(fh, " double vcutSq;~%"), + printf(fh, " double vcutsq_n;~%"), 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 +102,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(fullratsimp((vparUp + vparLo)/2))), + printf(fh, " double dv = ~a;~%", float(fullratsimp(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,9 +123,27 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfMu_c : gcfac(fullratsimp(calcInnerProdList([vpar], 1, bVpar, fSurfMu_q[nodeIdx]))), fSurfMu_e : doExpand(fSurfMu_c, bVpar), + /* Get vcutsq at this mu node from the DG field. */ + printf(fh, " vcutsq_n = ~a; ~%", gcfac(float(fullratsimp(vcutsq_q[nodeIdx])))), + printf(fh, "~%"), + + /* Check regime at this mu node. */ + printf(fh, " if (vcutsq_n <= vparAbsSqLo) { // absorb (no reflection) ~%"), + printf(fh, "~%"), + writeCExprsWithZeros1(fReflSurfMu[j-1], makelist(0,i,1,length(bVpar))), + printf(fh, "~%"), + + printf(fh, " } else if (vcutsq_n > 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;~%"), + printf(fh, " double vcut_l = 2.*(sqrt(vcutsq_n)-wv)/dv;~%"), printf(fh, "~%"), intLims : [[-1,vcut_l]], @@ -166,7 +154,7 @@ genGkSheathKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( printf(fh, " } else {~%"), printf(fh, " // vcut in logical space.~%"), - printf(fh, " double vcut_l = 2.*(-sqrt(vcutSq)-wv)/dv;~%"), + printf(fh, " double vcut_l = 2.*(-sqrt(vcutsq_n)-wv)/dv;~%"), printf(fh, "~%"), intLims : [[vcut_l,1]], @@ -193,6 +181,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 +199,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..ffc406d4 100644 --- a/maxima/g0/gk-sheath/sheath-2x2v.mac +++ b/maxima/g0/gk-sheath/sheath-2x2v.mac @@ -11,9 +11,9 @@ fpprec : 24$ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,surfVars,surfVarsC, 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, + 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,vcutsq_e,nodesVparMu,tempVars1,tempVars2], pdim : cdim+vdim, @@ -36,12 +36,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 +75,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^2 on (x,mu) basis (bXMu). */ + vcutsq_e : doExpand1(vcutsq, bXMu), + + /* Evaluate vcutsq at (x,mu)_ij nodes. */ + vcutsq_q : evAtNodes(vcutsq_e, nodesXMu, zvSurfVars), + /* Generate separate kernels for lower/upper boundaries. */ edge : ["lower", "upper"], surfPerpVal : [-1, 1], @@ -80,17 +88,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), - - /* 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. */ - phiSheath_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), - phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phiWall, bC)), - deltaPhi_e : phiSheath_e - phiWall_e, - - /* Evaluate vcutSq = vcut^2 at x nodes. */ - deltaPhi_q : evAtNodes(deltaPhi_e, nodesX, surfVarsC), - vcutSq_q : gcfac(float(fullratsimp(-q2Dm*deltaPhi_q))), + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double *vcutsq, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), /* Evaluate f at the z surface. */ fSurf_c : calcInnerProdList(surfVars, 1, bSurf, subst(surfPerpVar=surfPerpVal[sI], f_e)), @@ -101,7 +99,7 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfX_q : evAtNodes(fSurf_e, nodesX, surfVarsC), fSurfXMu_q : evAtNodes(fSurf_e, nodesXMu, zvSurfVars), - printf(fh, " double vcutSq;~%"), + printf(fh, " double vcutsq_n;~%"), printf(fh, " double fReflSurfX[~a][~a] = {0.}; ~%", numNodesX, length(bV)), /* Reflected f at (x)_i nodes. */ @@ -117,51 +115,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(fullratsimp((vparUp + vparLo)/2))), + printf(fh, " dv = ~a;~%", float(fullratsimp(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,9 +144,27 @@ genGkSheathKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfXMu_c : gcfac(fullratsimp(calcInnerProdList([vpar], 1, bVpar, fSurfXMu_q[nodeIdx]))), fSurfXMu_e : doExpand(fSurfXMu_c, bVpar), + /* Get vcut^2 at this mu node. */ + printf(fh, " vcutsq_n = ~a; ~%", gcfac(float(fullratsimp(vcutsq_q[nodeIdx])))), + printf(fh, "~%"), + + /* Check regime at this (x,mu) node. */ + printf(fh, " if (vcutsq_n <= vparAbsSqLo) { // absorb (no reflection) ~%"), + printf(fh, "~%"), + writeCExprsWithZeros1(fReflSurfXMu[j-1], makelist(0,k,1,length(bVpar))), + printf(fh, "~%"), + + printf(fh, " } else if (vcutsq_n > 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;~%"), + printf(fh, " double vcut_l = 2.*(sqrt(vcutsq_n)-wv)/dv;~%"), printf(fh, "~%"), intLims : [[-1,vcut_l]], @@ -179,10 +172,10 @@ 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;~%"), + printf(fh, " double vcut_l = 2.*(-sqrt(vcutsq_n)-wv)/dv;~%"), printf(fh, "~%"), intLims : [[vcut_l,1]], @@ -198,15 +191,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 +220,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..0e3422be 100644 --- a/maxima/g0/gk-sheath/sheath-3x2v.mac +++ b/maxima/g0/gk-sheath/sheath-3x2v.mac @@ -12,9 +12,10 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,surfVars,surfVarsC, bSurf,bV,bVpar,bMu,bXY,numBsurf,zvSurfVars,nodesXY,nodesXYMu,nodesMu,numNodesXY,numNodesXYMu, 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, + 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, + vcutsq_vpar0_e,vcutsq_e,nodesVparMu,vcutsq_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 vcutsq on (x,y,mu) basis (bXYMu). */ + vcutsq_e : doExpand1(vcutsq, bXYMu), + + /* Evaluate vcutsq at (x,y,mu)_ij nodes. */ + vcutsq_q : evAtNodes(vcutsq_e, nodesXYMu, zvSurfVars), + /* Generate separate kernels for lower/upper boundaries. */ edge : ["lower", "upper"], surfPerpVal : [-1, 1], @@ -81,17 +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), - - /* 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. */ - phiSheath_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), - phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phiWall, bC)), - deltaPhi_e : phiSheath_e - phiWall_e, - - /* Evaluate vcutSq = vcut^2 at x,y nodes. */ - deltaPhi_q : evAtNodes(deltaPhi_e, nodesXY, surfVarsC), - vcutSq_q : gcfac(float(fullratsimp(-q2Dm*deltaPhi_q))), + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_reflectedf_~a_~ax~av_~a_p~a(const double *vmap, const double *vcutsq, const double *f, double *fRefl) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), /* Evaluate f at the z surface. */ fSurf_c : calcInnerProdList(surfVars, 1, bSurf, subst(surfPerpVar=surfPerpVal[sI], f_e)), @@ -102,7 +101,7 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfXY_q : evAtNodes(fSurf_e, nodesXY, surfVarsC), fSurfXYMu_q : evAtNodes(fSurf_e, nodesXYMu, zvSurfVars), - printf(fh, " double vcutSq;~%"), + printf(fh, " double vcutsq_n;~%"), printf(fh, " double fReflSurfXY[~a][~a] = {0.}; ~%", numNodesXY, length(bV)), /* Reflected f at (x,y)_i nodes. */ @@ -118,51 +117,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(fullratsimp((vparUp + vparLo)/2))), + printf(fh, " dv = ~a;~%", float(fullratsimp(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,9 +144,27 @@ genGkSheathKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( fSurfXYMu_c : gcfac(fullratsimp(calcInnerProdList([vpar], 1, bVpar, fSurfXYMu_q[nodeIdx]))), fSurfXYMu_e : doExpand(fSurfXYMu_c, bVpar), + /* Get vcutsq at this (x,y,mu) node from the DG field. */ + printf(fh, " vcutsq_n = ~a; ~%", gcfac(float(fullratsimp(vcutsq_q[nodeIdx])))), + printf(fh, "~%"), + + /* Check regime at this (x,y,mu) node. */ + printf(fh, " if (vcutsq_n <= vparAbsSqLo) { // absorb (no reflection) ~%"), + printf(fh, "~%"), + writeCExprsWithZeros1(fReflSurfXYMu[j-1], makelist(0,k,1,length(bVpar))), + printf(fh, "~%"), + + printf(fh, " } else if (vcutsq_n > 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;~%"), + printf(fh, " double vcut_l = 2.*(sqrt(vcutsq_n)-wv)/dv;~%"), printf(fh, "~%"), intLims : [[-1,vcut_l]], @@ -180,10 +172,10 @@ 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;~%"), + printf(fh, " double vcut_l = 2.*(-sqrt(vcutsq_n)-wv)/dv;~%"), printf(fh, "~%"), intLims : [[vcut_l,1]], @@ -194,20 +186,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 +220,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) diff --git a/maxima/g0/gk-sheath/sheath_surrogate-1x2v.mac b/maxima/g0/gk-sheath/sheath_surrogate-1x2v.mac new file mode 100644 index 00000000..636c2bad --- /dev/null +++ b/maxima/g0/gk-sheath/sheath_surrogate-1x2v.mac @@ -0,0 +1,149 @@ +/* + Gyrokinetic sheath function for 1x2v kernel. +*/ +load("modal-basis")$ +load("out-scripts")$ +load("nodal_operations/nodal_functions")$ +load("utilities_gyrokinetic")$ +load(stringproc)$ +fpprec : 24$ + +genGkSheathVcutCalcKer1x2v(fh, cdim, vdim, basisFun, polyOrder, useSurrogate) := block( + [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,d,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,surfVarsC, + bMu,bCperpMu,zvSurfVars,nodesMu,numNodesMu,j,bNMu,edge,surfPerpVal,sI,surfVarsCmu, + phi_e,phiWall_e,density_e,temperature_e,bmag_e,bimpactAngle_e,vcut_e,vcut_c], + + pdim : cdim+vdim, + + /* Get desired basis. */ + [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder), + NP : length(bP), + varsV : copylist(varsP), for d : 1 thru cdim do (varsV : delete(varsC[d],varsV)), + + [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP), + + /* Get name of last config space dimension, which is always assumed to be + the direction parallel to the magnetic field (z). */ + surfPerpVar : varsC[cdim], + surfVarsC : delete(surfPerpVar, varsC), /* = [] */ + surfVarsCmu : append(surfVarsC, [mu]), /* = mu. */ + + /* Set up mu basis and perp conf + mu basis. */ + bMu : basisFromVars("ser",[mu],polyOrder), + bCperpMu : basisFromVars("ser",surfVarsCmu,polyOrder), + + zvSurfVars : append(surfVarsC, [mu]), /* = mu. */ + /* Set up surface nodes. */ + /* nodesMu : gaussOrdGkHyb(polyOrder+1, [], [mu]) */ + nodesMu : gaussOrd(polyOrder+1, length([mu])), + numNodesMu : length(nodesMu), + + /* Get nodal basis sets. */ + bNMu : getVarsNodalBasisWithNodesHyb("gkhyb", 0, 1, [mu], nodesMu), + + /* Generate separate kernels for lower/upper boundaries. */ + edge : ["lower", "upper"], + surfPerpVal : [-1, 1], + + for sI : 1 thru length(edge) do ( + printf(fh, "~%"), + if useSurrogate = 1 then ( + /* vcut kernel: takes pre-computed NN output nn_out (SRGRZ_N_MU values), evaluates vcut at mu nodes via interpolation, and projects to DG. */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_vcutsq_surr_~a_~ax~av_~a_p~a(const double *vmap, const float *nn_out, int n_out, const double *temperature, const double *bmag, double *vcutsq_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + /* Only temperature and bmag are needed to compute mu_ref = T/B. */ + temperature_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(temperature, bC)), + bmag_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(bmag, bC)), + printf(fh, " double temp_n = ~a; ~%", gcfac(float(fullratsimp(temperature_e)))), + printf(fh, " double bmag_n = ~a; ~%", gcfac(float(fullratsimp(bmag_e)))), + printf(fh, " double mu_ref = temp_n/bmag_n;~%"), + printf(fh, " double vthSq = temp_n/GKYL_ELECTRON_MASS;~%"), + printf(fh, " double mu;~%"), + printf(fh, " double vcut;~%") + ) else ( + /* vcut kernel: calculates vcut using conducting sheath model (no surrogate). */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_vcutsq_const_~a_~ax~av_~a_p~a(const double *phi, const double *phi_wall, double q2Dm, double *vcutsq_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + phi_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), + phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi_wall, bC)), + deltaPhi_e : phi_e - phiWall_e, + deltaPhi_q : deltaPhi_e, + printf(fh, " double deltaphi_n = ~a; ~%", gcfac(float(fullratsimp(deltaPhi_q)))) + ), + printf(fh, " double vcutsq_n[~a] = {0.}; ~%", numNodesMu), + printf(fh, "~%"), + + for j : 1 thru numNodesMu do ( + printf(fh, " // node (mu)_~a ~%", j-1), + if useSurrogate = 1 then ( + mu_e : subst(mu=nodesMu[j][1], vmap_e[2]), + printf(fh, " mu = ~a; ~%", float(expand(mu_e))), + printf(fh, " bc_sheath_gyrokinetic_surr_interpf(nn_out+0*n_out, nn_out+0*n_out+n_out/2, n_out/2, &mu, 1, mu_ref, &vcut);~%"), + printf(fh, " vcutsq_n[~a] = vcut*vcut * vthSq;~%", j-1) + ) else ( + printf(fh, " vcutsq_n[~a] = -q2Dm * deltaphi_n;~%", j-1) + ), + printf(fh, "~%") + ), + + /* Reconstruct the mu modal representation via a nodal-modal conversion. */ + vcut_c : calcInnerProdList([mu], 1, bMu, doExpand1(vcutsq_n, bNMu)), + + /* Write coefficients. */ + writeCExprsWithZerosNoExpand1(vcutsq_out, vcut_c), + + printf(fh, "}~%") + ) +)$ + +/* Generate the input kernel: constructs NN input features at the sheath surface + and runs gkyl_kann_net_apply, writing SRGRZ_N_MU raw outputs to *out. + For 1x the perpendicular space is 0D so there is a single inference call. */ +genGkSheathInputKer1x2v(fh, cdim, vdim, basisFun, polyOrder) := block( + [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,surfVarsC, + phi_e,phiWall_e,density_e,temperature_e,bmag_e,bimpactAngle_e, + edge,surfPerpVal,sI], + + pdim : cdim+vdim, + + [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder), + NP : length(bP), + varsV : copylist(varsP), for d : 1 thru cdim do (varsV : delete(varsC[d],varsV)), + + [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP), + + surfPerpVar : varsC[cdim], + surfVarsC : delete(surfPerpVar, varsC), /* = [] */ + + edge : ["lower", "upper"], + surfPerpVal : [-1, 1], + + for sI : 1 thru length(edge) do ( + printf(fh, "~%"), + /* Input kernel: constructs NN input features at the sheath surface, + calls gkyl_kann_net_apply, stores SRGRZ_N_MU outputs to *out. */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_build_input_~a_~ax~av_~a_p~a(const double *phi, const double *phi_wall, const double *density, const double *temperature, const double *bmag, const double *bimpact_angle, int n_inp, float *nn_inp_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + + phi_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), + phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi_wall, bC)), + density_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(density, bC)), + temperature_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(temperature, bC)), + bmag_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(bmag, bC)), + bimpactAngle_e : doExpand1(bimpact_angle, [1]), /* 0D surface for 1x */ + + printf(fh, " double phi_n = ~a; ~%", gcfac(float(fullratsimp(phi_e)))), + printf(fh, " double phi_wall_n = ~a; ~%", gcfac(float(fullratsimp(phiWall_e)))), + printf(fh, " double dens_n = ~a; ~%", gcfac(float(fullratsimp(density_e)))), + printf(fh, " double temp_n = ~a; ~%", gcfac(float(fullratsimp(temperature_e)))), + printf(fh, " double bmag_n = ~a; ~%", gcfac(float(fullratsimp(bmag_e)))), + printf(fh, " double angle_n = ~a; ~%", gcfac(float(fullratsimp(bimpactAngle_e)))), + printf(fh, "~%"), + printf(fh, " nn_inp_out[0+n_inp*0] = (float)(angle_n*(180.0/GKYL_PI));~%"), + printf(fh, " nn_inp_out[1+n_inp*0] = (float)((1.0/bmag_n)*sqrt(GKYL_ELECTRON_MASS*fmax(0.0, dens_n)/GKYL_EPSILON0));~%"), + printf(fh, " nn_inp_out[2+n_inp*0] = (float)(fmax(0.0, (GKYL_ELEMENTARY_CHARGE*(phi_n-phi_wall_n))/temp_n));~%"), + /* Note: In the future we can easily extend this to more parameters by adding more nn_inp_out entries, like: + nn_inp_out[3+n_inp*0] = (float)(tperp/tpar); // temperature anisotropy + */ + printf(fh, "~%"), + printf(fh, "}~%") + ) +)$ + diff --git a/maxima/g0/gk-sheath/sheath_surrogate-2x2v.mac b/maxima/g0/gk-sheath/sheath_surrogate-2x2v.mac new file mode 100644 index 00000000..2fd68a31 --- /dev/null +++ b/maxima/g0/gk-sheath/sheath_surrogate-2x2v.mac @@ -0,0 +1,198 @@ +/* + Gyrokinetic sheath function for 1x2v kernel. +*/ +load("modal-basis")$ +load("out-scripts")$ +load("nodal_operations/nodal_functions")$ +load("utilities_gyrokinetic")$ +load(stringproc)$ +fpprec : 24$ + +genGkSheathVcutCalcKer2x2v(fh, cdim, vdim, basisFun, polyOrder, useSurrogate) := block( + [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,d,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,VarsX, + bMu,bX,bXMu,VarsXMu,nodesMu,nodesX,numNodesMu,numNodesX,j,i,bNMu,bNX,edge,surfPerpVal,sI,surfVarsCmu, + phi_e,phiWall_e,density_e,temperature_e,bmag_e,bimpactAngle_e,vcut_e,vcut_c, + phi_q,phiWall_q,density_q,temperature_q,bmag_q,bimpactAngle_q,vcut,vcutX_e], + + pdim : cdim+vdim, + + /* Get desired basis. */ + [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder), + NP : length(bP), + varsV : copylist(varsP), for d : 1 thru cdim do (varsV : delete(varsC[d],varsV)), + + [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP), + + /* Get name of last config space dimension, which is always assumed to be + the direction parallel to the magnetic field (z). */ + surfPerpVar : varsC[cdim], + VarsX : delete(surfPerpVar, varsC), /* = x. */ + surfVarsCmu : append(VarsX, [mu]), /* = x,mu. */ + + /* Set up mu basis and perp conf + mu basis. */ + bMu : basisFromVars("ser",[mu],polyOrder), + bX : basisFromVars("ser",[x],polyOrder), + bXMu : basisFromVars("ser",surfVarsCmu,polyOrder), + + VarsXMu : append(VarsX, [mu]), /* = x,mu. */ + /* Set up surface nodes. */ + nodesMu : gaussOrdGkHyb(polyOrder+1, [], [mu]), + nodesX : gaussOrdGkHyb(polyOrder+1, VarsX, []), + + numNodesMu : length(nodesMu), + numNodesX : length(nodesX), + + /* Get nodal basis sets. */ + bNMu : getVarsNodalBasisWithNodesHyb("gkhyb", 0, 1, [mu], nodesMu), + bNX : getVarsNodalBasisWithNodesHyb("gkhyb", length(VarsX), polyOrder, VarsX, nodesX), + + /* Generate separate kernels for lower/upper boundaries. */ + edge : ["lower", "upper"], + surfPerpVal : [-1, 1], + + for sI : 1 thru length(edge) do ( + printf(fh, "~%"), + if useSurrogate = 1 then ( + /* vcut kernel: takes pre-computed NN output nn_out (numNodesX*SRGRZ_N_MU values), + evaluates vcut at mu nodes via interpolation per x-node, projects to DG. */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_vcutsq_surr_~a_~ax~av_~a_p~a(const double *vmap, const float *nn_out, int n_out, const double *temperature, const double *bmag, double *vcutsq_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + /* Only temperature and bmag are needed: mu_ref = T/B varies per x-node. */ + temperature_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(temperature, bC)), + bmag_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(bmag, bC)), + temperature_q : evAtNodes(temperature_e, nodesX, VarsX), + bmag_q : evAtNodes(bmag_e, nodesX, VarsX), + printf(fh, " double temp_n; ~%"), + printf(fh, " double bmag_n; ~%"), + printf(fh, " double mu_ref; ~%"), + printf(fh, " double vthSq; ~%"), + printf(fh, " double mu;~%"), + printf(fh, " double vcut;~%") + ) else ( + /* vcut kernel: calculates vcut using conducting sheath model (no surrogate). */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_vcutsq_const_~a_~ax~av_~a_p~a(const double *phi, const double *phi_wall, double q2Dm, double *vcutsq_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + phi_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), + phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi_wall, bC)), + deltaPhi_e : phi_e - phiWall_e, + deltaPhi_q : evAtNodes(deltaPhi_e, nodesX, VarsX), + printf(fh, " double deltaphi_n; ~%") + ), + for i : 1 thru numNodesX do ( + printf(fh, " double vcutsq_n_~a[~a] = {0.}; ~%", i-1, numNodesMu) + ), + printf(fh, "~%"), + + vcutX_e : makelist(0,i,1,numNodesX), + + for i : 1 thru numNodesX do ( + printf(fh, " // node (x)_~a ~%", i-1), + if useSurrogate = 1 then ( + printf(fh, " temp_n = ~a; ~%", gcfac(float(fullratsimp(temperature_q[i])))), + printf(fh, " bmag_n = ~a; ~%", gcfac(float(fullratsimp(bmag_q[i])))), + printf(fh, " mu_ref = temp_n/bmag_n;~%"), + printf(fh, " vthSq = temp_n/GKYL_ELECTRON_MASS;~%") + ) else ( + printf(fh, " deltaphi_n = ~a; ~%", gcfac(float(fullratsimp(deltaPhi_q[i])))) + ), + printf(fh, "~%"), + + for j : 1 thru numNodesMu do ( + printf(fh, " // node (mu)_~a ~%", j-1), + if useSurrogate = 1 then ( + mu_e : subst(mu=nodesMu[j][1], vmap_e[2]), + printf(fh, " mu = ~a; ~%", float(expand(mu_e))), + printf(fh, " bc_sheath_gyrokinetic_surr_interpf(nn_out+~a*n_out, nn_out+~a*n_out+n_out/2, n_out/2, &mu, 1, mu_ref, &vcut);~%", i-1, i-1), + printf(fh, " vcutsq_n_~a[~a] = vcut*vcut * vthSq;~%", i-1, j-1) + ) else ( + printf(fh, " vcutsq_n_~a[~a] = -q2Dm * deltaphi_n;~%", i-1, j-1) + ), + printf(fh, "~%") + ), + + vcutX_e[i] : doExpand(makelist(concat(vcutsq_n_, i-1)[j-1], j, 1, numNodesMu), bNMu) + ), + + /* We project the mu expansion at each x node onto the x,mu basis. */ + vcut_c : gcfac(fullratsimp(calcInnerProdList(VarsX, 1, bX, doExpand(vcutX_e, bNX)))), + vcut_e : doExpand(vcut_c, bX), + + /* Project expansion onto confperp + mu basis */ + vcut_c : gcfac(fullratsimp(calcInnerProdList(VarsXMu, 1, bXMu, vcut_e))), + + /* Write coefficients. */ + writeCExprsWithZerosNoExpand1(vcutsq_out, vcut_c), + + printf(fh, "}~%") + ) +)$ + +/* Generate the infer kernel for 2x2v: loops over x-nodes, evaluates NN features + at each node, calls gkyl_kann_net_apply, writes SRGRZ_N_MU outputs per node + to out[node*SRGRZ_N_MU .. (node+1)*SRGRZ_N_MU - 1]. */ +genGkSheathInputKer2x2v(fh, cdim, vdim, basisFun, polyOrder) := block( + [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,VarsX, + bX,nodesX,numNodesX,i,edge,surfPerpVal,sI, + phi_e,phiWall_e,density_e,temperature_e,bmag_e,bimpactAngle_e, + phi_q,phiWall_q,density_q,temperature_q,bmag_q,bimpactAngle_q], + + pdim : cdim+vdim, + + [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder), + NP : length(bP), + varsV : copylist(varsP), for d : 1 thru cdim do (varsV : delete(varsC[d],varsV)), + + [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP), + + surfPerpVar : varsC[cdim], + VarsX : delete(surfPerpVar, varsC), /* = x. */ + + bX : basisFromVars("ser",[x],polyOrder), + nodesX : gaussOrdGkHyb(polyOrder+1, VarsX, []), + numNodesX : length(nodesX), + + edge : ["lower", "upper"], + surfPerpVal : [-1, 1], + + for sI : 1 thru length(edge) do ( + printf(fh, "~%"), + /* Input kernel: constructs NN input features at the sheath surface, + calls gkyl_kann_net_apply, stores SRGRZ_N_MU outputs to *out. */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_build_input_~a_~ax~av_~a_p~a(const double *phi, const double *phi_wall, const double *density, const double *temperature, const double *bmag, const double *bimpact_angle, int n_inp, float *nn_inp_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + + phi_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), + phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi_wall, bC)), + density_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(density, bC)), + temperature_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(temperature, bC)), + bmag_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(bmag, bC)), + bimpactAngle_e : doExpand1(bimpact_angle, bX), /* bimpact_angle lives on x basis */ + + phi_q : evAtNodes(phi_e, nodesX, VarsX), + phiWall_q : evAtNodes(phiWall_e, nodesX, VarsX), + density_q : evAtNodes(density_e, nodesX, VarsX), + temperature_q : evAtNodes(temperature_e, nodesX, VarsX), + bmag_q : evAtNodes(bmag_e, nodesX, VarsX), + bimpactAngle_q : evAtNodes(bimpactAngle_e, nodesX, VarsX), + + printf(fh, " double phi_n, phi_wall_n, dens_n, temp_n, bmag_n, angle_n;~%"), + printf(fh, "~%"), + + for i : 1 thru numNodesX do ( + printf(fh, " // node (x)_~a ~%", i-1), + printf(fh, " phi_n = ~a; ~%", gcfac(float(fullratsimp(phi_q[i])))), + printf(fh, " phi_wall_n = ~a; ~%", gcfac(float(fullratsimp(phiWall_q[i])))), + printf(fh, " dens_n = ~a; ~%", gcfac(float(fullratsimp(density_q[i])))), + printf(fh, " temp_n = ~a; ~%", gcfac(float(fullratsimp(temperature_q[i])))), + printf(fh, " bmag_n = ~a; ~%", gcfac(float(fullratsimp(bmag_q[i])))), + printf(fh, " angle_n = ~a; ~%", gcfac(float(fullratsimp(bimpactAngle_q[i])))), + printf(fh, " nn_inp_out[0+n_inp*~a] = (float)(angle_n*(180.0/GKYL_PI));~%", i-1), + printf(fh, " nn_inp_out[1+n_inp*~a] = (float)((1.0/bmag_n)*sqrt(GKYL_ELECTRON_MASS*fmax(0.0, dens_n)/GKYL_EPSILON0));~%", i-1), + printf(fh, " nn_inp_out[2+n_inp*~a] = (float)(fmax(0.0, (GKYL_ELEMENTARY_CHARGE*(phi_n-phi_wall_n))/temp_n));~%", i-1), + /* Note: In the future we can easily extend this to more parameters by adding more nn_inp_out entries, like: + nn_inp_out[3+n_inp*~a] = (float)(tperp/tpar); // temperature anisotropy + */ + printf(fh, "~%") + ), + + printf(fh, "}~%") + ) +)$ + diff --git a/maxima/g0/gk-sheath/sheath_surrogate-3x2v.mac b/maxima/g0/gk-sheath/sheath_surrogate-3x2v.mac new file mode 100644 index 00000000..a849d2ed --- /dev/null +++ b/maxima/g0/gk-sheath/sheath_surrogate-3x2v.mac @@ -0,0 +1,198 @@ +/* + Gyrokinetic sheath function for 3x2v kernel. +*/ +load("modal-basis")$ +load("out-scripts")$ +load("nodal_operations/nodal_functions")$ +load("utilities_gyrokinetic")$ +load(stringproc)$ +fpprec : 24$ + +genGkSheathVcutCalcKer3x2v(fh, cdim, vdim, basisFun, polyOrder, useSurrogate) := block( + [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,d,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,VarsXY, + bMu,bXY,bXYMu,VarsXYMu,nodesMu,nodesXY,numNodesMu,numNodesXY,j,i,bNMu,bNXY,edge,surfPerpVal,sI,surfVarsCmu, + phi_e,phiWall_e,density_e,temperature_e,bmag_e,bimpactAngle_e,vcut_e,vcut_c, + phi_q,phiWall_q,density_q,temperature_q,bmag_q,bimpactAngle_q,vcut,vcutXY_e], + + pdim : cdim+vdim, + + /* Get desired basis. */ + [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder), + NP : length(bP), + varsV : copylist(varsP), for d : 1 thru cdim do (varsV : delete(varsC[d],varsV)), + + [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP), + + /* Get name of last config space dimension, which is always assumed to be + the direction parallel to the magnetic field (z). */ + surfPerpVar : varsC[cdim], + VarsXY : delete(surfPerpVar, varsC), /* = [x, y]. */ + surfVarsCmu : append(VarsXY, [mu]), /* = [x, y, mu]. */ + + /* Set up mu basis and perp conf + mu basis. */ + bMu : basisFromVars("ser",[mu],polyOrder), + bXY : basisFromVars("ser",VarsXY,polyOrder), + bXYMu : basisFromVars("ser",surfVarsCmu,polyOrder), + + VarsXYMu : append(VarsXY, [mu]), /* = [x, y, mu]. */ + /* Set up surface nodes. */ + nodesMu : gaussOrdGkHyb(polyOrder+1, [], [mu]), + nodesXY : gaussOrdGkHyb(polyOrder+1, VarsXY, []), + + numNodesMu : length(nodesMu), + numNodesXY : length(nodesXY), + + /* Get nodal basis sets. */ + bNMu : getVarsNodalBasisWithNodesHyb("gkhyb", 0, 1, [mu], nodesMu), + bNXY : getVarsNodalBasisWithNodesHyb("gkhyb", length(VarsXY), polyOrder, VarsXY, nodesXY), + + /* Generate separate kernels for lower/upper boundaries. */ + edge : ["lower", "upper"], + surfPerpVal : [-1, 1], + + for sI : 1 thru length(edge) do ( + printf(fh, "~%"), + if useSurrogate = 1 then ( + /* vcut kernel: takes pre-computed NN output nn_out (numNodesXY*SRGRZ_N_MU values), + evaluates vcut at mu nodes via interpolation per (x,y)-node, projects to DG. */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_vcutsq_surr_~a_~ax~av_~a_p~a(const double *vmap, const float *nn_out, int n_out, const double *temperature, const double *bmag, double *vcutsq_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + /* Only temperature and bmag are needed: mu_ref = T/B varies per (x,y)-node. */ + temperature_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(temperature, bC)), + bmag_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(bmag, bC)), + temperature_q : evAtNodes(temperature_e, nodesXY, VarsXY), + bmag_q : evAtNodes(bmag_e, nodesXY, VarsXY), + printf(fh, " double temp_n; ~%"), + printf(fh, " double bmag_n; ~%"), + printf(fh, " double mu_ref; ~%"), + printf(fh, " double mu;~%"), + printf(fh, " double vthSq;~%"), + printf(fh, " double vcut;~%") + ) else ( + /* vcut kernel: calculates vcut using conducting sheath model (no surrogate). */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_vcutsq_const_~a_~ax~av_~a_p~a(const double *phi, const double *phi_wall, double q2Dm, double *vcutsq_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + phi_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), + phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi_wall, bC)), + deltaPhi_e : phi_e - phiWall_e, + deltaphi_q : evAtNodes(deltaPhi_e, nodesXY, VarsXY), + printf(fh, " double deltaphi_n; ~%") + ), + + for i : 1 thru numNodesXY do ( + printf(fh, " double vcutsq_n_~a[~a] = {0.}; ~%", i-1, numNodesMu) + ), + printf(fh, "~%"), + + vcutXY_e : makelist(0,i,1,numNodesXY), + + for i : 1 thru numNodesXY do ( + printf(fh, " // node (x,y)_~a ~%", i-1), + if useSurrogate = 1 then ( + printf(fh, " temp_n = ~a; ~%", gcfac(float(fullratsimp(temperature_q[i])))), + printf(fh, " bmag_n = ~a; ~%", gcfac(float(fullratsimp(bmag_q[i])))), + printf(fh, " mu_ref = temp_n/bmag_n;~%"), + printf(fh, " vthSq = temp_n/GKYL_ELECTRON_MASS;~%") + ) else ( + printf(fh, " deltaphi_n = ~a; ~%", gcfac(float(fullratsimp(deltaphi_q[i])))) + ), + printf(fh, "~%"), + + for j : 1 thru numNodesMu do ( + printf(fh, " // node (mu)_~a ~%", j-1), + if useSurrogate = 1 then ( + mu_e : subst(mu=nodesMu[j][1], vmap_e[2]), + printf(fh, " mu = ~a; ~%", float(expand(mu_e))), + printf(fh, " bc_sheath_gyrokinetic_surr_interpf(nn_out+~a*n_out, nn_out+~a*n_out+n_out/2, n_out/2, &mu, 1, mu_ref, &vcut);~%", i-1, i-1), + printf(fh, " vcutsq_n_~a[~a] = vcut*vcut * vthSq;~%", i-1, j-1) + ) else ( + printf(fh, " vcutsq_n_~a[~a] = -q2Dm * deltaphi_n;~%", i-1, j-1) + ), + printf(fh, "~%") + ), + + vcutXY_e[i] : doExpand(makelist(concat(vcutsq_n_, i-1)[j-1], j, 1, numNodesMu), bNMu) + ), + + /* We project the mu expansion at each perpendicular config space node onto the perp config + mu basis. */ + vcut_c : gcfac(fullratsimp(calcInnerProdList(VarsXY, 1, bXY, doExpand(vcutXY_e, bNXY)))), + vcut_e : doExpand(vcut_c, bXY), + + /* Project expansion onto confperp + mu basis */ + vcut_c : gcfac(fullratsimp(calcInnerProdList(VarsXYMu, 1, bXYMu, vcut_e))), + + /* Write coefficients. */ + writeCExprsWithZerosNoExpand1(vcutsq_out, vcut_c), + + printf(fh, "}~%") + ) +)$ + +/* Generate the infer kernel for 3x2v: loops over (x,y)-nodes, evaluates NN features + at each node, calls gkyl_kann_net_apply, writes SRGRZ_N_MU outputs per node + to out[node*SRGRZ_N_MU .. (node+1)*SRGRZ_N_MU - 1]. */ +genGkSheathInputKer3x2v(fh, cdim, vdim, basisFun, polyOrder) := block( + [pdim,varsC,bC,varsP,bP,vSub,NP,varsV,d,vmap_e,vmapSq_e,vmap_prime_e,surfPerpVar,VarsXY, + bXY,nodesXY,numNodesXY,i,edge,surfPerpVal,sI, + phi_e,phiWall_e,density_e,temperature_e,bmag_e,bimpactAngle_e, + phi_q,phiWall_q,density_q,temperature_q,bmag_q,bimpactAngle_q], + + pdim : cdim+vdim, + + [varsC,bC,varsP,bP,vSub] : loadGkBasis(basisFun, cdim, vdim, polyOrder), + NP : length(bP), + varsV : copylist(varsP), for d : 1 thru cdim do (varsV : delete(varsC[d],varsV)), + + [vmap_e,vmapSq_e,vmap_prime_e] : expandVmapFields(varsP), + + surfPerpVar : varsC[cdim], + VarsXY : delete(surfPerpVar, varsC), /* = [x, y]. */ + + bXY : basisFromVars("ser",VarsXY,polyOrder), + nodesXY : gaussOrdGkHyb(polyOrder+1, VarsXY, []), + numNodesXY : length(nodesXY), + + edge : ["lower", "upper"], + surfPerpVal : [-1, 1], + + for sI : 1 thru length(edge) do ( + printf(fh, "~%"), + /* Input kernel: constructs NN input features at the sheath surface, + calls gkyl_kann_net_apply, stores SRGRZ_N_MU outputs to *out. */ + printf(fh, "GKYL_CU_DH void bc_sheath_gyrokinetic_build_input_~a_~ax~av_~a_p~a(const double *phi, const double *phi_wall, const double *density, const double *temperature, const double *bmag, const double *bimpact_angle, int n_inp, float *nn_inp_out) ~%{ ~%", edge[sI], cdim, vdim, basisFun, polyOrder), + + phi_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi, bC)), + phiWall_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(phi_wall, bC)), + density_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(density, bC)), + temperature_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(temperature, bC)), + bmag_e : subst(surfPerpVar=surfPerpVal[sI],doExpand1(bmag, bC)), + bimpactAngle_e : doExpand1(bimpact_angle, bXY), /* bimpact_angle lives on (x,y) basis */ + + phi_q : evAtNodes(phi_e, nodesXY, VarsXY), + phiWall_q : evAtNodes(phiWall_e, nodesXY, VarsXY), + density_q : evAtNodes(density_e, nodesXY, VarsXY), + temperature_q : evAtNodes(temperature_e, nodesXY, VarsXY), + bmag_q : evAtNodes(bmag_e, nodesXY, VarsXY), + bimpactAngle_q : evAtNodes(bimpactAngle_e, nodesXY, VarsXY), + + printf(fh, " double phi_n, phi_wall_n, dens_n, temp_n, bmag_n, angle_n;~%"), + printf(fh, "~%"), + + for i : 1 thru numNodesXY do ( + printf(fh, " // node (x,y)_~a ~%", i-1), + printf(fh, " phi_n = ~a; ~%", gcfac(float(fullratsimp(phi_q[i])))), + printf(fh, " phi_wall_n = ~a; ~%", gcfac(float(fullratsimp(phiWall_q[i])))), + printf(fh, " dens_n = ~a; ~%", gcfac(float(fullratsimp(density_q[i])))), + printf(fh, " temp_n = ~a; ~%", gcfac(float(fullratsimp(temperature_q[i])))), + printf(fh, " bmag_n = ~a; ~%", gcfac(float(fullratsimp(bmag_q[i])))), + printf(fh, " angle_n = ~a; ~%", gcfac(float(fullratsimp(bimpactAngle_q[i])))), + printf(fh, " nn_inp_out[0+n_inp*~a] = (float)(angle_n*(180.0/GKYL_PI));~%", i-1), + printf(fh, " nn_inp_out[1+n_inp*~a] = (float)((1.0/bmag_n)*sqrt(GKYL_ELECTRON_MASS*fmax(0.0, dens_n)/GKYL_EPSILON0));~%", i-1), + printf(fh, " nn_inp_out[2+n_inp*~a] = (float)(fmax(0.0, (GKYL_ELEMENTARY_CHARGE*(phi_n-phi_wall_n))/temp_n));~%", i-1), + /* Note: In the future we can easily extend this to more parameters by adding more nn_inp_out entries, like: + nn_inp_out[3+n_inp*~a] = (float)(tperp/tpar); // temperature anisotropy + */ + printf(fh, "~%") + ), + + printf(fh, "}~%") + ) +)$ diff --git a/maxima/g0/twist_shift_calc/twistShift-calc.mac b/maxima/g0/twist_shift_calc/twistShift-calc.mac index 097320c4..868f546d 100644 --- a/maxima/g0/twist_shift_calc/twistShift-calc.mac +++ b/maxima/g0/twist_shift_calc/twistShift-calc.mac @@ -41,7 +41,7 @@ for bInd : 1 thru length(bName) do ( print("pOrder = ",pOrder), vStr : "", if (v>0) then (vStr: sconcat(v,"v")), - fname : sconcat("/home/akash/max-out/bc_twistshift_gyrokinetic_", bName[bInd], "_", c, "x", vStr, "_p", pOrder, ".c"), + fname : sconcat("~/max-out/bc_twistshift_gyrokinetic_", bName[bInd], "_", c, "x", vStr, "_p", pOrder, ".c"), fh : openw(fname), disp(printf(false,sconcat("Creating ~ax", vStr, "P~a ", bName[bInd]),c,pOrder)),