Merge pull request #98 from sintefmath/dev

Coarsening features, lower default report level

Project.toml

@@ -1,7 +1,7 @@
 name = "Jutul"
 uuid = "2b460a1a-8a2b-45b2-b125-b5c536396eb9"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
-version = "0.2.39"
+version = "0.2.40"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -69,7 +69,7 @@ ExprTools = "0.1.8"
 ForwardDiff = "0.10.30"
 GraphMakie = "0.5.3"
 Graphs = "1.8.0"
-HYPRE = "1.6.0"
+HYPRE = "1.6.0, 1.7"
 ILUZero = "0.2.0"
 JLD2 = "0.5"
 KaHyPar = "0.3.0"

ext/JutulMakieExt/interactive_3d.jl

@@ -133,8 +133,7 @@ function plot_interactive_impl(grid, states;
 
             for s in states
                 di = s[k]
-                mv = min(minimum(x -> isnan(x) ? Inf : x, di), mv)
-                Mv = max(maximum(x -> isnan(x) ? -Inf : x, di), Mv)
+                mv, Mv = my_minmax(di, mv, Mv)
             end
             if mv == Mv
                 Mv = 1.01*mv + 1e-12
@@ -151,8 +150,7 @@ function plot_interactive_impl(grid, states;
                         Mv = -Inf
                         for s in states
                             di = view(s[k], row, :)
-                            mv = min(minimum(x -> isnan(x) ? Inf : x, di), mv)
-                            Mv = max(maximum(x -> isnan(x) ? -Inf : x, di), Mv)
+                            mv, Mv = my_minmax(di, mv, Mv)
                         end
                         row_limits["$k"][row] = (mv, Mv)
                     end
@@ -551,13 +549,17 @@ function plot_interactive_impl(grid, states;
 
     # Actual plotting call
     if plot_type == :mesh
-        tri = primitives.triangulation
-        scat = Makie.mesh!(ax, pts, tri; color = ys,
-                                        colorrange = lims,
-                                        backlight = 1,
-                                        colormap = cmap,
-                                        transparency = transparency,
-                                        kwarg...)
+        # TODO: Not sure if this speeds things up
+        tri_c = Makie.to_triangles(primitives.triangulation)
+        pts_c = Makie.to_vertices(pts)
+        scat = Makie.mesh!(ax, pts_c, tri_c;
+            color = ys,
+            colorrange = lims,
+            backlight = 1,
+            colormap = cmap,
+            transparency = transparency,
+            kwarg...
+        )
         if !isnothing(edge_color)
             eplt = Jutul.plot_mesh_edges!(ax, grid; color = edge_color, edge_arg...)
             connect!(eplt.visible, edge_toggle.active)
@@ -912,3 +914,14 @@ function commercial_colormap()
     end
     return cmap
 end
+
+function my_minmax(di, mv, Mv)
+    for v in di
+        if !isfinite(v)
+            continue
+        end
+        mv = min(mv, v)
+        Mv = max(Mv, v)
+    end
+    return (mv, Mv)
+end

src/Jutul.jl

@@ -95,6 +95,9 @@ module Jutul
     include("meshes/meshes.jl")
     include("discretization/discretization.jl")
 
+    # Coarsening utilities
+    include("coarsening.jl")
+
     # Extensions' interfaces
     include("ext/extensions.jl")
 

src/coarsening.jl

@@ -0,0 +1,102 @@
+abstract type AbstractCoarseningFunction end
+
+struct CoarsenByVolumeAverage <: AbstractCoarseningFunction end
+
+function inner_apply_coarsening_function(finevals, fine_indices, op::CoarsenByVolumeAverage, coarse, fine, row, name, entity)
+    subvols = fine[:volumes][fine_indices]
+    return sum(finevals.*subvols)/sum(subvols)
+end
+
+struct CoarsenByHarmonicAverage <: AbstractCoarseningFunction end
+
+function inner_apply_coarsening_function(finevals, fine_indices, op::CoarsenByHarmonicAverage, coarse, fine, row, name, entity)
+    invvals = 1.0./finevals
+    return length(invvals)/sum(invvals)
+end
+
+struct CoarsenByArithemticAverage <: AbstractCoarseningFunction end
+
+function inner_apply_coarsening_function(finevals, fine_indices, op::CoarsenByArithemticAverage, coarse, fine, row, name, entity)
+    return sum(finevals)/length(finevals)
+end
+
+struct CoarsenByFirstValue <: AbstractCoarseningFunction end
+
+function inner_apply_coarsening_function(finevals, fine_indices, op::CoarsenByFirstValue, coarse, fine, row, name, entity)
+    return finevals[1]
+end
+
+function apply_coarsening_function!(coarsevals, finevals, op, coarse::DataDomain, fine::DataDomain, name, entity::Union{Cells, Faces})
+    CG = physical_representation(coarse)
+    function block_indices(CG, block, ::Cells)
+        return findall(isequal(block), CG.partition)
+    end
+    function block_indices(CG, block, ::Faces)
+        return CG.coarse_faces_to_fine[block]
+    end
+    function block_indices(CG, block, ::BoundaryFaces)
+        return CG.coarse_boundary_to_fine[block]
+    end
+    ncoarse = count_entities(coarse, entity)
+    if finevals isa AbstractVector
+        for block in 1:ncoarse
+            ix = block_indices(CG, block, entity)
+            coarsevals[block] = inner_apply_coarsening_function(view(finevals, ix), ix, op, coarse, fine, 1, name, entity)
+        end
+    else
+        for block in 1:ncoarse
+            ix = block_indices(CG, block, entity)
+            for j in axes(coarsevals, 1)
+                coarsevals[j, block] = inner_apply_coarsening_function(view(finevals, j, ix), ix, op, coarse, fine, j, name, entity)
+            end
+        end
+    end
+    return coarsevals
+end
+
+function coarsen_data_domain(D::DataDomain, partition;
+        functions = Dict(),
+        default = CoarsenByArithemticAverage(),
+        default_other = CoarsenByFirstValue(),
+        kwarg...
+    )
+    for (k, v) in pairs(kwarg)
+        functions[k] = v
+    end
+    g = physical_representation(D)
+    cg = CoarseMesh(g, partition)
+    cD = DataDomain(cg)
+    for name in keys(D)
+        if !haskey(cD, name)
+            val = D[name]
+            e = Jutul.associated_entity(D, name)
+            Te = eltype(val)
+            if !(e in (Cells(), Faces(), BoundaryFaces(), nothing))
+                # Other entities are not supported yet.
+                continue
+            end
+            if isnothing(e)
+                # No idea about coarse dims, just copy
+                coarseval = deepcopy(val)
+            elseif val isa AbstractVecOrMat
+                ne = count_entities(cg, e)
+                if val isa AbstractVector
+                    coarseval = zeros(Te, ne)
+                else
+                    coarseval = zeros(Te, size(val, 1), ne)
+                end
+                if eltype(Te)<:AbstractFloat
+                    f = get(functions, name, default)
+                else
+                    f = get(functions, name, default_other)
+                end
+                coarseval = apply_coarsening_function!(coarseval, val, f, cD, D, name, e)
+            else
+                # Don't know what's going on
+                coarseval = deepcopy(val)
+            end
+            cD[name, e] = coarseval
+        end
+    end
+    return cD
+end

src/dd/submodels.jl

@@ -148,7 +148,7 @@ Get subvariable of Jutul variable
 """
 function subvariable(var, map)
     if hasproperty(var, :regions)
-        @warn "Default subvariable called for $(typeof(var)) that contains regions property. Potential missing interface specialization."
+        @warn "Default subvariable called for $(typeof(var)) that contains regions property. Potential missing interface specialization. Map was: $(typeof(map))"
     end
     return var
 end

src/meshes/cart.jl

@@ -70,6 +70,13 @@ struct CartesianMesh{D, Δ, O, T} <: FiniteVolumeMesh
 end
 Base.show(io::IO, g::CartesianMesh) = print(io, "CartesianMesh ($(dim(g))D) with $(join(grid_dims_ijk(g), "x"))=$(number_of_cells(g)) cells")
 
+function CartesianMesh(v::Int, d = 1.0; kwarg...)
+    if d isa Float64
+        d = (d, )
+    end
+    return CartesianMesh((v, ), d; kwarg...)
+end
+
 dim(t::CartesianMesh) = length(t.dims)
 number_of_cells(t::CartesianMesh) = prod(t.dims)
 function number_of_faces(t::CartesianMesh)

src/meshes/coarse.jl

@@ -99,6 +99,10 @@ function dim(G::CoarseMesh)
     return dim(G.parent)
 end
 
+function mesh_z_is_depth(G::CoarseMesh)
+    return mesh_z_is_depth(G.parent)
+end
+
 function number_of_cells(G::CoarseMesh)
     return length(G.partition_lookup)
 end
@@ -188,6 +192,49 @@ function triangulate_mesh(m::CoarseMesh; kwarg...)
         Cells = (cell_data) -> cell_data[cell_ix],
         Faces = (face_data) -> face_data[face_index],
         indices = (Cells = cell_ix, Faces = face_index)
-      )
+    )
     return (points = points, triangulation = triangulation, mapper = mapper)
 end
+
+function plot_primitives(mesh::CoarseMesh, plot_type; kwarg...)
+    if plot_type == :mesh
+        out = triangulate_mesh(mesh; kwarg...)
+    else
+        out = nothing
+    end
+    return out
+end
+
+function cell_dims(cg::CoarseMesh, pos)
+    # TODO: This function is inefficient
+    index = cell_index(cg, pos)
+    g = cg.parent
+    T = eltype(g.node_points)
+    minv = Inf .+ zero(T)
+    maxv = -Inf .+ zero(T)
+
+    for (face, lr) in enumerate(cg.face_neighbors)
+        l, r = lr
+        if l == index || r == index
+            pt, = compute_centroid_and_measure(cg, Faces(), face)
+            if any(isnan.(pt))
+                continue
+            end
+            minv = min.(pt, minv)
+            maxv = max.(pt, maxv)
+        end
+    end
+    for (bface, bcell) in enumerate(cg.boundary_cells)
+        if bcell == index
+            pt, = compute_centroid_and_measure(cg, BoundaryFaces(), bface)
+            if any(isnan.(pt))
+                continue
+            end
+            minv = min.(pt, minv)
+            maxv = max.(pt, maxv)
+        end
+    end
+    Δ = maxv - minv
+    @assert all(x -> x > 0, Δ) "Cell dimensions were zero? Computed $Δ = $maxv - $minv for cell $index."
+    return Tuple(Δ)
+end