Flux reconstruction and visualization for the Laplace operator
https://github.com/j-fu/VoronoiFVM.jl/blob/master/pluto-examples/outflow
We demonstrate the reconstruction of the gradient vector field from the solution of the Laplace operator and its visualization.
begin
import Pkg as _Pkg
haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"])
using Revise
using Test
using SimplexGridFactory, Triangulate, ExtendableGrids, VoronoiFVM
using PlutoUI, GridVisualize
using CairoMakie
CairoMakie.activate!(; type = "svg", visible = false)
GridVisualize.default_plotter!(CairoMakie)
end;
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 = 1:(n - 1)
facet!(builder, points[i], points[i + 1])
end
facet!(builder, points[end], points[1])
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 = 1:length(radii)
facetregion!(builder, i + 1)
circlehole!(builder, holes[:, i], radii[i])
end
simplexgrid(builder)
end
swiss_cheese_2d (generic function with 1 method)
grid = swiss_cheese_2d()
ExtendableGrids.ExtendableGrid{Float64, Int32}; dim: 2 nodes: 1434 cells: 2443 bfaces: 459
System + solution
mutable struct Params
val11::Float64
end
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)
end
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}, Matrix{Float64}}(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)
Test Passed
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)
Test Passed
vis = GridVisualizer(; dim = 2, resolution = (400, 400))
GridVisualizer(Plotter=CairoMakie)
$v_{11}:$
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) = y[1] = src(node[1])
source1d (generic function with 1 method)
flux1d(y, u, edge) = y[1] = u[1, 1]^2 - u[1, 2]^2
flux1d (generic function with 1 method)
function bc1d(y, u, bnode)
boundary_dirichlet!(y, u, bnode; region = 1, value = 0.01)
boundary_dirichlet!(y, u, bnode; region = 2, value = 0.01)
end
bc1d (generic function with 1 method)
grid1d = simplexgrid(-1:0.01:1)
ExtendableGrids.ExtendableGrid{Float64, Int32}; dim: 1 nodes: 201 cells: 200 bfaces: 2
sys1d = VoronoiFVM.System(grid1d;
flux = flux1d,
bcondition = bc1d,
source = source1d,
species = [1],)
VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}, Matrix{Float64}}(num_species=1)
sol1d = solve(sys1d; inival = 0.1)
1×201 Matrix{Float64}: 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
Test Passed
@test nf1d[1, 1, 1] ≈ -nf1d[1, 1, end]
Test Passed
Built with Julia 1.10.2 and
CairoMakie 0.11.5ExtendableGrids 1.2.3
GridVisualize 1.5.0
Pkg 1.10.0
PlutoUI 0.7.55
Revise 3.5.13
SimplexGridFactory 1.0.0
Triangulate 2.3.2
VoronoiFVM 1.16.0