StructuredColumns FixupHaloForVectors method interferes with SphericalVector interpolation

[odlomax]

What happened?

StructuredColumns contains a struct template FixupHaloForVectors which is applied to fields which have type set to vector in their metadata. This multiplies field values beyond the poles by -1 directly after a halo exchange.

The SphericalVector interpolation scheme also checks for type == vector before applying the vector interpolation method. However, the factor of -1 applied by FixupHaloForVectors corrupts the interpolation at the poles (see below).

The work-around, currently used in the atlas_test_interpolation_spherical_vector ctest is to make sure that the dirty value in the field metadata is set to false. This produces a clean interpolation at the poles (see below), but it's clearly not useful in the general case where halo-exchanges are required.

What are the steps to reproduce the bug?

The bug can be reproduced by commenting out this line from the interpolation test.

Version

develop

Platform (OS and architecture)

All builds

Relevant log output

No response

Accompanying data

No response

Organisation

Met Office

[odlomax]

Hi @wdeconinck

I've drafted this PR that demonstrates a workaround for the StructuredColumns halo exchange. Basically, you've got to make sure the type in the field metadata is not set to vector when the halo exchange is called.

I'm happy to add an option to the StructuredColumns config that turns off FixupHalos, but it seems best to check with you in case you have any other ideas.

[odlomax]

Through a bit of experimentation, I found that that using the halo "longitudes" in the complex weight generation almost completely cancels out the factor of minus one applied by FixupHalos, i.e.:

weight(lon_true, lat) ~= weight(lon_halo, lat) * exp(i pi) = - weight(lon_halo, lat)

The error on this approximation is of order [pi - abs(lat)]^2, which is about 1e-4 for source points about a degree off the pole, and reduces quadratically to zero at the pole.

tl;dr, if was okay with very tiny inaccuracies, we can simply use the StructuredColumns lonlat field and forget about FixupHalos entirely.

p.s. I'm surprised this tripped me up at all first time around...

[wdeconinck]

That's great @odlomax , could you show a zoomed in picture around the Pole to show the difference in error visually?

[odlomax]

That's great @odlomax , could you show a zoomed in picture around the Pole to show the difference in error visually?

I suspect I've gotten my maths slightly wrong. Either the error is machine precision or very tiny, as stated above. I might have to pull apart this equation with a pencil a paper to give a definitive answer.

[odlomax]

I did some maths!

Basically, an exact phase difference of pi sneaks into the weights when you use structured columns (lon, lat)_halo instead of (lon, lat)_true. This cancels the -1 from FixupHalos.

This caught me off guard, because (lon, lat)_halo and (lon, lat)_true both map to identical (x, y, z) coordinates.

tl;dr; there is no additional error (as far as I can tell). I got cocky with trigonometry and it bit me!

I can post some TeX here later if you want to see what happened.

[odlomax]

Following the conventions here, we can define:
$\phi_1 = \frac{\pi}{2} - \Delta\phi; \quad\phi_1' = \frac{\pi}{2} + \Delta\phi;\quad\lambda_1' = \lambda_1\pm\pi,$
where $\Delta\phi$ is an altitude angle descending from the North Pole. If $(\phi_1, \lambda_1)$ are the coordinates of an owned point near the pole, then $(\phi_1', \lambda_1')$ are the coordinates of the corresponding halo point near the pole.

Examining the great-circle course equation:
$\alpha_1 = \mathrm{atan2}(y_1, x_1);$
$\alpha_2 = \mathrm{atan2}(y_2, x_2),$
where,
$y_1 = \cos\phi_2\sin\lambda_{12};$
$x_1 = \cos\phi_1\sin\phi_2 - \sin\phi_1\cos\phi_2\cos\lambda_{12};$
$y_2 = \cos\phi_1\sin\lambda_{12};$
$x_2 = -\cos\phi_2\sin\phi_1 - \sin\phi_2\cos\phi_1\cos\lambda_{12}.$

If we recall the following trigonometric identities,
$\cos(\theta\pm\pi) = - \cos\theta;$
$\sin(\theta\pm\pi) = - \sin\theta;$
$\cos\left(\frac{\pi}{2}+\theta\right) = - \cos\left(\frac{\pi}{2}-\theta\right);$
$\sin\left(\frac{\pi}{2}+\theta\right) = \sin\left(\frac{\pi}{2}-\theta\right),$

and we substitute $(\phi_1',\lambda_1')$ into the great-circle course equation, we find that:
$y_1' = -y_1;\quad x_1' = -x_1;\quad y_2' = y_2;\quad x_2' = x_2,$
which means that $\alpha_1' = \alpha_1 \pm\pi$ and $\alpha_2' = \alpha_2$.

If our owned-point interpolation weight is $w_{ij} \exp(\mathrm{i}\alpha_{21})$, then the corresponding halo weight is $w_{ij} \exp\left(\mathrm{i}(\alpha_{21} \pm\pi)\right) = -w_{ij} \exp(\mathrm{i}\alpha_{21})$, which conveniently cancels out the modification from FixupHalos.

[wdeconinck]

Thanks for the extensive explanations! This is now merged.