Flux reconstruction and visualization for the Laplace operator

We demonstrate the reconstruction of the gradient vector field from the solution of the Laplace operator and its visualization.

# Convert this cell to markdown in order to enable Pluto's inbuilt package managerif isdefined(Main, :PlutoRunner)    using Pkg    docsdir = joinpath(@__DIR__, "..", "docs")    if isdir(docsdir)        Pkg.activate(docsdir)    endend

begin    using Revise    using Test    using SimplexGridFactory, Triangulate, ExtendableGrids, VoronoiFVM    using PlutoUI, GridVisualize    if isdefined(Main, :PlutoRunner)        using CairoMakie        CairoMakie.activate!(; type = "png", visible = false)        GridVisualize.default_plotter!(CairoMakie)    endend;

Grid

Define a "Swiss cheese domain" with punched-out holes, where each hole boundary corresponds to a different boundary condition.

function swiss_cheese_2d()    function circlehole!(builder, center, radius; n = 20)        points = [            point!(builder, center[1] + radius * sin(t), center[2] + radius * cos(t))                for t in range(0, 2π; length = n)        ]        for i in 1:(n - 1)            facet!(builder, points[i], points[i + 1])        end        facet!(builder, points[end], points[1])        return holepoint!(builder, center)    end    builder = SimplexGridBuilder(; Generator = Triangulate)    cellregion!(builder, 1)    maxvolume!(builder, 0.1)    regionpoint!(builder, 0.1, 0.1)    p1 = point!(builder, 0, 0)    p2 = point!(builder, 10, 0)    p3 = point!(builder, 10, 10)    p4 = point!(builder, 0, 10)    facetregion!(builder, 1)    facet!(builder, p1, p2)    facet!(builder, p2, p3)    facet!(builder, p3, p4)    facet!(builder, p4, p1)    holes = [        1.0 2.0        8.0 9.0        2.0 8.0        8.0 4.0        9.0 1.0        3.0 4.0        4.0 6.0        7.0 9.0        4.0 7.0        7.0 5.0        2.0 1.0        4.0 1.0        4.0 8.0        3.0 6.0        4.0 9.0        6.0 9.0        3.0 5.0        1.0 4.0    ]'    radii = [        0.15,        0.15,        0.1,        0.35,        0.2,        0.3,        0.1,        0.4,        0.1,        0.4,        0.2,        0.2,        0.2,        0.35,        0.15,        0.25,        0.15,        0.25,    ]    for i in 1:length(radii)        facetregion!(builder, i + 1)        circlehole!(builder, holes[:, i], radii[i])    end    return simplexgrid(builder)end

swiss_cheese_2d (generic function with 1 method)

grid = swiss_cheese_2d()

ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       2
   nnodes =    1434
   ncells =    2443
  nbfaces =     459

System + solution

mutable struct Params    val11::Float64end

params = Params(5)

Params(5.0)

Simple flux function for Laplace operator

flux(y, u, edge, data) = y[1] = u[1, 1] - u[1, 2];

At hole #11, the value will be bound to a slider defined below

function bc(y, u, bnode, data)    boundary_dirichlet!(y, u, bnode; region = 2, value = 10.0)    boundary_dirichlet!(y, u, bnode; region = 3, value = 0.0)    boundary_dirichlet!(y, u, bnode; region = 11, value = data.val11)    return nothingend

bc (generic function with 1 method)

Define a finite volume system with Dirichlet boundary conditions at some of the holes

system = VoronoiFVM.System(grid; flux = flux, species = 1, bcondition = bc, data = params)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=2, nnodes=1434, ncells=2443,
  nbfaces=459),  
  physics = Physics(data=Main.var"workspace#4".Params, flux=flux,
  storage=default_storage, breaction=bc, ),  
  num_species = 1)

Solve, and trigger solution upon boundary value change

begin    params.val11 = val11    sol = solve(system)end;

@test params.val11 != 5.0 || isapprox(sum(sol), 7842.2173682050525; rtol = 1.0e-12)

