# 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) endendbegin using Revise using Test using VoronoiFVM using GridVisualize using SimplexGridFactory, Triangulate using ExtendableGrids using PlutoUI, HypertextLiteral, UUIDs if isdefined(Main, :PlutoRunner) using CairoMakie CairoMakie.activate!(; type = "png", visible = false) GridVisualize.default_plotter!(CairoMakie) endendCairoMakie
Outflow boundary conditions
We show how to implement outflow boundary conditions when the velocities at the boundary are calculated by another equation in the system. A typical case is solute transport in porous media where fluid flow is calculated by Darcy's law which defines the convective velocity in the solute transport equation.
Regard the following system of equations in domain \(\Omega\subset \mathbb R^d\):
$$\begin{aligned} \nabla \cdot \vec v &=0\\ \vec v&=-k\nabla p\\ \nabla \cdot \vec j &=0\\ \vec j&= - D\nabla c - c\vec v \end{aligned}$$
The variable p can be seen as a the pressure of a fluid in porous medium. c is the concentration of a transported species.
We subdivide the boundary: \(\partial\Omega=Γ_{in}\cup Γ_{out}\cup Γ_{noflow}\) abs set
$$\begin{aligned} p=&1 \quad & c&=c_{in} & \text{on}\quad Γ_{in}\\ p=&0 \quad & \vec j\cdot \vec n &= c\vec v \cdot \vec n & \text{on}\quad Γ_{out}\\ \vec v\cdot \vec n &=0 \quad & \vec j\cdot \vec n &= 0 & \text{on}\quad Γ_{noflow}\\ \end{aligned}$$
Discretization data
Base.@kwdef struct FlowTransportData k = 1.0 v_in = 1.0 c_in = 0.5 D = 1.0 Γ_in = 1 Γ_out = 2 ip = 1 ic = 2endFlowTransportData
X = 0:0.1:10.0:0.1:1.0
darcyvelo(u, data) = data.k * (u[data.ip, 1] - u[data.ip, 2])darcyvelo (generic function with 1 method)
function flux(y, u, edge, data) vh = darcyvelo(u, data) y[data.ip] = vh bp, bm = fbernoulli_pm(vh / data.D) y[data.ic] = data.D * (bm * u[data.ic, 1] - bp * u[data.ic, 2]) return nothingendflux (generic function with 1 method)
function bcondition(y, u, bnode, data) boundary_neumann!(y, u, bnode; species = data.ip, region = data.Γ_in, value = data.v_in) boundary_dirichlet!(y, u, bnode; species = data.ip, region = data.Γ_out, value = 0) boundary_dirichlet!(y, u, bnode; species = data.ic, region = data.Γ_in, value = data.c_in) return nothingendbcondition (generic function with 1 method)
This function describes the outflow boundary condition. It is called on edges (including interior ones) which have at least one ode on one of the outflow boundaries. Within this function outflownode can be used to identify that node. There is some ambiguity in the case that both nodes are outflow nodes, in that case it is assumed that the contribution is zero. In the present case this is guaranteed by the constant Dirichlet boundary condition for the pressure.
function boutflow(y, u, edge, data) y[data.ic] = -darcyvelo(u, data) * u[data.ic, outflownode(edge)] return nothingendboutflow (generic function with 1 method)
function flowtransportsystem(grid; kwargs...) data = FlowTransportData(; kwargs...) return VoronoiFVM.System( grid; flux, bcondition, boutflow, data, outflowboundaries = [data.Γ_out], species = [1, 2], )endflowtransportsystem (generic function with 1 method)
function checkinout(sys, sol) data = sys.physics.data tfact = TestFunctionFactory(sys) tf_in = testfunction(tfact, [data.Γ_out], [data.Γ_in]) tf_out = testfunction(tfact, [data.Γ_in], [data.Γ_out]) return (; in = integrate(sys, tf_in, sol), out = integrate(sys, tf_out, sol))endcheckinout (generic function with 1 method)
1D Case
grid = simplexgrid(X)ExtendableGrids.ExtendableGrid{Float64, Int32}
dim = 1
nnodes = 11
ncells = 10
nbfaces = 2
sys1 = flowtransportsystem(grid);sol1 = solve(sys1; verbose = "n");
t1 = checkinout(sys1, sol1)(in = [1.0, 0.5000000000000003], out = [-1.0, -0.5000000000000003])
@test t1.in ≈ -t1.out�[32m�[1mTest Passed�[22m�[39m
@test maximum(sol1[2, :]) ≈ 0.5�[32m�[1mTest Passed�[22m�[39m
@test minimum(sol1[2, :]) ≈ 0.5�[32m�[1mTest Passed�[22m�[39m
2D Case
begin g2 = simplexgrid(X, X) bfacemask!(g2, [1, 0.3], [1, 0.7], 5)endExtendableGrids.ExtendableGrid{Float64, Int32}
dim = 2
nnodes = 121
ncells = 200
nbfaces = 40
gridplot(g2; size = (300, 300))
sys2 = flowtransportsystem(g2; Γ_in = 4, Γ_out = 5);sol2 = solve(sys2; verbose = "n")2×121 VoronoiFVM.DenseSolutionArray{Float64, 2}:
1.12521 1.02537 0.925877 0.826948 … 0.453027 0.377299 0.321519 0.298904
0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
let vis = GridVisualizer(; size = (700, 300), layout = (1, 2)) scalarplot!(vis[1, 1], g2, sol2[1, :]) scalarplot!(vis[1, 2], g2, sol2[2, :]; limits = (0, 1)) reveal(vis)end
t2 = checkinout(sys2, sol2)(in = [1.0000000000000009, 0.5000000000000002], out = [-1.0000000000000004, -0.5000000000000002])
@test t2.in ≈ -t2.out�[32m�[1mTest Passed�[22m�[39m
@test maximum(sol2[2, :]) ≈ 0.5�[32m�[1mTest Passed�[22m�[39m
@test minimum(sol2[2, :]) ≈ 0.5�[32m�[1mTest Passed�[22m�[39m
3D Case
begin g3 = simplexgrid(X, X, X) bfacemask!(g3, [0.3, 0.3, 0], [0.7, 0.7, 0], 7)endExtendableGrids.ExtendableGrid{Float64, Int32}
dim = 3
nnodes = 1331
ncells = 6000
nbfaces = 1200
gridplot(g3; size = (300, 300))
sys3 = flowtransportsystem(g3; Γ_in = 6, Γ_out = 7);sol3 = solve(sys3; verbose = "n")2×1331 VoronoiFVM.DenseSolutionArray{Float64, 2}:
0.547438 0.539229 0.516946 0.488884 … 1.32867 1.32914 1.32951 1.32966
0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
let vis = GridVisualizer(; size = (700, 300), layout = (1, 2)) scalarplot!(vis[1, 1], g3, sol3[1, :]) scalarplot!(vis[1, 2], g3, sol3[2, :]; limits = (0, 1)) reveal(vis)end
t3 = checkinout(sys3, sol3)(in = [1.0000000000000369, 0.5000000000000183], out = [-1.0000000000000375, -0.5000000000000188])
@test t3.in ≈ -t3.out�[32m�[1mTest Passed�[22m�[39m
@test maximum(sol3[2, :]) ≈ 0.5�[32m�[1mTest Passed�[22m�[39m
@test minimum(sol3[2, :]) ≈ 0.5�[32m�[1mTest Passed�[22m�[39m