Global mass conservation and diffusion equation

[JanSuntajs]

Hello, I am following the code in Example 106, since I am interested in implementing the nonlinear diffusion case in 1D for some benchmarking purposes. However, upon running the code, I've noticed the mass is not conserved globally in the example.

If I understand the tutorials and Github issues correctly, if no boundary conditions are specified, the default are the homogenous Neumann boundary conditions and so the global mass of the system should remain conserved at all times. Already slightly increasing the final time of the simulation in the example and summing the solution vectors shows this is not the case. I am appending a slightly modified example below in which I changed the initial condition (it is a step function now) and am also simulating the standard linear diffusion equation for simplicity.

module CheckMassConservation
using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize



function step(x, a, b, n)
    return (a < x < b) ? 1.0 : 0.0 / n
end

function main(; n = 200, m = 2, Plotter = nothing, verbose = false,
              unknown_storage = :sparse, tend = 1., tstep = 0.001, assembly = :edgewise)

    # 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)
        f[1] = u[1, 1] - u[1, 2]
    end

    # Storage term
    function storage!(f, u, node)
        f[1] = u[1]
    end

    # Create a physics structure
    physics = VoronoiFVM.Physics(; flux = flux!,
                                 storage = storage!)

    # 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])

    # boundary_neumann!(sys, 1, 1, 0.)
    # boundary_neumann!(sys, 1, 2, 0.)
    # Create a solution array
    inival = unknowns(sys)
    t0 = 0.001

    # Broadcast the initial value
    # inival[1, :] .= map(x -> step(x, ), X)
    inival[1, :] .= map(x -> step(x, -0.5, 0.5, length(X)), X)

    # Create solver control info for constant time step size
    control = VoronoiFVM.NewtonControl()
    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)

    return grid, tsol, sum.(tsol.u)
end

using Test

grid, tsol, sumsol = main()
println(sumsol)
end

Given that the sums of the solution vector increase from $\approx 462$ to $\approx 462$ in a rather short time span, something is clearly not right; am I implementing the boundary conditions incorrectly (the same also happens if I explicitly specify them using boundary_neumann!(...), is there an error in post-processing or perhaps I should impose stricter solver tolerances?

Apart from this, is it possible there is a bug in the code that the tests do not catch?

Regards, Jan

EDIT:

After some more tinkering and trying out different initial profiles, it looks as if I choose an initial profile that was too difficult for the solver to handle in the initial stages of relaxation. I will report back once I am convinced about that.