Oceananigans.jl: halo points behaving differently in master branch compared with 0.36.0
I am using the master branch for nutrient tracers with a 3rd DST advection scheme. I found that the halo points don't have the same values as the boundaries in a periodic domain. Below is the configuration I use and a slice of u velocity. Not sure it's a bug or you intended to do so...
Also, the fancy logger disappeared after I updated to 0.36.0
grid = RegularCartesianGrid(size=(32, 2, 32), extent=(2*32, 2*2, 2*32), halo = (2, 2, 2))
RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}
domain: x ∈ [0.0, 66.0], y ∈ [0.0, 6.0], z ∈ [-64.0, 2.0]
topology: (Periodic, Periodic, Bounded)
resolution (Nx, Ny, Nz): (32, 2, 32)
halo size (Hx, Hy, Hz): (2, 2, 2)
grid spacing (Δx, Δy, Δz): (2.0, 2.0, 2.0)
```julia
T_bcs = TracerBoundaryConditions(Ogrid, top = BoundaryCondition(Flux, Qáµ€),
bottom = BoundaryCondition(Gradient, ∂T∂z))
model = IncompressibleModel(
architecture = CPU(),
grid = grid,
coriolis = FPlane(f=f),
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=α, β=β)),
closure = AnisotropicMinimumDissipation(),
boundary_conditions = (T=T_bcs,)
)
```julia
model.velocities.u.data.parent[:,:,3]
36×6 Array{Float64,2}:
-0.0085938 -0.00844015 -0.0085938 -0.00844015 -0.0085938 -0.00844015
-0.010009 -0.0104645 -0.010009 -0.0104645 -0.010009 -0.0104645
-0.00859588 -0.00912747 -0.00829848 -0.00823567 -0.00859588 -0.00912747
-0.00742535 -0.00655195 -0.00594664 -0.00595502 -0.00742535 -0.00655195
-0.00536066 -0.0038463 -0.00437428 -0.00422142 -0.00536066 -0.0038463
-0.00241687 -0.00325761 -0.00216266 -0.00211736 -0.00241687 -0.00325761
0.00258328 0.00188971 0.00159199 0.00141804 0.00258328 0.00188971
0.00289388 0.00378436 0.00257106 0.00259916 0.00289388 0.00378436
0.00363445 0.00407903 0.00342151 0.00360496 0.00363445 0.00407903
0.00575018 0.00480254 0.00481157 0.00447708 0.00575018 0.00480254
0.00667745 0.00577658 0.00495125 0.00499728 0.00667745 0.00577658
0.00464067 0.00487611 0.00543163 0.00506227 0.00464067 0.00487611
0.00418414 0.00342895 0.00512322 0.0045781 0.00418414 0.00342895
0.00527862 0.00577499 0.00466384 0.00464566 0.00527862 0.00577499
0.00437337 0.00429671 0.00476728 0.00471958 0.00437337 0.00429671
0.003937 0.0035391 0.0039966 0.00382132 0.003937 0.0035391
0.0040974 0.0042016 0.00400427 0.00382213 0.0040974 0.0042016
0.00160725 0.00113037 0.00263925 0.00234604 0.00160725 0.00113037
0.00163343 0.000980412 0.00162656 0.00150644 0.00163343 0.000980412
0.000711606 -0.000278507 0.000465185 0.000347766 0.000711606 -0.000278507
-9.73378e-5 0.00143122 0.00054061 0.00055209 -9.73378e-5 0.00143122
0.0023124 0.00251638 0.00250684 0.00271011 0.0023124 0.00251638
0.00363353 0.00380998 0.00309503 0.00301228 0.00363353 0.00380998
0.00192508 0.00112694 0.00231717 0.00216197 0.00192508 0.00112694
0.000759942 0.000315939 8.66382e-5 9.44978e-5 0.000759942 0.000315939
-0.00305414 -0.00279308 -0.00245104 -0.00225708 -0.00305414 -0.00279308
-0.00400213 -0.00332518 -0.00415305 -0.00403566 -0.00400213 -0.00332518
-0.00682336 -0.00530647 -0.00655179 -0.0062438 -0.00682336 -0.00530647
-0.00732349 -0.00843704 -0.00802866 -0.00780196 -0.00732349 -0.00843704
-0.00769086 -0.00736913 -0.00865324 -0.0082686 -0.00769086 -0.00736913
-0.00893387 -0.00923972 -0.00800158 -0.00804158 -0.00893387 -0.00923972
-0.0106417 -0.0106265 -0.00994619 -0.0101376 -0.0106417 -0.0106265
-0.0085938 -0.00844015 -0.0102716 -0.0103836 -0.0085938 -0.00844015
-0.010009 -0.0104645 -0.0101071 -0.0100802 -0.010009 -0.0104645
-0.00859588 -0.00912747 -0.00859588 -0.00912747 -0.00859588 -0.00912747
-0.00742535 -0.00655195 -0.00742535 -0.00655195 -0.00742535 -0.00655195
zhenwu0728
All 16 comments
Will reply to the rest in a bit!
Also, the fancy logger disappeared after I updated to 0.36.0
Ah yeah sorry. It was setting the global logger to OceananigansLogger after using Oceananigans but then I had some issues when using Oceananigans with PyPlot that seemed related to the logger so I stopped hijacking the global logger with Oceananigans.
If you want the Oceananigans logger now you should add
using Logging
Logging.global_logger(OceananigansLogger())
to the top of your script.
ali-ramadhan
on 16 Sep 2020
Hmmm, yeah indeed the halos aren't exactly correct. For example, rows 1-2 (j=-1,0)
-0.0085938 -0.00844015 -0.0085938 -0.00844015 -0.0085938 -0.00844015
-0.010009 -0.0104645 -0.010009 -0.0104645 -0.010009 -0.0104645
should exactly match rows 33-34 (j=Ny-1,Ny) but don't
-0.0085938 -0.00844015 -0.0102716 -0.0103836 -0.0085938 -0.00844015
-0.010009 -0.0104645 -0.0101071 -0.0100802 -0.010009 -0.0104645
I don't think this is a problem for time-stepping as the Oceananigans algorithm only fills halo regions right before it is needed (to improve performance).
If you need the halo regions to be correct then you could manually call fill_halo_regions!(model.velocities.u, CPU()).
A cleaner alternative if this is a desirable feature might be to add a fill_halo_regions keyword argument to time_step! (and Simulation) to ensure all halos are filled and the fields are consistent at the end of each time step.
ali-ramadhan
on 16 Sep 2020
Thanks!
If you need the halo regions to be correct then you could manually call
fill_halo_regions!(model.velocities.u, CPU()).
Yes, I tried to fill halo regions myself and it worked.
A cleaner alternative if this is a desirable feature might be to add a
fill_halo_regionskeyword argument totime_step!(andSimulation) to ensure all halos are filled and the fields are consistent at the end of each time step.
I don't know if someone else also needs this feature.
But I can handle this myself in the 3rd order DST scheme I have.
zhenwu0728
on 16 Sep 2020
I did a quick test (see below) and indeed it does seem like the halos look different between Oceananigans v0.36.0 and master.
I believe this is because @glwagner modified how halos are filled in PRs #904 and #929 but I don't think it should affect your simulations. But yeah, if you need good halo regions then I might manually fill them in at the end of each time step (or just before when you actually need them).
using Oceananigans
grid = RegularCartesianGrid(size=(3, 3, 3), extent=(1, 1, 1))
model = IncompressibleModel(grid=grid)
set!(model, u = (x, y, z) -> x)
time_step!(model, 1)
model.velocities.u.data.parent
Oceananigans v0.36.0
5×5×5 Array{Float64,3}:
[:, :, 1] =
1.33332 1.33332 1.33332 1.33332 1.33332
0.400015 0.333333 0.333333 0.333333 0.400015
-0.733333 0.333333 0.333333 0.333333 -0.733333
1.33332 0.333333 0.333333 0.333333 1.33332
0.400015 0.400015 0.400015 0.400015 0.400015
[:, :, 2] =
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
[:, :, 3] =
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
[:, :, 4] =
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
0.333333 0.333333 0.333333 0.333333 0.333333
[:, :, 5] =
1.33332 1.33332 1.33332 1.33332 1.33332
0.400015 0.333333 0.333333 0.333333 0.400015
-0.733333 0.333333 0.333333 0.333333 -0.733333
1.33332 0.333333 0.333333 0.333333 1.33332
0.400015 0.400015 0.400015 0.400015 0.400015
Oceananigans#master
5×5×5 Array{Float64,3}:
[:, :, 1] =
0.666667 0.666667 0.666667 0.666667 0.666667
0.0 0.400015 0.400015 0.400015 0.0
0.333333 -0.733333 -0.733333 -0.733333 0.333333
0.666667 1.33332 1.33332 1.33332 0.666667
0.0 0.0 0.0 0.0 0.0
[:, :, 2] =
1.33332 1.33332 1.33332 1.33332 1.33332
0.400015 0.333333 0.333333 0.333333 0.400015
-0.733333 0.333333 0.333333 0.333333 -0.733333
1.33332 0.333333 0.333333 0.333333 1.33332
0.400015 0.400015 0.400015 0.400015 0.400015
[:, :, 3] =
1.33332 1.33332 1.33332 1.33332 1.33332
0.400015 0.333333 0.333333 0.333333 0.400015
-0.733333 0.333333 0.333333 0.333333 -0.733333
1.33332 0.333333 0.333333 0.333333 1.33332
0.400015 0.400015 0.400015 0.400015 0.400015
[:, :, 4] =
1.33332 1.33332 1.33332 1.33332 1.33332
0.400015 0.333333 0.333333 0.333333 0.400015
-0.733333 0.333333 0.333333 0.333333 -0.733333
1.33332 0.333333 0.333333 0.333333 1.33332
0.400015 0.400015 0.400015 0.400015 0.400015
[:, :, 5] =
0.666667 0.666667 0.666667 0.666667 0.666667
0.0 0.400015 0.400015 0.400015 0.0
0.333333 -0.733333 -0.733333 -0.733333 0.333333
0.666667 1.33332 1.33332 1.33332 0.666667
0.0 0.0 0.0 0.0 0.0
I am using the master branch for nutrient tracers with a 3rd DST advection scheme.
That's awesome!
I know we ended up closing your last PR #183 since we didn't have a clean way of adding new advection schemes but we've recently gotten serious about adding higher-order advection schemes (see #815 and #592).
If you have time it would be awesome if you could open a PR with the 3rd order DST advection scheme and we can work together on merging it into Oceananigans so other people could use it.
ali-ramadhan
on 16 Sep 2020
If you have time it would be awesome if you could open a PR with the 3rd order DST advection scheme and we can work together on merging it into Oceananigans so other people could use it.
Sure! I'll open a PR about that. I'm also using KernelAbstractions.jl for the scheme, so it would not be hard to implement.
zhenwu0728
on 16 Sep 2020
Apologies for not understanding --- what is incorrect? I couldn't find a problem with the arrays originally posted.
Halo filling for flux boundary conditions was changed. This could affect corner points?
I'm not sure how an advection scheme could be implemented external to the time_step! function but I think some care should be taken. It'd be better (more accurate, less prone to bugs) to change the native advection scheme. We also have RK3 time stepping now, which could improve the accuracy of tracer advection substantially when using a higher-order scheme.
glwagner
on 16 Sep 2020
As for @ali-ramadhan's suggestion --- I think the solution here and elsewhere should be to compute things immediately prior to use.
Now that I think about it, halo regions would never be correct for the velocity field after time step is complete, either nor or in prior versions of Oceananigans. At the end of a time step, the velocity field is advanced to new values and halos would have to be refilled for them to be correct. We fill halos prior to making computations, which ensures they are correct during the computation in question.
glwagner
on 16 Sep 2020
Apologies for not understanding --- what is incorrect? I couldn't find a problem with the arrays originally posted.
I don't think anything is incorrect. I think @zhenwu0728 noticed a change and assumed it may be a bug, or at least a regression. But both are fine, halo filling was just done a little differently.
ali-ramadhan
on 16 Sep 2020
What change was observed?
glwagner
on 17 Sep 2020
It'd be great to have more time-stepping schemes!
glwagner
on 17 Sep 2020
What change was observed?
As you see in my original message, u.parent[3,3,3] = -0.00829848 should be the same as u.parent[35,3,3] = -0.00859588, but it's not. It's because the halo regions are not filled after each time step.
I noticed this difference because I use halo regions for the advection of a periodic domain and the mass was not conserved after I updated to the master branch.
Now I know why the halo regions are not correct and I can fill the halo regions myself.
Problem solved!
zhenwu0728
on 17 Sep 2020
Ah ok! Was the behavior different in prior versions of Oceananigans? This is what I would like to understand. We did not fill halo regions after time stepping in prior versions of Oceananigans, either, as far as I know.
Do you want to contribute an advection scheme to Oceananigans? A lot of people might benefit from having the kind of advection scheme that you are using! It may not be very much work, since we have an interface for adding new advection schemes.
Do you evolve the tracer fully outside of the Oceananigans infrastructure? Or is the tracer "included" in the list of tracers for your IncompressibleModel, but you somehow negate tracer advection and add it outside of time_step!, somehow? It would be interesting to see this method, and figure out how to make it easier for users to solve problems like the one you're encountering in the future.
glwagner
on 17 Sep 2020
Do you want to contribute an advection scheme to Oceananigans? A lot of people might benefit from having the kind of advection scheme that you are using! It may not be very much work, since we have an interface for adding new advection schemes.
Sure! I'll submit a PR soon.
Do you evolve the tracer fully outside of the
Oceananigansinfrastructure?
Yes, it's fully outside of Oceananigans. I just use the velocities from Oceananigans.
zhenwu0728
on 17 Sep 2020
Should we resolve this issue?
glwagner
on 18 Sep 2020
Should we resolve this issue?
Yes, I think so. I'll submit a PR of the advection scheme soon.
zhenwu0728
on 18 Sep 2020
Related issues
glwagner
·
3Comments
glwagner
·
6Comments
ali-ramadhan
·
7Comments
freemin7
·
5Comments
ali-ramadhan
·
5Comments