Regarding projections

[abarret]

After regridding, the interpolated velocities will in general no longer satisfy the co-incompressibiloty condition. While this does not appear to affect the order of accuracy, we should do something about this.

One possibility is to project the volume averaged velocity onto an incompressible field. However, it's unclear what the correct phase velocities should be in this case.

An alternative is to use divergence preserving interpolation schemes. This would minimize cost and hopefully be just as effective. See https://doi.org/10.1016/j.jcp.2022.111500

[bindi-nagda]

Python code as a first attempt to reproduce the results from the paper: https://github.com/bindi-nagda/div-free-interpolation

[abarret]

Python code as a first attempt to reproduce the results from the paper: https://github.com/bindi-nagda/div-free-interpolation

The questions to answer are:

1. Can you reproduce the results from the paper?
2. For velocity fields that either are divergence free or not:
   • If you interpolate to a grid that's twice as fine, does that preserve divergence better than linear interpolation?
   • If you interpolate to a grid that's half as coarse, does that preserve divergence better than averaging?

[bindi-nagda]

Can you reproduce the results from the paper?

Yes I can for schemes C^0 and C^0i if I follow their instructions.

But I am confused with their instructions. In the paper in Section 3.1, they state:

We perform our tests on a 16^D grid over the region [0, 1]^D (for D = 2, 3). We seed the domain with one million uniformly random test locations, which are held fixed across all tests, and reject and resample locations that are within a distance of δ from a cell face.

Following these instructions exactly yields the promised divergence-free results. However, it is not clear to me the reason behind: a) why do they use δ in their equation to approximate divergence instead of dx, dy? b) why do they reject and resample locations that are within a distance of δ from a cell face.

To investigate b) further, if I pick test locations that fall on the input (coarser) grid, the divergence is computed to be non-zero. I suspect this is why they reject those locations, although I don't know the mathematical reason behind this. Also, I think this poses a problem...if we want to interpolate to, say a 128 x 128 grid, from a coarser grid of size 64 x 64 (or 32 x 32, or 16 x 16), then we will have grid points on our refined 128 x 128 grid that will fall on the coarser 64 x 64 grid. And then using their divergence formula, we will get non-zero divergence at those grid points.

[abarret]

a) why do they use δ in their equation to approximate divergence instead of dx, dy?

They construct a global interpolant $\hat{u}^n$, which they can evaluate at any spatial location. To compute a discrete divergence, they require some notion of a "spatial distance." They use $\delta$ to show that their results are independent of grid spacing.

b) why do they reject and resample locations that are within a distance of δ from a cell face.

If you are within $\delta$ of a grid face, then using their discrete divergence formula, the interpolant would be evaluated in different cells. I think a key result here is "$\hat{u}^n$ is $C^{n-1}$ continuous," so if you use $n=1$, the interpolant does not have continuous derivatives across the boundary, which sounds bad for divergence interpolation. What happens if you use the $C^1$ interpolant?

[bindi-nagda]

If you are within $\delta$ of a grid face, then using their discrete divergence formula, the interpolant would be evaluated in different cells. I think a key result here is "$\hat{u}^n$ is continuous," so if you use , the interpolant does not have continuous derivatives across the boundary, which sounds bad for divergence interpolation.

Okay, this makes sense. So if we want to use their scheme, we simply map the coarse grid values directly to the fine grid at those points falling on a cell face shared by both coarse and fine grids?

What happens if you use the $C^1$ interpolant?

When I use the $C^1$ interpolant (coming from Sec 2.2), the divergence is also non-zero for points that are within $\delta$ of a grid face.

[abarret]

When I use the C1 interpolant (coming from Sec 2.2), the divergence is also non-zero for points that are within δ of a grid face.

That's a bit disappointing but maybe salvageable. How large is the divergence? If you refine the original grid, does the divergence appear to converge? What happens if you "fudge" the piecewise definition of $\hat{u}$ so that when you evaluate the divergence, you always use the same "piece" of the piecewise interpolant?

[bindi-nagda]

FYI - this is the surface plot of the x-velocity field interpolated from an 4 x 4 grid to a 128 x 128 using $\hat{u}^2$. Those spikes in the velocity are appearing at points on the shared cell faces.

How large is the divergence? If you refine the original grid, does the divergence appear to converge?

The magnitude of the divergence across cell faces depends on the input grid size. Here are some snippets of the divergence values if I interpolate onto a 128 x 128 grid. As you can see, anywhere the points fall on the shared cell edge, the divergence value is "bad."

