510: Mixture
Test mixture diffusion flux. The problem is here that in the flux function we need to solve a linear system of equations which calculates the fluxes from the gradients.# To do so without (heap) allocations can be achieved using StrideArrays, together with the possibility to have static (compile time) information about the size of the local arrays to be allocated.
$u_i$ are the species partial pressures, $\vec N_i$ are the species fluxes. $D_i^K$ are the Knudsen diffusion coefficients, and $D^B_{ij}$ are the binary diffusion coefficients.
\[ -\nabla \cdot \vec N_i =0 \quad (i=1\dots n)\\ \frac{\vec N_i}{D^K_i} + \sum_{j\neq i}\frac{u_j \vec N_i - u_i \vec N_j}{D^B_{ij}} = -\vec \nabla u_i \quad (i=1\dots n)\]
From this representation, we can derive the matrix $M=(m_{ij})$ with
\[m_{ii}= \frac{1}{D^K_i} + \sum_{j\neq i} \frac{u_j}{D_ij}\\ m_{ij}= -\sum_{j\neq i} \frac{u_i}{D_ij}\]
such that
\[ M\begin{pmatrix} \vec N_1\\ \vdots\\ \vec N_n \end{pmatrix} = \begin{pmatrix} \vec \nabla u_1\\ \vdots\\ \vec \nabla u_n \end{pmatrix}\]
In the two point flux finite volume discretization, this results into a corresponding linear system which calculates the discrete edge fluxes from the discrete gradients. Here we demonstrate how to implement this in a fast, (heap) allocation free way.
For this purpose, intermediate arrays need to be allocated on the stack with via StrideArrays.jl or MArray from StaticArrays.jl They need to have the same element type as the unknowns passed to the flux function (which could be Float64 or some dual number). Size information must be static, e.g. a global constant, or, as demonstrated here, a type parameter. Broadcasting probably should be avoided.
As documented in StrideArrays.jl, use @gc_preserve when passing a StrideArray as a function parameter.
Alternatively, we can try to avoid StrideArrays.jl and use MArrays together with inlined code. In this case, one should be aware of this issue, which requires @inbounds e.g. with reverse order loops.
See also this Discourse thread.
module Example510_Mixture
using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize
using LinearAlgebra
using AMGCLWrap
using Random
using StrideArraysCore: @gc_preserve, StrideArray, StaticInt, PtrArray
using LinearSolve, ExtendableSparse
using StaticArrays
using ExtendableSparse
struct MyData{NSPec}
DBinary::Symmetric{Float64, Matrix{Float64}}
DKnudsen::Vector{Float64}
diribc::Vector{Int}
end
nspec(::MyData{NSpec}) where {NSpec} = NSpec
function flux_strided(f, u, edge, data)
T = eltype(u)
M = StrideArray{T}(undef, StaticInt(nspec(data)), StaticInt(nspec(data)))
au = StrideArray{T}(undef, StaticInt(nspec(data)))
du = StrideArray{T}(undef, StaticInt(nspec(data)))
ipiv = StrideArray{Int}(undef, StaticInt(nspec(data)))
for ispec = 1:nspec(data)
M[ispec, ispec] = 1.0 / data.DKnudsen[ispec]
du[ispec] = u[ispec, 1] - u[ispec, 2]
au[ispec] = 0.5 * (u[ispec, 1] + u[ispec, 2])
end
for ispec = 1:nspec(data)
for jspec = 1:nspec(data)
if ispec != jspec
M[ispec, ispec] += au[jspec] / data.DBinary[ispec, jspec]
M[ispec, jspec] = -au[ispec] / data.DBinary[ispec, jspec]
end
end
end
if VERSION >= v"1.9-rc0"
# Pivoting linear system solution via RecursiveFactorizations.jl
@gc_preserve inplace_linsolve!(M, du, ipiv)
else
# Non-pivoting implementation currently implemented in vfvm_functions.jl
@gc_preserve inplace_linsolve!(M, du)
end
for ispec = 1:nspec(data)
f[ispec] = du[ispec]
end
end
function flux_marray(f, u, edge, data)
T = eltype(u)
n = nspec(data)
M = MMatrix{nspec(data), nspec(data), T}(undef)
au = MVector{nspec(data), T}(undef)
du = MVector{nspec(data), T}(undef)
ipiv = MVector{nspec(data), Int}(undef)
for ispec = 1:nspec(data)
M[ispec, ispec] = 1.0 / data.DKnudsen[ispec]
du[ispec] = u[ispec, 1] - u[ispec, 2]
au[ispec] = 0.5 * (u[ispec, 1] + u[ispec, 2])
end
for ispec = 1:nspec(data)
for jspec = 1:nspec(data)
if ispec != jspec
M[ispec, ispec] += au[jspec] / data.DBinary[ispec, jspec]
M[ispec, jspec] = -au[ispec] / data.DBinary[ispec, jspec]
end
end
end
# Here, we also could use @gc_preserve.
# As this function is inlined one can avoid StrideArrays.jl
# Starting with Julia 1.8 one also can use callsite @inline.
inplace_linsolve!(M, du)
for ispec = 1:nspec(data)
f[ispec] = du[ispec]
end
end
function bcondition(f, u, node, data)
for species = 1:nspec(data)
boundary_dirichlet!(f, u, node; species, region = data.diribc[1],
value = species % 2)
boundary_dirichlet!(f, u, node; species, region = data.diribc[2],
value = 1 - species % 2)
end
end
function main(; n = 11, nspec = 5,
dim = 2,
Plotter = nothing,
verbose = false,
unknown_storage = :dense,
flux = :flux_strided,
strategy = nothing,
assembly = :cellwise)
h = 1.0 / convert(Float64, n - 1)
X = collect(0.0:h:1.0)
DBinary = Symmetric(fill(0.1, nspec, nspec))
for ispec = 1:nspec
DBinary[ispec, ispec] = 0
end
DKnudsen = fill(1.0, nspec)
if dim == 1
grid = simplexgrid(X)
diribc = [1, 2]
elseif dim == 2
grid = simplexgrid(X, X)
diribc = [4, 2]
else
grid = simplexgrid(X, X, X)
diribc = [4, 2]
end
function storage(f, u, node, data)
f .= u
end
_flux = flux == :flux_strided ? flux_strided : flux_marray
data = MyData{nspec}(DBinary, DKnudsen, diribc)
sys = VoronoiFVM.System(grid; flux = _flux, storage, bcondition, species = 1:nspec, data, assembly,unknown_storage)
@info "Strategy: $(strategy)"
control = SolverControl(strategy, sys)
control.maxiters = 500
@info control.method_linear
u = solve(sys; verbose, control, log = true)
@show norm(u)
norm(u)
end
using Test
function runtests()
strategies = [DirectSolver(UMFPACKFactorization()),
GMRESIteration(UMFPACKFactorization()),
GMRESIteration(UMFPACKFactorization(), EquationBlock()),
GMRESIteration(AMGCL_AMGPreconditioner(),EquationBlock()),
BICGstabIteration(AMGCL_AMGPreconditioner(),EquationBlock()),
GMRESIteration(ILUZeroPreconditioner()), GMRESIteration(ILUZeroPreconditioner(), EquationBlock()),
GMRESIteration(ILUZeroPreconditioner(), PointBlock()) ]
val1D = 4.788926530387466
val2D = 15.883072449873742
val3D = 52.67819183426213
@test main(; dim = 1, assembly = :cellwise) ≈ val1D
@test main(; dim = 2, assembly = :cellwise) ≈ val2D
@test main(; dim = 3, assembly = :cellwise) ≈ val3D
@test main(; dim = 1, flux = :flux_marray, assembly = :cellwise) ≈ val1D
@test main(; dim = 2, flux = :flux_marray, assembly = :cellwise) ≈ val2D
@test main(; dim = 3, flux = :flux_marray, assembly = :cellwise) ≈ val3D
@test all(map(strategy -> main(; dim = 2, flux = :flux_marray, strategy) ≈ val2D, strategies))
end
endThis page was generated using Literate.jl.