Merge branch 'main' into multiscale

AlgebraicJulia · Oct 14, 2024 · 61cce2d · 61cce2d
2 parents 609f330 + d0a896f
commit 61cce2d
Show file tree

Hide file tree

Showing 11 changed files with 657 additions and 774 deletions.
diff --git a/.buildkite/jobscript.sh b/.buildkite/jobscript.sh
@@ -8,4 +8,5 @@ echo "Running Tests..."
 julia --project -t 16 -e 'using Pkg; Pkg.status(); Pkg.test()'
 
 echo "Building Documentation..."
+
 julia --project=docs -t 16 -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.status(); Pkg.instantiate(); include("docs/make.jl")'
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -1,24 +1,21 @@
 env:
   JULIA_VERSION: "1.10.2"
-  GATAS_HOME: "~/.gatas/buildkite/agents/$BUILDKITE_AGENT_NAME"
+  JULIA_DEPOT_PATH: "$DEPOT"
 
 steps:
 
-  - label: ":hammer: Build Project"
-    env:
-      JULIA_DEPOT_PATH: "$GATAS_HOME"
-    command: 
-      - "module load julia"
-      - "julia --project=docs --color=yes -e 'using Pkg; Pkg.update(); Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate(); Pkg.precompile()'"
-
-  - wait 
+  - label: ":arrow_down: Load AlgebraicJulia pipeline"
+    command: |
+      curl -s https://raw.githubusercontent.com/AlgebraicJulia/.github/main/buildkite/pipeline.yml | buildkite-agent pipeline upload
 
-  - label: ":scroll: Build docs and run tests"
-    env:
-      JULIA_DEPOT_PATH: "$GATAS_HOME"
-      JULIA_PROJECT: "docs/"
+  - wait
+
+  - label: ":racehorse: Build docs and run tests on A100 GPU partition"
     command:
-      - "srun --cpus-per-task=16 --mem=64G --time=1:00:00 --output=.buildkite/log_%j.log --unbuffered .buildkite/jobscript.sh"
+      - "srun --cpus-per-task=16 --mem=16G  --gres=gpu:a100:1 --time=1:00:00 -p gpu --output=.buildkite/log_gpu_a100_%j.log --unbuffered .buildkite/jobscript.sh"
 
-  - wait
+  - label: ":racehorse: Build docs and run tests on Quadro GPU partition"
+    command:
+      - "srun --cpus-per-task=16 --mem=16G  --gres=gpu:geforce:1 --time=1:00:00 -p gpu --output=.buildkite/log_gpu_quadro_%j.log --unbuffered .buildkite/jobscript.sh"
 
+  - wait
diff --git a/Project.toml b/Project.toml
@@ -16,6 +16,7 @@ Debugger = "31a5f54b-26ea-5ae9-a837-f05ce5417438"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -45,6 +46,8 @@ DataMigrations = "0.0.3, 0.1"
 FileIO = "^1"
 GeometryBasics = "0.3, 0.4"
 JSON = "0.21.4"
+KernelAbstractions = "0.9"
+Krylov = "0.9.6"
 LazyArrays = "0.20, 0.21, 0.22, 1.0, 2"
 LinearAlgebra = "1.9"
 Makie = "0.21"

diff --git a/ext/CombinatorialSpacesCUDAExt.jl b/ext/CombinatorialSpacesCUDAExt.jl
@@ -1,187 +1,94 @@
 module CombinatorialSpacesCUDAExt
 
 using CombinatorialSpaces
-using GeometryBasics
+using CombinatorialSpaces.DiscreteExteriorCalculus: DiscreteHodge
 using CUDA
 using CUDA.CUSPARSE
-using LinearAlgebra
-using SparseArrays
+using KernelAbstractions
 using Krylov
-Point2D = Point2{Float64}
-Point3D = Point3{Float64}
-import CombinatorialSpaces: dec_wedge_product, dec_c_wedge_product, dec_c_wedge_product!, dec_p_wedge_product, 
-dec_boundary, dec_differential, dec_dual_derivative, dec_hodge_star, dec_inv_hodge_star
+import CombinatorialSpaces: cache_wedge,
+  dec_boundary, dec_differential, dec_dual_derivative,
+  dec_hodge_star, dec_inv_hodge_star
 
