From 9a9694cde251ebea663306804217fbf0121202a2 Mon Sep 17 00:00:00 2001 From: Ludovic Raess Date: Tue, 5 Mar 2024 23:52:38 +0100 Subject: [PATCH] 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