Fixup

.JuliaFormatter.toml

@@ -1,5 +1,5 @@
 style = "yas"
-margin = 140
+margin = 165
 align_assignment = true
 align_conditional = true
 whitespace_ops_in_indices = false

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

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)

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
 

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
+end

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
+    return
+end

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
-

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