OrdinaryDiffEq.jl changing mass matrix

begin
    import Pkg as _Pkg
    haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"])
    using Test	
    using Revise
    using Printf
    using VoronoiFVM
    using SciMLBase: NoInit
    using OrdinaryDiffEqBDF
    using OrdinaryDiffEqRosenbrock
    using OrdinaryDiffEqSDIRK
    using LinearAlgebra
    using PlutoUI, HypertextLiteral,UUIDs
    using DataStructures
    using GridVisualize,CairoMakie
    CairoMakie.activate!(type="svg")
end

1D Nonlinear Storage

TableOfContents(aside=false)

This equation comes from the transformation of the nonlinear diffusion equation

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

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

barenblatt (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-2
end

0.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=VoronoiFVM.Grid(X)
end

ExtendableGrids.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]
    end

flux! (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)
    end

storage! (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 = 4.27, tasm = 2.62, tlinsolve = 1.63, steps = 731, iters = 2920, maxabsnorm = 1.73e-12, maxrelnorm = 1.75e-11, maxroundoff = 9.06e-15, iters_per_step = 4.0, facts_per_step = 0.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]
end

dae_storage! (generic function with 1 method)

function dae_reaction!(y,u,node,data)
    y[2]= u[2]^m-u[1]
end

dae_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 = 2.92, tasm = 2.41, tlinsolve = 0.504, steps = 732, iters = 2205, maxabsnorm = 9.72e-11, maxrelnorm = 9.82e-10, maxroundoff = 6.99e-13, iters_per_step = 3.02, facts_per_step = 0.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=0.0019999

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

Test Passed

@test norm(dae_tsol[1,:,end]-exact,Inf)<0.1

Test Passed

 @test norm(de_tsol[1,:,end]-exact,Inf)<0.05

Test Passed

myaside (generic function with 1 method)

function plotsolutions()
    vis=GridVisualizer(resolution=(380,200),dim=1,Plotter=CairoMakie,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")
    reveal(vis)
end

plotsolutions (generic function with 1 method)

Built with Julia 1.11.1 and