107: 1D Nonlinear Storage

(source code)

This equation comes from the transformation of the nonlinear diffuision equation.

\[\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 previous example.

module Example107_NonlinearStorage1Dusing Printfusing VoronoiFVMusing ExtendableGridsusing GridVisualize# Problem data structure to avoid global variablesmutable struct ProblemData    m::Float64    # Nonlinearity parameter    ϵ::Float64    # Regularization 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.0, Plotter = nothing, verbose = false,        unknown_storage = :sparse, tend = 0.01, tstep = 0.0001, assembly = :edgewise    )    # Create problem data structure    problem_data = ProblemData(m, 1.0e-10)    # 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] - u[1, 2]        return nothing    end    # Storage term    # This needs to be regularized as its derivative    # at 0 is infinity    function storage!(f, u, node, data)        f[1] = (data.ϵ + u[1])^(1.0 / data.m)        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)    solution = unknowns(sys)    t0 = 0.001    # Broadcast the initial value    inival[1, :] .= map(x -> barenblatt(x, t0, problem_data.m)^problem_data.m, X)    # Create solver control info    control = VoronoiFVM.SolverControl()    control.verbose = verbose    control.Δu_opt = 0.1    control.force_first_step = true    tsol = solve(sys; inival, times = [t0, tend], control)    if Plotter != nothing        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"            )            scalarplot!(                p[1, 1], grid, map(x -> barenblatt(x, time, problem_data.m)^problem_data.m, grid); clear = false,                color = :green, label = "exact"            )            reveal(p)            sleep(1.0e-2)        end    end    return sum(tsol.u[end])endusing Testfunction runtests()    testval = 174.72418935404414    @test main(; unknown_storage = :sparse, assembly = :edgewise) ≈ testval rtol = 1.0e-5    @test main(; unknown_storage = :dense, assembly = :edgewise) ≈ testval rtol = 1.0e-5    @test main(; unknown_storage = :sparse, assembly = :cellwise) ≈ testval rtol = 1.0e-5    @test main(; unknown_storage = :dense, assembly = :cellwise) ≈ testval rtol = 1.0e-5    return nothingendend

This page was generated using Literate.jl.