Construct grid parameters with `BigFloat` to avoid round-off error

[glwagner]

This PR cherry-picks a commit from #2253 . Before this PR:

julia> grid = RectilinearGrid(size=(128, 128), extent=(2π, 2π), topology=(Periodic, Periodic, Flat))
128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── Periodic x ∈ [-7.51279e-18, 6.28319) regularly spaced with Δx=0.0490874
├── Periodic y ∈ [-7.51279e-18, 6.28319) regularly spaced with Δy=0.0490874
└── Flat z

after:

julia> grid = RectilinearGrid(size=(128, 128), extent=(2π, 2π), topology=(Periodic, Periodic, Flat))
128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── Periodic x ∈ [0.0, 6.28319) regularly spaced with Δx=0.0490874
├── Periodic y ∈ [0.0, 6.28319) regularly spaced with Δy=0.0490874
└── Flat z

Basically, this PR computes the location of the first halo point more accurate, which in turn means that the location of the right interior boundary is computed more accurately, without round-off error (which pollutes main at the moment).

I think our method could be slightly improved. Here we first truncate to Float32, and then expand to BigFloat for subsequent arithmetic. The truncation is important, but also means that the input is limited to single-precision accurate. Is that inevitable, or is there another way to achieve this outcome that doesn't lose precision?

PS we may need to do something similar with precision / round-off error for align_time_step.

[simone-silvestri]

It seems like with this method the grid's length is less accurate...

I ll take a look

[francispoulin]

I don't know BigFloat but I do like the the result from this PR.

[glwagner]

Yeah @simone-silvestri thats my prime concern esp when we use stuff like "2*pi"

[simone-silvestri]

By removing the conversion to Float32 worked for me...

[glwagner]

@simone-silvestri now I get

julia> grid = RectilinearGrid(size=(128, 128), extent=(2π, 2π), topology=(Periodic, Periodic, Flat))
128×128×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── Periodic x ∈ [-7.51279e-18, 6.28319) regularly spaced with Δx=0.0490874
├── Periodic y ∈ [-7.51279e-18, 6.28319) regularly spaced with Δy=0.0490874
└── Flat z

do you get something different?

[glwagner]

@simonbyrne do you have any advice for how to avoid these annoying round-off errors? Basically I'm wondering if theres a robust way to get

├── Periodic x ∈ [0.0, 6.28319) regularly spaced with Δx=0.0490874

here we are given the length of the domain (2pi), but then we have to construct a range of nodes which extends left of 0, so that it includes our halo region. The problem is when we do that, the first face (which is then something like range[4] if we have 3 halo points) has round-off error (or at least that's how I'm interpreting it... maybe there's no way around it).

[simonbyrne]

You want it to be a Range and not a Vector?

[simonbyrne]

So you can make use of a feature I put in there that lets you specify ranges from an arbitrary point:

julia> n = 128
128

julia> r = range(0.0,2pi,length=n+1) # original range
0.0:0.04908738521234052:6.283185307179586

julia> h = 3
3

julia> r = StepRangeLen(r.ref, r.step, r.len + 2*h, r.offset + h) # add halo to range
-0.14726215563702155:0.04908738521234052:6.430447462816608

julia> r[4]
0.0

[glwagner]

You want it to be a Range and not a Vector?

Well... we could convert to Vector now. Originally, it was nice that we didn't have to worry about changing architecture betweeen CPU/GPU (but now we have features for that so it wouldn't matter).

But even if we convert to Vector it would be convenient to build the vector from a range, right?

So you can make use of a feature I put in there that lets you specify ranges from an arbitrary point:

julia> n = 128
128

julia> r = range(0.0,2pi,length=n+1) # original range
0.0:0.04908738521234052:6.283185307179586

julia> h = 3
3

julia> r = StepRangeLen(r.ref, r.step, r.len + 2*h, r.offset + h) # add halo to range
-0.14726215563702155:0.04908738521234052:6.430447462816608

julia> r[4]
0.0

nice!

[simone-silvestri]

Thanks @simonbyrne, that works beautifully!

[glwagner]

@simone-silvestri i guess this was superceded by #2253 maybe?

[simone-silvestri]

Not really, I just used BigFloat, it's better but still not perfect (now the round-off is down to 1e-18). The range implementation is much better, we should get that in

[navidcy]

@glwagner, merge main and resolve conflicts?

[glwagner]

This is stale and not too urgent so I'm closing