[DR] Excision BC - Githubissues

The excision BC is presently implemented in the gk-g0-app_excision branch (created before the DR process was instituted), simply as an updater in /zero (bc_excision), which gets called by the GK app just like bc_basic or bc_sheath.

For illustration purposes consider a domain such as

        ghost           skin
  |...............|...............|
  |     g3        |     s3        |
  |...............|...............|
  |     g2        |     s2        |
  |...............|...............|
  |     g1        |     s1        |
  |...............|...............|
  |     g0        |     s0        |
  |...............|...............|
x=x_lo-dx       x_lo          x_lo+dx

y^
  x->

The excision BC works by

Assume we have a domain with a whole in the middle, i.e. the lower boundary in one direction wraps onto itself, and the direction in which it wraps is periodic. In the sample diagram the hole is below x_lo in the x direction, and y is periodic.
Take J*f from a skin cell along this boundary, e.g. cell s2 above.
Copy it to the ghost cell on the opposite side of the hole, e.g. cell g1 above.
Flip the sign of the DG coefficients corresponding to odd-ordered monomials in the direction perpendicular to the boundary.

In order to avoid complicated interpolations or flux-matching between overlapping cells, this is to be limited to cases in which the number of cells along this boundary is even (e.g. 4 in the above example). That way a skin cell is matched to a single "ghost cell on the opposite side of the hole".

There's a question of whether one should copy Jf or f, and whether we should copy the volume expansion, or some surface expansion. I'll address these two: A) For the present use case, which is with cylindrical/circular holes (e.g. LAPD, mirrors, black holes), the Jacobian is constant along the boundary (in the direction in which it wraps), i.e. along the x=x_lo line above. Therefore, Jf in s2 evaluated at x_lo will be equal to the result of applying this BC, J*f in g1, at the x_lo surface. B) For the case of J being constant along the x_lo boundary, we can do these operations with volume DG coefficients or surface DG coefficients. That is, in Maxima, these three options are all equivalent (assuming p=1 in x):

fGhost_c : calcInnerProdList(varsP,1,basis,subst(x=-1,fSkin_e));

fGhostSurf_c : calcInnerProdList(surfIntVars,1,bSurf,subst(x=-1,fSkin_e));
fGhost_c : calcInnerProdList(vars,1,basis,doExpand(fGhostSurf_c,bSurf));

fGhost_c : makelist(fin[i-1],i,1,numB)$
for k : 1 thru numB do (
    if not freeof(x,basis[k]) then (
      fGhost_c[k] : -fGhost_c[k]
    )
)$

because ultimately what matters, and all these three guarantee it, is that J, f, and J*f at the s2 x_lo boundary is the same as that at the g1 x_lo boundary.

Two other relevant observatons: i) As mentioned, this is constrained to the case in which J is constant along the boundary line (e.g. x_lo above). It's also OK if J changes along z as long as the change in z is the same in s2 as it is in s1 (as it would be in LAPD and mirrors). However for tokamaks this would not be the case. But although there are lessons here relevant to tokamaks, I don't think these BCs can be used for them because the banana orbits are more complicated and dependent on things other than the shape of the hole. We'll have to create different BCs to handle that case. ii) One may think that what we want is flux out of s2 goes into s1. However, this isn't easy to do because:

The flux depends on the ExB drift, which can be of opposite signs on the s2 and s1 x_lo surfaces, i.e. you could have a instance of the flux being outward at both of these surfaces, in which case it doesn't make sense to try to take flux from s2 and put it into s1.
Our code (GK/Vlasov kernels and hyper_dg) are not set up to take fluxes, they take f in the neighboring (ghost) cell.

ammarhakim / gkylzero

[DR] Excision BC #419