With a 16 x 16 input grid: Divergence non-zero at (0.1640625, 0.25). Divergence = 7.106181388910615 Divergence-free at cell (0.1640625, 0.2578125). Divergence = 5.342943865116467e-10 Divergence-free at cell (0.1640625, 0.265625). Divergence = 2.914344321425233e-10 Divergence-free at cell (0.1640625, 0.2734375). Divergence = 1.7173817923321621e-10 Divergence-free at cell (0.1640625, 0.28125). Divergence = 2.879634308783352e-10 Divergence-free at cell (0.1640625, 0.2890625). Divergence = 4.198028591417824e-10 Divergence-free at cell (0.1640625, 0.296875). Divergence = 6.661338147750939e-10 Divergence-free at cell (0.1640625, 0.3046875). Divergence = 1.0321610233177125e-09 Divergence non-zero at (0.1640625, 0.3125). Divergence = 2.2146394310189335

With a 64 x 64 input grid: Divergence-free at cell (0.1640625, 0.25). Divergence = -22.98608008567682 Divergence-free at cell (0.1640625, 0.2578125). Divergence = -7.216449660063518e-10 Divergence-free at cell (0.1640625, 0.265625). Divergence = -29.1199890298105 Divergence-free at cell (0.1640625, 0.2734375). Divergence = -7.771561172376096e-10 Divergence-free at cell (0.1640625, 0.28125). Divergence = -28.070176252989466 Divergence-free at cell (0.1640625, 0.2890625). Divergence = -7.07768066376957e-10 Divergence-free at cell (0.1640625, 0.296875). Divergence = -20.095624093705027 Divergence-free at cell (0.1640625, 0.3046875). Divergence = -2.532694054480089e-10 Divergence-free at cell (0.1640625, 0.3125). Divergence = -7.163605393292055

[abarret]

That plot does not look not continuous (which even the C^0 reconstruction should be). You should double check that there's not a bug in your code.

[bindi-nagda]

Okay, I think I have fixed the code. The image below shows the interpolated function on a 64 x 64 grid, where the original grid was 32 x 32. It's hard to tell from the image if the surface is C1 continuous. The paper says that the C1 interpolant is not an interpolating polynomial (property 9 in Sec 2.2). Computing the divergence shows the interpolation is divergence-free.

When I calculate the divergence using the test methodology outlined in the paper, I am able to reproduce their results in Sec 3.1 (to the same order of magnitude) for the velocity fields $u{2a}, u{2b}, u{2b}, u{2d}$.

The image below shows the same function interpolated using bilinear interpolation from a 32 x 32 grid to a 64 x 64 grid. This looks smoother than the first plot. But computing the divergence shows the interpolation is not divergence-free.

[abarret]

The image below shows the same function interpolated using bilinear interpolation from a 32 x 32 grid to a 64 x 64 grid. This looks smoother than the first plot. But computing the divergence shows the interpolation is not divergence-free.

What is not divergence free: the bilinear interpolant or the C^1 approximation? How do the divergences compare?

[bindi-nagda]

The bilinear interpolation is not divergence-free. The divergence values are on the order of 1e-1 to 1e1.

The C^1 approximation is divergence-free, with the divergence values falling between 1e-12 to 1e-9.

P.S: The plots show the x-component of the vector field: $[\sin(6x+2)\sin(6y+4), \cos(6x+2)\cos(6y+4)]$

[abarret]

The bilinear interpolation is not divergence-free. The divergence values are on the order of 1e-1 to 1e1.

The C^1 approximation is divergence-free, with the divergence values falling between 1e-12 to 1e-9.

P.S: The plots show the x-component of the vector field: [sin⁡(6x+2)sin⁡(6y+4),cos⁡(6x+2)cos⁡(6y+4)]

Even on the boundary of a cell? What about a vector field that is not divergence-free?

[bindi-nagda]

Even on the boundary of a cell?

No, neither are divergence-free for test locations on the boundary of a cell. The divergence values for bilinear and C^1 are on the same order of magnitude.

What about a vector field that is not divergence-free?

For such a vector field, with test locations that aren't within $\delta$ of a cell edge, the C^1 approximation is not divergence-free. This vector field used to test this was $u_{2d}$ from the paper.

[bindi-nagda]

Btw, by boundary of a cell, I’m referring to the boundaries of the original coarse grid. Test locations on cell boundaries of the coarse grid give divergence values that are non-zero.