# 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 Test using Revise using Printf using ExtendableGrids using VoronoiFVM using SciMLBase: NoInit using OrdinaryDiffEqBDF using OrdinaryDiffEqRosenbrock using OrdinaryDiffEqSDIRK using LinearAlgebra using PlutoUI, HypertextLiteral, UUIDs using DataStructures using GridVisualize if isdefined(Main, :PlutoRunner) using CairoMakie default_plotter!(CairoMakie) CairoMakie.activate!(type = "png") endend1D Nonlinear Storage
TableOfContents(aside = false)This equation comes from the transformation of the nonlinear diffusion equation
$$\partial_t v - \Delta v^m = 0$$
to
$$\partial_t u^\frac{1}{m} -\Delta u = 0$$
in \(\Omega=(-1,1)\) with homogeneous Neumann boundary conditions. We can derive an exact solution from the Barenblatt solution of the equation for u.
function barenblatt(x, t, m) tx = t^(-1.0 / (m + 1.0)) xx = x * tx xx = xx * xx xx = 1 - xx * (m - 1) / (2.0 * m * (m + 1)) if xx < 0.0 xx = 0.0 end return tx * xx^(1.0 / (m - 1.0))endbarenblatt (generic function with 1 method)
begin const m = 2 const ε = 1.0e-10 const n = 50 const t0 = 1.0e-3 const tend = 1.0e-2end0.01
X = collect(-1:(2.0 / n):1)51-element Vector{Float64}:
-1.0
-0.96
-0.92
-0.88
-0.84
-0.8
-0.76
⋮
0.8
0.84
0.88
0.92
0.96
1.0
u0 = map(x -> barenblatt(x, t0, m)^m, X)51-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.0
0.0
0.0
⋮
0.0
0.0
0.0
0.0
0.0
0.0
begin grid = simplexgrid(X)endExtendableGrids.ExtendableGrid{Float64, Int32}
dim = 1
nnodes = 51
ncells = 50
nbfaces = 2
Direct implementation with VoronoiFVM
function flux!(f, u, edge, data) f[1] = u[1, 1] - u[1, 2] return nothingendflux! (generic function with 1 method)
Storage term needs to be regularized as its derivative at 0 is infinity:
function storage!(f, u, node, data) f[1] = (ε + u[1])^(1.0 / m) return nothingendstorage! (generic function with 1 method)
begin physics = VoronoiFVM.Physics( flux = flux!, storage = storage! ) sys = VoronoiFVM.System(grid, physics, species = 1) inival = unknowns(sys) inival[1, :] = u0 control = VoronoiFVM.SolverControl() tsol = VoronoiFVM.solve(sys; inival, times = (t0, tend), Δt_min = 1.0e-4, Δt = 1.0e-4, Δu_opt = 0.1, force_first_step = true, log = true) summary(tsol.history)end(seconds = 0.311, tasm = 0.226, tlinsolve = 0.0825, steps = 731, iters = 2920, maxabsnorm = 1.73e-12, maxrelnorm = 1.75e-11, maxroundoff = 9.06e-15, iters_per_step = 4.0, facts_per_step = 1.0, liniters_per_step = 0.0)
Implementation as DAE
If we want to solve the problem with DifferentialEquations.jl solvers, we see that the problem structure does not fit into the setting of that package due to the nonlinearity under the time derivative. Here we propose a reformulation to a DAE as a way to achieve this possibility:
$$\begin{cases} \partial_t w -\Delta u &= 0\\ w^m - u &=0 \end{cases}$$
function dae_storage!(y, u, node, data) y[1] = u[2] return nothingenddae_storage! (generic function with 1 method)
function dae_reaction!(y, u, node, data) y[2] = u[2]^m - u[1] return nothingenddae_reaction! (generic function with 1 method)
First, we test this with the implicit Euler method of VoronoiFVM
begin dae_physics = VoronoiFVM.Physics( flux = flux!, storage = dae_storage!, reaction = dae_reaction! ) dae_sys = VoronoiFVM.System(grid, dae_physics, species = [1, 2]) dae_inival = unknowns(dae_sys) dae_inival[1, :] .= u0 dae_inival[2, :] .= u0 .^ (1 / m) dae_control = VoronoiFVM.SolverControl() dae_tsol = VoronoiFVM.solve(dae_sys; inival = dae_inival, times = (t0, tend), Δt_min = 1.0e-4, Δt = 1.0e-4, Δu_opt = 0.1, force_first_step = true, log = true) summary(dae_tsol.history)end(seconds = 0.475, tasm = 0.175, tlinsolve = 0.298, steps = 732, iters = 2205, maxabsnorm = 9.72e-11, maxrelnorm = 9.82e-10, maxroundoff = 6.99e-13, iters_per_step = 3.02, facts_per_step = 1.0, liniters_per_step = 0.0)
Implementation via OrdinaryDiffEq.jl
OrderedDict{String, UnionAll} with 5 entries:
"QNDF2 (Like matlab's ode15s)" => QNDF2
"Rodas5" => Rodas5
"Rosenbrock23 (Rosenbrock)" => Rosenbrock23
"FBDF" => FBDF
"Implicit Euler" => ImplicitEuler
method:
begin de_sys = VoronoiFVM.System(grid, dae_physics, species = [1, 2]) problem = ODEProblem(de_sys, dae_inival, (t0, tend)) de_odesol = solve( problem, diffeqmethods[method](), adaptive = true, reltol = 1.0e-3, abstol = 1.0e-3, initializealg = NoInit() ) de_tsol = reshape(de_odesol, de_sys)end;
t=
exact = map(x -> barenblatt(x, tend, m) .^ m, X)51-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.0
0.0
0.0
⋮
0.0
0.0
0.0
0.0
0.0
0.0
@test norm(tsol[1, :, end] - exact, Inf) < 0.1�[32m�[1mTest Passed�[22m�[39m
@test norm(dae_tsol[1, :, end] - exact, Inf) < 0.1�[32m�[1mTest Passed�[22m�[39m
@test norm(de_tsol[1, :, end] - exact, Inf) < 0.05�[32m�[1mTest Passed�[22m�[39m
myaside (generic function with 1 method)
function plotsolutions() vis = GridVisualizer(resolution = (380, 200), dim = 1, legend = :lt) u = tsol(t) u_dae = dae_tsol(t) u_de = de_tsol(t) scalarplot!(vis, X, map(x -> barenblatt(x, t, m) .^ m, X), clear = true, color = :red, linestyle = :solid, flimits = (0, 100), label = "exact") scalarplot!(vis, grid, u_dae[1, :], clear = false, color = :green, linestyle = :solid, label = "vfvm_dae") scalarplot!(vis, grid, u_de[1, :], clear = false, color = :blue, markershape = :cross, linestyle = :dot, label = "vfvm_diffeq") scalarplot!(vis, grid, u[1, :], clear = false, color = :black, markershape = :none, linestyle = :dash, title = "t=$(t)", label = "vfvm_default") return reveal(vis)endplotsolutions (generic function with 1 method)