From 626ab8ceb29a8cadfdd211344df3937c816e1e1a Mon Sep 17 00:00:00 2001 From: fquarenghi Date: Tue, 5 Mar 2024 13:59:54 +0100 Subject: [PATCH 01/32] level sets and volume fractions with Chmy backend --- scripts_future_API/Project.toml | 4 +- .../test_levelsets_volumefractions.jl | 109 ++++++++++++++ src/LevelSets/LevelSets.jl | 2 + src/LevelSets/compute_level_sets.jl | 27 ++-- src/LevelSets/compute_volume_fractions.jl | 142 ++++++++++++++++++ src/LevelSets/signed_distances.jl | 21 ++- src/LevelSets/volume_fractions_kernels.jl | 70 +++++++++ 7 files changed, 357 insertions(+), 18 deletions(-) create mode 100644 scripts_future_API/test_levelsets_volumefractions.jl create mode 100644 src/LevelSets/compute_volume_fractions.jl create mode 100644 src/LevelSets/volume_fractions_kernels.jl diff --git a/scripts_future_API/Project.toml b/scripts_future_API/Project.toml index 9b70abea..c6735f55 100644 --- a/scripts_future_API/Project.toml +++ b/scripts_future_API/Project.toml @@ -1,8 +1,8 @@ [deps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" FastIce = "e0de9f13-a007-490e-b696-b07d031015ca" -FastIceTools = "14532042-4700-46e0-8dec-55d517bc1ae0" -FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl new file mode 100644 index 00000000..1de7b789 --- /dev/null +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -0,0 +1,109 @@ +using CUDA +using JLD2 +using FastIce.LevelSets +using Chmy.Grids +using Chmy.Fields +using Chmy.Architectures + + +# Data +# vavilov_path = "../data/vavilov.jld2" +# vavilov_path = "../data/Vavilov/Vavilov_80m.jld2" +vavilov_path = "../data/Vavilov/Vavilov_500m.jld2" +# synthetic_data = "../data/Synthetic/dome_glacier.jld2" +# synthetic_data = "../data/Synthetic/low_res_dome_glacier.jld2" + +# Select backend +backend = CPU() +# backend = CUDABackend() +arch = Arch(backend) + +""" + load_dem_on_GPU(backend, arch::Architecture, path::String) + +Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, +initiate grids, copy data on gpu. +""" +function load_dem_on_GPU(backend, arch::Architecture, path::String) + data = load(path) + z_surf = data["DataElevation"].z_surf + z_bed = data["DataElevation"].z_bed + x = data["DataElevation"].x + y = data["DataElevation"].y + nx = length(x) - 1 + ny = length(y) - 1 + nz = 100 + # TODO: choose nz accordingly + surf_lz = maximum(z_surf) - minimum(z_surf) + bed_lz = maximum(z_bed) - minimum(z_bed) + surf_Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_surf)), extent=(surf_lz, surf_lz, surf_lz), dims=(nx, ny, nz)) + bed_Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_bed)), extent=(bed_lz, bed_lz, bed_lz), dims=(nx, ny, nz)) + surf_dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(surf_lz, surf_lz), dims=(nx, ny)) + bed_dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(bed_lz, bed_lz), dims=(nx, ny)) + surf_field = Field(backend, surf_dem_grid, Vertex()) + bed_field = Field(backend, bed_dem_grid, Vertex()) + set!(surf_field, z_surf) + set!(bed_field, z_bed) + return surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid +end + + +""" + load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::String) + +Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, +initiate grids, copy data on gpu. +""" +function load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::String) + data = load(synthetic_data) + z_surf = data["SyntheticElevation"].z_surf + z_bed = data["SyntheticElevation"].z_bed + nx = size(z_surf)[1] - 1 + ny = size(z_surf)[2] - 1 + nz = 10 + lz = maximum(z_surf) - minimum(z_surf) + Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_surf)), extent=(lz, lz, lz), dims=(nx, ny, nz)) + # Ψ = Field(backend, Ψ_grid, Vertex()) + dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(lz, lz), dims=(nx, ny)) + dem_bed = Field(backend, dem_grid, Vertex()) + dem_surf = Field(backend, dem_grid, Vertex()) + set!(dem_bed, z_bed) + set!(dem_surf, z_surf) + return dem_surf, dem_bed, dem_grid, Ψ_grid +end + +# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem_on_GPU(backend, arch, synthetic_data); +surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem_on_GPU(backend, arch, vavilov_path); + +Ψ = ( + na=Field(backend, surf_Ψ_grid, Vertex()), + ns=Field(backend, bed_Ψ_grid, Vertex()), +) + +compute_level_set_from_dem!(backend, Ψ[1], surf_field, surf_dem_grid, surf_Ψ_grid) +compute_level_set_from_dem!(backend, Ψ[2], bed_field, bed_dem_grid, bed_Ψ_grid) + +# for phase in eachindex(Ψ) +# compute_level_set_from_dem!(backend, Ψ[phase], dem_surf, dem_grid, Ψ_grid) +# end + +wt = ( + na=volfrac_field(backend, surf_Ψ_grid), + ns=volfrac_field(backend, bed_Ψ_grid), +) + +compute_volume_fractions_from_level_set!(backend, wt[1], Ψ[1], surf_Ψ_grid) +compute_volume_fractions_from_level_set!(backend, wt[2], Ψ[2], bed_Ψ_grid) + + +# for phase in eachindex(Ψ) +# compute_volume_fractions_from_level_set!(backend, wt[phase], Ψ[phase], Ψ_grid) +# end + +# Save +jldopen(vavilov_path, "a+") do file + file["level_sets_na"] = Array(interior(Ψ[1])) + file["level_sets_ns"] = Array(interior(Ψ[2])) + file["volume_frac_na"] = Array.(interior.(Tuple(wt[1]))) + file["volume_frac_ns"] = Array.(interior.(Tuple(wt[2]))) +end diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index 5eb8428c..37ad3edf 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -1,6 +1,7 @@ module LevelSets export compute_level_set_from_dem! +export volfrac_field, compute_volume_fractions_from_level_set! using Chmy using Chmy.Grids @@ -11,5 +12,6 @@ using LinearAlgebra, GeometryBasics include("signed_distances.jl") include("compute_level_sets.jl") +include("compute_volume_fractions.jl") end diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index 552dc932..3b6b6b36 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -1,20 +1,25 @@ -@kernel inbounds = true function _init_level_set!(Ψ::Field, dem::Field, dem_grid::StructuredGrid, Ψ_grid::StructuredGrid, cutoff, R) - I = @index(Global, Cartesian) - x, y, z = coord(Ψ_grid, location(Ψ), I) +""" + _init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff::AbstractFloat, R::AbstractMatrix) + +Initialize level sets. +""" +@kernel function init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R) + I = @index(Global, NTuple) + x, y, z = coord(Ψ_grid, location(Ψ), I...) P = R * Point3(x, y, z) ud, sgn = sd_dem(P, cutoff, dem, dem_grid) - Ψ[I] = ud * sgn + Ψ[I...] = ud * sgn end + """ - compute_level_set_from_dem!(arch, Ψ, dem, dem_grid, Ψ_grid, R[=LinearAlgebra.I]) + compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) -Compute a level set from a given dem. +Compute level sets from dem. """ -function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::StructuredGrid, Ψ_grid::StructuredGrid, - R=LinearAlgebra.I) - cutoff = 4maximum(spacing(Ψ_grid)) - kernel = _init_level_set!(backend(arch), 256, size(Ψ)) +function compute_level_set_from_dem!(backend, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) + kernel = init_level_set!(backend, 256, size(Ψ)) + cutoff = 4maximum(spacing(Ψ_grid, Center(), 1, 1, 1)) kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) return -end +end \ No newline at end of file diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl new file mode 100644 index 00000000..37c8dbb5 --- /dev/null +++ b/src/LevelSets/compute_volume_fractions.jl @@ -0,0 +1,142 @@ +export volfrac + +@inline perturb(ϕ) = abs(ϕ) > 1e-20 ? ϕ : (ϕ > 0 ? 1e-20 : -1e-20) + +@inline trivol(v1, v2, v3) = 0.5 * abs(cross(v3 - v1, v2 - v1)) + +function volfrac(tri, ϕ::Vec3{T})::T where {T} + v1, v2, v3 = tri + if ϕ[1] < 0 && ϕ[2] < 0 && ϕ[3] < 0 # --- + return trivol(v1, v2, v3) + elseif ϕ[1] > 0 && ϕ[2] > 0 && ϕ[3] > 0 # +++ + return 0.0 + end + @inline vij(i, j) = tri[j] * (ϕ[i] / (ϕ[i] - ϕ[j])) - tri[i] * (ϕ[j] / (ϕ[i] - ϕ[j])) + v12, v13, v23 = vij(1, 2), vij(1, 3), vij(2, 3) + if ϕ[1] < 0 + if ϕ[2] < 0 + trivol(v1, v23, v13) + trivol(v1, v2, v23) # --+ + else + if ϕ[3] < 0 + trivol(v3, v12, v23) + trivol(v3, v1, v12) # -+- + else + trivol(v1, v12, v13) # -++ + end + end + else + if ϕ[2] < 0 + if ϕ[3] < 0 + trivol(v2, v13, v12) + trivol(v2, v3, v13) # +-- + else + trivol(v12, v2, v23) # +-+ + end + else + trivol(v13, v23, v3) # ++- + end + end +end + +function volfrac(rect::Rect2{T}, ϕ::Vec4{T}) where {T} + or, ws = Grids.origin(rect), widths(rect) + v1, v2, v3, v4 = or, or + Vec(ws[1], 0.0), or + ws, or + Vec(0.0, ws[2]) + ϕ1, ϕ2, ϕ3, ϕ4 = perturb.(ϕ) + return volfrac(Vec(v1, v2, v3), Vec3{T}(ϕ1, ϕ2, ϕ3)) + + volfrac(Vec(v1, v3, v4), Vec3{T}(ϕ1, ϕ3, ϕ4)) +end + +@inline tetvol(v1, v2, v3, v4) = abs(det([v2 - v1 v3 - v1 v4 - v1])) / 6.0 + +function volfrac(tet, ϕ::Vec4) + v1, v2, v3, v4 = tet + @inline vij(i, j) = tet[j] * (ϕ[i] / (ϕ[i] - ϕ[j])) - tet[i] * (ϕ[j] / (ϕ[i] - ϕ[j])) + nneg = count(ϕ .< 0) + if nneg == 0 # ++++ + return 0.0 + elseif nneg == 1 # -+++ + if ϕ[1] < 0 + return tetvol(v1, vij(1, 2), vij(1, 3), vij(1, 4)) + elseif ϕ[2] < 0 + return tetvol(v2, vij(2, 1), vij(2, 3), vij(2, 4)) + elseif ϕ[3] < 0 + return tetvol(v3, vij(3, 1), vij(3, 2), vij(3, 4)) + else # ϕ[4] < 0 + return tetvol(v4, vij(4, 1), vij(4, 2), vij(4, 3)) + end + elseif nneg == 2 # --++ + if ϕ[1] < 0 && ϕ[2] < 0 + return tetvol(v1, v2, vij(1, 3), vij(2, 4)) + + tetvol(vij(2, 3), v2, vij(1, 3), vij(2, 4)) + + tetvol(v1, vij(1, 4), vij(1, 3), vij(2, 4)) + elseif ϕ[1] < 0 && ϕ[3] < 0 + return tetvol(v1, v3, vij(1, 4), vij(3, 2)) + + tetvol(vij(3, 4), v3, vij(1, 4), vij(3, 2)) + + tetvol(v1, vij(1, 2), vij(1, 4), vij(3, 2)) + elseif ϕ[1] < 0 && ϕ[4] < 0 + return tetvol(v1, v4, vij(1, 2), vij(4, 3)) + + tetvol(vij(4, 2), v4, vij(1, 2), vij(4, 3)) + + tetvol(v1, vij(1, 3), vij(1, 2), vij(4, 3)) + elseif ϕ[2] < 0 && ϕ[3] < 0 + return tetvol(v3, v2, vij(3, 1), vij(2, 4)) + + tetvol(vij(2, 1), v2, vij(3, 1), vij(2, 4)) + + tetvol(v3, vij(3, 4), vij(3, 1), vij(2, 4)) + elseif ϕ[2] < 0 && ϕ[4] < 0 + return tetvol(v4, v2, vij(4, 1), vij(2, 3)) + + tetvol(vij(2, 1), v2, vij(4, 1), vij(2, 3)) + + tetvol(v4, vij(4, 3), vij(4, 1), vij(2, 3)) + else # ϕ[3] < 0 && ϕ[4] < 0 + return tetvol(v3, v4, vij(3, 1), vij(4, 2)) + + tetvol(vij(4, 1), v4, vij(3, 1), vij(4, 2)) + + tetvol(v3, vij(3, 2), vij(3, 1), vij(4, 2)) + end + elseif nneg == 3 # ---+ + vol_tot = tetvol(v1, v2, v3, v4) + if ϕ[1] >= 0 + return vol_tot - tetvol(v1, vij(1, 2), vij(1, 3), vij(1, 4)) + elseif ϕ[2] >= 0 + return vol_tot - tetvol(v2, vij(2, 1), vij(2, 3), vij(2, 4)) + elseif ϕ[3] >= 0 + return vol_tot - tetvol(v3, vij(3, 1), vij(3, 2), vij(3, 4)) + else # ϕ[4] >= 0 + return vol_tot - tetvol(v4, vij(4, 1), vij(4, 2), vij(4, 3)) + end + else # ---- + return tetvol(v1, v2, v3, v4) + end +end + +function volfrac(rect::Rect3, ϕ::Vec{8}) + or, ws = GeometryBasics.origin(rect), widths(rect) + v000, v001, v100, v101 = or, or + Vec(ws[1], 0.0, 0.0), or + Vec(0.0, ws[2], 0.0), or + Vec(ws[1], ws[2], 0.0) + v010, v011, v110, v111 = or + Vec(0.0, 0.0, ws[3]), or + Vec(ws[1], 0.0, ws[3]), or + Vec(0.0, ws[2], ws[3]), + or + Vec(ws[1], ws[2], ws[3]) + ϕ = perturb.(ϕ) + return volfrac(Vec(v000, v100, v010, v001), Vec(ϕ[1], ϕ[5], ϕ[3], ϕ[2])) + + volfrac(Vec(v110, v100, v010, v111), Vec(ϕ[7], ϕ[5], ϕ[3], ϕ[7])) + + volfrac(Vec(v101, v100, v111, v001), Vec(ϕ[6], ϕ[5], ϕ[7], ϕ[2])) + + volfrac(Vec(v011, v111, v010, v001), Vec(ϕ[4], ϕ[7], ϕ[3], ϕ[2])) + + volfrac(Vec(v111, v100, v010, v001), Vec(ϕ[7], ϕ[5], ϕ[3], ϕ[2])) +end + +function volfrac_field(backend::Backend, grid::UniformGrid{2}, args...; kwargs...) + (c=Field(backend, grid, Center()), + x=Field(backend, grid, (Vertex(), Center())), + y=Field(backend, grid, (Center(), Vertex())), + xy=Field(backend, grid, Vertex())) +end + +function volfrac_field(backend::Backend, grid::UniformGrid{3}, args...; kwargs...) + (c=Field(backend, grid, Center()), + x=Field(backend, grid, (Vertex(), Center(), Center())), + y=Field(backend, grid, (Center(), Vertex(), Center())), + z=Field(backend, grid, (Center(), Center(), Vertex())), + xy=Field(backend, grid, (Vertex(), Vertex(), Center())), + xz=Field(backend, grid, (Vertex(), Center(), Vertex())), + yz=Field(backend, grid, (Center(), Vertex(), Vertex()))) +end + +include("volume_fractions_kernels.jl") + +function compute_volume_fractions_from_level_set!(backend::Backend, wt, Ψ, grid::UniformGrid) + kernel = kernel_compute_volume_fractions_from_level_set!(backend, 256, size(Ψ)) + kernel(wt, Ψ, grid) +end \ No newline at end of file diff --git a/src/LevelSets/signed_distances.jl b/src/LevelSets/signed_distances.jl index 5af7c964..36d97ee3 100644 --- a/src/LevelSets/signed_distances.jl +++ b/src/LevelSets/signed_distances.jl @@ -1,5 +1,8 @@ +export sd_dem + @inline S(x) = x == zero(x) ? oneunit(x) : sign(x) -@inline sign_triangle(p, a, b, c) = S(dot(p - a, cross(b - a, c - a))) +sign_triangle(p, a, b, c) = S(dot(p - a, cross(b - a, c - a))) + @inline function ud_triangle(p, a, b, c) dot2(v) = dot(v, v) @@ -23,18 +26,22 @@ dot(nor, pa) * dot(nor, pa) / dot2(nor)) end + @inline function closest_vertex_index(P, grid) - Δ = spacing(grid) - O = Grids.origin(grid) - I = @. Int(fld(P - O, Δ)) + 1 + Δs = spacing(grid, Vertex(), 1, 1) + O = Grids.origin(grid, Vertex()) + X = @. fld(P - O, Δs) + I = @. Int(X) + 1 I1 = 1 I2 = size(grid, Vertex()) return clamp.(I, I1, I2) |> CartesianIndex end + @inline inc(I, dim) = Base.setindex(I, I[dim] + 1, dim) @inline inc(I) = I + oneunit(I) + @inline function triangle_pair(Iv, dem, grid) @inline function sample_dem(I) @inbounds x, y = coord(grid, location(dem), I) @@ -45,13 +52,16 @@ end return T_BL, T_TR end + @inline function distance_to_triangle_pair(P, Iv, dem, grid) T_BL, T_TR = triangle_pair(Iv, dem, grid) ud = min(ud_triangle(P, T_BL...), ud_triangle(P, T_TR...)) return ud, sign_triangle(P, T_BL...) end -Base.clamp(p::NTuple{N}, grid::StructuredGrid{N}) where {N} = clamp.(p, Grids.origin(grid), Grids.origin(grid) .+ extent(grid)) + +Base.clamp(p::NTuple{N}, grid::UniformGrid{N}) where {N} = clamp.(p, Grids.origin(grid, Vertex()), Grids.origin(grid, Vertex()) .+ extent(grid, Vertex())) + function sd_dem(P, cutoff, dem, grid) @inbounds Pp = clamp((P[1], P[2]), grid) @@ -68,3 +78,4 @@ function sd_dem(P, cutoff, dem, grid) end return ud, sgn end + diff --git a/src/LevelSets/volume_fractions_kernels.jl b/src/LevelSets/volume_fractions_kernels.jl new file mode 100644 index 00000000..dd16cbf9 --- /dev/null +++ b/src/LevelSets/volume_fractions_kernels.jl @@ -0,0 +1,70 @@ +# Volume fraction kernel (2D) +@kernel function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid) + ix, iy = @index(Global, NTuple) + dx, dy = spacing(grid, Center(), 1, 1) + cell = Rect(Vec(0.0, 0.0), Vec(dx, dy)) + ω = GeometryBasics.volume(cell) + Ψ_ax(dix, diy) = 0.5 * (Ψ[ix+dix, iy+diy] + Ψ[ix+dix+1, iy+diy]) + Ψ_ay(dix, diy) = 0.5 * (Ψ[ix+dix, iy+diy] + Ψ[ix+dix, iy+diy+1]) + Ψ_axy(dix, diy) = 0.25 * (Ψ[ix+dix, iy+diy] + Ψ[ix+dix+1, iy+diy] + Ψ[ix+dix, iy+diy+1] + Ψ[ix+dix+1, iy+diy+1]) + # cell centers + Ψs = Vec{4}(Ψ[ix, iy], Ψ[ix+1, iy], Ψ[ix+1, iy+1], Ψ[ix, iy+1]) + wt.c[ix, iy] = volfrac(cell, Ψs) / ω + # x faces + Ψs = Vec{4}(Ψ_ax(0, 0), Ψ_ax(1, 0), Ψ_ax(1, 1), Ψ_ax(0, 1)) + wt.x[ix, iy] = volfrac(cell, Ψs) / ω + # y faces + Ψs = Vec{4}(Ψ_ay(0, 0), Ψ_ay(1, 0), Ψ_ay(1, 1), Ψ_ay(0, 1)) + wt.y[ix, iy] = volfrac(cell, Ψs) / ω + # xy edges + Ψs = Vec{4}(Ψ_axy(0, 0), Ψ_axy(1, 0), Ψ_axy(1, 1), Ψ_axy(0, 1)) + wt.xy[ix, iy] = volfrac(cell, Ψs) / ω +end + +# Volume fraction kernel (3D) +@kernel function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid) + ix, iy, iz = @index(Global, NTuple) + dx, dy, dz = spacing(grid, Center(), 1, 1, 1) + cell = Rect(Vec(0.0, 0.0, 0.0), Vec(dx, dy, dz)) + ω = GeometryBasics.volume(cell) + Ψ_ax(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix+1, iy+diy, iz+diz]) + Ψ_ay(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix, iy+diy+1, iz+diz]) + Ψ_az(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix, iy+diy, iz+diz+1]) + function Ψ_axy(dix, diy, diz) + 0.25 * + (Ψ[ix+dix, iy+diy, iz+diz+1] + Ψ[ix+dix+1, iy+diy, iz+diz+1] + Ψ[ix+dix, iy+diy+1, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) + end + function Ψ_axz(dix, diy, diz) + 0.25 * + (Ψ[ix+dix, iy+diy+1, iz+diz] + Ψ[ix+dix+1, iy+diy+1, iz+diz] + Ψ[ix+dix, iy+diy+1, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) + end + function Ψ_ayz(dix, diy, diz) + 0.25 * + (Ψ[ix+dix+1, iy+diy, iz+diz] + Ψ[ix+dix+1, iy+diy+1, iz+diz] + Ψ[ix+dix+1, iy+diy, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) + end + # cell centers + Ψs = Vec{8}(Ψ[ix, iy, iz], Ψ[ix+1, iy, iz], Ψ[ix, iy+1, iz], Ψ[ix+1, iy+1, iz], Ψ[ix, iy, iz+1], Ψ[ix+1, iy, iz+1], Ψ[ix, iy+1, iz+1], + Ψ[ix+1, iy+1, iz+1]) + wt.c[ix, iy, iz] = volfrac(cell, Ψs) / ω + # x faces + Ψs = Vec{8}(Ψ_ax(0, 0, 0), Ψ_ax(1, 0, 0), Ψ_ax(0, 1, 0), Ψ_ax(1, 1, 0), Ψ_ax(0, 0, 1), Ψ_ax(1, 0, 1), Ψ_ax(0, 1, 1), Ψ_ax(1, 1, 1)) + wt.x[ix, iy, iz] = volfrac(cell, Ψs) / ω + # y faces + Ψs = Vec{8}(Ψ_ay(0, 0, 0), Ψ_ay(1, 0, 0), Ψ_ay(0, 1, 0), Ψ_ay(1, 1, 0), Ψ_ay(0, 0, 1), Ψ_ay(1, 0, 1), Ψ_ay(0, 1, 1), Ψ_ay(1, 1, 1)) + wt.y[ix, iy, iz] = volfrac(cell, Ψs) / ω + # z faces + Ψs = Vec{8}(Ψ_az(0, 0, 0), Ψ_az(1, 0, 0), Ψ_az(0, 1, 0), Ψ_az(1, 1, 0), Ψ_az(0, 0, 1), Ψ_az(1, 0, 1), Ψ_az(0, 1, 1), Ψ_az(1, 1, 1)) + wt.z[ix, iy, iz] = volfrac(cell, Ψs) / ω + # xy edges + Ψs = Vec{8}(Ψ_axy(0, 0, 0), Ψ_axy(1, 0, 0), Ψ_axy(0, 1, 0), Ψ_axy(1, 1, 0), Ψ_axy(0, 0, 1), Ψ_axy(1, 0, 1), Ψ_axy(0, 1, 1), + Ψ_axy(1, 1, 1)) + wt.xy[ix, iy, iz] = volfrac(cell, Ψs) / ω + # xz edges + Ψs = Vec{8}(Ψ_axz(0, 0, 0), Ψ_axz(1, 0, 0), Ψ_axz(0, 1, 0), Ψ_axz(1, 1, 0), Ψ_axz(0, 0, 1), Ψ_axz(1, 0, 1), Ψ_axz(0, 1, 1), + Ψ_axz(1, 1, 1)) + wt.xz[ix, iy, iz] = volfrac(cell, Ψs) / ω + # yz edges + Ψs = Vec{8}(Ψ_ayz(0, 0, 0), Ψ_ayz(1, 0, 0), Ψ_ayz(0, 1, 0), Ψ_ayz(1, 1, 0), Ψ_ayz(0, 0, 1), Ψ_ayz(1, 0, 1), Ψ_ayz(0, 1, 1), + Ψ_ayz(1, 1, 1)) + wt.yz[ix, iy, iz] = volfrac(cell, Ψs) / ω +end From ca353e008969df51d67e0cec4255b2c219a0d8ef Mon Sep 17 00:00:00 2001 From: Filippo Quarenghi Date: Tue, 5 Mar 2024 16:28:06 +0100 Subject: [PATCH 02/32] Update src/LevelSets/compute_level_sets.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Ludovic Räss <61313342+luraess@users.noreply.github.com> --- src/LevelSets/compute_level_sets.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index 3b6b6b36..e486cedd 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -19,7 +19,7 @@ Compute level sets from dem. """ function compute_level_set_from_dem!(backend, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) kernel = init_level_set!(backend, 256, size(Ψ)) - cutoff = 4maximum(spacing(Ψ_grid, Center(), 1, 1, 1)) + cutoff = 4maximum(spacing(Ψ_grid, Center())) kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) return end \ No newline at end of file From 766136bf215156698c37700de2a81268e708685a Mon Sep 17 00:00:00 2001 From: fquarenghi Date: Tue, 5 Mar 2024 16:31:45 +0100 Subject: [PATCH 03/32] fixed names and toml file --- scripts_future_API/Project.toml | 3 --- .../test_levelsets_volumefractions.jl | 16 +++++++--------- 2 files changed, 7 insertions(+), 12 deletions(-) diff --git a/scripts_future_API/Project.toml b/scripts_future_API/Project.toml index c6735f55..94563595 100644 --- a/scripts_future_API/Project.toml +++ b/scripts_future_API/Project.toml @@ -1,9 +1,6 @@ [deps] -CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" FastIce = "e0de9f13-a007-490e-b696-b07d031015ca" -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl index 1de7b789..965fa1a7 100644 --- a/scripts_future_API/test_levelsets_volumefractions.jl +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -4,7 +4,7 @@ using FastIce.LevelSets using Chmy.Grids using Chmy.Fields using Chmy.Architectures - +using KernelAbstractions # Data # vavilov_path = "../data/vavilov.jld2" @@ -19,12 +19,12 @@ backend = CPU() arch = Arch(backend) """ - load_dem_on_GPU(backend, arch::Architecture, path::String) + load_dem(backend, arch::Architecture, path::String) Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, initiate grids, copy data on gpu. """ -function load_dem_on_GPU(backend, arch::Architecture, path::String) +function load_dem(backend, arch::Architecture, path::String) data = load(path) z_surf = data["DataElevation"].z_surf z_bed = data["DataElevation"].z_bed @@ -47,14 +47,13 @@ function load_dem_on_GPU(backend, arch::Architecture, path::String) return surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid end - """ - load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::String) + load_synth_dem(backend, arch::Architecture, synthetic_data::String) Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, initiate grids, copy data on gpu. """ -function load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::String) +function load_synth_dem(backend, arch::Architecture, synthetic_data::String) data = load(synthetic_data) z_surf = data["SyntheticElevation"].z_surf z_bed = data["SyntheticElevation"].z_bed @@ -63,7 +62,6 @@ function load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::Stri nz = 10 lz = maximum(z_surf) - minimum(z_surf) Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_surf)), extent=(lz, lz, lz), dims=(nx, ny, nz)) - # Ψ = Field(backend, Ψ_grid, Vertex()) dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(lz, lz), dims=(nx, ny)) dem_bed = Field(backend, dem_grid, Vertex()) dem_surf = Field(backend, dem_grid, Vertex()) @@ -72,8 +70,8 @@ function load_synth_dem_on_GPU(backend, arch::Architecture, synthetic_data::Stri return dem_surf, dem_bed, dem_grid, Ψ_grid end -# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem_on_GPU(backend, arch, synthetic_data); -surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem_on_GPU(backend, arch, vavilov_path); +# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem(backend, arch, synthetic_data); +surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem(backend, arch, vavilov_path); Ψ = ( na=Field(backend, surf_Ψ_grid, Vertex()), From e0d803c3bf696b2eda0016bdcb6a8e8c9cd50e67 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Tue, 5 Mar 2024 23:12:28 +0100 Subject: [PATCH 04/32] Update --- scripts_future_API/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts_future_API/Project.toml b/scripts_future_API/Project.toml index 9b70abea..57a5026b 100644 --- a/scripts_future_API/Project.toml +++ b/scripts_future_API/Project.toml @@ -3,7 +3,6 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" FastIce = "e0de9f13-a007-490e-b696-b07d031015ca" FastIceTools = "14532042-4700-46e0-8dec-55d517bc1ae0" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" From 9a9694cde251ebea663306804217fbf0121202a2 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Tue, 5 Mar 2024 23:52:38 +0100 Subject: [PATCH 05/32] Fixup --- .JuliaFormatter.toml | 2 +- scripts_future_API/test_levelsets.jl | 136 +++++++++--------- .../test_levelsets_volumefractions.jl | 15 +- src/LevelSets/LevelSets.jl | 5 +- src/LevelSets/compute_level_sets.jl | 5 +- src/LevelSets/compute_volume_fractions.jl | 3 +- src/LevelSets/signed_distances.jl | 32 ++--- src/LevelSets/volume_fractions_kernels.jl | 35 ++--- 8 files changed, 108 insertions(+), 125 deletions(-) diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 19024d0c..82423681 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -1,5 +1,5 @@ style = "yas" -margin = 140 +margin = 165 align_assignment = true align_conditional = true whitespace_ops_in_indices = false diff --git a/scripts_future_API/test_levelsets.jl b/scripts_future_API/test_levelsets.jl index 94f9833b..b1eca39e 100644 --- a/scripts_future_API/test_levelsets.jl +++ b/scripts_future_API/test_levelsets.jl @@ -1,78 +1,78 @@ -using CUDA -using FileIO -using FastIceTools -using JLD2 -using KernelAbstractions -using HDF5 +# using CUDA +# using FileIO +# using FastIceTools +# using JLD2 +# using KernelAbstractions +# using HDF5 -using FastIce.Grids -using FastIce.Fields -using FastIce.LevelSets -using FastIce.Architectures -using FastIce.Writers +# using FastIce.Grids +# using FastIce.Fields +# using FastIce.LevelSets +# using FastIce.Architectures +# using FastIce.Writers -vavilov_path = "../data/Vavilov/vavilov.jld2" -synthetic_data = "../data/synthetic.jld2" +# vavilov_path = "../data/Vavilov/vavilov.jld2" +# synthetic_data = "../data/synthetic.jld2" -# Select backend (CPU(), CUDABackend()) -backend = CUDABackend() -arch = Architecture(backend) +# # Select backend (CPU(), CUDABackend()) +# backend = CUDABackend() +# arch = Architecture(backend) -""" - load_dem_on_GPU(path::String, arch::Architecture) +# """ +# load_dem_on_GPU(path::String, arch::Architecture) -Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, -initiate grids, copy data on gpu. -""" -function load_dem_on_GPU(path::String, arch::Architecture) - data = load(path) - x = data["DataElevation"].x - y = data["DataElevation"].y - R = data["DataElevation"].rotation - z_surf = permutedims(data["DataElevation"].z_surf) - z_bed = permutedims(data["DataElevation"].z_bed) - domain = data["DataElevation"].domain - nx = length(x) - 1 - ny = length(y) - 1 - nz = 10 - z = LinRange(domain.zmin, domain.zmax, nz) - Ψ_grid = CartesianGrid(; origin=(0.0, 0.0, 0.0), extent=(1.0, 1.0, 1.0), size=(nx, ny, nz)) - Ψ = Field(arch, Ψ_grid, Vertex()) - dem_grid = CartesianGrid(; origin=(0.0, 0.0), extent=(1.0, 1.0), size=(nx, ny)) - dem_bed = Field(arch, dem_grid, Vertex()) - dem_surf = Field(arch, dem_grid, Vertex()) - set!(dem_bed, z_bed) - set!(dem_surf, z_surf) - return Ψ, dem_surf, dem_bed, dem_grid, Ψ_grid -end +# Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, +# initiate grids, copy data on gpu. +# """ +# function load_dem_on_GPU(path::String, arch::Architecture) +# data = load(path) +# x = data["DataElevation"].x +# y = data["DataElevation"].y +# R = data["DataElevation"].rotation +# z_surf = permutedims(data["DataElevation"].z_surf) +# z_bed = permutedims(data["DataElevation"].z_bed) +# domain = data["DataElevation"].domain +# nx = length(x) - 1 +# ny = length(y) - 1 +# nz = 10 +# z = LinRange(domain.zmin, domain.zmax, nz) +# Ψ_grid = CartesianGrid(; origin=(0.0, 0.0, 0.0), extent=(1.0, 1.0, 1.0), size=(nx, ny, nz)) +# Ψ = Field(arch, Ψ_grid, Vertex()) +# dem_grid = CartesianGrid(; origin=(0.0, 0.0), extent=(1.0, 1.0), size=(nx, ny)) +# dem_bed = Field(arch, dem_grid, Vertex()) +# dem_surf = Field(arch, dem_grid, Vertex()) +# set!(dem_bed, z_bed) +# set!(dem_surf, z_surf) +# return Ψ, dem_surf, dem_bed, dem_grid, Ψ_grid +# end -function load_synth_dem_on_GPU(path::String, arch::Architecture) - data = load(path) - z_surf = data["z_surf"] - z_bed = data["z_bed"] - nx = size(z_bed)[1] - 1 - ny = size(z_bed)[2] - 1 - nz = 100 - lz = maximum(z_bed) - minimum(z_bed) - Ψ_grid = CartesianGrid(; origin=(0.0, 0.0, minimum(z_bed)), extent=(lz, lz, lz), size=(nx, ny, nz)) - Ψ = Field(arch, Ψ_grid, Vertex()) - dem_grid = CartesianGrid(; origin=(0.0, 0.0), extent=(lz, lz), size=(nx, ny)) - dem_bed = Field(arch, dem_grid, Vertex()) - dem_surf = Field(arch, dem_grid, Vertex()) - set!(dem_bed, z_bed) - set!(dem_surf, z_surf) - return Ψ, dem_surf, dem_bed, dem_grid, Ψ_grid -end +# function load_synth_dem_on_GPU(path::String, arch::Architecture) +# data = load(path) +# z_surf = data["z_surf"] +# z_bed = data["z_bed"] +# nx = size(z_bed)[1] - 1 +# ny = size(z_bed)[2] - 1 +# nz = 100 +# lz = maximum(z_bed) - minimum(z_bed) +# Ψ_grid = CartesianGrid(; origin=(0.0, 0.0, minimum(z_bed)), extent=(lz, lz, lz), size=(nx, ny, nz)) +# Ψ = Field(arch, Ψ_grid, Vertex()) +# dem_grid = CartesianGrid(; origin=(0.0, 0.0), extent=(lz, lz), size=(nx, ny)) +# dem_bed = Field(arch, dem_grid, Vertex()) +# dem_surf = Field(arch, dem_grid, Vertex()) +# set!(dem_bed, z_bed) +# set!(dem_surf, z_surf) +# return Ψ, dem_surf, dem_bed, dem_grid, Ψ_grid +# end -# Ψ, dem_surf, dem_bed, dem_grid, Ψ_grid = load_dem_on_GPU(vavilov_path, arch); -Ψ, dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem_on_GPU(synthetic_data, arch); +# # Ψ, dem_surf, dem_bed, dem_grid, Ψ_grid = load_dem_on_GPU(vavilov_path, arch); +# Ψ, dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem_on_GPU(synthetic_data, arch); -compute_level_set_from_dem!(arch, Ψ, dem_bed, dem_grid, Ψ_grid) +# compute_level_set_from_dem!(arch, Ψ, dem_bed, dem_grid, Ψ_grid) -jldopen(synthetic_data, "a+") do file - file["level_sets"] = Array(interior(Ψ)) -end +# jldopen(synthetic_data, "a+") do file +# file["level_sets"] = Array(interior(Ψ)) +# end -jldopen(vavilov_path, "a+") do file - file["level_sets"] = Array(interior(Ψ)) -end +# jldopen(vavilov_path, "a+") do file +# file["level_sets"] = Array(interior(Ψ)) +# end diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl index 965fa1a7..03aa1901 100644 --- a/scripts_future_API/test_levelsets_volumefractions.jl +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -1,11 +1,14 @@ -using CUDA -using JLD2 -using FastIce.LevelSets +using Chmy.Architectures using Chmy.Grids using Chmy.Fields -using Chmy.Architectures + +using FastIce.LevelSets + using KernelAbstractions +using CUDA +using JLD2 + # Data # vavilov_path = "../data/vavilov.jld2" # vavilov_path = "../data/Vavilov/Vavilov_80m.jld2" @@ -21,7 +24,7 @@ arch = Arch(backend) """ load_dem(backend, arch::Architecture, path::String) -Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, +Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, initiate grids, copy data on gpu. """ function load_dem(backend, arch::Architecture, path::String) @@ -50,7 +53,7 @@ end """ load_synth_dem(backend, arch::Architecture, synthetic_data::String) -Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, +Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, initiate grids, copy data on gpu. """ function load_synth_dem(backend, arch::Architecture, synthetic_data::String) diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index 37ad3edf..0391eeed 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -3,10 +3,11 @@ module LevelSets export compute_level_set_from_dem! export volfrac_field, compute_volume_fractions_from_level_set! -using Chmy -using Chmy.Grids +# using Chmy using Chmy.Architectures +using Chmy.Grids using Chmy.Fields + using KernelAbstractions using LinearAlgebra, GeometryBasics diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index e486cedd..6a41fe64 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -1,5 +1,5 @@ """ - _init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff::AbstractFloat, R::AbstractMatrix) + init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R) Initialize level sets. """ @@ -11,7 +11,6 @@ Initialize level sets. Ψ[I...] = ud * sgn end - """ compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) @@ -22,4 +21,4 @@ function compute_level_set_from_dem!(backend, Ψ::Field, dem::Field, dem_grid::U cutoff = 4maximum(spacing(Ψ_grid, Center())) kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) return -end \ No newline at end of file +end diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl index 37c8dbb5..7ecbf270 100644 --- a/src/LevelSets/compute_volume_fractions.jl +++ b/src/LevelSets/compute_volume_fractions.jl @@ -139,4 +139,5 @@ include("volume_fractions_kernels.jl") function compute_volume_fractions_from_level_set!(backend::Backend, wt, Ψ, grid::UniformGrid) kernel = kernel_compute_volume_fractions_from_level_set!(backend, 256, size(Ψ)) kernel(wt, Ψ, grid) -end \ No newline at end of file + return +end diff --git a/src/LevelSets/signed_distances.jl b/src/LevelSets/signed_distances.jl index 36d97ee3..95befcd2 100644 --- a/src/LevelSets/signed_distances.jl +++ b/src/LevelSets/signed_distances.jl @@ -3,7 +3,6 @@ export sd_dem @inline S(x) = x == zero(x) ? oneunit(x) : sign(x) sign_triangle(p, a, b, c) = S(dot(p - a, cross(b - a, c - a))) - @inline function ud_triangle(p, a, b, c) dot2(v) = dot(v, v) ba = b - a @@ -13,20 +12,17 @@ sign_triangle(p, a, b, c) = S(dot(p - a, cross(b - a, c - a))) ac = a - c pc = p - c nor = cross(ba, ac) - return sqrt( - (sign(dot(cross(ba, nor), pa)) + - sign(dot(cross(cb, nor), pb)) + - sign(dot(cross(ac, nor), pc)) < 2) - ? - min( - dot2(ba * clamp(dot(ba, pa) / dot2(ba), 0, 1) - pa), - dot2(cb * clamp(dot(cb, pb) / dot2(cb), 0, 1) - pb), - dot2(ac * clamp(dot(ac, pc) / dot2(ac), 0, 1) - pc)) - : - dot(nor, pa) * dot(nor, pa) / dot2(nor)) + return sqrt((sign(dot(cross(ba, nor), pa)) + + sign(dot(cross(cb, nor), pb)) + + sign(dot(cross(ac, nor), pc)) < 2) + ? + min(dot2(ba * clamp(dot(ba, pa) / dot2(ba), 0, 1) - pa), + dot2(cb * clamp(dot(cb, pb) / dot2(cb), 0, 1) - pb), + dot2(ac * clamp(dot(ac, pc) / dot2(ac), 0, 1) - pc)) + : + dot(nor, pa) * dot(nor, pa) / dot2(nor)) end - @inline function closest_vertex_index(P, grid) Δs = spacing(grid, Vertex(), 1, 1) O = Grids.origin(grid, Vertex()) @@ -37,11 +33,9 @@ end return clamp.(I, I1, I2) |> CartesianIndex end - @inline inc(I, dim) = Base.setindex(I, I[dim] + 1, dim) @inline inc(I) = I + oneunit(I) - @inline function triangle_pair(Iv, dem, grid) @inline function sample_dem(I) @inbounds x, y = coord(grid, location(dem), I) @@ -52,16 +46,15 @@ end return T_BL, T_TR end - @inline function distance_to_triangle_pair(P, Iv, dem, grid) T_BL, T_TR = triangle_pair(Iv, dem, grid) ud = min(ud_triangle(P, T_BL...), ud_triangle(P, T_TR...)) return ud, sign_triangle(P, T_BL...) end - -Base.clamp(p::NTuple{N}, grid::UniformGrid{N}) where {N} = clamp.(p, Grids.origin(grid, Vertex()), Grids.origin(grid, Vertex()) .+ extent(grid, Vertex())) - +function Base.clamp(p::NTuple{N}, grid::UniformGrid{N}) where {N} + clamp.(p, Grids.origin(grid, Vertex()), Grids.origin(grid, Vertex()) .+ extent(grid, Vertex())) +end function sd_dem(P, cutoff, dem, grid) @inbounds Pp = clamp((P[1], P[2]), grid) @@ -78,4 +71,3 @@ function sd_dem(P, cutoff, dem, grid) end return ud, sgn end - diff --git a/src/LevelSets/volume_fractions_kernels.jl b/src/LevelSets/volume_fractions_kernels.jl index dd16cbf9..0b197e69 100644 --- a/src/LevelSets/volume_fractions_kernels.jl +++ b/src/LevelSets/volume_fractions_kernels.jl @@ -1,7 +1,7 @@ # Volume fraction kernel (2D) -@kernel function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid) +@kernel inbounds = true function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid{2}) ix, iy = @index(Global, NTuple) - dx, dy = spacing(grid, Center(), 1, 1) + dx, dy = spacing(grid, Center()) cell = Rect(Vec(0.0, 0.0), Vec(dx, dy)) ω = GeometryBasics.volume(cell) Ψ_ax(dix, diy) = 0.5 * (Ψ[ix+dix, iy+diy] + Ψ[ix+dix+1, iy+diy]) @@ -22,29 +22,19 @@ end # Volume fraction kernel (3D) -@kernel function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid) +@kernel inbounds = true function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid{3}) ix, iy, iz = @index(Global, NTuple) - dx, dy, dz = spacing(grid, Center(), 1, 1, 1) + dx, dy, dz = spacing(grid, Center()) cell = Rect(Vec(0.0, 0.0, 0.0), Vec(dx, dy, dz)) ω = GeometryBasics.volume(cell) Ψ_ax(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix+1, iy+diy, iz+diz]) Ψ_ay(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix, iy+diy+1, iz+diz]) Ψ_az(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix, iy+diy, iz+diz+1]) - function Ψ_axy(dix, diy, diz) - 0.25 * - (Ψ[ix+dix, iy+diy, iz+diz+1] + Ψ[ix+dix+1, iy+diy, iz+diz+1] + Ψ[ix+dix, iy+diy+1, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) - end - function Ψ_axz(dix, diy, diz) - 0.25 * - (Ψ[ix+dix, iy+diy+1, iz+diz] + Ψ[ix+dix+1, iy+diy+1, iz+diz] + Ψ[ix+dix, iy+diy+1, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) - end - function Ψ_ayz(dix, diy, diz) - 0.25 * - (Ψ[ix+dix+1, iy+diy, iz+diz] + Ψ[ix+dix+1, iy+diy+1, iz+diz] + Ψ[ix+dix+1, iy+diy, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) - end + Ψ_axy(dix, diy, diz) = 0.25 * (Ψ[ix+dix, iy+diy, iz+diz+1] + Ψ[ix+dix+1, iy+diy, iz+diz+1] + Ψ[ix+dix, iy+diy+1, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) + Ψ_axz(dix, diy, diz) = 0.25 * (Ψ[ix+dix, iy+diy+1, iz+diz] + Ψ[ix+dix+1, iy+diy+1, iz+diz] + Ψ[ix+dix, iy+diy+1, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) + Ψ_ayz(dix, diy, diz) = 0.25 * (Ψ[ix+dix+1, iy+diy, iz+diz] + Ψ[ix+dix+1, iy+diy+1, iz+diz] + Ψ[ix+dix+1, iy+diy, iz+diz+1] + Ψ[ix+dix+1, iy+diy+1, iz+diz+1]) # cell centers - Ψs = Vec{8}(Ψ[ix, iy, iz], Ψ[ix+1, iy, iz], Ψ[ix, iy+1, iz], Ψ[ix+1, iy+1, iz], Ψ[ix, iy, iz+1], Ψ[ix+1, iy, iz+1], Ψ[ix, iy+1, iz+1], - Ψ[ix+1, iy+1, iz+1]) + Ψs = Vec{8}(Ψ[ix, iy, iz], Ψ[ix+1, iy, iz], Ψ[ix, iy+1, iz], Ψ[ix+1, iy+1, iz], Ψ[ix, iy, iz+1], Ψ[ix+1, iy, iz+1], Ψ[ix, iy+1, iz+1], Ψ[ix+1, iy+1, iz+1]) wt.c[ix, iy, iz] = volfrac(cell, Ψs) / ω # x faces Ψs = Vec{8}(Ψ_ax(0, 0, 0), Ψ_ax(1, 0, 0), Ψ_ax(0, 1, 0), Ψ_ax(1, 1, 0), Ψ_ax(0, 0, 1), Ψ_ax(1, 0, 1), Ψ_ax(0, 1, 1), Ψ_ax(1, 1, 1)) @@ -56,15 +46,12 @@ end Ψs = Vec{8}(Ψ_az(0, 0, 0), Ψ_az(1, 0, 0), Ψ_az(0, 1, 0), Ψ_az(1, 1, 0), Ψ_az(0, 0, 1), Ψ_az(1, 0, 1), Ψ_az(0, 1, 1), Ψ_az(1, 1, 1)) wt.z[ix, iy, iz] = volfrac(cell, Ψs) / ω # xy edges - Ψs = Vec{8}(Ψ_axy(0, 0, 0), Ψ_axy(1, 0, 0), Ψ_axy(0, 1, 0), Ψ_axy(1, 1, 0), Ψ_axy(0, 0, 1), Ψ_axy(1, 0, 1), Ψ_axy(0, 1, 1), - Ψ_axy(1, 1, 1)) + Ψs = Vec{8}(Ψ_axy(0, 0, 0), Ψ_axy(1, 0, 0), Ψ_axy(0, 1, 0), Ψ_axy(1, 1, 0), Ψ_axy(0, 0, 1), Ψ_axy(1, 0, 1), Ψ_axy(0, 1, 1), Ψ_axy(1, 1, 1)) wt.xy[ix, iy, iz] = volfrac(cell, Ψs) / ω # xz edges - Ψs = Vec{8}(Ψ_axz(0, 0, 0), Ψ_axz(1, 0, 0), Ψ_axz(0, 1, 0), Ψ_axz(1, 1, 0), Ψ_axz(0, 0, 1), Ψ_axz(1, 0, 1), Ψ_axz(0, 1, 1), - Ψ_axz(1, 1, 1)) + Ψs = Vec{8}(Ψ_axz(0, 0, 0), Ψ_axz(1, 0, 0), Ψ_axz(0, 1, 0), Ψ_axz(1, 1, 0), Ψ_axz(0, 0, 1), Ψ_axz(1, 0, 1), Ψ_axz(0, 1, 1), Ψ_axz(1, 1, 1)) wt.xz[ix, iy, iz] = volfrac(cell, Ψs) / ω # yz edges - Ψs = Vec{8}(Ψ_ayz(0, 0, 0), Ψ_ayz(1, 0, 0), Ψ_ayz(0, 1, 0), Ψ_ayz(1, 1, 0), Ψ_ayz(0, 0, 1), Ψ_ayz(1, 0, 1), Ψ_ayz(0, 1, 1), - Ψ_ayz(1, 1, 1)) + Ψs = Vec{8}(Ψ_ayz(0, 0, 0), Ψ_ayz(1, 0, 0), Ψ_ayz(0, 1, 0), Ψ_ayz(1, 1, 0), Ψ_ayz(0, 0, 1), Ψ_ayz(1, 0, 1), Ψ_ayz(0, 1, 1), Ψ_ayz(1, 1, 1)) wt.yz[ix, iy, iz] = volfrac(cell, Ψs) / ω end From 739c201ac64524c3b1296e6556864977dd59ab80 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 6 Mar 2024 00:09:10 +0100 Subject: [PATCH 06/32] Fixup --- src/LevelSets/LevelSets.jl | 3 +-- src/LevelSets/compute_level_sets.jl | 5 +++-- src/LevelSets/compute_volume_fractions.jl | 10 ++++++++-- src/LevelSets/volume_fractions_kernels.jl | 4 ++-- 4 files changed, 14 insertions(+), 8 deletions(-) diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index 0391eeed..a4ad2ba8 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -1,9 +1,8 @@ module LevelSets export compute_level_set_from_dem! -export volfrac_field, compute_volume_fractions_from_level_set! +export volfrac_field, compute_volfrac_from_level_set! -# using Chmy using Chmy.Architectures using Chmy.Grids using Chmy.Fields diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index 6a41fe64..985d9c08 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -3,7 +3,7 @@ Initialize level sets. """ -@kernel function init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R) +@kernel inbounds = true function init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R) I = @index(Global, NTuple) x, y, z = coord(Ψ_grid, location(Ψ), I...) P = R * Point3(x, y, z) @@ -16,7 +16,8 @@ end Compute level sets from dem. """ -function compute_level_set_from_dem!(backend, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) +function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) + backend = Architectures.get_backend(arch) kernel = init_level_set!(backend, 256, size(Ψ)) cutoff = 4maximum(spacing(Ψ_grid, Center())) kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl index 7ecbf270..e9f2f4e5 100644 --- a/src/LevelSets/compute_volume_fractions.jl +++ b/src/LevelSets/compute_volume_fractions.jl @@ -136,8 +136,14 @@ end include("volume_fractions_kernels.jl") -function compute_volume_fractions_from_level_set!(backend::Backend, wt, Ψ, grid::UniformGrid) - kernel = kernel_compute_volume_fractions_from_level_set!(backend, 256, size(Ψ)) +""" + compute_volfrac_from_level_set!(arch::Architecture, wt::Field, Ψ::Field, grid::UniformGrid) + +Compute volume fractions from level sets. +""" +function compute_volfrac_from_level_set!(arch::Architecture, wt::Field, Ψ::Field, grid::UniformGrid) + backend = Architectures.get_backend(arch) + kernel = compute_volfrac_from_level_set!(backend, 256, size(Ψ)) kernel(wt, Ψ, grid) return end diff --git a/src/LevelSets/volume_fractions_kernels.jl b/src/LevelSets/volume_fractions_kernels.jl index 0b197e69..86fb0ce3 100644 --- a/src/LevelSets/volume_fractions_kernels.jl +++ b/src/LevelSets/volume_fractions_kernels.jl @@ -1,5 +1,5 @@ # Volume fraction kernel (2D) -@kernel inbounds = true function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid{2}) +@kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{2}) ix, iy = @index(Global, NTuple) dx, dy = spacing(grid, Center()) cell = Rect(Vec(0.0, 0.0), Vec(dx, dy)) @@ -22,7 +22,7 @@ end # Volume fraction kernel (3D) -@kernel inbounds = true function kernel_compute_volume_fractions_from_level_set!(wt, Ψ, grid::UniformGrid{3}) +@kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{3}) ix, iy, iz = @index(Global, NTuple) dx, dy, dz = spacing(grid, Center()) cell = Rect(Vec(0.0, 0.0, 0.0), Vec(dx, dy, dz)) From c22f2284c65c63604748fe4c425d2dfff30b77ae Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 6 Mar 2024 00:21:17 +0100 Subject: [PATCH 07/32] More fixup --- .../test_levelsets_volumefractions.jl | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl index 03aa1901..32ee8144 100644 --- a/scripts_future_API/test_levelsets_volumefractions.jl +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -22,12 +22,12 @@ backend = CPU() arch = Arch(backend) """ - load_dem(backend, arch::Architecture, path::String) + load_dem(arch::Architecture, path::String) -Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, -initiate grids, copy data on gpu. +Load digital elevation model of surface and bedrock from (JLD2) file and initialise the grid. """ -function load_dem(backend, arch::Architecture, path::String) +function load_dem(arch::Architecture, path::String) + backend = arch.backend data = load(path) z_surf = data["DataElevation"].z_surf z_bed = data["DataElevation"].z_bed @@ -51,12 +51,12 @@ function load_dem(backend, arch::Architecture, path::String) end """ - load_synth_dem(backend, arch::Architecture, synthetic_data::String) + load_synth_dem(arch::Architecture, synthetic_data::String) -Load digital elevation map of surface and bedrock from (jld2) file, set dimensions of simulation, -initiate grids, copy data on gpu. +Load synthetic elevation of surface and bedrock from (JLD2) file and initialise the grid. """ -function load_synth_dem(backend, arch::Architecture, synthetic_data::String) +function load_synth_dem(arch::Architecture, synthetic_data::String) + backend = arch.backend data = load(synthetic_data) z_surf = data["SyntheticElevation"].z_surf z_bed = data["SyntheticElevation"].z_bed @@ -73,19 +73,19 @@ function load_synth_dem(backend, arch::Architecture, synthetic_data::String) return dem_surf, dem_bed, dem_grid, Ψ_grid end -# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem(backend, arch, synthetic_data); -surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem(backend, arch, vavilov_path); +# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem(arch, synthetic_data); +surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem(arch, vavilov_path); Ψ = ( na=Field(backend, surf_Ψ_grid, Vertex()), ns=Field(backend, bed_Ψ_grid, Vertex()), ) -compute_level_set_from_dem!(backend, Ψ[1], surf_field, surf_dem_grid, surf_Ψ_grid) -compute_level_set_from_dem!(backend, Ψ[2], bed_field, bed_dem_grid, bed_Ψ_grid) +compute_level_set_from_dem!(arch, Ψ[1], surf_field, surf_dem_grid, surf_Ψ_grid) +compute_level_set_from_dem!(arch, Ψ[2], bed_field, bed_dem_grid, bed_Ψ_grid) # for phase in eachindex(Ψ) -# compute_level_set_from_dem!(backend, Ψ[phase], dem_surf, dem_grid, Ψ_grid) +# compute_level_set_from_dem!(arch, Ψ[phase], dem_surf, dem_grid, Ψ_grid) # end wt = ( @@ -93,12 +93,12 @@ wt = ( ns=volfrac_field(backend, bed_Ψ_grid), ) -compute_volume_fractions_from_level_set!(backend, wt[1], Ψ[1], surf_Ψ_grid) -compute_volume_fractions_from_level_set!(backend, wt[2], Ψ[2], bed_Ψ_grid) +compute_volfrac_from_level_set!(arch, wt[1], Ψ[1], surf_Ψ_grid) +compute_volfrac_from_level_set!(arch, wt[2], Ψ[2], bed_Ψ_grid) # for phase in eachindex(Ψ) -# compute_volume_fractions_from_level_set!(backend, wt[phase], Ψ[phase], Ψ_grid) +# compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], Ψ_grid) # end # Save From 7e52a2a78d6db09ce173fee3883e75b436641445 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 6 Mar 2024 09:31:20 +0100 Subject: [PATCH 08/32] Fixup --- src/LevelSets/compute_level_sets.jl | 2 +- src/LevelSets/compute_volume_fractions.jl | 8 +++++--- src/LevelSets/volume_fractions_kernels.jl | 4 ++-- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index 985d9c08..b9b3244f 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -19,7 +19,7 @@ Compute level sets from dem. function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) backend = Architectures.get_backend(arch) kernel = init_level_set!(backend, 256, size(Ψ)) - cutoff = 4maximum(spacing(Ψ_grid, Center())) + cutoff = 4maximum(spacing(Ψ_grid)) kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) return end diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl index e9f2f4e5..ce8ee647 100644 --- a/src/LevelSets/compute_volume_fractions.jl +++ b/src/LevelSets/compute_volume_fractions.jl @@ -46,7 +46,7 @@ end @inline tetvol(v1, v2, v3, v4) = abs(det([v2 - v1 v3 - v1 v4 - v1])) / 6.0 -function volfrac(tet, ϕ::Vec4) +function volfrac(tet, ϕ::Vec4{T})::T where {T} v1, v2, v3, v4 = tet @inline vij(i, j) = tet[j] * (ϕ[i] / (ϕ[i] - ϕ[j])) - tet[i] * (ϕ[j] / (ϕ[i] - ϕ[j])) nneg = count(ϕ .< 0) @@ -137,13 +137,15 @@ end include("volume_fractions_kernels.jl") """ - compute_volfrac_from_level_set!(arch::Architecture, wt::Field, Ψ::Field, grid::UniformGrid) + compute_volfrac_from_level_set!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) Compute volume fractions from level sets. """ -function compute_volfrac_from_level_set!(arch::Architecture, wt::Field, Ψ::Field, grid::UniformGrid) +function compute_volfrac_from_level_set!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) backend = Architectures.get_backend(arch) kernel = compute_volfrac_from_level_set!(backend, 256, size(Ψ)) kernel(wt, Ψ, grid) + # BC Neumann for x,y in 2D + # BC Neumann for x,y,z in 3D return end diff --git a/src/LevelSets/volume_fractions_kernels.jl b/src/LevelSets/volume_fractions_kernels.jl index 86fb0ce3..d6a464aa 100644 --- a/src/LevelSets/volume_fractions_kernels.jl +++ b/src/LevelSets/volume_fractions_kernels.jl @@ -1,7 +1,7 @@ # Volume fraction kernel (2D) @kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{2}) ix, iy = @index(Global, NTuple) - dx, dy = spacing(grid, Center()) + dx, dy = spacing(grid) cell = Rect(Vec(0.0, 0.0), Vec(dx, dy)) ω = GeometryBasics.volume(cell) Ψ_ax(dix, diy) = 0.5 * (Ψ[ix+dix, iy+diy] + Ψ[ix+dix+1, iy+diy]) @@ -24,7 +24,7 @@ end # Volume fraction kernel (3D) @kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{3}) ix, iy, iz = @index(Global, NTuple) - dx, dy, dz = spacing(grid, Center()) + dx, dy, dz = spacing(grid) cell = Rect(Vec(0.0, 0.0, 0.0), Vec(dx, dy, dz)) ω = GeometryBasics.volume(cell) Ψ_ax(dix, diy, diz) = 0.5 * (Ψ[ix+dix, iy+diy, iz+diz] + Ψ[ix+dix+1, iy+diy, iz+diz]) From d6508348e8eb0fe253b45d0690d6ebcaf62d00e5 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 6 Mar 2024 14:23:29 +0100 Subject: [PATCH 09/32] Export Geometry for FastIceTools --- src/FastIce.jl | 5 +++-- src/Geometry.jl | 27 +++++++++++++++++++++++++++ 2 files changed, 30 insertions(+), 2 deletions(-) create mode 100644 src/Geometry.jl diff --git a/src/FastIce.jl b/src/FastIce.jl index b8408213..69537a2b 100644 --- a/src/FastIce.jl +++ b/src/FastIce.jl @@ -25,11 +25,12 @@ https://github.com/PTsolvers/FastIce.jl greet(; kwargs...) = printstyled(GREETING; kwargs...) greet_fast(; kwargs...) = printstyled(GREETING_FAST; kwargs...) -# core modules +# core modules (included in alphabetical order) +include("Geometry.jl") +include("LevelSets/LevelSets.jl") include("Logging.jl") include("Physics.jl") include("Writers.jl") -include("LevelSets/LevelSets.jl") # ice flow models # include("Models/Models.jl") diff --git a/src/Geometry.jl b/src/Geometry.jl new file mode 100644 index 00000000..ecff2cab --- /dev/null +++ b/src/Geometry.jl @@ -0,0 +1,27 @@ +module Geometry + +export AbstractElevation, DataElevation, AABB + +abstract type AbstractElevation{T<:Real} end + +struct DataElevation{T,V<:AbstractArray{T},R<:AbstractRange{T},M<:AbstractMatrix{T}} <: AbstractElevation{T} + x::R + y::R + offsets::V + z_bed::M + z_surf::M + rotation::M + domain::AABB{T} + rotated_domain::AABB{T} +end + +struct AABB{T<:Union{Real}} + xmin::T + xmax::T + ymin::T + ymax::T + zmin::T + zmax::T +end + +end From 4cc66370a99bafc36b8478e40179324564767ae9 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 6 Mar 2024 14:23:44 +0100 Subject: [PATCH 10/32] Use launcher --- src/LevelSets/LevelSets.jl | 1 + src/LevelSets/compute_level_sets.jl | 11 +++++------ src/LevelSets/compute_volume_fractions.jl | 7 ++----- 3 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index a4ad2ba8..255f34ee 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -6,6 +6,7 @@ export volfrac_field, compute_volfrac_from_level_set! using Chmy.Architectures using Chmy.Grids using Chmy.Fields +using Chmy.KernelLaunch using KernelAbstractions using LinearAlgebra, GeometryBasics diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index b9b3244f..f92c77a4 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -12,14 +12,13 @@ Initialize level sets. end """ - compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) + compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I) Compute level sets from dem. """ -function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, R=LinearAlgebra.I) - backend = Architectures.get_backend(arch) - kernel = init_level_set!(backend, 256, size(Ψ)) - cutoff = 4maximum(spacing(Ψ_grid)) - kernel(Ψ, dem, dem_grid, Ψ_grid, cutoff, R) +function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::Field, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I; outer_width=(16, 8, 4)) + cutoff = 4maximum(spacing(grid)) + launch = Launcher(arch, grid; outer_width=outer_width) + launch(arch, grid, init_level_set! => (Ψ, dem, dem_grid2D, grid, cutoff, R); bc=batch(grid, Ψ => Neumann(); exchange=Ψ)) return end diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl index ce8ee647..bd83e7ba 100644 --- a/src/LevelSets/compute_volume_fractions.jl +++ b/src/LevelSets/compute_volume_fractions.jl @@ -142,10 +142,7 @@ include("volume_fractions_kernels.jl") Compute volume fractions from level sets. """ function compute_volfrac_from_level_set!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) - backend = Architectures.get_backend(arch) - kernel = compute_volfrac_from_level_set!(backend, 256, size(Ψ)) - kernel(wt, Ψ, grid) - # BC Neumann for x,y in 2D - # BC Neumann for x,y,z in 3D + launch = Launcher(arch, grid) + launch(arch, grid, compute_volfrac_from_level_set! => (wt, Ψ, grid)) return end From 5212609966c4c3c79f4f4be019b1f52aa0b24c11 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 6 Mar 2024 14:31:22 +0100 Subject: [PATCH 11/32] Update --- src/Geometry.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/Geometry.jl b/src/Geometry.jl index ecff2cab..8fa17000 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -4,6 +4,15 @@ export AbstractElevation, DataElevation, AABB abstract type AbstractElevation{T<:Real} end +struct AABB{T<:Union{Real}} + xmin::T + xmax::T + ymin::T + ymax::T + zmin::T + zmax::T +end + struct DataElevation{T,V<:AbstractArray{T},R<:AbstractRange{T},M<:AbstractMatrix{T}} <: AbstractElevation{T} x::R y::R @@ -15,13 +24,4 @@ struct DataElevation{T,V<:AbstractArray{T},R<:AbstractRange{T},M<:AbstractMatrix rotated_domain::AABB{T} end -struct AABB{T<:Union{Real}} - xmin::T - xmax::T - ymin::T - ymax::T - zmin::T - zmax::T -end - end From 60ee0fecb7a14b3305e7a315c59f8db12d69eced Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 6 Mar 2024 14:55:28 +0100 Subject: [PATCH 12/32] Lower deps version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index afc47649..4ed6522a 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Adapt = "4" ElasticArrays = "1" GeometryBasics = "0.4" -HDF5 = "0.17" +HDF5 = "0.16, 0.17" KernelAbstractions = "0.9" LightXML = "0.9" MPI = "0.20" From 8335ff86e247de717dd55de93a679de278eacc2d Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 6 Mar 2024 15:53:59 +0100 Subject: [PATCH 13/32] Fixup and update --- scripts_future_API/test_levelsets_volumefractions.jl | 10 +++++++--- src/LevelSets/LevelSets.jl | 1 + src/LevelSets/compute_level_sets.jl | 3 ++- src/LevelSets/volume_fractions_kernels.jl | 10 ++++++---- 4 files changed, 16 insertions(+), 8 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl index 32ee8144..25f4fb75 100644 --- a/scripts_future_API/test_levelsets_volumefractions.jl +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -3,22 +3,25 @@ using Chmy.Grids using Chmy.Fields using FastIce.LevelSets +using FastIce.Geometry using KernelAbstractions - -using CUDA using JLD2 +# using CUDA +# using AMDGPU + # Data # vavilov_path = "../data/vavilov.jld2" # vavilov_path = "../data/Vavilov/Vavilov_80m.jld2" -vavilov_path = "../data/Vavilov/Vavilov_500m.jld2" +vavilov_path = "./Vavilov_500m.jld2" # synthetic_data = "../data/Synthetic/dome_glacier.jld2" # synthetic_data = "../data/Synthetic/low_res_dome_glacier.jld2" # Select backend backend = CPU() # backend = CUDABackend() +# backend = ROCBackend() arch = Arch(backend) """ @@ -102,6 +105,7 @@ compute_volfrac_from_level_set!(arch, wt[2], Ψ[2], bed_Ψ_grid) # end # Save +# save_path = "Vavilov_500m_out.jld2" jldopen(vavilov_path, "a+") do file file["level_sets_na"] = Array(interior(Ψ[1])) file["level_sets_ns"] = Array(interior(Ψ[2])) diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index 255f34ee..7117fc4a 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -4,6 +4,7 @@ export compute_level_set_from_dem! export volfrac_field, compute_volfrac_from_level_set! using Chmy.Architectures +using Chmy.BoundaryConditions using Chmy.Grids using Chmy.Fields using Chmy.KernelLaunch diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl index f92c77a4..11130d98 100644 --- a/src/LevelSets/compute_level_sets.jl +++ b/src/LevelSets/compute_level_sets.jl @@ -3,8 +3,9 @@ Initialize level sets. """ -@kernel inbounds = true function init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R) +@kernel inbounds = true function init_level_set!(Ψ::Field, dem::Field, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R, O) I = @index(Global, NTuple) + I = I + O x, y, z = coord(Ψ_grid, location(Ψ), I...) P = R * Point3(x, y, z) ud, sgn = sd_dem(P, cutoff, dem, dem_grid) diff --git a/src/LevelSets/volume_fractions_kernels.jl b/src/LevelSets/volume_fractions_kernels.jl index d6a464aa..472f8f0b 100644 --- a/src/LevelSets/volume_fractions_kernels.jl +++ b/src/LevelSets/volume_fractions_kernels.jl @@ -1,6 +1,7 @@ # Volume fraction kernel (2D) -@kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{2}) - ix, iy = @index(Global, NTuple) +@kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{2}, O) + I = @index(Global, NTuple) + ix, iy = I + O dx, dy = spacing(grid) cell = Rect(Vec(0.0, 0.0), Vec(dx, dy)) ω = GeometryBasics.volume(cell) @@ -22,8 +23,9 @@ end # Volume fraction kernel (3D) -@kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{3}) - ix, iy, iz = @index(Global, NTuple) +@kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{3}, O) + I = @index(Global, NTuple) + ix, iy, iz = I + O dx, dy, dz = spacing(grid) cell = Rect(Vec(0.0, 0.0, 0.0), Vec(dx, dy, dz)) ω = GeometryBasics.volume(cell) From 0d5cd5b4b005e1d34a577740ae2c9ce351a3e3d0 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Thu, 7 Mar 2024 11:42:43 +0100 Subject: [PATCH 14/32] Update Geometry --- src/Geometry.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/Geometry.jl b/src/Geometry.jl index 8fa17000..590b8f4d 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -1,6 +1,6 @@ module Geometry -export AbstractElevation, DataElevation, AABB +export AbstractElevation, DataElevation, SyntheticElevation, AABB abstract type AbstractElevation{T<:Real} end @@ -24,4 +24,10 @@ struct DataElevation{T,V<:AbstractArray{T},R<:AbstractRange{T},M<:AbstractMatrix rotated_domain::AABB{T} end +struct SyntheticElevation{T,B,S} <: AbstractElevation{T} + z_bed::B + z_surf::S + domain::AABB{T} +end + end From 03398a95d495950c9ce7fbd12fc57066b13959ff Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Thu, 7 Mar 2024 23:00:54 +0100 Subject: [PATCH 15/32] Cleanup --- .../test_levelsets_volumefractions.jl | 101 ++++++------------ 1 file changed, 33 insertions(+), 68 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl index 25f4fb75..db2e9fbe 100644 --- a/scripts_future_API/test_levelsets_volumefractions.jl +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -12,100 +12,65 @@ using JLD2 # using AMDGPU # Data -# vavilov_path = "../data/vavilov.jld2" -# vavilov_path = "../data/Vavilov/Vavilov_80m.jld2" -vavilov_path = "./Vavilov_500m.jld2" -# synthetic_data = "../data/Synthetic/dome_glacier.jld2" -# synthetic_data = "../data/Synthetic/low_res_dome_glacier.jld2" +data_path = "./low_res_dome_rough_glacier.jld2" # Select backend backend = CPU() # backend = CUDABackend() # backend = ROCBackend() -arch = Arch(backend) -""" - load_dem(arch::Architecture, path::String) +arch = Arch(backend) -Load digital elevation model of surface and bedrock from (JLD2) file and initialise the grid. -""" -function load_dem(arch::Architecture, path::String) - backend = arch.backend - data = load(path) - z_surf = data["DataElevation"].z_surf - z_bed = data["DataElevation"].z_bed - x = data["DataElevation"].x - y = data["DataElevation"].y - nx = length(x) - 1 - ny = length(y) - 1 - nz = 100 - # TODO: choose nz accordingly - surf_lz = maximum(z_surf) - minimum(z_surf) - bed_lz = maximum(z_bed) - minimum(z_bed) - surf_Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_surf)), extent=(surf_lz, surf_lz, surf_lz), dims=(nx, ny, nz)) - bed_Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_bed)), extent=(bed_lz, bed_lz, bed_lz), dims=(nx, ny, nz)) - surf_dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(surf_lz, surf_lz), dims=(nx, ny)) - bed_dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(bed_lz, bed_lz), dims=(nx, ny)) - surf_field = Field(backend, surf_dem_grid, Vertex()) - bed_field = Field(backend, bed_dem_grid, Vertex()) - set!(surf_field, z_surf) - set!(bed_field, z_bed) - return surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid +function load_data(data_path) + dat = load(data_path) + return first(keys(dat)), dat end + """ - load_synth_dem(arch::Architecture, synthetic_data::String) + extract_dem(arch::Architecture, data_path::String, z_resol=nothing) -Load synthetic elevation of surface and bedrock from (JLD2) file and initialise the grid. +Load digital elevation model of surface and bedrock from (JLD2) file and initialise the grid. """ -function load_synth_dem(arch::Architecture, synthetic_data::String) - backend = arch.backend - data = load(synthetic_data) - z_surf = data["SyntheticElevation"].z_surf - z_bed = data["SyntheticElevation"].z_bed - nx = size(z_surf)[1] - 1 - ny = size(z_surf)[2] - 1 - nz = 10 - lz = maximum(z_surf) - minimum(z_surf) - Ψ_grid = UniformGrid(arch; origin=(0.0, 0.0, minimum(z_surf)), extent=(lz, lz, lz), dims=(nx, ny, nz)) - dem_grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(lz, lz), dims=(nx, ny)) - dem_bed = Field(backend, dem_grid, Vertex()) - dem_surf = Field(backend, dem_grid, Vertex()) - set!(dem_bed, z_bed) +function extract_dem(arch::Architecture, data_path::String, z_resol=nothing) + dtype, data = load_data(data_path) + z_surf = data[dtype].z_surf + z_bed = data[dtype].z_bed + dm = data[dtype].domain + lx, ly, lz = dm.xmax - dm.xmin, dm.ymax - dm.ymin, dm.zmax - dm.zmin + nx, ny = size(z_surf) .- 1 + nz = isnothing(z_resol) ? ceil(Int, lz / lx * nx) : z_resol + Ψ_grid = UniformGrid(arch; origin=(dm.xmin, dm.ymin, dm.zmin), extent=(lx, ly, lz), dims=(nx, ny, nz)) + dem_grid = UniformGrid(arch; origin=(dm.xmin, dm.ymin), extent=(lx, ly), dims=(nx, ny)) + dem_surf = Field(arch.backend, dem_grid, Vertex()) + dem_bed = Field(arch.backend, dem_grid, Vertex()) set!(dem_surf, z_surf) + set!(dem_bed, z_bed) return dem_surf, dem_bed, dem_grid, Ψ_grid end -# dem_surf, dem_bed, dem_grid, Ψ_grid = load_synth_dem(arch, synthetic_data); -surf_field, bed_field, surf_dem_grid, bed_dem_grid, surf_Ψ_grid, bed_Ψ_grid = load_dem(arch, vavilov_path); +dem_surf, dem_bed, dem_grid, Ψ_grid = extract_dem(arch, data_path) Ψ = ( - na=Field(backend, surf_Ψ_grid, Vertex()), - ns=Field(backend, bed_Ψ_grid, Vertex()), + na=Field(backend, Ψ_grid, Vertex()), + ns=Field(backend, Ψ_grid, Vertex()), ) -compute_level_set_from_dem!(arch, Ψ[1], surf_field, surf_dem_grid, surf_Ψ_grid) -compute_level_set_from_dem!(arch, Ψ[2], bed_field, bed_dem_grid, bed_Ψ_grid) - -# for phase in eachindex(Ψ) -# compute_level_set_from_dem!(arch, Ψ[phase], dem_surf, dem_grid, Ψ_grid) -# end +for phase in eachindex(Ψ) + compute_level_set_from_dem!(arch, Ψ[phase], dem_surf, dem_grid, Ψ_grid) +end wt = ( - na=volfrac_field(backend, surf_Ψ_grid), - ns=volfrac_field(backend, bed_Ψ_grid), + na=volfrac_field(backend, Ψ_grid), + ns=volfrac_field(backend, Ψ_grid), ) -compute_volfrac_from_level_set!(arch, wt[1], Ψ[1], surf_Ψ_grid) -compute_volfrac_from_level_set!(arch, wt[2], Ψ[2], bed_Ψ_grid) - - -# for phase in eachindex(Ψ) -# compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], Ψ_grid) -# end +for phase in eachindex(Ψ) + compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], Ψ_grid) +end # Save -# save_path = "Vavilov_500m_out.jld2" +save_path = "Vavilov_500m_out.jld2" jldopen(vavilov_path, "a+") do file file["level_sets_na"] = Array(interior(Ψ[1])) file["level_sets_ns"] = Array(interior(Ψ[2])) From 97407236e3c3edde76530876fa2caa1d986e72bb Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Fri, 8 Mar 2024 12:06:20 +0100 Subject: [PATCH 16/32] Add geom --- src/Geometry.jl | 33 ++++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/src/Geometry.jl b/src/Geometry.jl index 590b8f4d..397be35a 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -1,6 +1,7 @@ module Geometry -export AbstractElevation, DataElevation, SyntheticElevation, AABB +export AbstractElevation, AABB, DataElevation, SyntheticElevation +export domain, rotated_domain, rotation, generate_elevation, extents, center, dilate abstract type AbstractElevation{T<:Real} end @@ -30,4 +31,34 @@ struct SyntheticElevation{T,B,S} <: AbstractElevation{T} domain::AABB{T} end +domain(dem::SyntheticElevation) = dem.domain +rotated_domain(dem::AbstractElevation) = domain(dem) +rotation(::AbstractElevation) = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] + +"Construct AABB from coordinates" +AABB(xs, ys, zs) = AABB(extrema(xs)..., extrema(ys)..., extrema(zs)...) + +"Generate SyntheticElevation data." +function generate_elevation(lx, ly, zmin, zmax, z_bed, z_surf) + domain = AABB(-lx / 2, lx / 2, -ly / 2, ly / 2, zmin, zmax) + return SyntheticElevation(z_bed, z_surf, domain) end + +"AABB extents" +function extents(box::AABB) + return box.xmax - box.xmin, box.ymax - box.ymin, box.zmax - box.zmin +end + +"AABB center" +function center(box::AABB{T}) where {T} + half = convert(T, 0.5) + return half * (box.xmin + box.xmax), half * (box.ymin + box.ymax), half * (box.zmin + box.zmax) +end + +"Dilate AABB by extending its limits around the center by certain fraction in each dimension" +function dilate(box::AABB, fractions) + Δx, Δy, Δz = extents(box) .* fractions + return AABB(box.xmin - Δx, box.xmax + Δx, box.ymin - Δy, box.ymax + Δy, box.zmin - Δz, box.zmax + Δz) +end + +end # module From f394f81f1d60924ea35f35221ee72692df017068 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Fri, 8 Mar 2024 13:27:44 +0100 Subject: [PATCH 17/32] Try export --- src/FastIce.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/FastIce.jl b/src/FastIce.jl index 69537a2b..927429e4 100644 --- a/src/FastIce.jl +++ b/src/FastIce.jl @@ -25,6 +25,8 @@ https://github.com/PTsolvers/FastIce.jl greet(; kwargs...) = printstyled(GREETING; kwargs...) greet_fast(; kwargs...) = printstyled(GREETING_FAST; kwargs...) +export Geometry + # core modules (included in alphabetical order) include("Geometry.jl") include("LevelSets/LevelSets.jl") From e202c258228a4a572726f8cf1705f6e58c0d15ab Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Fri, 8 Mar 2024 13:29:44 +0100 Subject: [PATCH 18/32] revert --- src/FastIce.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FastIce.jl b/src/FastIce.jl index 927429e4..d9fc4571 100644 --- a/src/FastIce.jl +++ b/src/FastIce.jl @@ -25,7 +25,7 @@ https://github.com/PTsolvers/FastIce.jl greet(; kwargs...) = printstyled(GREETING; kwargs...) greet_fast(; kwargs...) = printstyled(GREETING_FAST; kwargs...) -export Geometry +# export Geometry # core modules (included in alphabetical order) include("Geometry.jl") From 17023f92c68d38e1f38b296e6fe9f6b5686a0200 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Sun, 10 Mar 2024 22:25:00 +0100 Subject: [PATCH 19/32] Rework Geometry --- src/FastIce.jl | 2 -- src/Geometry.jl | 59 ++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 48 insertions(+), 13 deletions(-) diff --git a/src/FastIce.jl b/src/FastIce.jl index d9fc4571..69537a2b 100644 --- a/src/FastIce.jl +++ b/src/FastIce.jl @@ -25,8 +25,6 @@ https://github.com/PTsolvers/FastIce.jl greet(; kwargs...) = printstyled(GREETING; kwargs...) greet_fast(; kwargs...) = printstyled(GREETING_FAST; kwargs...) -# export Geometry - # core modules (included in alphabetical order) include("Geometry.jl") include("LevelSets/LevelSets.jl") diff --git a/src/Geometry.jl b/src/Geometry.jl index 397be35a..e383886c 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -1,7 +1,7 @@ module Geometry -export AbstractElevation, AABB, DataElevation, SyntheticElevation -export domain, rotated_domain, rotation, generate_elevation, extents, center, dilate +# export AbstractElevation, AABB, DataElevation, SyntheticElevation +# export domain, rotated_domain, rotation, generate_elevation, extents, center, dilate abstract type AbstractElevation{T<:Real} end @@ -35,19 +35,15 @@ domain(dem::SyntheticElevation) = dem.domain rotated_domain(dem::AbstractElevation) = domain(dem) rotation(::AbstractElevation) = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] +domain(dem::DataElevation) = dem.domain +rotated_domain(dem::DataElevation) = dem.rotated_domain +rotation(dem::DataElevation) = dem.rotation + "Construct AABB from coordinates" AABB(xs, ys, zs) = AABB(extrema(xs)..., extrema(ys)..., extrema(zs)...) -"Generate SyntheticElevation data." -function generate_elevation(lx, ly, zmin, zmax, z_bed, z_surf) - domain = AABB(-lx / 2, lx / 2, -ly / 2, ly / 2, zmin, zmax) - return SyntheticElevation(z_bed, z_surf, domain) -end - "AABB extents" -function extents(box::AABB) - return box.xmax - box.xmin, box.ymax - box.ymin, box.zmax - box.zmin -end +extents(box::AABB) = box.xmax - box.xmin, box.ymax - box.ymin, box.zmax - box.zmin "AABB center" function center(box::AABB{T}) where {T} @@ -61,4 +57,45 @@ function dilate(box::AABB, fractions) return AABB(box.xmin - Δx, box.xmax + Δx, box.ymin - Δy, box.ymax + Δy, box.zmin - Δz, box.zmax + Δz) end +"Filter NaNs." +filtered(X) = filter(_x -> !isnan(_x), X) + +"Create AABB enclosing both box1 and box2" +function union(box1::AABB, box2::AABB) + return AABB(min(box1.xmin, box2.xmin), max(box1.xmax, box2.xmax), + min(box1.ymin, box2.ymin), max(box1.ymax, box2.ymax), + min(box1.zmin, box2.zmin), max(box1.zmax, box2.zmax)) +end + +"Rotate field `X`, `Y`, `Z` with rotation matrix `R`." +function rotate(X, Y, Z, R) + xrot = R[1, 1] .* X .+ R[1, 2] .* Y .+ R[1, 3] .* Z + yrot = R[2, 1] .* X .+ R[2, 2] .* Y .+ R[2, 3] .* Z + zrot = R[3, 1] .* X .+ R[3, 2] .* Y .+ R[3, 3] .* Z + return xrot, yrot, zrot +end + +"Rotate field `X`, `Y`, `Z` with rotation matrix `R` and return extents." +rotate_minmax(X, Y, Z, R) = rotate(collect(extrema(X)), collect(extrema(Y)), collect(extrema(filtered(Z))), R) + +"Generate SyntheticElevation data." +function generate_elevation(lx, ly, zmin, zmax, z_bed, z_surf) + domain = AABB(-lx / 2, lx / 2, -ly / 2, ly / 2, zmin, zmax) + return SyntheticElevation(z_bed, z_surf, domain) +end + +"Generate DataElevation data." +function DataElevation(x, y, offsets, z_bed, z_surf, R) + # get non-rotated domain + domain = AABB(extrema(x)..., extrema(y)..., Float64(minimum([minimum(filtered(z_bed)), minimum(filtered(z_surf))])), + Float64(maximum([maximum(filtered(z_bed)), maximum(filtered(z_surf))]))) + # rotate bed and surface + bed_extents = AABB(rotate_minmax(x, y, z_bed, R)...) + surf_extents = AABB(rotate_minmax(x, y, z_surf, R)...) + # get rotated domain + rotated_domain = union(bed_extents, surf_extents) + dem = DataElevation(x, y, offsets, z_bed, z_surf, R, domain, rotated_domain) + return dem +end + end # module From c2f402252098de83ecdbf4eed0db9a60d0490a1b Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Sun, 10 Mar 2024 22:56:53 +0100 Subject: [PATCH 20/32] Rework synthetic geometries --- scripts_future_API/generate_synthetic.jl | 29 +++++++ src/FastIce.jl | 2 +- src/Geometries/Geometries.jl | 11 +++ src/Geometries/elevation_data.jl | 46 +++++++++++ src/Geometries/geometry_helpers.jl | 55 ++++++++++++ src/Geometries/synthetic_geometries.jl | 89 ++++++++++++++++++++ src/Geometry.jl | 101 ----------------------- 7 files changed, 231 insertions(+), 102 deletions(-) create mode 100644 scripts_future_API/generate_synthetic.jl create mode 100644 src/Geometries/Geometries.jl create mode 100644 src/Geometries/elevation_data.jl create mode 100644 src/Geometries/geometry_helpers.jl create mode 100644 src/Geometries/synthetic_geometries.jl delete mode 100644 src/Geometry.jl diff --git a/scripts_future_API/generate_synthetic.jl b/scripts_future_API/generate_synthetic.jl new file mode 100644 index 00000000..c18b7427 --- /dev/null +++ b/scripts_future_API/generate_synthetic.jl @@ -0,0 +1,29 @@ +using FastIce.Geometries + +nx, ny = 128, 128 +lx, ly = 5.0, 5.0 +zmin, zmax = 0.0, 1.0 +synth_dem = make_synthetic(nx, ny, lx, ly, zmin, zmax, :turtle) + +@views function visme(synth_dem) + x = LinRange(synth_dem.domain.xmin, synth_dem.domain.xmax, nx + 1) + y = LinRange(synth_dem.domain.ymin, synth_dem.domain.ymax, ny + 1) + bed = synth_dem.z_bed + surf = synth_dem.z_surf + + surf2 = copy(surf) + surf2[surf. !isnan(_x), X) + +"Create AABB enclosing both box1 and box2" +function union(box1::AABB, box2::AABB) + return AABB(min(box1.xmin, box2.xmin), max(box1.xmax, box2.xmax), + min(box1.ymin, box2.ymin), max(box1.ymax, box2.ymax), + min(box1.zmin, box2.zmin), max(box1.zmax, box2.zmax)) +end + +"Rotate field `X`, `Y`, `Z` with rotation matrix `R`." +function rotate(X, Y, Z, R) + xrot = R[1, 1] .* X .+ R[1, 2] .* Y .+ R[1, 3] .* Z + yrot = R[2, 1] .* X .+ R[2, 2] .* Y .+ R[2, 3] .* Z + zrot = R[3, 1] .* X .+ R[3, 2] .* Y .+ R[3, 3] .* Z + return xrot, yrot, zrot +end + +"Rotate field `X`, `Y`, `Z` with rotation matrix `R` and return extents." +rotate_minmax(X, Y, Z, R) = rotate(collect(extrema(X)), collect(extrema(Y)), collect(extrema(filtered(Z))), R) diff --git a/src/Geometries/synthetic_geometries.jl b/src/Geometries/synthetic_geometries.jl new file mode 100644 index 00000000..58c41956 --- /dev/null +++ b/src/Geometries/synthetic_geometries.jl @@ -0,0 +1,89 @@ +function generate_turtle(nx, ny, lx, ly, zmin, zmax) + @info "expects lx, ly, lz = 5.0, 5.0, 1.0" + # Numerics + amp = 0.1 + ω = 10π + tanβ = tan(-π/6) + el = 0.35 + gl = 0.9 + x = LinRange(-lx / 2, lx / 2, nx + 1) + y = LinRange(-lx / 2, ly / 2, ny + 1) + # Functions + generate_bed(x, y) = amp * sin(ω * x) * sin(ω * y) + tanβ * x + el + (1.5 * y)^2 + generate_surf(x, y) = 1.0 - ((1.7 * x + 0.22)^2 + (0.5 * y)^2) + # Generate + bed = [(zmax - zmin) * generate_bed(x / lx, y / ly) for x in x, y in y] + surf = [gl * generate_surf(x / lx, y / ly) for x in x, y in y] + return bed, surf +end + +function generate_dome(nx, ny, lx, ly) + # Numerics + h_0 = 0.8 + x = LinRange(-lx / 2, lx / 2, nx + 1) + y = LinRange(-lx / 2, ly / 2, ny + 1) + # Functions + generate_bed(x, y) = 0.0 + generate_surf(x, y) = h_0 - (x^2 + y^2) + # Compute + bed = [generate_bed(_x, _y) for _x in x, _y in y] + surf = [generate_surf(_x, _y) for _x in x, _y in y] + surf[surf. !isnan(_x), X) - -"Create AABB enclosing both box1 and box2" -function union(box1::AABB, box2::AABB) - return AABB(min(box1.xmin, box2.xmin), max(box1.xmax, box2.xmax), - min(box1.ymin, box2.ymin), max(box1.ymax, box2.ymax), - min(box1.zmin, box2.zmin), max(box1.zmax, box2.zmax)) -end - -"Rotate field `X`, `Y`, `Z` with rotation matrix `R`." -function rotate(X, Y, Z, R) - xrot = R[1, 1] .* X .+ R[1, 2] .* Y .+ R[1, 3] .* Z - yrot = R[2, 1] .* X .+ R[2, 2] .* Y .+ R[2, 3] .* Z - zrot = R[3, 1] .* X .+ R[3, 2] .* Y .+ R[3, 3] .* Z - return xrot, yrot, zrot -end - -"Rotate field `X`, `Y`, `Z` with rotation matrix `R` and return extents." -rotate_minmax(X, Y, Z, R) = rotate(collect(extrema(X)), collect(extrema(Y)), collect(extrema(filtered(Z))), R) - -"Generate SyntheticElevation data." -function generate_elevation(lx, ly, zmin, zmax, z_bed, z_surf) - domain = AABB(-lx / 2, lx / 2, -ly / 2, ly / 2, zmin, zmax) - return SyntheticElevation(z_bed, z_surf, domain) -end - -"Generate DataElevation data." -function DataElevation(x, y, offsets, z_bed, z_surf, R) - # get non-rotated domain - domain = AABB(extrema(x)..., extrema(y)..., Float64(minimum([minimum(filtered(z_bed)), minimum(filtered(z_surf))])), - Float64(maximum([maximum(filtered(z_bed)), maximum(filtered(z_surf))]))) - # rotate bed and surface - bed_extents = AABB(rotate_minmax(x, y, z_bed, R)...) - surf_extents = AABB(rotate_minmax(x, y, z_surf, R)...) - # get rotated domain - rotated_domain = union(bed_extents, surf_extents) - dem = DataElevation(x, y, offsets, z_bed, z_surf, R, domain, rotated_domain) - return dem -end - -end # module From ebb15f714994d5838ceb13304ac9bba3c85d0471 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Sun, 10 Mar 2024 23:08:37 +0100 Subject: [PATCH 21/32] Fixup --- src/Geometries/elevation_data.jl | 16 ++++++++++++++++ src/Geometries/geometry_helpers.jl | 16 ---------------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/Geometries/elevation_data.jl b/src/Geometries/elevation_data.jl index 4b73c0c8..847bb050 100644 --- a/src/Geometries/elevation_data.jl +++ b/src/Geometries/elevation_data.jl @@ -1,3 +1,19 @@ +struct AABB{T<:Union{Real}} + xmin::T + xmax::T + ymin::T + ymax::T + zmin::T + zmax::T +end + +""" + AABB(xs, ys, zs) + +Construct an axis aligend bounding box `AABB` from the coordinates `xs`, `ys`, and `zs`. +""" +AABB(xs, ys, zs) = AABB(extrema(xs)..., extrema(ys)..., extrema(zs)...) + abstract type AbstractElevation{T<:Real} end struct DataElevation{T,V<:AbstractArray{T},R<:AbstractRange{T},M<:AbstractMatrix{T}} <: AbstractElevation{T} diff --git a/src/Geometries/geometry_helpers.jl b/src/Geometries/geometry_helpers.jl index 232d279a..2c4e7033 100644 --- a/src/Geometries/geometry_helpers.jl +++ b/src/Geometries/geometry_helpers.jl @@ -1,19 +1,3 @@ -struct AABB{T<:Union{Real}} - xmin::T - xmax::T - ymin::T - ymax::T - zmin::T - zmax::T -end - -""" - AABB(xs, ys, zs) - -Construct an axis aligend bounding box `AABB` from the coordinates `xs`, `ys`, and `zs`. -""" -AABB(xs, ys, zs) = AABB(extrema(xs)..., extrema(ys)..., extrema(zs)...) - """ extents(box::AABB) From 842d5de660122cc497ab6a4e453a7839602dbbe7 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Mon, 11 Mar 2024 15:52:58 +0100 Subject: [PATCH 22/32] Non-MPI workflow --- .../test_levelsets_volumefractions.jl | 80 +++++++++++++++---- src/Geometries/Geometries.jl | 2 + 2 files changed, 66 insertions(+), 16 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl index db2e9fbe..f522a351 100644 --- a/scripts_future_API/test_levelsets_volumefractions.jl +++ b/scripts_future_API/test_levelsets_volumefractions.jl @@ -3,17 +3,13 @@ using Chmy.Grids using Chmy.Fields using FastIce.LevelSets -using FastIce.Geometry +using FastIce.Geometries using KernelAbstractions -using JLD2 # using CUDA # using AMDGPU -# Data -data_path = "./low_res_dome_rough_glacier.jld2" - # Select backend backend = CPU() # backend = CUDABackend() @@ -49,16 +45,40 @@ function extract_dem(arch::Architecture, data_path::String, z_resol=nothing) return dem_surf, dem_bed, dem_grid, Ψ_grid end -dem_surf, dem_bed, dem_grid, Ψ_grid = extract_dem(arch, data_path) +function extract_dem(arch::Architecture, data::SyntheticElevation, z_resol=nothing) + z_surf = data.z_surf + z_bed = data.z_bed + dm = data.domain + lx, ly, lz = extents(dm) + nx, ny = size(z_surf) .- 1 + nz = isnothing(z_resol) ? ceil(Int, lz / lx * nx) : z_resol + Ψ_grid = UniformGrid(arch; origin=(dm.xmin, dm.ymin, dm.zmin), extent=(lx, ly, lz), dims=(nx, ny, nz)) + dem_grid = UniformGrid(arch; origin=(dm.xmin, dm.ymin), extent=(lx, ly), dims=(nx, ny)) + dem_surf = Field(arch.backend, dem_grid, Vertex()) + dem_bed = Field(arch.backend, dem_grid, Vertex()) + set!(dem_surf, z_surf) + set!(dem_bed, z_bed) + dem = (surf=dem_surf, bed=dem_bed) + return dem, dem_grid, Ψ_grid +end + +nx, ny = 128, 128 +lx, ly = 5.0, 5.0 +zmin, zmax = 0.0, 1.0 +type = :turtle + +synth_dem = make_synthetic(nx, ny, lx, ly, zmin, zmax, type) +dem, dem_grid, Ψ_grid = extract_dem(arch, synth_dem) Ψ = ( na=Field(backend, Ψ_grid, Vertex()), ns=Field(backend, Ψ_grid, Vertex()), ) -for phase in eachindex(Ψ) - compute_level_set_from_dem!(arch, Ψ[phase], dem_surf, dem_grid, Ψ_grid) -end +compute_level_set_from_dem!(arch, Ψ.na, dem.surf, dem_grid, Ψ_grid) +compute_level_set_from_dem!(arch, Ψ.ns, dem.bed, dem_grid, Ψ_grid) +# invert level set to set what's below the DEM surface as inside +@. Ψ.ns *= -1.0 wt = ( na=volfrac_field(backend, Ψ_grid), @@ -69,11 +89,39 @@ for phase in eachindex(Ψ) compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], Ψ_grid) end -# Save -save_path = "Vavilov_500m_out.jld2" -jldopen(vavilov_path, "a+") do file - file["level_sets_na"] = Array(interior(Ψ[1])) - file["level_sets_ns"] = Array(interior(Ψ[2])) - file["volume_frac_na"] = Array.(interior.(Tuple(wt[1]))) - file["volume_frac_ns"] = Array.(interior.(Tuple(wt[2]))) + +@views function visme(synth_dem, Ψ, wt) + nx, ny, nz = size(Ψ.ns) + x = LinRange(synth_dem.domain.xmin, synth_dem.domain.xmax, nx) + y = LinRange(synth_dem.domain.ymin, synth_dem.domain.ymax, ny) + z = LinRange(synth_dem.domain.zmin, synth_dem.domain.zmax, nz) + bed = synth_dem.z_bed + surf = synth_dem.z_surf + + surf2 = copy(surf) + surf2[surf. Date: Mon, 11 Mar 2024 15:53:21 +0100 Subject: [PATCH 23/32] WIP - MPI (distributed) workflow --- .../test_levelsets_volumefractions_mpi.jl | 104 ++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 scripts_future_API/test_levelsets_volumefractions_mpi.jl diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl new file mode 100644 index 00000000..3ef49cf1 --- /dev/null +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -0,0 +1,104 @@ +using Chmy.Architectures +using Chmy.Grids +using Chmy.Fields + +using FastIce.LevelSets +using FastIce.Geometries + +using KernelAbstractions + +# using CUDA +# using AMDGPU + +# Select backend +backend = CPU() +# backend = CUDABackend() +# backend = ROCBackend() + +using Chmy.Distributed +using MPI +MPI.Init() + +# function main(backend=CPU(); nxyz_l=(126, 126)) + backend=CPU() + nxyz_l=(126, 126) + arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) + arch2 = Arch(backend, MPI.COMM_WORLD, (0, 0)) + topo = topology(arch) + me = global_rank(topo) + # geometry + zmin, zmax = 0.0, 1.0 + lx, ly, lz = 5.0, 5.0, zmax - zmin + nx, ny = nxyz_l[1:2] + nz = (length(nxyz_l) < 3) ? ceil(Int, lz / lx * nxyz_l[1]) : nxyz_l[3] + dims_l = (nx, ny, nz) + dims_g = dims_l .* dims(topo) + grid = UniformGrid(arch; origin=(-lx/2, -ly/2, -lz/2), extent=(lx, ly, lz), dims=dims_g) + grid_dem = UniformGrid(arch2; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=dims_g[1:2]) + # synthetic topo + amp = 0.1 + ω = 10π + tanβ = tan(-π/6) + el = 0.35 + gl = 0.9 + # bed and surf functions + # type = :turtle + generate_surf(x, y, gl, lx, ly) = gl * (1.0 - ((1.7 * x / lx + 0.22)^2 + (0.5 * y / ly)^2)) + generate_bed(x, y, amp, ω, tanβ, el, lx, ly, lz) = lz * (amp * sin(ω * x / lx) * sin(ω * y /ly) + tanβ * x / lx + el + (1.5 * y / ly)^2) + # init fields + dem_surf = Field(backend, grid_dem, Vertex()) + dem_bed = Field(backend, grid_dem, Vertex()) + Ψ = (na=Field(backend, grid, Vertex()), + ns=Field(backend, grid, Vertex())) + wt = (na=volfrac_field(backend, grid), + ns=volfrac_field(backend, grid)) + # generate 2D synthetic dem + set!(dem_surf, grid_dem, generate_surf; parameters=(gl=gl, lx=lx, ly=ly)) + set!(dem_bed, grid_dem, generate_bed; parameters=(amp=amp, ω=ω, tanβ=tanβ, el=el, lx=lx, ly=ly, lz=lz)) + # comput level set + compute_level_set_from_dem!(arch, Ψ.na, dem_surf, grid_dem, grid) + compute_level_set_from_dem!(arch, Ψ.ns, dem_bed, grid_dem, grid) + # invert level set to set what's below the DEM surface as inside + @. Ψ.ns *= -1.0 + # volume fractions + for phase in eachindex(Ψ) + compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], grid) + end + # gather + dem_bed_v = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(dem_bed)) .* dims(topo)[1:2]) : nothing + dem_surf_v = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(dem_surf)) .* dims(topo)[1:2]) : nothing + wt_na_c = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.na.c)) .* dims(topo)) : nothing + wt_ns_c = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.ns.c)) .* dims(topo)) : nothing + gather!(arch2, dem_bed_v, dem_bed) + gather!(arch2, dem_surf_v, dem_surf) + gather!(arch, wt_na_c, wt.na.c) + gather!(arch, wt_ns_c, wt.ns.c) + # visualise + dem_surf_v[dem_surf_v. Date: Mon, 11 Mar 2024 17:31:03 +0100 Subject: [PATCH 24/32] Fixup turtle --- .../test_levelsets_volumefractions_mpi.jl | 68 ++++++++++--------- 1 file changed, 35 insertions(+), 33 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl index 3ef49cf1..79268d3a 100644 --- a/scripts_future_API/test_levelsets_volumefractions_mpi.jl +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -6,7 +6,7 @@ using FastIce.LevelSets using FastIce.Geometries using KernelAbstractions - +using CairoMakie # using CUDA # using AMDGPU @@ -22,8 +22,8 @@ MPI.Init() # function main(backend=CPU(); nxyz_l=(126, 126)) backend=CPU() nxyz_l=(126, 126) - arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) - arch2 = Arch(backend, MPI.COMM_WORLD, (0, 0)) + arch = Arch(backend, MPI.COMM_WORLD, (2, 2, 2)) + arch2 = Arch(backend, MPI.COMM_WORLD, (2, 2)) topo = topology(arch) me = global_rank(topo) # geometry @@ -33,7 +33,7 @@ MPI.Init() nz = (length(nxyz_l) < 3) ? ceil(Int, lz / lx * nxyz_l[1]) : nxyz_l[3] dims_l = (nx, ny, nz) dims_g = dims_l .* dims(topo) - grid = UniformGrid(arch; origin=(-lx/2, -ly/2, -lz/2), extent=(lx, ly, lz), dims=dims_g) + grid = UniformGrid(arch; origin=(-lx/2, -ly/2, 0.0), extent=(lx, ly, lz), dims=dims_g) grid_dem = UniformGrid(arch2; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=dims_g[1:2]) # synthetic topo amp = 0.1 @@ -65,40 +65,42 @@ MPI.Init() compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], grid) end # gather - dem_bed_v = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(dem_bed)) .* dims(topo)[1:2]) : nothing - dem_surf_v = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(dem_surf)) .* dims(topo)[1:2]) : nothing - wt_na_c = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.na.c)) .* dims(topo)) : nothing - wt_ns_c = (me==0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.ns.c)) .* dims(topo)) : nothing + dem_bed_v = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(dem_bed)) .* dims(topo)[1:2]) : nothing + dem_surf_v = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(dem_surf)) .* dims(topo)[1:2]) : nothing + wt_na_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.na.c)) .* dims(topo)) : nothing + wt_ns_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.ns.c)) .* dims(topo)) : nothing gather!(arch2, dem_bed_v, dem_bed) gather!(arch2, dem_surf_v, dem_surf) gather!(arch, wt_na_c, wt.na.c) gather!(arch, wt_ns_c, wt.ns.c) # visualise - dem_surf_v[dem_surf_v. Date: Tue, 12 Mar 2024 13:32:09 +0100 Subject: [PATCH 25/32] Update --- .../test_levelsets_volumefractions_mpi.jl | 76 ++++++++++--------- src/LevelSets/compute_level_sets.jl | 8 +- 2 files changed, 43 insertions(+), 41 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl index 79268d3a..59517aff 100644 --- a/scripts_future_API/test_levelsets_volumefractions_mpi.jl +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -3,7 +3,6 @@ using Chmy.Grids using Chmy.Fields using FastIce.LevelSets -using FastIce.Geometries using KernelAbstractions using CairoMakie @@ -19,68 +18,67 @@ using Chmy.Distributed using MPI MPI.Init() -# function main(backend=CPU(); nxyz_l=(126, 126)) - backend=CPU() - nxyz_l=(126, 126) - arch = Arch(backend, MPI.COMM_WORLD, (2, 2, 2)) - arch2 = Arch(backend, MPI.COMM_WORLD, (2, 2)) - topo = topology(arch) - me = global_rank(topo) +function main(backend=CPU(); nxyz_l=(126, 126)) + arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) + topo = topology(arch) + me = global_rank(topo) # geometry zmin, zmax = 0.0, 1.0 lx, ly, lz = 5.0, 5.0, zmax - zmin - nx, ny = nxyz_l[1:2] - nz = (length(nxyz_l) < 3) ? ceil(Int, lz / lx * nxyz_l[1]) : nxyz_l[3] + nx, ny = nxyz_l[1:2] + nz = (length(nxyz_l) < 3) ? ceil(Int, lz / lx * nxyz_l[1]) : nxyz_l[3] dims_l = (nx, ny, nz) dims_g = dims_l .* dims(topo) grid = UniformGrid(arch; origin=(-lx/2, -ly/2, 0.0), extent=(lx, ly, lz), dims=dims_g) - grid_dem = UniformGrid(arch2; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=dims_g[1:2]) + + # synthetic topo ############### + arch_2d = Arch(backend) + grid_2d = UniformGrid(arch_2d; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=dims_g[1:2]) # synthetic topo amp = 0.1 ω = 10π tanβ = tan(-π/6) el = 0.35 gl = 0.9 - # bed and surf functions + # synthetic bed and surf # type = :turtle generate_surf(x, y, gl, lx, ly) = gl * (1.0 - ((1.7 * x / lx + 0.22)^2 + (0.5 * y / ly)^2)) generate_bed(x, y, amp, ω, tanβ, el, lx, ly, lz) = lz * (amp * sin(ω * x / lx) * sin(ω * y /ly) + tanβ * x / lx + el + (1.5 * y / ly)^2) + dem_surf = FunctionField(generate_surf, grid_2d, Vertex(); parameters=(gl=gl, lx=lx, ly=ly)) + dem_bed = FunctionField(generate_bed, grid_2d, Vertex(); parameters=(amp=amp, ω=ω, tanβ=tanβ, el=el, lx=lx, ly=ly, lz=lz)) + # synthetic topo ############### + # init fields - dem_surf = Field(backend, grid_dem, Vertex()) - dem_bed = Field(backend, grid_dem, Vertex()) Ψ = (na=Field(backend, grid, Vertex()), ns=Field(backend, grid, Vertex())) wt = (na=volfrac_field(backend, grid), ns=volfrac_field(backend, grid)) - # generate 2D synthetic dem - set!(dem_surf, grid_dem, generate_surf; parameters=(gl=gl, lx=lx, ly=ly)) - set!(dem_bed, grid_dem, generate_bed; parameters=(amp=amp, ω=ω, tanβ=tanβ, el=el, lx=lx, ly=ly, lz=lz)) # comput level set - compute_level_set_from_dem!(arch, Ψ.na, dem_surf, grid_dem, grid) - compute_level_set_from_dem!(arch, Ψ.ns, dem_bed, grid_dem, grid) - # invert level set to set what's below the DEM surface as inside - @. Ψ.ns *= -1.0 + compute_level_set_from_dem!(arch, Ψ.na, dem_surf, grid_2d, grid) + compute_level_set_from_dem!(arch, Ψ.ns, dem_bed, grid_2d, grid) + @. Ψ.ns *= -1.0 # invert level set to set what's below the DEM surface as inside # volume fractions for phase in eachindex(Ψ) compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], grid) end + + # compute physics or else + + # postprocessing # gather - dem_bed_v = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(dem_bed)) .* dims(topo)[1:2]) : nothing - dem_surf_v = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(dem_surf)) .* dims(topo)[1:2]) : nothing - wt_na_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.na.c)) .* dims(topo)) : nothing - wt_ns_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.ns.c)) .* dims(topo)) : nothing - gather!(arch2, dem_bed_v, dem_bed) - gather!(arch2, dem_surf_v, dem_surf) + wt_na_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.na.c)) .* dims(topo)) : nothing + wt_ns_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.ns.c)) .* dims(topo)) : nothing gather!(arch, wt_na_c, wt.na.c) gather!(arch, wt_ns_c, wt.ns.c) # visualise if me == 0 - dem_surf_v[dem_surf_v. (Ψ, dem, dem_grid2D, grid, cutoff, R); bc=batch(grid, Ψ => Neumann(); exchange=Ψ)) From 1400897203fea81e5d0bd2fc620ece4a46558c70 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Tue, 12 Mar 2024 16:38:50 +0100 Subject: [PATCH 26/32] Fixup --- .../test_levelsets_volumefractions_mpi.jl | 4 ++-- src/LevelSets/compute_volume_fractions.jl | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl index 59517aff..e4263a65 100644 --- a/scripts_future_API/test_levelsets_volumefractions_mpi.jl +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -44,8 +44,8 @@ function main(backend=CPU(); nxyz_l=(126, 126)) # type = :turtle generate_surf(x, y, gl, lx, ly) = gl * (1.0 - ((1.7 * x / lx + 0.22)^2 + (0.5 * y / ly)^2)) generate_bed(x, y, amp, ω, tanβ, el, lx, ly, lz) = lz * (amp * sin(ω * x / lx) * sin(ω * y /ly) + tanβ * x / lx + el + (1.5 * y / ly)^2) - dem_surf = FunctionField(generate_surf, grid_2d, Vertex(); parameters=(gl=gl, lx=lx, ly=ly)) - dem_bed = FunctionField(generate_bed, grid_2d, Vertex(); parameters=(amp=amp, ω=ω, tanβ=tanβ, el=el, lx=lx, ly=ly, lz=lz)) + dem_surf = FunctionField(generate_surf, grid_2d, Vertex(); parameters=(; gl, lx, ly)) + dem_bed = FunctionField(generate_bed, grid_2d, Vertex(); parameters=(; amp, ω, tanβ, el, lx, ly, lz)) # synthetic topo ############### # init fields diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl index bd83e7ba..17b2b03d 100644 --- a/src/LevelSets/compute_volume_fractions.jl +++ b/src/LevelSets/compute_volume_fractions.jl @@ -118,17 +118,17 @@ function volfrac(rect::Rect3, ϕ::Vec{8}) end function volfrac_field(backend::Backend, grid::UniformGrid{2}, args...; kwargs...) - (c=Field(backend, grid, Center()), - x=Field(backend, grid, (Vertex(), Center())), - y=Field(backend, grid, (Center(), Vertex())), + (c =Field(backend, grid, Center()), + x =Field(backend, grid, (Vertex(), Center())), + y =Field(backend, grid, (Center(), Vertex())), xy=Field(backend, grid, Vertex())) end function volfrac_field(backend::Backend, grid::UniformGrid{3}, args...; kwargs...) - (c=Field(backend, grid, Center()), - x=Field(backend, grid, (Vertex(), Center(), Center())), - y=Field(backend, grid, (Center(), Vertex(), Center())), - z=Field(backend, grid, (Center(), Center(), Vertex())), + (c =Field(backend, grid, Center()), + x =Field(backend, grid, (Vertex(), Center(), Center())), + y =Field(backend, grid, (Center(), Vertex(), Center())), + z =Field(backend, grid, (Center(), Center(), Vertex())), xy=Field(backend, grid, (Vertex(), Vertex(), Center())), xz=Field(backend, grid, (Vertex(), Center(), Vertex())), yz=Field(backend, grid, (Center(), Vertex(), Vertex()))) From b0259487f7b8ebc82bbe6290af6846558c21f94c Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Tue, 12 Mar 2024 17:16:17 +0100 Subject: [PATCH 27/32] Fixup GPU --- .../test_levelsets_volumefractions_mpi.jl | 25 +++++------ src/LevelSets/LevelSets.jl | 6 +-- src/LevelSets/compute_level_sets.jl | 25 ----------- src/LevelSets/compute_levelsets.jl | 42 +++++++++++++++++++ src/LevelSets/compute_volume_fractions.jl | 6 +-- src/LevelSets/volume_fractions_kernels.jl | 4 +- 6 files changed, 63 insertions(+), 45 deletions(-) delete mode 100644 src/LevelSets/compute_level_sets.jl create mode 100644 src/LevelSets/compute_levelsets.jl diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl index e4263a65..77ae2c5f 100644 --- a/scripts_future_API/test_levelsets_volumefractions_mpi.jl +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -7,7 +7,7 @@ using FastIce.LevelSets using KernelAbstractions using CairoMakie # using CUDA -# using AMDGPU +using AMDGPU # Select backend backend = CPU() @@ -54,12 +54,13 @@ function main(backend=CPU(); nxyz_l=(126, 126)) wt = (na=volfrac_field(backend, grid), ns=volfrac_field(backend, grid)) # comput level set - compute_level_set_from_dem!(arch, Ψ.na, dem_surf, grid_2d, grid) - compute_level_set_from_dem!(arch, Ψ.ns, dem_bed, grid_2d, grid) - @. Ψ.ns *= -1.0 # invert level set to set what's below the DEM surface as inside + compute_levelset_from_dem!(arch, Ψ.na, dem_surf, grid_2d, grid) + compute_levelset_from_dem!(arch, Ψ.ns, dem_bed, grid_2d, grid) + invert_levelset!(arch, Ψ.ns, grid) + # @. Ψ.ns *= -1.0 # invert level set to set what's below the DEM surface as inside # volume fractions for phase in eachindex(Ψ) - compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], grid) + compute_volfrac_from_levelset!(arch, wt[phase], Ψ[phase], grid) end # compute physics or else @@ -77,8 +78,8 @@ function main(backend=CPU(); nxyz_l=(126, 126)) dem_surf_[dem_surf. Array; colormap=:turbo), + p1_ = surface!(axs.ax1, x_g, y_g, interior(dem_surf_) |> Array; colormap=:turbo), + p2 = plot!(axs.ax2, x_g, interior(dem_bed)[:, sly] |> Array), + p2_ = plot!(axs.ax2, x_g, interior(dem_surf_)[:, sly] |> Array), p3 = heatmap!(axs.ax3, wt_na_c[:, sly, :]; colormap=:turbo), p4 = heatmap!(axs.ax4, wt_ns_c[:, sly, :]; colormap=:turbo), p5 = heatmap!(axs.ax5, wt_na_c[slx, :, :]; colormap=:turbo), @@ -103,6 +104,6 @@ function main(backend=CPU(); nxyz_l=(126, 126)) return end -main() +main(backend; nxyz_l=(126, 126)) MPI.Finalize() diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index 7117fc4a..d4ff52fc 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -1,7 +1,7 @@ module LevelSets -export compute_level_set_from_dem! -export volfrac_field, compute_volfrac_from_level_set! +export compute_levelset_from_dem!, invert_levelset! +export volfrac_field, compute_volfrac_from_levelset! using Chmy.Architectures using Chmy.BoundaryConditions @@ -13,7 +13,7 @@ using KernelAbstractions using LinearAlgebra, GeometryBasics include("signed_distances.jl") -include("compute_level_sets.jl") +include("compute_levelsets.jl") include("compute_volume_fractions.jl") end diff --git a/src/LevelSets/compute_level_sets.jl b/src/LevelSets/compute_level_sets.jl deleted file mode 100644 index cfe99653..00000000 --- a/src/LevelSets/compute_level_sets.jl +++ /dev/null @@ -1,25 +0,0 @@ -""" - init_level_set!(Ψ::Field, dem::AbstractField, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R) - -Initialize level sets. -""" -@kernel inbounds = true function init_level_set!(Ψ::Field, dem::AbstractField, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R, O) - I = @index(Global, NTuple) - I = I + O - x, y, z = coord(Ψ_grid, location(Ψ), I...) - P = R * Point3(x, y, z) - ud, sgn = sd_dem(P, cutoff, dem, dem_grid) - Ψ[I...] = ud * sgn -end - -""" - compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I) - -Compute level sets from dem. -""" -function compute_level_set_from_dem!(arch::Architecture, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I; outer_width=(16, 8, 4)) - cutoff = 4maximum(spacing(grid)) - launch = Launcher(arch, grid; outer_width=outer_width) - launch(arch, grid, init_level_set! => (Ψ, dem, dem_grid2D, grid, cutoff, R); bc=batch(grid, Ψ => Neumann(); exchange=Ψ)) - return -end diff --git a/src/LevelSets/compute_levelsets.jl b/src/LevelSets/compute_levelsets.jl new file mode 100644 index 00000000..aa72a891 --- /dev/null +++ b/src/LevelSets/compute_levelsets.jl @@ -0,0 +1,42 @@ +""" + init_levelset!(Ψ::Field, dem::AbstractField, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R) + +Initialize level sets. +""" +@kernel inbounds = true function init_levelset!(Ψ::Field, dem::AbstractField, dem_grid::UniformGrid, Ψ_grid::UniformGrid, cutoff, R, O) + I = @index(Global, NTuple) + I = I + O + x, y, z = coord(Ψ_grid, location(Ψ), I...) + P = R * Point3(x, y, z) + ud, sgn = sd_dem(P, cutoff, dem, dem_grid) + Ψ[I...] = ud * sgn +end + +@kernel inbounds = true function invert_levelset!(Ψ::Field, O) + I = @index(Global, NTuple) + I = I + O + Ψ[I...] *= -1.0 +end + +""" + compute_levelset_from_dem!(arch::Architecture, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I) + +Compute level sets from dem. +""" +function compute_levelset_from_dem!(arch::Architecture, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I; outer_width=(16, 8, 4)) + cutoff = 4maximum(spacing(grid)) + launch = Launcher(arch, grid; outer_width=outer_width) + launch(arch, grid, init_levelset! => (Ψ, dem, dem_grid2D, grid, cutoff, R); bc=batch(grid, Ψ => Neumann(); exchange=Ψ)) + return +end + +""" + invert_levelset!(arch::Architecture, Ψ::Field, grid::UniformGrid) + +Invert level set `Ψ` to set what's below the surface as "inside". +""" +function invert_levelset!(arch::Architecture, Ψ::Field, grid::UniformGrid) + launch = Launcher(arch, grid) + launch(arch, grid, invert_levelset! => (Ψ,)) + return +end diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl index 17b2b03d..0eac3ac8 100644 --- a/src/LevelSets/compute_volume_fractions.jl +++ b/src/LevelSets/compute_volume_fractions.jl @@ -137,12 +137,12 @@ end include("volume_fractions_kernels.jl") """ - compute_volfrac_from_level_set!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) + compute_volfrac_from_levelset!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) Compute volume fractions from level sets. """ -function compute_volfrac_from_level_set!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) +function compute_volfrac_from_levelset!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) launch = Launcher(arch, grid) - launch(arch, grid, compute_volfrac_from_level_set! => (wt, Ψ, grid)) + launch(arch, grid, compute_volfrac_from_levelset! => (wt, Ψ, grid)) return end diff --git a/src/LevelSets/volume_fractions_kernels.jl b/src/LevelSets/volume_fractions_kernels.jl index 472f8f0b..e8782d87 100644 --- a/src/LevelSets/volume_fractions_kernels.jl +++ b/src/LevelSets/volume_fractions_kernels.jl @@ -1,5 +1,5 @@ # Volume fraction kernel (2D) -@kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{2}, O) +@kernel inbounds = true function compute_volfrac_from_levelset!(wt, Ψ, grid::UniformGrid{2}, O) I = @index(Global, NTuple) ix, iy = I + O dx, dy = spacing(grid) @@ -23,7 +23,7 @@ end # Volume fraction kernel (3D) -@kernel inbounds = true function compute_volfrac_from_level_set!(wt, Ψ, grid::UniformGrid{3}, O) +@kernel inbounds = true function compute_volfrac_from_levelset!(wt, Ψ, grid::UniformGrid{3}, O) I = @index(Global, NTuple) ix, iy, iz = I + O dx, dy, dz = spacing(grid) From fa8f0b1788a177974c1163c39de2ee21a39af1b5 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 13 Mar 2024 00:36:38 +0100 Subject: [PATCH 28/32] Fixup --- .../test_levelsets_volumefractions_mpi.jl | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl index 77ae2c5f..fc5ffadb 100644 --- a/scripts_future_API/test_levelsets_volumefractions_mpi.jl +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -6,12 +6,12 @@ using FastIce.LevelSets using KernelAbstractions using CairoMakie -# using CUDA -using AMDGPU +using CUDA +# using AMDGPU # Select backend -backend = CPU() -# backend = CUDABackend() +# backend = CPU() +backend = CUDABackend() # backend = ROCBackend() using Chmy.Distributed @@ -73,13 +73,13 @@ function main(backend=CPU(); nxyz_l=(126, 126)) gather!(arch, wt_ns_c, wt.ns.c) # visualise if me == 0 - dem_surf_ = Field(backend, grid_2d, Vertex()) - dem_surf_ .= dem_surf # materialise function field for NaN filtering - dem_surf_[dem_surf. Array; colormap=:turbo), - p1_ = surface!(axs.ax1, x_g, y_g, interior(dem_surf_) |> Array; colormap=:turbo), - p2 = plot!(axs.ax2, x_g, interior(dem_bed)[:, sly] |> Array), - p2_ = plot!(axs.ax2, x_g, interior(dem_surf_)[:, sly] |> Array), + plt = (p1 = surface!(axs.ax1, x_g, y_g, dem_bed; colormap=:turbo), + p1_ = surface!(axs.ax1, x_g, y_g, dem_surf; colormap=:turbo), + p2 = plot!(axs.ax2, x_g, dem_bed[:, sly] |> Array), + p2_ = plot!(axs.ax2, x_g, dem_surf[:, sly] |> Array), p3 = heatmap!(axs.ax3, wt_na_c[:, sly, :]; colormap=:turbo), p4 = heatmap!(axs.ax4, wt_ns_c[:, sly, :]; colormap=:turbo), p5 = heatmap!(axs.ax5, wt_na_c[slx, :, :]; colormap=:turbo), From 5c13d19c9d9ac1d68b4a8c8a419eb3087c3089e9 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 13 Mar 2024 08:47:44 +0100 Subject: [PATCH 29/32] Update --- scripts_future_API/test_levelsets_volumefractions_mpi.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl index fc5ffadb..da568773 100644 --- a/scripts_future_API/test_levelsets_volumefractions_mpi.jl +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -10,8 +10,8 @@ using CUDA # using AMDGPU # Select backend -# backend = CPU() -backend = CUDABackend() +backend = CPU() +# backend = CUDABackend() # backend = ROCBackend() using Chmy.Distributed From dccaed3f9a1add9dc84f275d0c238a8335227da3 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 13 Mar 2024 13:58:33 +0100 Subject: [PATCH 30/32] Update --- .../test_levelsets_volumefractions_mpi.jl | 41 +++++++++++-------- src/LevelSets/LevelSets.jl | 1 - src/LevelSets/compute_levelsets.jl | 10 ++--- src/LevelSets/compute_volume_fractions.jl | 5 +-- 4 files changed, 29 insertions(+), 28 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl index da568773..aeff7c29 100644 --- a/scripts_future_API/test_levelsets_volumefractions_mpi.jl +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -1,18 +1,20 @@ using Chmy.Architectures using Chmy.Grids using Chmy.Fields +using Chmy.KernelLaunch using FastIce.LevelSets using KernelAbstractions using CairoMakie -using CUDA -# using AMDGPU +# using CUDA +using AMDGPU +AMDGPU.allowscalar(false) # Select backend -backend = CPU() +# backend = CPU() # backend = CUDABackend() -# backend = ROCBackend() +backend = ROCBackend() using Chmy.Distributed using MPI @@ -23,13 +25,15 @@ function main(backend=CPU(); nxyz_l=(126, 126)) topo = topology(arch) me = global_rank(topo) # geometry - zmin, zmax = 0.0, 1.0 - lx, ly, lz = 5.0, 5.0, zmax - zmin - nx, ny = nxyz_l[1:2] - nz = (length(nxyz_l) < 3) ? ceil(Int, lz / lx * nxyz_l[1]) : nxyz_l[3] - dims_l = (nx, ny, nz) - dims_g = dims_l .* dims(topo) - grid = UniformGrid(arch; origin=(-lx/2, -ly/2, 0.0), extent=(lx, ly, lz), dims=dims_g) + zmin, zmax = 0.0, 1.0 + lx, ly, lz = 5.0, 5.0, zmax - zmin + nx, ny = nxyz_l[1:2] + nz = (length(nxyz_l) < 3) ? ceil(Int, lz / lx * nxyz_l[1]) : nxyz_l[3] + @info "nz = $nz (lz = $lz)" + dims_l = (nx, ny, nz) + dims_g = dims_l .* dims(topo) + grid = UniformGrid(arch; origin=(-lx/2, -ly/2, 0.0), extent=(lx, ly, lz), dims=dims_g) + launch = Launcher(arch, grid; outer_width=(16, 8, 4)) # synthetic topo ############### arch_2d = Arch(backend) @@ -54,13 +58,13 @@ function main(backend=CPU(); nxyz_l=(126, 126)) wt = (na=volfrac_field(backend, grid), ns=volfrac_field(backend, grid)) # comput level set - compute_levelset_from_dem!(arch, Ψ.na, dem_surf, grid_2d, grid) - compute_levelset_from_dem!(arch, Ψ.ns, dem_bed, grid_2d, grid) - invert_levelset!(arch, Ψ.ns, grid) + compute_levelset_from_dem!(arch, launch, Ψ.na, dem_surf, grid_2d, grid) + compute_levelset_from_dem!(arch, launch, Ψ.ns, dem_bed, grid_2d, grid) + invert_levelset!(arch, launch, Ψ.ns, grid) # @. Ψ.ns *= -1.0 # invert level set to set what's below the DEM surface as inside # volume fractions for phase in eachindex(Ψ) - compute_volfrac_from_levelset!(arch, wt[phase], Ψ[phase], grid) + compute_volfrac_from_levelset!(arch, launch, wt[phase], Ψ[phase], grid) end # compute physics or else @@ -73,7 +77,7 @@ function main(backend=CPU(); nxyz_l=(126, 126)) gather!(arch, wt_ns_c, wt.ns.c) # visualise if me == 0 - dem_bed = Array(interior(dem_bed)) + dem_bed = Array(interior(dem_bed)) dem_surf = Array(interior(dem_surf)) dem_surf[dem_surf .< dem_bed] .= NaN slx = ceil(Int, size(wt_na_c, 1) / 2) # for visu @@ -99,11 +103,12 @@ function main(backend=CPU(); nxyz_l=(126, 126)) Colorbar(fig[1, 1][1, 2], plt.p1) Colorbar(fig[2, 2][1, 2], plt.p4) # display(fig) - save("levset.png", fig) + save("levset_$nx.png", fig) end return end -main(backend; nxyz_l=(126, 126)) +n = 126 +main(backend; nxyz_l=(n, n)) MPI.Finalize() diff --git a/src/LevelSets/LevelSets.jl b/src/LevelSets/LevelSets.jl index d4ff52fc..98936c71 100644 --- a/src/LevelSets/LevelSets.jl +++ b/src/LevelSets/LevelSets.jl @@ -7,7 +7,6 @@ using Chmy.Architectures using Chmy.BoundaryConditions using Chmy.Grids using Chmy.Fields -using Chmy.KernelLaunch using KernelAbstractions using LinearAlgebra, GeometryBasics diff --git a/src/LevelSets/compute_levelsets.jl b/src/LevelSets/compute_levelsets.jl index aa72a891..0023e8d1 100644 --- a/src/LevelSets/compute_levelsets.jl +++ b/src/LevelSets/compute_levelsets.jl @@ -19,24 +19,22 @@ end end """ - compute_levelset_from_dem!(arch::Architecture, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I) + compute_levelset_from_dem!(arch::Architecture, launch, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I) Compute level sets from dem. """ -function compute_levelset_from_dem!(arch::Architecture, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I; outer_width=(16, 8, 4)) +function compute_levelset_from_dem!(arch::Architecture, launch, Ψ::Field, dem::AbstractField, dem_grid2D::UniformGrid, grid::UniformGrid, R=LinearAlgebra.I) cutoff = 4maximum(spacing(grid)) - launch = Launcher(arch, grid; outer_width=outer_width) launch(arch, grid, init_levelset! => (Ψ, dem, dem_grid2D, grid, cutoff, R); bc=batch(grid, Ψ => Neumann(); exchange=Ψ)) return end """ - invert_levelset!(arch::Architecture, Ψ::Field, grid::UniformGrid) + invert_levelset!(arch::Architecture, launch, Ψ::Field, grid::UniformGrid) Invert level set `Ψ` to set what's below the surface as "inside". """ -function invert_levelset!(arch::Architecture, Ψ::Field, grid::UniformGrid) - launch = Launcher(arch, grid) +function invert_levelset!(arch::Architecture, launch, Ψ::Field, grid::UniformGrid) launch(arch, grid, invert_levelset! => (Ψ,)) return end diff --git a/src/LevelSets/compute_volume_fractions.jl b/src/LevelSets/compute_volume_fractions.jl index 0eac3ac8..1b6b44b5 100644 --- a/src/LevelSets/compute_volume_fractions.jl +++ b/src/LevelSets/compute_volume_fractions.jl @@ -137,12 +137,11 @@ end include("volume_fractions_kernels.jl") """ - compute_volfrac_from_levelset!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) + compute_volfrac_from_levelset!(arch::Architecture, launch, wt, Ψ::Field, grid::UniformGrid) Compute volume fractions from level sets. """ -function compute_volfrac_from_levelset!(arch::Architecture, wt, Ψ::Field, grid::UniformGrid) - launch = Launcher(arch, grid) +function compute_volfrac_from_levelset!(arch::Architecture, launch, wt, Ψ::Field, grid::UniformGrid) launch(arch, grid, compute_volfrac_from_levelset! => (wt, Ψ, grid)) return end From a326b2d3f25dbe8238f6c10462cb1adacb7877c7 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Wed, 13 Mar 2024 22:44:11 +0100 Subject: [PATCH 31/32] rework workflow --- .../test_levelsets_volumefractions_mpi.jl | 89 +++++++++++-------- 1 file changed, 50 insertions(+), 39 deletions(-) diff --git a/scripts_future_API/test_levelsets_volumefractions_mpi.jl b/scripts_future_API/test_levelsets_volumefractions_mpi.jl index aeff7c29..f3491ec8 100644 --- a/scripts_future_API/test_levelsets_volumefractions_mpi.jl +++ b/scripts_future_API/test_levelsets_volumefractions_mpi.jl @@ -8,60 +8,73 @@ using FastIce.LevelSets using KernelAbstractions using CairoMakie # using CUDA -using AMDGPU -AMDGPU.allowscalar(false) +# using AMDGPU +# AMDGPU.allowscalar(false) # Select backend -# backend = CPU() +backend = CPU() # backend = CUDABackend() -backend = ROCBackend() +# backend = ROCBackend() using Chmy.Distributed using MPI MPI.Init() -function main(backend=CPU(); nxyz_l=(126, 126)) - arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) - topo = topology(arch) - me = global_rank(topo) - # geometry - zmin, zmax = 0.0, 1.0 - lx, ly, lz = 5.0, 5.0, zmax - zmin - nx, ny = nxyz_l[1:2] - nz = (length(nxyz_l) < 3) ? ceil(Int, lz / lx * nxyz_l[1]) : nxyz_l[3] - @info "nz = $nz (lz = $lz)" - dims_l = (nx, ny, nz) - dims_g = dims_l .* dims(topo) - grid = UniformGrid(arch; origin=(-lx/2, -ly/2, 0.0), extent=(lx, ly, lz), dims=dims_g) - launch = Launcher(arch, grid; outer_width=(16, 8, 4)) - - # synthetic topo ############### +@views function make_dem(backend=CPU(); nx, ny, lx, ly, zmin, zmax, amp, ω, tanβ, el, gl) arch_2d = Arch(backend) - grid_2d = UniformGrid(arch_2d; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=dims_g[1:2]) + grid_2d = UniformGrid(arch_2d; origin=(-lx/2, -ly/2), extent=(lx, ly), dims=(nx, ny)) + + lz = zmax - zmin + + # type = :turtle + generate_surf(x, y, gl, lx, ly) = gl * (1.0 - ((1.7 * x / lx + 0.22)^2 + (0.5 * y / ly)^2)) + generate_bed(x, y, amp, ω, tanβ, el, lx, ly, lz) = lz * (amp * sin(ω * x / lx) * sin(ω * y /ly) + tanβ * x / lx + el + (1.5 * y / ly)^2) + + surf = FunctionField(generate_surf, grid_2d, Vertex(); parameters=(; gl, lx, ly)) + bed = FunctionField(generate_bed, grid_2d, Vertex(); parameters=(; amp, ω, tanβ, el, lx, ly, lz)) + + return (; bed, surf, grid_2d, lz) +end + +@views function main(backend=CPU()) + lx, ly = 5.0, 5.0 + zmin, zmax = 0.0, 1.0 + nx = ny = 126 # synthetic topo amp = 0.1 ω = 10π tanβ = tan(-π/6) el = 0.35 gl = 0.9 - # synthetic bed and surf - # type = :turtle - generate_surf(x, y, gl, lx, ly) = gl * (1.0 - ((1.7 * x / lx + 0.22)^2 + (0.5 * y / ly)^2)) - generate_bed(x, y, amp, ω, tanβ, el, lx, ly, lz) = lz * (amp * sin(ω * x / lx) * sin(ω * y /ly) + tanβ * x / lx + el + (1.5 * y / ly)^2) - dem_surf = FunctionField(generate_surf, grid_2d, Vertex(); parameters=(; gl, lx, ly)) - dem_bed = FunctionField(generate_bed, grid_2d, Vertex(); parameters=(; amp, ω, tanβ, el, lx, ly, lz)) - # synthetic topo ############### + dem_data = make_dem(backend; nx, ny, lx, ly, zmin, zmax, amp, ω, tanβ, el, gl) + + run_simulation(backend; nxyz_l=(nx, ny), dem_data) + return +end + +@views function run_simulation(backend=CPU(); nxyz_l=(126, 126), dem_data) + arch = Arch(backend, MPI.COMM_WORLD, (0, 0, 0)) + topo = topology(arch) + me = global_rank(topo) + # geometry + lx, ly, lz = extent(dem_data.grid_2d, Vertex())..., dem_data.lz + nx, ny = nxyz_l[1:2] + nz = (length(nxyz_l) < 3) ? ceil(Int, lz / lx * nxyz_l[1]) : nxyz_l[3] + @info "nz = $nz (lz = $lz)" + dims_l = (nx, ny, nz) + dims_g = dims_l .* dims(topo) + grid = UniformGrid(arch; origin=(-lx/2, -ly/2, 0.0), extent=(lx, ly, lz), dims=dims_g) + launch = Launcher(arch, grid; outer_width=(16, 8, 4)) # init fields Ψ = (na=Field(backend, grid, Vertex()), ns=Field(backend, grid, Vertex())) wt = (na=volfrac_field(backend, grid), ns=volfrac_field(backend, grid)) # comput level set - compute_levelset_from_dem!(arch, launch, Ψ.na, dem_surf, grid_2d, grid) - compute_levelset_from_dem!(arch, launch, Ψ.ns, dem_bed, grid_2d, grid) - invert_levelset!(arch, launch, Ψ.ns, grid) - # @. Ψ.ns *= -1.0 # invert level set to set what's below the DEM surface as inside + compute_levelset_from_dem!(arch, launch, Ψ.na, dem_data.surf, dem_data.grid_2d, grid) + compute_levelset_from_dem!(arch, launch, Ψ.ns, dem_data.bed, dem_data.grid_2d, grid) + invert_levelset!(arch, launch, Ψ.ns, grid) # invert level set to set what's below the DEM surface as inside # volume fractions for phase in eachindex(Ψ) compute_volfrac_from_levelset!(arch, launch, wt[phase], Ψ[phase], grid) @@ -70,15 +83,14 @@ function main(backend=CPU(); nxyz_l=(126, 126)) # compute physics or else # postprocessing - # gather wt_na_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.na.c)) .* dims(topo)) : nothing wt_ns_c = (me == 0) ? KernelAbstractions.zeros(CPU(), Float64, size(interior(wt.ns.c)) .* dims(topo)) : nothing gather!(arch, wt_na_c, wt.na.c) gather!(arch, wt_ns_c, wt.ns.c) # visualise if me == 0 - dem_bed = Array(interior(dem_bed)) - dem_surf = Array(interior(dem_surf)) + dem_bed = Array(interior(dem_data.bed)) + dem_surf = Array(interior(dem_data.surf)) dem_surf[dem_surf .< dem_bed] .= NaN slx = ceil(Int, size(wt_na_c, 1) / 2) # for visu sly = ceil(Int, size(wt_na_c, 2) / 2) # for visu @@ -102,13 +114,12 @@ function main(backend=CPU(); nxyz_l=(126, 126)) p6 = heatmap!(axs.ax6, wt_ns_c[slx, :, :]; colormap=:turbo)) Colorbar(fig[1, 1][1, 2], plt.p1) Colorbar(fig[2, 2][1, 2], plt.p4) - # display(fig) - save("levset_$nx.png", fig) + display(fig) + # save("levset_$nx.png", fig) end return end -n = 126 -main(backend; nxyz_l=(n, n)) +main(backend) MPI.Finalize() From 51614290825b10933b209175309a51702af9fa5e Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Tue, 19 Mar 2024 15:02:23 +0100 Subject: [PATCH 32/32] Cleanup --- .../test_levelsets_volumefractions.jl | 127 ------------------ 1 file changed, 127 deletions(-) delete mode 100644 scripts_future_API/test_levelsets_volumefractions.jl diff --git a/scripts_future_API/test_levelsets_volumefractions.jl b/scripts_future_API/test_levelsets_volumefractions.jl deleted file mode 100644 index f522a351..00000000 --- a/scripts_future_API/test_levelsets_volumefractions.jl +++ /dev/null @@ -1,127 +0,0 @@ -using Chmy.Architectures -using Chmy.Grids -using Chmy.Fields - -using FastIce.LevelSets -using FastIce.Geometries - -using KernelAbstractions - -# using CUDA -# using AMDGPU - -# Select backend -backend = CPU() -# backend = CUDABackend() -# backend = ROCBackend() - -arch = Arch(backend) - -function load_data(data_path) - dat = load(data_path) - return first(keys(dat)), dat -end - - -""" - extract_dem(arch::Architecture, data_path::String, z_resol=nothing) - -Load digital elevation model of surface and bedrock from (JLD2) file and initialise the grid. -""" -function extract_dem(arch::Architecture, data_path::String, z_resol=nothing) - dtype, data = load_data(data_path) - z_surf = data[dtype].z_surf - z_bed = data[dtype].z_bed - dm = data[dtype].domain - lx, ly, lz = dm.xmax - dm.xmin, dm.ymax - dm.ymin, dm.zmax - dm.zmin - nx, ny = size(z_surf) .- 1 - nz = isnothing(z_resol) ? ceil(Int, lz / lx * nx) : z_resol - Ψ_grid = UniformGrid(arch; origin=(dm.xmin, dm.ymin, dm.zmin), extent=(lx, ly, lz), dims=(nx, ny, nz)) - dem_grid = UniformGrid(arch; origin=(dm.xmin, dm.ymin), extent=(lx, ly), dims=(nx, ny)) - dem_surf = Field(arch.backend, dem_grid, Vertex()) - dem_bed = Field(arch.backend, dem_grid, Vertex()) - set!(dem_surf, z_surf) - set!(dem_bed, z_bed) - return dem_surf, dem_bed, dem_grid, Ψ_grid -end - -function extract_dem(arch::Architecture, data::SyntheticElevation, z_resol=nothing) - z_surf = data.z_surf - z_bed = data.z_bed - dm = data.domain - lx, ly, lz = extents(dm) - nx, ny = size(z_surf) .- 1 - nz = isnothing(z_resol) ? ceil(Int, lz / lx * nx) : z_resol - Ψ_grid = UniformGrid(arch; origin=(dm.xmin, dm.ymin, dm.zmin), extent=(lx, ly, lz), dims=(nx, ny, nz)) - dem_grid = UniformGrid(arch; origin=(dm.xmin, dm.ymin), extent=(lx, ly), dims=(nx, ny)) - dem_surf = Field(arch.backend, dem_grid, Vertex()) - dem_bed = Field(arch.backend, dem_grid, Vertex()) - set!(dem_surf, z_surf) - set!(dem_bed, z_bed) - dem = (surf=dem_surf, bed=dem_bed) - return dem, dem_grid, Ψ_grid -end - -nx, ny = 128, 128 -lx, ly = 5.0, 5.0 -zmin, zmax = 0.0, 1.0 -type = :turtle - -synth_dem = make_synthetic(nx, ny, lx, ly, zmin, zmax, type) -dem, dem_grid, Ψ_grid = extract_dem(arch, synth_dem) - -Ψ = ( - na=Field(backend, Ψ_grid, Vertex()), - ns=Field(backend, Ψ_grid, Vertex()), -) - -compute_level_set_from_dem!(arch, Ψ.na, dem.surf, dem_grid, Ψ_grid) -compute_level_set_from_dem!(arch, Ψ.ns, dem.bed, dem_grid, Ψ_grid) -# invert level set to set what's below the DEM surface as inside -@. Ψ.ns *= -1.0 - -wt = ( - na=volfrac_field(backend, Ψ_grid), - ns=volfrac_field(backend, Ψ_grid), -) - -for phase in eachindex(Ψ) - compute_volfrac_from_level_set!(arch, wt[phase], Ψ[phase], Ψ_grid) -end - - -@views function visme(synth_dem, Ψ, wt) - nx, ny, nz = size(Ψ.ns) - x = LinRange(synth_dem.domain.xmin, synth_dem.domain.xmax, nx) - y = LinRange(synth_dem.domain.ymin, synth_dem.domain.ymax, ny) - z = LinRange(synth_dem.domain.zmin, synth_dem.domain.zmax, nz) - bed = synth_dem.z_bed - surf = synth_dem.z_surf - - surf2 = copy(surf) - surf2[surf.