�[32m�[1mTest Passed�[22m�[39m

Flux reconstruction

Reconstruct the node flux. It is a \(d\times n_{spec}\times n_{nodes}\) tensor. nf[:,ispec,:] then is a vector function representing the flux density of species ispec in each node of the domain. This readily can be fed into GridVisualize.vectorplot.

The reconstruction is based on the "magic formula" R. Eymard, T. Gallouet, R. Herbin, IMA Journal of Numerical Analysis (2006) 26, 326−353 (Arxive version), Lemma 2.4 .

nf = nodeflux(system, sol)

2×1×1434 Array{Float64, 3}:
[:, :, 1] =
  0.05468367090633706
 -0.05468367090633706

[:, :, 2] =
 0.01852257911924369
 0.014818063295393813

[:, :, 3] =
 0.0
 0.0

;;; … 

[:, :, 1432] =
 0.020523328431052094
 0.036034112502081016

[:, :, 1433] =
 0.028906869669800477
 0.012529162794698063

[:, :, 1434] =
 0.03545781788403497
 0.0342143068249023

@test params.val11 != 5.0 || isapprox(sum(nf), 978.000534849034; rtol = 1.0e-14)

�[32m�[1mTest Passed�[22m�[39m

vis = GridVisualizer(; dim = 2, resolution = (400, 400))

GridVisualizer(Plotter=CairoMakie)

\(v_{11}:\)5.0

Joint plot of solution and flux reconstruction

begin    scalarplot!(vis, grid, sol[1, :]; levels = 9, colormap = :summer, clear = true)    vectorplot!(vis, grid, nf[:, 1, :]; clear = false, vscale = 1.5)    reveal(vis)end

The 1D case

src(x) = exp(-x^2 / 0.01)

src (generic function with 1 method)

source1d(y, node, data) = y[1] = src(node[1])

source1d (generic function with 1 method)

flux1d(y, u, edge, data) = y[1] = u[1, 1]^2 - u[1, 2]^2

flux1d (generic function with 1 method)

function bc1d(y, u, bnode, data)    boundary_dirichlet!(y, u, bnode; region = 1, value = 0.01)    boundary_dirichlet!(y, u, bnode; region = 2, value = 0.01)    return nothingend

bc1d (generic function with 1 method)

grid1d = simplexgrid(-1:0.01:1)

ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       1
   nnodes =     201
   ncells =     200
  nbfaces =       2

sys1d = VoronoiFVM.System(    grid1d;    flux = flux1d,    bcondition = bc1d,    source = source1d,    species = [1],)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=201, ncells=200,
  nbfaces=2),  
  physics = Physics(flux=flux1d, storage=default_storage, source=source1d,
  breaction=bc1d, ),  
  num_species = 1)

sol1d = solve(sys1d; inival = 0.1)

1×201 VoronoiFVM.DenseSolutionArray{Float64, 2}:
 0.01  0.0314043  0.0432719  0.0525231  …  0.0525231  0.0432719  0.0314043  0.01

nf1d = nodeflux(sys1d, sol1d)

1×1×201 Array{Float64, 3}:
[:, :, 1] =
 -0.08862269254527583

[:, :, 2] =
 -0.08862269254527581

[:, :, 3] =
 -0.08862269254527581

;;; … 

[:, :, 199] =
 0.08862269254527581

[:, :, 200] =
 0.08862269254527581

[:, :, 201] =
 0.08862269254527583

let    vis1d = GridVisualizer(; dim = 1, resolution = (500, 250), legend = :lt)    scalarplot!(vis1d, grid1d, map(src, grid1d); label = "rhs", color = :blue)    scalarplot!(vis1d, grid1d, sol1d[1, :]; label = "solution", color = :red, clear = false)    vectorplot!(vis1d, grid1d, nf1d[:, 1, :]; label = "flux", clear = false, color = :green)    reveal(vis1d)end

@test nf1d[1, 1, 101] ≈ 0.0

�[32m�[1mTest Passed�[22m�[39m

@test nf1d[1, 1, 1] ≈ -nf1d[1, 1, end]

�[32m�[1mTest Passed�[22m�[39m