Get em back

maxima/g0/gk_collisionless/dg_gk-surf.mac

@@ -26,7 +26,7 @@ calcGKSurfUpdateInDir(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, vars
    BstarXdBmag_e,BstarYdBmag_e,BstarZdBmag_e,BstardBmag_e,
    hamil_e,alphaSurfL_e,alphaSurfR_e,
    fl_e,fc_e,fr_e,fUpL_e,fUpR_e,GhatL_c,GhatR_c,GhatL_e,GhatR_e,incrL_c,incrR_c,pOrderCFL,
-   fnodal_l_e, fnodal_r_e, fmodproj_e],
+   fnodal_l_e,fnodal_r_e,fmodproj_e,numC,vmap_prime_fac_l,vmap_prime_fac_c,vmap_prime_fac_r,vmap_prime_e,basisNodal,i],
 
   kill(varsC,varsP,bC,bP),
   pDim : cdim+vdim,
@@ -143,7 +143,8 @@ calcGKBoundarySurfUpdateInDir(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrd
    rdx2vec,rdv2vec,rdSurfVar2,bmagBasis,ignoreVars,inFlds_e,cmag_e,b_x_e,b_y_e,b_z_e,jacobTotInv_e,
    BstarXdBmag_e,BstarYdBmag_e,BstarZdBmag_e,BstardBmag_e,
    hamil_e,alphaUpL_e,alphaSurfL_e,alphaUpSurfL_e,alphaUpR_e,alphaSurfR_e,alphaUpSurfR_e,
-   fEdge_e,fSkin_e,fUpL_e,fUpR_e,GhatL_c,GhatR_c,GhatL_e,GhatR_e,incrL_c,incrR_c,pOrderCFL],
+   fEdge_e,fSkin_e,fUpL_e,fUpR_e,GhatL_c,GhatR_c,GhatL_e,GhatR_e,incrL_c,incrR_c,pOrderCFL,
+   numC,vmap_prime_fac_edge,vmap_prime_fac_skin,vmap_prime_e,basisNodal,i,fnodal_l_e,fnodal_r_e,fmodproj_e],
 
   kill(varsC,varsP,bC,bP),
   pDim : cdim+vdim,

maxima/g0/gk_collisionless/gk_collisionless_flux-surf-conf.mac

@@ -6,13 +6,15 @@ load("utilities_gyrokinetic")$
 load("nodal_operations/nodal_functions")$
 fpprec : 24$
 
-buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no_by, edge, mb_bound) := block(
+buildGKFluxConfKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no_by, edge, mb_bound, scheme) := block(
   [pDim,varsC,bC,varsP,bP,vSub,numC,numP,surfVar,varLabel,dirLabel,surfIntVars,surf_cvars,surf_vvars,
    surfNodes,nodeVars,bSurf,basisNodal,surfConfigNodes,numSurfNodes,numSurfConfigNodes,numVelNodes,
-   numMuNodes,numVparNodes,d,rdx2vec,rdv2vec,rdSurfVar2,bmagBasis,phi_e,bmagSurf_e,vmap_e,vmapSq_e,
+   numMuNodes,numVparNodes,d,rdx2vec,rdv2vec,rdSurfVar2,bmagBasis,phi_e,apar_e,
+   bmagSurf_e,vmap_e,vmapSq_e,
    vmap_prime_e,evPoint,hamil_e,hamil_c,replaceList,hamilNoZero_c,JfL_e,JfR_e,JfL_c,JfR_c,
    jacobgeo_rat_surfR_e,jacobgeo_rat_surfL_e,JfL_nodes,JfR_nodes,vmap_prime_nodes,vpardim,
-   dH_dz_nodes,mvpar_nodes,di3,i,j,j0index,j1index,vparindex,vpar0index,pOrderCFL
+   dH_dz_nodes,mvpar_nodes,i,j,j0index,j1index,vparindex,vpar0index,pOrderCFL,
+   surfIntVarsC,bSurfC,hamilCvar,aparCvar,apar_nodes,apar_c,dA_dx_nodes
    ],
 
   kill(varsC,varsP,bC,bP),
@@ -60,7 +62,8 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
     const double *w, const double *dxv,
     const double *vmap, const double *vmapSq, const double q_, const double m_,
     const struct gkyl_dg_surf_geom *dgs, const struct gkyl_gk_dg_surf_geom *gkdgs, 
-    const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR, const double *phi,
+    const double *bmag, const double *jacobgeo_rat_surfL, const double *jacobgeo_rat_surfR, 
+    const double *phi, const double *apar, const double *apardot,
     const double *JfL, const double *JfR, double* GKYL_RESTRICT flux_surf) ~%{ ~%", funcNm),
   printf(fh, "  // w[NDIM]: cell-center.~%"),
   printf(fh, "  // dxv[NDIM]: cell length.~%"),
@@ -73,6 +76,8 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   printf(fh, "  // jacobgeo_rat_surfL: Ratio of surface conf-space Jacobians in left cell.~%"),
   printf(fh, "  // jacobgeo_rat_surfR: Ratio of surface conf-space Jacobians in right cell.~%"),
   printf(fh, "  // phi: electrostatic potential.~%"),
+  printf(fh, "  // apar: parallel component of vector potential.~%"),
+  printf(fh, "  // apardot: time derivative of parallel component of vector potential.~%"),
   printf(fh, "  // JfL: distribution times total jacobian in left cell.~%"),
   printf(fh, "  // JfR: distribution times total jacobian in right cell.~%"),
   printf(fh, "  // flux_surf: output surface phase space flux in each direction (cdim + 1 components).~%"),
@@ -120,7 +125,7 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   flush_output(fh),
   hamilNoZero_c : makelistNoZeros1(hamil_c, hamilCvar),
   /* Expand projected Hamiltonian on basis. */
-  hamil_e : hamilNoZero_c . bSurf,
+  hamil_e : doExpand(hamilNoZero_c, bSurf),
 
   /* fl and fr */
   JfL_e : doExpand1(JfL, bP),
@@ -173,16 +178,28 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
     )
   ),
 
