Cut-cell algorithm used in Linear Solvers

[ankithadas]

Dear AMReX team,

What is the algorithm used for calculating cut-cell gradients in MLEBABecLap solver? Any paper/references detailing the method would be much appreciated.

Kind regards,
Ankith

[ankithadas]

It looks likes the gradients at the cut faces are calculated using Johansen-Colella framework (https://doi.org/10.1006/jcph.1998.5965). However, it is not clear how the gradients at the cut boundaries are calculated.

I have attached the code snippet that calculates the boundary EB gradients (extracted from EB Poisson Solver)

if (is_dirichlet) {
    Real dapx = (apxm-apxp)/dxinv[1];
    Real dapy = (apym-apyp)/dxinv[0];
    Real anorm = std::hypot(dapx,dapy);
    Real anorminv = Real(1.0)/anorm;
    Real anrmx = dapx * anorminv;
    Real anrmy = dapy * anorminv;
    Real phib = is_inhomog ? phieb(i,j,k,n) : Real(0.0);
    Real bctx = bc(i,j,k,0);
    Real bcty = bc(i,j,k,1);
    Real dx_eb = get_dx_eb(kappa);
    Real dg, gx, gy, sx, sy;
    if (std::abs(anrmx) > std::abs(anrmy)) {
        dg = dx_eb / std::abs(anrmx);
    } else {
        dg = dx_eb / std::abs(anrmy);
    }
    gx = (bctx - dg*anrmx);
    gy = (bcty - dg*anrmy);
    sx = std::copysign(Real(1.0),anrmx);
    sy = std::copysign(Real(1.0),anrmy);
    int ii = i - static_cast<int>(sx);
    int jj = j - static_cast<int>(sy);
    Real phig = (Real(1.0) + gx*sx + gy*sy + gx*gy*sx*sy) * x(i ,j ,k,n)
        +       (    - gx*sx         - gx*gy*sx*sy) * x(ii,j ,k,n)
        +       (            - gy*sy - gx*gy*sx*sy) * x(i ,jj,k,n)
        +       (                    + gx*gy*sx*sy) * x(ii,jj,k,n) ;
    Real dphidn = (phib-phig) / dg;
    feb = dphidn * ba(i,j,k) * beb(i,j,k,n);
}

Any suggestions why get_dx_eb() is used or the mathematics behind it would be greatly appreciated.

[WeiqunZhang]

If it's 1D, it would be very easy to compute the gradient on the EB. The formula used in get_dx_eb is what we would use if the EB face is perpendicular to one of the axis. Admittedly, this is not the case in general. But it simplifies the stencil a lot and in practice it appears to work well. The factor of 0.3 in get_dx_eb came from testing.

[ankithadas]

Thank you for your response. Is the analysis explained in any paper that I could look into ? Could you please explain how phig is calculated.

[WeiqunZhang]

We don't have a paper on this. phig is computed with linear interpolation (extrapolation). Suppose this is 1D. Suppose there are two cells and the cell on the left is a cut cell. Suppose the EB is at x=(1-kappa) where kappa is volume fraction, the two cell centers are 0.5 and 1.5. Given the three points, one can draw a parabola and compute the gradient at x=1-kappa. This is equivalent to using the two cell centers to interpolate or extrapolation to x=get_dx_eb(kappa) x=(1-kappa)+get_dx_eb(kappa) (I guess? Anyway you should be able to get the idea and ignore whatever mistakes I might have made), and then compute the gradient with (phig-phieb)/dx_eb.

[ankithadas]

Thank you for your response. Your response does clear a lot of things up. Could you please be able to give me some insight into the technique for 2D ? Any insight about the math performed would be much appreciated