-""" dec_wedge_product(::Type{Tuple{0,0}}, sd::HasDeltaSet, ::Type{Val{:CUDA}})
+# Wedge Product
+#--------------
 
-Wedge Primal 0, 0 for CUDA
-"""
-function dec_wedge_product(::Type{Tuple{0,0}}, sd::HasDeltaSet, ::Type{Val{:CUDA}})
-  (f, g) -> f .* g
+function cache_wedge(::Type{Tuple{m,n}}, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p}, ::Type{Val{:CUDA}}, arr_cons=CuArray, cast_float=nothing) where {float_type,_p,m,n}
+  cache_wedge(m, n, sd, float_type, arr_cons, cast_float)
 end
-
-""" function dec_wedge_product(::Type{Tuple{k,0}}, sd::HasDeltaSet, ::Type{Val{:CUDA}})
-
-Wedge Primal K, 0 for CUDA
-"""
-function dec_wedge_product(::Type{Tuple{k,0}}, sd::HasDeltaSet, ::Type{Val{:CUDA}}) where {k}
-  val_pack = dec_p_wedge_product(Tuple{0,k}, sd, Val{:CUDA})
-  (α, g) -> dec_c_wedge_product(Tuple{0,k}, g, α, val_pack, Val{:CUDA})
+function cache_wedge(::Type{Tuple{m,n}}, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p}, ::Type{Val{:CUDA}}, arr_cons=CuArray, cast_float=nothing) where {float_type,_p,m,n}
+  cache_wedge(m, n, sd, float_type, arr_cons, cast_float)
 end
 
