OrdinaryDiffEq.jl nonlinear diffusion

# 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 Test    using Revise    using Printf    using ExtendableGrids    using VoronoiFVM    using OrdinaryDiffEqBDF    using OrdinaryDiffEqRosenbrock    using OrdinaryDiffEqSDIRK    using LinearAlgebra    using PlutoUI    using DataStructures    using GridVisualize    if isdefined(Main, :PlutoRunner)        using CairoMakie        default_plotter!(CairoMakie)    endend

CairoMakie

Solve the nonlinear diffusion equation

$$\partial_t u -\Delta u^m = 0$$

in $\Omega=(-1,1)$ with homogeneous Neumann boundary conditions using the implicit Euler method.

This equation is also called "porous medium equation". The Barenblatt solution

$$b(x,t)=\max\left(0,t^{-\alpha}\left(1-\frac{\alpha(m-1)r^2}{2dmt^{\frac{2\alpha}{d}}}\right)^{\frac{1}{m-1}}\right)$$

is an exact solution of this problem which for m>1 has a finite support. We initialize this problem with the exact solution for $t=t_0=0.001$.

(see Barenblatt, G. I. "On nonsteady motions of gas and fluid in porous medium." Appl. Math. and Mech.(PMM) 16.1 (1952): 67-78.)

Here, we compare the implicit Euler approach in VoronoiFVM with the ODE solvers in DifferentialEquations.jl and demonstrate the possibility to use VoronoiFVM to define differential operators compatible with its ODEFunction interface.

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))end

barenblatt (generic function with 1 method)

function create_porous_medium_problem(n, m)    h = 1.0 / convert(Float64, n / 2)    X = collect(-1:h:1)    grid = simplexgrid(X)    function flux!(f, u, edge, data)        f[1] = u[1, 1]^m - u[1, 2]^m        return nothing    end    storage!(f, u, node, data) = f[1] = u[1]    sys = VoronoiFVM.System(grid, flux = flux!, storage = storage!, species = 1)    return sys, Xend

create_porous_medium_problem (generic function with 1 method)

begin    function run_vfvm(; n = 20, m = 2, t0 = 0.001, tend = 0.01, tstep = 1.0e-6)        sys, X = create_porous_medium_problem(n, m)        inival = unknowns(sys)        inival[1, :] .= map(x -> barenblatt(x, t0, m), X)        sol = VoronoiFVM.solve(sys; inival, times = (t0, tend), Δt = tstep, Δu_opt = 0.01, Δt_min = tstep, store_all = true, log = true, reltol = 1.0e-3)        err = norm(sol[1, :, end] - map(x -> barenblatt(x, tend, m), X))        return sol, sys, err    end    run_vfvm(m = 2, n = 10) # "Precompile"end;

begin    function run_diffeq(; n = 20, m = 2, t0 = 0.001, tend = 0.01, solver = nothing)        sys, X = create_porous_medium_problem(n, m)        inival = unknowns(sys)        inival[1, :] .= map(x -> barenblatt(x, t0, m), X)        state = VoronoiFVM.SystemState(sys)        problem = ODEProblem(state, inival, (t0, tend))        odesol = solve(problem, solver)        sol = reshape(odesol, sys; state)        err = norm(sol[1, :, end] - map(x -> barenblatt(x, tend, m), X))        return sol, sys, err    end    for method in diffeqmethods        run_diffeq(m = 2, n = 10, solver = method.second()) # "Precompile"    endend;

OrderedDict{String, UnionAll} with 4 entries:
  "Rosenbrock23 (Rosenbrock)"    => Rosenbrock23
  "QNDF2 (Like matlab's ode15s)" => QNDF2
  "FBDF"                         => FBDF
  "Implicit Euler"               => ImplicitEuler

t1 = @elapsed sol1, sys1, err1 = run_vfvm(m = m, n = n);history_summary(sol1)

(seconds = 0.349, tasm = 0.123, tlinsolve = 0.224, steps = 832, iters = 1662, maxabsnorm = 2.65e-6, maxrelnorm = 0.000263, maxroundoff = 2.46e-16, iters_per_step = 2.0, facts_per_step = 1.0, liniters_per_step = 0.0)

method:

m = 2; n = 50;

t2 = @elapsed sol2, sys2, err2 = run_diffeq(m = m, n = n, solver = diffeqmethods[method]());history_summary(sol2)

(nd = 165, njac = 82, nf = 248)

Left: VoronoiFVM implicit Euler: 351 ms e=5.17e-02

Right: Rosenbrock23 (Rosenbrock): 14 ms, e=4.62e-02

@test err2 < err1

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