Possible correctness bug in PowerLawConvectionTerm

[opensourcefan314]

PowerLawConvectionTerm seems to handle constraints at boundary faces incorrectly, or at least differently than I would expect.

This gist contains a minimal example, solving one step of an equation with a Grid1D with just 2 cells:
https://gist.github.com/opensourcefan314/c7fa68bf44e069406f617439a9417518

I would expect it to print [2/3, 2/3] but it prints [1/2, 1/2].

I think I'm doing the math correctly because I've been able to do the math to reproduce FiPy exactly for diffusion terms and for slightly different convection problems. I think the problem is specifically with the boundary conditions because I always match FiPy exactly if I set convective strength to 0 at the outermost faces.

The attached PDF contains my derivation of the [2/3, 2/3] expected result.

minfail.pdf

[guyer]

There are no ghost cells in FiPy. In your derivation, $u_{0-,1} \equiv u_{0-} = 1$.
As a result

$$ \begin{aligned} u_{0,1} &= u_{0-} - u_{0,1} + u_{0, 0} \\ &= \frac{u_{0-} + u_{0,0}}{2} \\ u_{1,1} &= u_{0,1} + u_{1,0} \end{aligned} $$

For $u_{i, 0} = [0, 0]$, this gives $u_{i,1} = [u_{0-} / 2, u_{0-} / 2] = [1/2, 1/2]$.

You might find this notebook helpful to understand some potentially non-intuitive things about FiPy's boundaries.

[opensourcefan314]

Thanks. I agree this test case can be explained by using $u_{0-,1} \equiv u_{0-} = 1$ to calculate the face values.

However, that process doesn't seem to be consistently applied to other convection boundary conditions.

If we take the same code and change v_face = np.array([1.0, 1.0, 0.0]) to v_face = np.array([0.0, 1.0, 1.0]) to exercise the right/east boundary condition, we get behavior that is not consistent with $u_{1+,1} \equiv u_{1+} = 1$.

In that case, FiPy returns [0, 0].
I can get [0, 0] using alpha-weighted faces values and a ghost cell (and there are presumably other ways to derive it), but can’t derive [0, 0] using $u_{1+,1} \equiv u_{1+} = 1$.

I’ve attached my derivation for the v=[0,1,1] case: minfail2.pdf

Any idea why the boundary conditions seem to be handled differently for v=[1, 1, 0] than for v=[0, 1, 1]? Any idea what is the correct process to follow (instead of ghost cell + alpha weighting) to derive [0, 0] in the v=[0, 1, 1] case? Thanks in advance.

I've looked at the linked notebook and didn't find a connection to this issue. The linked notebook is all about diffusion, and I've been able to successfully predict the behavior of diffusion term boundary conditions. I'm pretty sure the issue I'm describing here is specific to convection, and I'm guessing related to the way that convection uses alpha weighting when calculating face values.

[guyer]

What I said earlier was incorrect. FiPy does not use ghost cells, but the boundary constraints for a ConvectionTerm are handled semi-implicitly:

$$ \begin{aligned} u_{i+,j} &= \alpha_{i+} u_{i,j} + \left(1 - \alpha_{i+}\right) u_{i+} \\ u_{i-,j} &= \alpha_{i-} u_{i,j} + \left(1 - \alpha_{i-}\right) u_{i-} \end{aligned} $$

so there is an $\alpha$ contribution to the matrix diagonal and $(1 - \alpha) BC$ contribution to the RHS-vector.

$v = [0, 1, 1]$ and $v = [1, 1, 0]$ are different because, in the first case, there is no convection of the left boundary condition into the domain because $v_{0-} = 0$ and the convection within the domain is away from the boundary. In the second case, there is likewise no convection of the right boundary condition into the domain because $v_{1+} = 0$, but in this case the convection within the domain is toward the boundary, so some $\alpha$ balancing occurs.

[opensourcefan314]

Thanks! I understand it now. I confirm the semi-implicit approach explains both of my test cases.

[guyer]

Great. I'm going to convert this to a "Discussion".