-""" function dec_wedge_product(::Type{Tuple{0,k}}, sd::HasDeltaSet, ::Type{Val{:CUDA}})
+# Boundary and Co-boundary
+#-------------------------
 
-Wedge Primal 0, K for CUDA
-"""
-function dec_wedge_product(::Type{Tuple{0,k}}, sd::HasDeltaSet, ::Type{Val{:CUDA}}) where {k}
-  val_pack = dec_p_wedge_product(Tuple{0,k}, sd, Val{:CUDA})
-  (f, β) -> dec_c_wedge_product(Tuple{0,k}, f, β, val_pack, Val{:CUDA})
-end
+"""    dec_boundary(n::Int, sd::HasDeltaSet, ::Type{Val{:CUDA}})
 
-""" function dec_wedge_product(::Type{Tuple{1,1}}, sd::HasDeltaSet2D, ::Type{Val{:CUDA}})
-
-Wedge Primal 1, 1 for CUDA
+Compute a boundary matrix as a sparse CUDA matrix.
 """
-function dec_wedge_product(::Type{Tuple{1,1}}, sd::HasDeltaSet2D, ::Type{Val{:CUDA}})
-  val_pack = dec_p_wedge_product(Tuple{1,1}, sd, Val{:CUDA})
-  (α, β) -> dec_c_wedge_product(Tuple{1,1}, α, β, val_pack, Val{:CUDA})
-end
+dec_boundary(n::Int, sd::HasDeltaSet, ::Type{Val{:CUDA}}) =
+  CuSparseMatrixCSC(dec_boundary(n, sd))
 
-""" function dec_p_wedge_product(::Type{Tuple{m,n}}, sd, ::Type{Val{:CUDA}})
+dec_boundary(n::Int, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type =
+  CuSparseMatrixCSC{float_type}(dec_boundary(n, sd))
+dec_boundary(n::Int, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type =
+  CuSparseMatrixCSC{float_type}(dec_boundary(n, sd))
 
-Preallocation for CUDA Wedges
-"""
-function dec_p_wedge_product(::Type{Tuple{m,n}}, sd, ::Type{Val{:CUDA}}) where {m,n}
-  val_pack = dec_p_wedge_product(Tuple{m,n}, sd)
-  (CuArray.(val_pack[1:end-1])..., val_pack[end]) 
-end
+"""    dec_dual_derivative(n::Int, sd::EmbeddedDeltaDualComplex1D, ::Type{Val{:CUDA}})
 
-""" function dec_c_wedge_product(::Type{Tuple{m,n}}, α, β, val_pack, ::Type{Val{:CUDA}})
-
-Computation for CUDA Wedges
+Compute a dual derivative matrix as a sparse CUDA matrix.
 """
-function dec_c_wedge_product(::Type{Tuple{m,n}}, α, β, val_pack, ::Type{Val{:CUDA}}) where {m,n}
-  wedge_terms = CUDA.zeros(eltype(α), last(last(val_pack)))
-  return dec_c_wedge_product!(Tuple{m,n}, wedge_terms, α, β, val_pack, Val{:CUDA})
-end
-
-function dec_c_wedge_product!(::Type{Tuple{0,1}}, wedge_terms, f, α, val_pack, ::Type{Val{:CUDA}})
-  num_threads = CUDA.max_block_size.x
-  num_blocks = min(ceil(Int, length(wedge_terms) / num_threads), CUDA.max_grid_size.x)
-
-  @cuda threads=num_threads blocks=num_blocks dec_cu_ker_c_wedge_product_01!(wedge_terms, f, α, val_pack[1])
-  return wedge_terms
-end
-
-function dec_c_wedge_product!(::Type{Tuple{0,2}}, wedge_terms, f, α, val_pack, ::Type{Val{:CUDA}})
-  num_threads = CUDA.max_block_size.x
-  num_blocks = min(ceil(Int, length(wedge_terms) / num_threads), CUDA.max_grid_size.x)
+dec_dual_derivative(n::Int, sd::HasDeltaSet, ::Type{Val{:CUDA}}) =
+  CuSparseMatrixCSC(dec_dual_derivative(n, sd))
 
-  @cuda threads=num_threads blocks=num_blocks dec_cu_ker_c_wedge_product_02!(wedge_terms, f, α, val_pack[1], val_pack[2])
-  return wedge_terms
-end
-
-function dec_c_wedge_product!(::Type{Tuple{1,1}}, wedge_terms, f, α, val_pack, ::Type{Val{:CUDA}})
-  num_threads = CUDA.max_block_size.x
-  num_blocks = min(ceil(Int, length(wedge_terms) / num_threads), CUDA.max_grid_size.x)
-
-  @cuda threads=num_threads blocks=num_blocks dec_cu_ker_c_wedge_product_11!(wedge_terms, f, α, val_pack[1], val_pack[2])
-  return wedge_terms
-end
-
-function dec_cu_ker_c_wedge_product_01!(wedge_terms::CuDeviceArray{T}, f, α, primal_vertices) where T
-  index = (blockIdx().x - Int32(1)) * blockDim().x + threadIdx().x   
-  stride = gridDim().x * blockDim().x
-  i = index
-  @inbounds while i <= Int32(length(wedge_terms))
-    wedge_terms[i] = T(0.5) * α[i] * (f[primal_vertices[i, Int32(1)]] + f[primal_vertices[i, Int32(2)]])
-    i += stride
-  end
-  return nothing
-end
-
-function dec_cu_ker_c_wedge_product_02!(wedge_terms::CuDeviceArray{T}, f, α, pv, coeffs) where T
-  index = (blockIdx().x - Int32(1)) * blockDim().x + threadIdx().x   
-  stride = gridDim().x * blockDim().x
-  i = index
-
-  @inbounds while i <= Int32(length(wedge_terms))
-    wedge_terms[i] = α[i] * (coeffs[Int32(1), i] * f[pv[Int32(1), i]] + coeffs[Int32(2), i] * f[pv[Int32(2), i]]
-                                + coeffs[Int32(3), i] * f[pv[Int32(3), i]] + coeffs[Int32(4), i] * f[pv[Int32(4), i]]
-                                + coeffs[Int32(5), i] * f[pv[Int32(5), i]] + coeffs[Int32(6), i] * f[pv[Int32(6), i]])
-    i += stride
-  end
-  return nothing
-end
+dec_dual_derivative(n::Int, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type =
+  CuSparseMatrixCSC{float_type}(dec_dual_derivative(n, sd))
+dec_dual_derivative(n::Int, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type =
+  CuSparseMatrixCSC{float_type}(dec_dual_derivative(n, sd))
 
-function dec_cu_ker_c_wedge_product_11!(wedge_terms, α, β, e, coeffs)
-  index = (blockIdx().x - Int32(1)) * blockDim().x + threadIdx().x   
-  stride = gridDim().x * blockDim().x
-  i = index
 
-  @inbounds while i <= Int32(length(wedge_terms))
-    e0, e1, e2 = e[Int32(1), i], e[Int32(2), i], e[Int32(3), i]
-    ae0, ae1, ae2 = α[e0], α[e1], α[e2]
-    be0, be1, be2 = β[e0], β[e1], β[e2]
+"""    dec_differential(n::Int, sd::HasDeltaSet, ::Type{Val{:CUDA}})
 
-    c1, c2, c3 = coeffs[Int32(1), i], coeffs[Int32(2), i], coeffs[Int32(3), i]
-
-    wedge_terms[i] = (c1 * (ae2 * be1 - ae1 * be2)
-                    + c2 * (ae2 * be0 - ae0 * be2)
-                    + c3 * (ae1 * be0 - ae0 * be1))
-    i += stride
-  end
-
-  return nothing
-end
-
-""" dec_boundary(n::Int, sd::EmbeddedDeltaDualComplex, ::Type{Val{:CUDA}})
-
-Boundary matrix as a sparse CUDA matrix
+Compute an exterior derivative matrix as a sparse CUDA matrix.
 """
-dec_boundary(n::Int, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type = CuSparseMatrixCSC{float_type}(dec_boundary(n, sd))
+dec_differential(n::Int, sd::HasDeltaSet, ::Type{Val{:CUDA}}) =
+  CuSparseMatrixCSC(dec_differential(n, sd))
 
-""" dec_dual_derivative(n::Int, sd::EmbeddedDeltaDualComplex, ::Type{Val{:CUDA}})
+dec_differential(n::Int, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type =
+  CuSparseMatrixCSC{float_type}(dec_differential(n, sd))
+dec_differential(n::Int, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type =
+  CuSparseMatrixCSC{float_type}(dec_differential(n, sd))
 
-Dual derivative matrix as a sparse CUDA matrix
-"""
-dec_dual_derivative(n::Int, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type = CuSparseMatrixCSC{float_type}(dec_dual_derivative(n, sd))
+# Hodge Star
+#-----------
 
-""" dec_differential(n::Int, sd::EmbeddedDeltaDualComplex, ::Type{Val{:CUDA}})
+"""    dec_hodge_star(n::Int, sd::HasDeltaSet, h::DiscreteHodge, ::Type{Val{:CUDA}})
 
-Exterior derivative matrix as a sparse CUDA matrix
+Compute a Hodge star as a diagonal or generic sparse CUDA matrix.
 """
-dec_differential(n::Int, sd::EmbeddedDeltaDualComplex1D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type = CuSparseMatrixCSC{float_type}(dec_differential(n, sd))
+dec_hodge_star(n::Int, sd::HasDeltaSet, h::DiscreteHodge, ::Type{Val{:CUDA}}) = 
+  dec_hodge_star(Val{n}, sd, h, Val{:CUDA})
 
-dec_boundary(n::Int, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type = CuSparseMatrixCSC{float_type}(dec_boundary(n, sd))
-dec_dual_derivative(n::Int, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type = CuSparseMatrixCSC{float_type}(dec_dual_derivative(n, sd))
-dec_differential(n::Int, sd::EmbeddedDeltaDualComplex2D{Bool, float_type, _p} where _p, ::Type{Val{:CUDA}}) where float_type = CuSparseMatrixCSC{float_type}(dec_differential(n, sd))
+dec_hodge_star(::Type{Val{n}}, sd::HasDeltaSet, h::DiscreteHodge, ::Type{Val{:CUDA}}) where n =
+  CuArray(dec_hodge_star(Val{n}, sd, h))
 
-dec_boundary(n::Int, sd::HasDeltaSet, ::Type{Val{:CUDA}}) = CuSparseMatrixCSC(dec_boundary(n, sd))
-dec_dual_derivative(n::Int, sd::HasDeltaSet, ::Type{Val{:CUDA}}) = CuSparseMatrixCSC(dec_dual_derivative(n, sd))
-dec_differential(n::Int, sd::HasDeltaSet, ::Type{Val{:CUDA}}) = CuSparseMatrixCSC(dec_differential(n, sd))
+dec_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D, h::GeometricHodge, ::Type{Val{:CUDA}}) =
+  CuSparseMatrixCSC(dec_hodge_star(Val{1}, sd, h))
 
-""" dec_hodge_star(n::Int, sd::HasDeltaSet, ::HodgeType, ::Type{Val{:CUDA}})
+"""    dec_inv_hodge_star(n::Int, sd::HasDeltaSet, h::DiscreteHodge, ::Type{Val{:CUDA}})
 
-Hodge matrices as appropriate diagonal or sparse CUDA matrices
+Compute an inverse Hodge star matrix as a diagonal CUDA matrix.
 """
-dec_hodge_star(n::Int, sd::HasDeltaSet, ::DiagonalHodge, ::Type{Val{:CUDA}}) = CuArray(dec_hodge_star(n, sd, DiagonalHodge()))
-dec_hodge_star(n::Int, sd::HasDeltaSet, ::GeometricHodge, ::Type{Val{:CUDA}}) = dec_hodge_star(Val{n}, sd, GeometricHodge(), Val{:CUDA})
-dec_hodge_star(::Type{Val{n}}, sd::HasDeltaSet, ::GeometricHodge, ::Type{Val{:CUDA}}) where n = CuArray(dec_hodge_star(Val{n}, sd, GeometricHodge()))
-dec_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge, ::Type{Val{:CUDA}}) = CuSparseMatrixCSC(dec_hodge_star(Val{1}, sd, GeometricHodge()))
+dec_inv_hodge_star(n::Int, sd::HasDeltaSet, h::DiscreteHodge, ::Type{Val{:CUDA}}) =
+  dec_inv_hodge_star(Val{n}, sd, h, Val{:CUDA})
 
-""" dec_inv_hodge_star(n::Int, sd::HasDeltaSet, ::HodgeType, ::Type{Val{:CUDA}})
-
-Inverse Hodge matrices as diagonal matrices
-"""
-dec_inv_hodge_star(n::Int, sd::HasDeltaSet, ::DiagonalHodge, ::Type{Val{:CUDA}}) = CuArray(dec_inv_hodge_star(n, sd, DiagonalHodge()))
-dec_inv_hodge_star(n::Int, sd::HasDeltaSet, ::GeometricHodge, ::Type{Val{:CUDA}}) = dec_inv_hodge_star(Val{n}, sd, GeometricHodge(), Val{:CUDA})
-dec_inv_hodge_star(::Type{Val{n}}, sd::HasDeltaSet, ::GeometricHodge, ::Type{Val{:CUDA}}) where n = CuArray(dec_inv_hodge_star(Val{n}, sd, GeometricHodge()))
+dec_inv_hodge_star(::Type{Val{n}}, sd::HasDeltaSet, h::DiscreteHodge, ::Type{Val{:CUDA}}) where n =
+  CuArray(dec_inv_hodge_star(Val{n}, sd, h))
 
-""" function dec_inv_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge, ::Type{Val{:CUDA}})
+"""    function dec_inv_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge, ::Type{Val{:CUDA}})
 
-Inverse Geometric Hodge for Primal 1 implemented as a GMRES solver. This returns a function that
-will compute the result.
+Return a function that computes the inverse geometric Hodge star of a primal 1-form via a GMRES solver.
 """
 function dec_inv_hodge_star(::Type{Val{1}}, sd::EmbeddedDeltaDualComplex2D, ::GeometricHodge, ::Type{Val{:CUDA}})
   hdg = -1 * dec_hodge_star(1, sd, GeometricHodge(), Val{:CUDA})

diff --git a/src/DiscreteExteriorCalculus.jl b/src/DiscreteExteriorCalculus.jl
@@ -1988,7 +1988,7 @@ end
 # XXX: This "left/right-hand-rule" trick only works when z=0.
 # XXX: So, do not use this function to test e.g. curved surfaces.
 function eval_constant_dual_form(s::EmbeddedDeltaDualComplex2D, α::SVector{3,Float64})
-  EForm(
+  DualForm{1}(
     hodge_star(1,s) *
       eval_constant_primal_form(s, SVector{3,Float64}(α[2], -α[1], α[3])))
 end