106: 1D Nonlinear Diffusion equation

(source code)

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

module Example106_NonlinearDiffusion1Dusing Printfusing VoronoiFVMusing ExtendableGridsusing GridVisualize# Problem data structure to avoid global variablesmutable struct ProblemData    m::Float64  # Nonlinearity parameterendfunction 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))endfunction main(;        n = 20, m = 2, Plotter = nothing, verbose = false,        unknown_storage = :sparse, tend = 0.01, tstep = 0.0001, assembly = :edgewise    )    # Create problem data structure    problem_data = ProblemData(m)    # Create a one-dimensional discretization    h = 1.0 / convert(Float64, n / 2)    X = collect(-1:h:1)    grid = simplexgrid(X)    # Flux function which describes the flux    # between neighboring control volumes    function flux!(f, u, edge, data)        f[1] = u[1, 1]^data.m - u[1, 2]^data.m        return nothing    end    # Storage term    function storage!(f, u, node, data)        f[1] = u[1]        return nothing    end    # Create a physics structure    physics = VoronoiFVM.Physics(;        flux = flux!,        storage = storage!,        data = problem_data    )    # Create a finite volume system - either    # in the dense or  the sparse version.    # The difference is in the way the solution object    # is stored - as dense or as sparse matrix    sys = VoronoiFVM.System(grid, physics; unknown_storage = unknown_storage, assembly = assembly)    # Add species 1 to region 1    enable_species!(sys, 1, [1])    # Create a solution array    inival = unknowns(sys)    t0 = 0.001    # Broadcast the initial value    inival[1, :] .= map(x -> barenblatt(x, t0, problem_data.m), X)    # Create solver control info for constant time step size    control = VoronoiFVM.SolverControl()    control.verbose = verbose    control.Δt_min = tstep    control.Δt_max = tstep    control.Δt = tstep    control.Δu_opt = 1    tsol = solve(sys; inival, times = [t0, tend], control)    p = GridVisualizer(; Plotter = Plotter, layout = (1, 1), fast = true)    for i in 1:length(tsol.t)        time = tsol.t[i]        scalarplot!(            p[1, 1], grid, tsol[1, :, i]; title = @sprintf("t=%.3g", time),            color = :red, label = "numerical",            markershape = :circle, markevery = 1        )        scalarplot!(            p[1, 1], grid, map(x -> barenblatt(x, time, problem_data.m), grid); clear = false,            color = :green,            label = "exact", markershape = :none        )        reveal(p)        sleep(1.0e-2)    end    return sum(tsol.u[end])endusing Testfunction runtests()    testval = 46.66666666647518    return @test main(; unknown_storage = :sparse, assembly = :edgewise) ≈ testval &&        main(; unknown_storage = :dense, assembly = :edgewise) ≈ testval &&        main(; unknown_storage = :sparse, assembly = :cellwise) ≈ testval &&        main(; unknown_storage = :dense, assembly = :cellwise) ≈ testval    nothingendend

This page was generated using Literate.jl.