Gk sheath bc with mu dependent vcut

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 <gkyl_bc_sheath_gyrokinetic_kernels.h> ~%"),
@@ -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 <gkyl_util.h>~%")$
@@ -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)
 
         )
       )

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. */

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,18 +77,17 @@ 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. */
     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. */
+    /* 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,54 +111,46 @@ 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],
 
       /* Project f at this (mu)_j node onto vpar basis (bVpar). */
       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),