Tests + less stringent types in WENO

Project.toml

@@ -1,7 +1,7 @@
 name = "Jutul"
 uuid = "2b460a1a-8a2b-45b2-b125-b5c536396eb9"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"

src/WENO/WENO.jl

@@ -66,7 +66,9 @@ module WENO
     function weno_discretize(domain::DataDomain; kwarg...)
         weno_cells = weno_discretize_cells(domain)
         g = physical_representation(domain)
-        g::UnstructuredMesh
+        if !(g isa UnstructuredMesh)
+            g = UnstructuredMesh(g)
+        end
         N = g.faces.neighbors
         D = dim(g)
         Point_t = SVector{D, Float64}
@@ -95,7 +97,7 @@ module WENO
     function weno_discretize_half_face(cell, wenodisc, fc::SVector{D, R}) where {D, R}
         function half_face_impl(planar_set, grad, V)
             cells = map(i -> wenodisc.stencil[i].cell, planar_set)
-            points = map(i -> wenodisc.stencil[i].point, planar_set)
+            points = map(i -> convert(SVector{D, R}, wenodisc.stencil[i].point), planar_set)
             # Now that we know the delta to the face centroid, we can collapse
             # the gradient basis to a single basis for the set of support cells.
             new_grad = grad[1]*V[1]
@@ -141,13 +143,16 @@ module WENO
         return WENOHalfFaceDiscretization(cell, V, vec(grad_disc))
     end
 
-    function area_from_points(points::NTuple{N, SVector{D, Float64}}) where {N, D}
+    function area_from_points(points)
+        N = length(points)
+        D = length(first(points))
         if D == 2
             @assert N == 3
             u, v, w = points
-            U = u - w
-            V = v - w
-            area = 0.5*norm(cross(U, V))
+            x1, y1 = u
+            x2, y2 = v
+            x3, y3 = w
+            area = 0.5*abs(x1*(y2 - y3) + x2*(y3 - y1) + x3*(y1 - y2))
         else
             @assert D == 3
             @assert N == 4
@@ -162,6 +167,9 @@ module WENO
 
     function weno_discretize_cells(domain::DataDomain)
         g = physical_representation(domain)
+        if !(g isa UnstructuredMesh)
+            g = UnstructuredMesh(g)
+        end
         D = dim(g)
         Point_t = SVector{D, Float64}
         cc = reinterpret(Point_t, domain[:cell_centroids])
@@ -224,32 +232,44 @@ module WENO
         pts = map(x -> x[:point], out)
         S = point_set_transformation_basis(pts)
         function convert_neighbor(i)
-            out[i][:point] = S*(out[i][:point])
+            pt = out[i][:point]
+            T = typeof(pt)
+            out[i][:point] = convert(T, S*pt)
             return (; pairs(out[i])...)
         end
         return (center = center, S = S, stencil = map(convert_neighbor, eachindex(out)))
     end
 
     function point_set_transformation_basis(pts::Vector{SVector{griddim, Float64}}) where griddim
         M = hcat(pts...)'
+        if griddim == 1
+            MT = SMatrix{1, 1, Float64, 1}
+        elseif griddim == 2
+            MT = SMatrix{2, 2, Float64, 4}
+        else
+            @assert griddim == 3
+            MT = SMatrix{3, 3, Float64, 9}
+        end
+        out = missing
         try
             UDV = svd(M)
             U = UDV.U
             D = UDV.S
             Vt = UDV.Vt
             Dbar = Diagonal(D[1:griddim])
-            return Dbar*inv(Vt')
+            out = Dbar*inv(Vt')
         catch
             @warn "Failure to decompose point set, using identity" M
             if griddim == 1
-                return @SMatrix [1.0]
+                out = @SMatrix [1.0]
             elseif griddim == 2
-                return @SMatrix [1.0 0.0; 0.0 1.0]
+                out = @SMatrix [1.0 0.0; 0.0 1.0]
             else
                 @assert griddim == 3
-                return @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
+                out = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
             end
         end
+        return convert(MT, out)
     end
 
     function cell_to_node_indices(g::UnstructuredMesh, c)

test/runtests.jl

@@ -9,3 +9,4 @@ include("partitioning.jl")
 include("units.jl")
 include("sparsity.jl")
 include("nfvm.jl")
+include("weno.jl")

test/weno.jl

@@ -0,0 +1,33 @@
+using Jutul, Test
+@testset "WENO setup" begin
+    mesh_2d = UnstructuredMesh(CartesianMesh((3, 3)))
+    @testset "2D" begin
+        call = Jutul.WENO.weno_discretize_cells(DataDomain(mesh_2d))
+        ix = 5
+        c = call[ix]
+        pset = c.planar_set
+        cells = map(x -> x.cell, c.stencil)
+
+        @test length(pset) == 4
+        for quad in pset
+            @test length(quad) == 3
+            i, j, k = quad
+            @test cells[i] == ix
+            # Neighbor set in 2D should TPFA pairs of 2, 4, 6, 8
+            c1 = cells[j]
+            c2 = cells[k]
+            if c1 > c2
+                c2, c1 = c1, c2
+            end
+            @test (c1, c2) in [(2, 4), (4, 8), (6, 8), (2, 6)]
+        end
+    end
+    mesh_3d = UnstructuredMesh(CartesianMesh((3, 3, 3)))
+
+    w2d = Jutul.WENO.weno_discretize(DataDomain(mesh_2d))
+    @test length(w2d) == number_of_faces(mesh_2d)
+    @test eltype(w2d) == Jutul.WENO.WENOFaceDiscretization{2, 3, Float64}
+    w3d = Jutul.WENO.weno_discretize(DataDomain(mesh_3d))
+    @test length(w3d) == number_of_faces(mesh_3d)
+    @test eltype(w3d) == Jutul.WENO.WENOFaceDiscretization{3, 4, Float64}
+end