-  mvpar_nodes : [],
+  dHdvpar_nodes : [],
   for i : 1 thru numVparNodes do (
-    mvpar_nodes : append(mvpar_nodes, [dH_dz_nodes[vpardim][i]])
+    dHdvpar_nodes : append(dHdvpar_nodes, [dH_dz_nodes[vpardim][i]])
   ),
 
-  if surfDir = cdim then(
-    di3 : true
-  )
-   else (
-    di3 : false
+  /* Expand Aparallel. */
+  apar_e : doExpand1(apar,bC),
+  apar_c : calcInnerProdList(surfIntVars, 1, bSurf, subst(surfVar=evPoint,apar_e)),
+  printf(fh, "  double apar_surf[~a] = {0.}; ~%", length(bSurf)),
+  aparSurfCvar : eval_string(sconcat("apar_surf")),
+  writeCExprsNoExpand1(aparSurfCvar, gcfac(float(expand(apar_c)))),
+  printf(fh, "~%"),
+  flush_output(fh),
+  apar_c : makelistNoZeros1(apar_c, aparSurfCvar),
+  /* Expand projected Apar on basis. */
+  apar_e : doExpand(apar_c, bSurf),
+  /* Eval Aparallel at nodes. */
+  apar_nodes : float(evAtNodes(apar_e,surfNodes,surfIntVars)),
+  /* Compute gradient of Aparallel */
+  dA_dx_nodes : makelist(0, i, 1, cdim),
+  for i : 1 thru cdim do (
+    dA_dx_nodes[i] : float(evAtNodes(diff(apar_e*rdx2vec[i],varsP[i]),surfNodes,surfIntVars))
   ),
 
   /* Now calculate flux at all quadrature nodes */
@@ -202,14 +219,14 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   printf(fh, "  double Jfavg_quad = 0.0; ~%"),
   printf(fh, "  double Jfjump_quad = 0.0; ~%"),
 
-  printf(fh, "  double mvpar_quad[3] = {0.0}; ~%"),
+  printf(fh, "  double dHdvpar_quad[3] = {0.0}; ~%"),
   for i : 1 thru numVparNodes do (
-    printf(fh, "  mvpar_quad[~a] = ~a; ~%", i-1, mvpar_nodes[i])
+    printf(fh, "  dHdvpar_quad[~a] = ~a; ~%", i-1, dHdvpar_nodes[i])
   ),
 
   printf(fh, "  double mvparsq_quad[3] = {0.0}; ~%"),
   for i : 1 thru numVparNodes do (
-    printf(fh, "  mvparsq_quad[~a] = mvpar_quad[~a]*mvpar_quad[~a]/m_; ~%", i-1, i-1,i-1)
+    printf(fh, "  mvparsq_quad[~a] = dHdvpar_quad[~a]*dHdvpar_quad[~a]/m_; ~%", i-1, i-1,i-1)
   ),
   printf(fh, "~%"),
 
@@ -229,64 +246,91 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
       vparindex : mod(j-1, numVparNodes) + 1,
       vpar0index : mod(j-1, numVparNodes),
       printf(fh, "~%"),
-      if no_by = true then (
-        if di3 = true then (
-          printf(fh, "  alpha_quad = (mvpar_quad[~a]*B3_quad/(m_*bmag_quad))*area_elem_quad/Jc_quad; ~%", vpar0index)
-        )
-        else (
-          printf(fh, "  alpha_quad = 0.0; ~%")
-        )
+
+      /* Electrostatic term */
+      if surfDir = cdim then (
+        /* Parallel streaming term */
+        printf(fh, "  alpha_quad = (dHdvpar_quad[~a]*B3_quad/(m_*bmag_quad)); ~%", vpar0index)
+      ) else (
+        printf(fh, "  alpha_quad = 0.0; ~%")
       ),
+
       if no_by = false then (
-        /*printf(fh, "  alpha_quad += mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) ;~%", vpar0index),*/
+        printf(fh, "  alpha_quad += mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_); ~%", vpar0index),
+        /* + 1/q b/B x ∇H . ∇ψ*/
         if cdim = 3 then (
            if surfDir = 1 then(
-             printf(fh, "  alpha_quad = (mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[1]*(~a) - bhat_quad[2]*(~a)))*area_elem_quad/Jc_quad; ~%", vpar0index, dH_dz_nodes[3][j1index], dH_dz_nodes[2][j1index])
+             printf(fh, "  alpha_quad += 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[1]*(~a) - bhat_quad[2]*(~a)); ~%", dH_dz_nodes[3][j1index], dH_dz_nodes[2][j1index])
            ),
            if surfDir = 2 then(
-             printf(fh, "  alpha_quad = (mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[2]*(~a) - bhat_quad[0]*(~a)))*area_elem_quad/Jc_quad; ~%", vpar0index, dH_dz_nodes[1][j1index], dH_dz_nodes[3][j1index])
+             printf(fh, "  alpha_quad += 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[2]*(~a) - bhat_quad[0]*(~a)); ~%", dH_dz_nodes[1][j1index], dH_dz_nodes[3][j1index])
            ),
            if surfDir = 3 then(
-             printf(fh, "  alpha_quad = (mvpar_quad[~a]*B3_quad/(m_*bmag_quad) + mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[0]*(~a) - bhat_quad[1]*(~a)))*area_elem_quad/Jc_quad; ~%", vpar0index, vpar0index, dH_dz_nodes[2][j1index], dH_dz_nodes[1][j1index])
+             printf(fh, "  alpha_quad += 1/(q_*bmag_quad*area_elem_quad) * (bhat_quad[0]*(~a) - bhat_quad[1]*(~a)); ~%", dH_dz_nodes[2][j1index], dH_dz_nodes[1][j1index])
            )
         ),
         if cdim = 2 then (
            if surfDir = 1 then(
-             printf(fh, "  alpha_quad = (mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * bhat_quad[1]*(~a))*area_elem_quad/Jc_quad; ~%", vpar0index, dH_dz_nodes[2][j1index])
+             printf(fh, "  alpha_quad += 1/(q_*bmag_quad*area_elem_quad) * bhat_quad[1]*(~a); ~%", dH_dz_nodes[2][j1index])
            ),
            if surfDir = 2 then(
-             printf(fh, "  alpha_quad = (mvpar_quad[~a]*B3_quad/(m_*bmag_quad) + mvparsq_quad[~a]*normcurlbhat_quad/(bmag_quad*q_) + 1/(q_*bmag_quad*area_elem_quad) * -bhat_quad[1]*(~a))*area_elem_quad/Jc_quad;~%", vpar0index, vpar0index, dH_dz_nodes[1][j1index])
+             printf(fh, "  alpha_quad -= 1/(q_*bmag_quad*area_elem_quad) * bhat_quad[1]*(~a); ~%", dH_dz_nodes[1][j1index])
            )
         ),
-        if cdim = 1 then (
-          printf(fh, "  alpha_quad = (mvpar_quad[~a]*B3_quad/(m_*bmag_quad))*area_elem_quad/Jc_quad; ~%", vpar0index)
+
+        /* Electromagnetic Apar contribution using ∇ x (A b) = A ∇ x b + ∇A x b */
+        /* No contribution for cdim = 1 */
+        if cdim > 1 then (
+          /* + 1/m A (∇ x b)/B . ∇ψ dH/dvpar */
+          printf(fh, "  alpha_quad += 1/m_ * (~a) * normcurlbhat_quad/bmag_quad * dHdvpar_quad[~a]; ~%", apar_nodes[j1index], vpar0index)
+        ),
+        /* + 1/m ∇A x b/B . ∇ψ dH/dvpar */
+        if cdim = 3 then (
+          if surfDir = 1 then(
+            printf(fh, "  alpha_quad += 1/(m_*bmag_quad*area_elem_quad) * ((~a) * bhat_quad[2] - (~a) * bhat_quad[1]) * dHdvpar_quad[~a]; ~%", dA_dx_nodes[2][j1index], dA_dx_nodes[3][j1index], vpar0index)
+          ),
+          if surfDir = 2 then(
+            printf(fh, "  alpha_quad += 1/(m_*bmag_quad*area_elem_quad) * ((~a) * bhat_quad[0] - (~a) * bhat_quad[2]) * dHdvpar_quad[~a]; ~%", dA_dx_nodes[3][j1index], dA_dx_nodes[1][j1index], vpar0index)
+          ),
+          if surfDir = 3 then(
+            printf(fh, "  alpha_quad += 1/(m_*bmag_quad*area_elem_quad) * ((~a) * bhat_quad[1] - (~a) * bhat_quad[0]) * dHdvpar_quad[~a]; ~%", dA_dx_nodes[1][j1index], dA_dx_nodes[2][j1index], vpar0index)
+          )
+        ),
+        /* in 2D, we take the first and last component of the 3D cross product, setting d/dy = 0.*/
+        if cdim = 2 then (
+          if surfDir = 1 then(
+            printf(fh, "  alpha_quad -= 1/(m_*bmag_quad*area_elem_quad) * (~a) * bhat_quad[1] * dHdvpar_quad[~a]; ~%", dA_dx_nodes[2][j1index], vpar0index)
+          ),
+          if surfDir = 2 then(
+            printf(fh, "  alpha_quad += 1/(m_*bmag_quad*area_elem_quad) * (~a) * bhat_quad[1] * dHdvpar_quad[~a]; ~%", dA_dx_nodes[1][j1index], vpar0index)
+          )
         )
+        /* Nothing in 1D */
       ),
 
+      /* Multiply by |e^i|, note: area_elem_quad = J_c |e^i| */
+      printf(fh, "  alpha_quad = alpha_quad * area_elem_quad/Jc_quad; ~%"),
+
       printf(fh, "~%"),
-      /*printf(fh, "  alpha_quad = alpha_quad*area_elem_quad/Jc_quad; ~%"),*/
       printf(fh, "  cfl = fmax(fabs(alpha_quad), fabs(cfl)); ~%"),
       printf(fh, "  JfL_quad = ~a; ~%", JfL_nodes[j1index]),
       printf(fh, "  JfR_quad = ~a; ~%", JfR_nodes[j1index]),
       printf(fh, "  Jfavg_quad = (JfL_quad + JfR_quad)/2.0; ~%"),
       printf(fh, "  Jfjump_quad = (JfR_quad - JfL_quad)/2.0; ~%"),
-      printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad - fabs(alpha_quad)*Jfjump_quad; ~%", j0index)
+      if scheme = "upwind" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad - fabs(alpha_quad)*Jfjump_quad; ~%", j0index)
+      ) else if scheme = "central" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad; ~%", j0index)
+      )       else if scheme = "downwind" then (
+        printf(fh, "  flux_surf_nodal[~a] = alpha_quad*Jfavg_quad + fabs(alpha_quad)*Jfjump_quad; ~%", j0index)
+      ) else (
+        /* Stop the script if an invalid scheme is provided. */
+        error("Invalid flux scheme provided. Options are: upwind, central, downwind.")
+      )
     ),
     printf(fh, "~%")
   ),
 
-  /* Do the quad nodal to modal ops directly here*/ 
-  /*printf(fh, "~%"),
-  printf(fh, "  double *fmodal = &flux_surf[~a]; ~%", length(bSurf)*(surfDir-1)),
-  flux_surf_nodal_e : doExpand1(flux_surf_nodal,basisNodal),
-  fmodproj_e : fullratsimp(calcInnerProdList(surfIntVars, 1, bSurf, flux_surf_nodal_e)),
-
-  for i : 1 thru length(fmodproj_e) do (
-    printf(fh, "  fmodal[~a] = ~a; ~%", i-1, float(expand(fmodproj_e[i])))
-  ),
-
-  printf(fh, "~%"),*/
-
   /*Calculate the cfl*/
   pOrderCFL : polyOrder,
   printf(fh, "~%"),
@@ -296,4 +340,4 @@ buildGKFluxConfESKernel(surfDir, fh, funcNm, cdim, vdim, basisFun, polyOrder, no
   flush_output(fh),
   printf(fh, "} ~%")
 
-)$
+)$