Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove Tullio dependency; take #2 and implement universally accurate CFL time-scale #3037

Merged
merged 16 commits into from
Apr 11, 2023
Merged
Show file tree
Hide file tree
Changes from 9 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 1 addition & 19 deletions Manifest.toml
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

julia_version = "1.8.5"
manifest_format = "2.0"
project_hash = "96eebad427fe7f986c00777f7738aa636a9e621d"
project_hash = "31dd08f0370dcfe9232b9df8c2131e894bcff3ca"

[[deps.AMGX]]
deps = ["AMGX_jll", "CEnum", "CUDA", "JSON", "Libdl", "SparseArrays"]
Expand Down Expand Up @@ -156,12 +156,6 @@ version = "1.0.0"
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DiffRules]]
deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "c5b6685d53f933c11404a3ae9822afe30d522494"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.12.2"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
Expand Down Expand Up @@ -431,12 +425,6 @@ git-tree-sha1 = "d3e32d2b0530d929a047ceab8799ae4204de6c88"
uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
version = "0.12.12"

[[deps.NaNMath]]
deps = ["OpenLibm_jll"]
git-tree-sha1 = "a7c3d1da1189a1c2fe843a3bfa04d18d20eb3211"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.1"

[[deps.NetCDF_jll]]
deps = ["Artifacts", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "072f8371f74c3b9e1b26679de7fbf059d45ea221"
Expand Down Expand Up @@ -703,12 +691,6 @@ git-tree-sha1 = "94f38103c984f89cf77c402f2a68dbd870f8165f"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.11"

[[deps.Tullio]]
deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"]
git-tree-sha1 = "7871a39eac745697ee512a87eeff06a048a7905b"
uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
version = "0.3.5"

[[deps.TupleTools]]
git-tree-sha1 = "3c712976c47707ff893cf6ba4354aa14db1d8938"
uuid = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
Expand Down
2 changes: 0 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,6 @@ SeawaterPolynomials = "d496a93d-167e-4197-9f49-d3af4ff8fe40"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"

[compat]
AMGX = "0.1.3"
Expand All @@ -62,7 +61,6 @@ PencilFFTs = "0.13.5, 0.14"
Rotations = "1.0"
SeawaterPolynomials = "0.3"
StructArrays = "0.4, 0.5, 0.6"
Tullio = "0.2, 0.3"
julia = "1.6"

[extras]
Expand Down
3 changes: 2 additions & 1 deletion src/AbstractOperations/AbstractOperations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,8 @@ include("show_abstract_operations.jl")
# Make some operators!

# Some operators:
import Base: sqrt, sin, cos, exp, tanh, -, +, /, ^, *
import Base: sqrt, sin, cos, exp, tanh, abs, -, +, /, ^, *
import Base: abs

@unary sqrt sin cos exp tanh abs
@unary -
Expand Down
7 changes: 3 additions & 4 deletions src/Advection/cell_advection_timescale.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,7 @@ end
Δy = Δyᶜᶠᶜ(i, j, k, grid)
Δz = Δzᶜᶜᶠ(i, j, k, grid)

return @inbounds min(Δx / abs(u[i, j, k]),
Δy / abs(v[i, j, k]),
Δz / abs(w[i, j, k]))
return @inbounds 1 / (abs(u[i, j, k]) / Δx +
abs(v[i, j, k]) / Δy +
abs(w[i, j, k]) / Δz)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is this the proper way to compute CFL? I'll trust you on that... That's like the geometric mean of the timescales in each direction?

Copy link
Collaborator Author

@simone-silvestri simone-silvestri Mar 29, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Imagine a 2D situation, we start from
$$\frac{\partial c}{\partial t} = u\frac{\partial c}{\partial x} + v\frac{\partial c}{\partial y}$$
We have advection in x and y, by extension from the 1D CFL condition we can approximate the CFL condition as
$$|u|\frac{\Delta t}{\Delta x} + |v|\frac{\Delta t}{\Delta y} < CFL (\text{ maximum 1 })$$
This means that
$$\Delta t / \tau < CFL$$
where
$$\tau = \left(\frac{|u|}{\Delta x} + \frac{|v|}{\Delta y} \right)^{-1}$$
The extension to 3D includes also the advection in z

This is just a choice, as to account for diagonal motion probably you can get away with a square root of 2 somewhere. But it is the usual definition of the CFL condition
https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good! Why don't we include all this in the docstring of the cell_advection?

simone-silvestri marked this conversation as resolved.
Show resolved Hide resolved
end

29 changes: 2 additions & 27 deletions src/Diagnostics/cfl.jl
Original file line number Diff line number Diff line change
Expand Up @@ -83,31 +83,6 @@ DiffusiveCFL(Δt) = CFL(Δt, cell_diffusion_timescale)
##### Accurate CFL via reduction
#####

using CUDA, CUDAKernels, KernelAbstractions, Tullio

using Oceananigans.Models
using Oceananigans.Grids: halo_size
using Oceananigans.Operators: Δxᶠᶜᶜ, Δyᶜᶠᶜ, Δzᶜᶜᶠ

accurate_cell_advection_timescale(model) = accurate_cell_advection_timescale(model.grid, model.velocities)

function accurate_cell_advection_timescale(grid, velocities)
Nx, Ny, Nz = size(grid)
Hx, Hy, Hz = halo_size(grid)

is = 1+Hx:Nx+Hx
js = 1+Hy:Ny+Hy
ks = 1+Hz:Nz+Hz

u = view(velocities.u.data.parent, is, js, ks)
v = view(velocities.v.data.parent, is, js, ks)
w = view(velocities.w.data.parent, is, js, ks)

min_timescale = minimum(
@tullio (min) timescale[k] := 1 / ( abs(u[i, j, k]) / Δxᶠᶜᶜ(i, j, k, grid)
+ abs(v[i, j, k]) / Δyᶜᶠᶜ(i, j, k, grid)
+ abs(w[i, j, k]) / Δzᶜᶜᶠ(i, j, k, grid))
)
using Oceananigans.Advection: cell_advection_timescale
navidcy marked this conversation as resolved.
Show resolved Hide resolved

return min_timescale
end
accurate_cell_advection_timescale(model) = cell_advection_timescale(model.grid, model.velocities)
Copy link
Member

@glwagner glwagner Apr 5, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why we would keep the function accurate_cell_advection_timescale as a duplicate of cell_advection_timescale? The word "accurate" is meaningless now.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agree!

glwagner marked this conversation as resolved.
Show resolved Hide resolved