Double precision

[delandmeterp]

This PR solves a numerical issue that we observe with Parcels.

In coastal flows, the particles tend to follow the coastlines, and can be very close to the coast. Impermeable boundary condition of the input data combined with a small time step (~15 minutes) should prevent the particles of beaching, but rounding errors result in beaching particles. By converting lon, lat, depth to double precision, we reduce this un-controlled process by ~95%.

So it can have a large effect on the particle general dynamics.

How could we solve this problem differently?

• Using unbeaching kernels (they have a variable effect depending on mesh resolution)
• Hacking the interpolation scheme (if UV is zero, and I am just on the edge of the beach, I put back the particle on the other side of the edge)
• By simply putting double precision as a flag of the particle set.

Any suggestion on the solution to adopt?

[delandmeterp]

About commit [694c8df]

It seems that we had to increase a lot the tolerance to pass the test. It is not really the case.
Indeed, we were comparing with a very small tolerance numbers stored as float, which resulted before in a difference of 0.

See attached code for a simple example: a is considered equal to b because the tolerance is smaller than float accuracy.

#include <stdio.h>
#include <math.h>

int main(){
  float a = 4.5;
  a += 1e-8;
  float b = 4.5;
  if (fabs(a-b) < 1e-12)
    printf("a == b\n");
  else
    printf("a !=  b\n");

  return 0;
}

[willirath]

There's considerable efforts in the modelling community to explore reduced precision, because this will yield much higher performance on, e.g., GPUs. Just defaulting to double is contrary to this tendency.

Would it make sense to consciously handle precision and formulate tolerances in terms of machine precision?

[delandmeterp]

Indeed and we try to avoid using double precision where it's not necessary. Then why this change?

It solves a numerical issue. In coastal flows, the particles tend to follow the coastlines, and can be very close to the coast. Impermeable boundary condition of the input data combined with a small time step (~15 minutes) should prevent the particles of beaching, but rounding errors result in beaching particles. By converting lon, lat, depth to double precision, we reduce this un-controlled process by ~95%.

So it can have a large effect on the particle general dynamics.

How could we solve this problem differently?

• Using unbeaching kernels (they have a variable effect depending on mesh resolution)
• Hacking the interpolation scheme (if UV is zero, and I am just on the edge of the beach, I put back the particle on the other side of the edge)
• By simply putting double precision as a flag of the particle set.

Any suggestion on the solution to adopt?

[erikvansebille]

Another option is to determine the precision based on the dtype of the lon and lat. Would that not be easier, and more backward-compatible?
The drawback would be that users may still by default use np.float32 for cgrid, but we could add to the nemo-tutorials that we recommend dtype=np.float64 for c-grids?
What do you think?

[delandmeterp]

User generally don't define a dtype to plon and plat. but the dtype is in pset constructor.

• our default way of defining it.

[erikvansebille]

So does this mean that all user code where dtype is explicitly set now breaks? Is this a reason to determine precision based on the lon/lat dtype?

[delandmeterp]

No. It was breaking only if we were comparing in python the value of a float with a double

[erikvansebille]

Shouldn't we then also rename the options to np.float32 and np.float64?

[delandmeterp]

Done

[VeckoTheGecko]

@erikvansebille this was the PR introducing the lon/lat/depth precision (related to our discussion today about removing this arg and making it a user warning).