C grid interpolation for VectorFields #2152
Conversation
Using correct mesh_type and conversion Also had to change the yi-indexing. To explore why!
parcels/tools/interpolation_utils.py
Outdated
| dxdeta = np.dot(px, dphideta) | ||
| dydxsi = np.dot(py, dphidxsi) | ||
| dydeta = np.dot(py, dphideta) | ||
| dxdxsi = np.dot(dphidxsi, px) |
There was a problem hiding this comment.
@erikvansebille, following up on our discussion about these four
dxdxsi_diag = np.einsum('ij,ji->i', dphidxsi, px)
dxdeta_diag = np.einsum('ij,ji->i', dphideta, px)
dydxsi_diag = np.einsum('ij,ji->i', dphidxsi, py)
dydeta_diag = np.einsum('ij,ji->i', dphideta, py)
jac_diag = dxdxsi_diag * dydeta_diag - dxdeta_diag * dydxsi_diag
return jac_diagUsing vectorized version of dot-product calculation as suggested by @michaeldenes in #2152 (comment)
Since np.where also calculates the parts that are not needed, these often throw a warning
VeckoTheGecko
left a comment
VeckoTheGecko
left a comment
There was a problem hiding this comment.
Still need to have a proper look at application_kernels/interpolation.py
There was a problem hiding this comment.
Do we have any tests that test C-grid vs A grid interpolation?
| with np.errstate(divide="ignore", invalid="ignore"): | ||
| det = np.where(det2 > 0, np.sqrt(det2), eta) | ||
|
|
||
| eta = np.where(abs(aa) < 1e-12, -cc / bb, np.where(det2 > 0, (-bb + det) / (2 * aa), eta)) | ||
|
|
There was a problem hiding this comment.
Which scenarios would produce the errors in particular are we ignoring here? And is the rest of the code robust to the nans/inf values that ignoring the errors would produce?
There was a problem hiding this comment.
The way that np.where works is that it always computes both the True and False results (as far as I understand).
That means that it will for example take the np.sqrt(det2) even if det2 < 0, leading to a lot of warnings (have you not seen them in your test results?).
The with np.errstate() filters these warnings out
So yes, the code is robust to NaNs and Infs, because these would not be warnings but errors
There was a problem hiding this comment.
one option would be to do:
mask = det2 > 0
det[mask] = np.sqrt(det2[mask])lets just merge for now and leave this 'unresolved' in the PR - minor thing
This PR implements C-grid interpolation for VectorFields (and Tracer Fields), building on #2122 (which will need to be merged into
v4-devfirst)mainfor v3 changes,v4-devfor v4 changes)