Provide higher order Interpolation

[isaac-blanc]

Update

For those new to this thread, a little context... This started out with me asking about 1st order interpolation. I thought that the FiPy behavior was wrong. I learned that FiPy is behaving as intended. Rather, I learned that what I had expected as 1st order interpolation is actually higher order interpolation. Hence, this turned from a question into an enhancement request.

Background

The CellVariable.__call__() method supports 0th and 1st order interpolation via its order argument. I have been doing some experiments and getting unexpected results, as demonstrated below.

MWE

All results are from the following script. I have been adjusting the two inputs cell_values and order.

import numpy as np
import matplotlib.pyplot as plt
from fipy import Grid1D, CellVariable

# Inputs
cell_values = [1., 2., 3., 4., 5.]
order = 1

# FiPy setup
mesh = Grid1D(dx=1, nx=5)
y = CellVariable(mesh)
y.setValue(cell_values)

# Get data
x_exact, = mesh.cellCenters
y_exact = y
x_inter = np.linspace(start=0, stop=5)
y_inter = []
for x in x_inter:
    y_inter.append(y((x,), order=order))

# Plot
plt.scatter(x_exact, y_exact, label='Cell centers')
plt.scatter(x_inter, y_inter, label='Interpolated', s=5)
plt.xlabel('x')
plt.ylabel('y')
plt.xticks([0, 1, 2, 3, 4, 5])
plt.grid(axis='x')
plt.legend()
plt.show()

Plots

First, lets trial order=0 (nearest neighbour) interpolation. The results are as expected.

Second, lets trial order=1 (linear) interpolation. The results are intuitively wrong for the cells at either end of the domain.

Finally, lets change the cell_values slightly. These results are wrong for even more of the cells.

[guyer]

Thank you for the nice illustration of the issue. This is actually working as designed.

• The 0th order interpolant is $\Phi(\vec{x}) \approx \phi_P$
• The 1st order interpolant is $\Phi(\vec{x}) \approx \phi_P + (\vec{x} - \vec{x}_P) \cdot \nabla_P\, \phi$

where $P$ is the cell closest to the point of interest $\vec{x}$, $\phi_P$ is the value of $\phi$ in cell $P$, $\vec{x}_P$ is the position of the cell center of cell $P$, and $\nabla_P \phi$ is the Gaussian gradient evaluated at the center of the cell $P$.

The issue you're seeing arises because $\nabla_P\, \phi \equiv {1/V_P} \sum_f \phi_f A_f \hat{n}_f$, where $V_P$ is the volume of cell $P$, $\sum_f$ is the sum over all faces of the cell $P$, $\phi_f$ is the face value on face $f$, $A_f$ is the area of face $f$, and $\hat{n}_f$ is the normal to face $f$.

For a smoothly varying function, this definition is fine, but with sharp changes and near boundaries, the results can be counterintuitive. Without constraints, the boundary gradients are zero, meaning that the value on a boundary face is equal to the value at the center of the nearest cell.

I've added the face values to your plot to hopefully illustrate why the interpolated line segments work out as they do.

This notebook from several years ago might be helpful.

Our own work has not called for interpolation of higher order, but if you want to make this a feature request, we'll consider it.

[Sorry for the delay in reply; we have been furloughed for the past six weeks due to the lapse in US Federal Government appropriations]

[isaac-blanc]

Sorry to hear that you were furloughed! I was aware of the US Government shutdown but had not appreciated it would affect you.

Thank you very much for the detailed reply. I need some time to think about this before responding. Then I will follow-up with more questions or mark this as answered.

[isaac-blanc]

Thanks again for your reply. Here is my response.

1st order interpolation

I understand your explanation that 1st order interpolation uses the gradient evaluated at the center of a single cell. I had wrongly assumed that it used the values at the centre of two adjacent cells.

Alternative interpolation

I would like to use a different interpolation scheme, but I am not sure what it would be called. Happy to research this myself, but asking you first in case there is an obvious answer. I want the following.

