Reducing resolution with `change_resolution` can lead to inconsistent boundary

[f0uriest]

eq = desc.examples.get("W7-X")
eq.change_resolution(3,3,3,6,6,6)

assert eq.M == 3
assert eq.surface.M == 3

R1 = eq.surface.R_lmn # these are coeffs from truncating boundary representation to M,N=3
R2 = eq.get_surface_at(rho=1.0).R_lmn # coeffs from evaluating truncate fourier-zernike series at rho=1
np.testing.assert_allclose(R1, R2) # fails

Not sure if this is really a bug, I think it's somewhat expected. IE, to start we have Rb[M,N] = sum R[L, M, N] over L This is fine if you only truncate in M, N since fourier modes are orthogonal, but if you also truncate in L (which is commonly the case when reducing resolution) then the sum has fewer terms so we don't expect them to match.

A few options

1. always do eq.surface = eq.get_surface_at after changing resolution. However, this means the surface will be different that one obtained just by truncating which may be unexpected
2. Throw a warning when reducing L resolution?
3. Adding note to the docstring?

I expect we'd see similar behavior in the axis when truncating in M, or in the xsection when truncating in N.