• Interpolated values are 'continuous' across the domain.
• Interpolated values match boundary conditions.
• Close to physical/analytical results for convection.
• Piecewise linear interpolation is good enough.

First, I'll apply this to my original example. This matches with all cell and face values, which I think means that its a good approximation of physical conditions, but I may have overlooked something.

Second, I'll apply this to a slightly modified version of your convection.exponential1D.mesh1D example. Code is attached here. To my naïve eyes, this looks like a better approximation than FiPy 1st order. Moreover, the results look more physically intuitive because they are continuous. Academically, I expect this doesn't matter, but in a commercial context there is value in producing 'intuitive' results. It saves me time explaining FVM to colleagues from a different background.

One problem I can see with this approach is that (in this case) it is biased towards underestimating the value. If I calculated the average of each interpolation scheme, this would be lower.

Face values

Without constraints, the boundary gradients are zero, meaning that the value on a boundary face is equal to the value at the center of the nearest cell.

As a final sanity check, am I right in thinking that you are referring to the boundaries of the whole domain? That is to say the left and right boundaries of the mesh; not the intermediate boundaries between cells.

[guyer]

Alternative interpolation

I would like to use a different interpolation scheme, but I am not sure what it would be called. Happy to research this myself, but asking you first in case there is an obvious answer. I want the following.

• Interpolated values are 'continuous' across the domain.
• Interpolated values match boundary conditions.
• Close to physical/analytical results for convection.
• Piecewise linear interpolation is good enough.

To my naïve eyes, this looks like a better approximation than FiPy 1st order. Moreover, the results look more physically intuitive because they are continuous. Academically, I expect this doesn't matter, but in a commercial context there is value in producing 'intuitive' results. It saves me time explaining FVM to colleagues from a different background.

I understand why this is preferable to you, but I have a few issues with it:

• Including face values in the set that you're interpolating is dubious. Face values aren't "real"; by and large, they're already interpolated from the cell values, and they're absolutely never solved for.
• What you propose works in 1D, but FiPy is not a 1D code. FiPy is a 3D code that can do things in 1D.
• it's absolutely a "better" approximation than FiPy 1st order, because it's not 1st order.

Better, I think, is to implement higher-order interpolation using something like a radial basis interpolator on the cell values with appropriate consideration of the boundary constraints. FiPy's 1st-order interpolation already accounts for Dirichlet constraints, but it does not deal with most other boundary conditions.

One problem I can see with this approach is that (in this case) it is biased towards underestimating the value. If I calculated the average of each interpolation scheme, this would be lower.

Not so much underestimating as biased toward the concave side. If the solution is convex-up, then this approach will overestimate.

Without constraints, the boundary gradients are zero, meaning that the value on a boundary face is equal to the value at the center of the nearest cell.

As a final sanity check, am I right in thinking that you are referring to the boundaries of the whole domain? That is to say the left and right boundaries of the mesh; not the intermediate boundaries between cells.

Correct

[isaac-blanc]

I may apply a simple interpolation of the kind I have described in the 1D part of my own model, but I recognise RBF interpolation would be better for FiPy, since it generalises to 3D and to unstructured grids. I can submit a feature request for this if you like, but I can't make a strong argument for it. I don't expect that better interpolation will make a meaningful difference to the overall accuracy of my model. At best, it will make some of my results look more physically intuitive. I would feel a bit guilty asking someone to implement RBF interpolation in all FiPy grids just so that a few of my results can look better.

[guyer]

A feature request would be welcome, and it's useful to know that it's not urgent for you.

If you like, you don't need to make a separate feature request; you can edit this issue and add an enhancement label.

[isaac-blanc]

Do you have a preference between those two options?

[guyer]

Do you have a preference between those two options?

I'd say convert this issue to an enhancement request

[isaac-blanc]

I'm afraid I don't have permission to edit labels in this repo so I can't do that.

[guyer]

OK. I've added the tag and edited the title. Feel free to either edit the original post or add comments to reflect any specifics you'd like to see