Adapode.jl

Adaptive multistep numerical ODE solver based on Grassmann.jl element assembly

This Julia project originally started as a FORTRAN 95 project called adapode.

using Grassmann, Adapode, Makie
basis"4"; x0 = 10.0v2+10.0v3+10.0v4
Lorenz(x::Chain{V}) where V = Chain{V,1}(
	1.0,
	10.0(x[3]-x[2]),
	x[2]*(28.0-x[4])-x[3],
	x[2]*x[3]-(8/3)*x[4])
lines(Point.((V(2,3,4)).(odesolve(Lorenz,x0))))

Supported ODE solvers include: explicit Euler, Heun's method (improved Euler), Midpoint 2nd order RK, Kutta's 3rd order RK, classical 4th order Runge-Kuta, adaptive Heun-Euler, adaptive Bogacki-Shampine RK23, adaptive Fehlberg RK45, adaptive Cash-Karp RK45, adaptive Dormand-Prince RK45, multistep Adams-Bashforth-Moulton 2nd,3rd,4th,5th order, adaptive multistep ABM 2nd,3rd,4th,5th order.

It is possible to work with L2 projection on a mesh with

L2Projector(t,f;args...) = mesh(t,color=\(assemblemassfunction(t,f)...);args...)
L2Projector(initmesh(0:1/5:1)[3],x->x[2]*sin(x[2]))
L2Projector(initmesh(0:1/5:1)[3],x->2x[2]*sin(2π*x[2])+3)

Partial differential equations can also be assembled with various additional methods:

PoissonSolver(p,e,t,c,f,κ,gD=1,gN=0) = mesh(t,color=solvepoisson(t,e,c,f,κ,gD,gN))
function solvepoisson(t,e,c,f,κ,gD=0,gN=0)
    m = volumes(t)
    b = assemblefunction(t,f,m)
    A = assemblestiffness(t,c,m)
    R,r = assemblerobin(e,κ,gD,gN)
    return (A+R)\(b+r)
end
function solveSD(t,e,c,f,δ,κ,gD=0,gN=0)
    m = volumes(t)
    g = gradienthat(t,m)
    A = assemblestiffness(t,c,m,g)
    b = means(t,f)
    C = assembleconvection(t,b,m,g)
    Sd = assembleSD(t,sqrt(δ)*b,m,g)
    R,r = assemblerobin(e,κ,gD,gN)
    return (A+R-C'+Sd)\r
end
function solvetransport(t,e,c,ϵ=0.1)
    m = volumes(t)
    g = gradienthat(t,m)
    A = assemblestiffness(t,ϵ,m,g)
    b = assembleload(t,m)
    C = assembleconvection(t,c,m,g)
    return solvedirichlet(A+C,b,e)
end

Such modular methods can work on input meshes of any dimension. The following examples are based on trivially generated 1 dimensional domains:

function BackwardEulerHeat1D()
    x,m = 0:0.01:1,100; p,e,t = initmesh(x)
    T = range(0,0.5,length=m+1) # time grid
    ξ = 0.5.-abs.(0.5.-x) # initial condition
    A = assemblestiffness(t) # assemble(t,1,2x)
    M,b = assemblemassfunction(t,2x).+assemblerobin(e,1e6,0,0)
    h = Float64(T.step); LHS = M+h*A # time step
    for l ∈ 1:m
        ξ = LHS\(M*ξ+h*b); l%10==0 && println(l*h)
    end
    mesh(t,color=ξ)
end
function PoissonAdaptive(g,p,e,t,c=1,a=0,f=1)
    ϵ = 1.0
    while ϵ > 5e-5 && length(t) < 10000
        m = volumes(t)
        h = gradienthat(t,m)
        A,M,b = assemble(t,c,a,f,m,h)
        ξ = solvedirichlet(A+M,b,e)
        η = jumps(t,c,a,f,ξ,m,h)
        display(mesh(t,color=ξ,shading=false))
        if typeof(g)<:AbstractRange
            scatter!(p,ξ,markersize=0.01)
        else
            wireframe!(t,color=(:red,0.6),linewidth=3)
        end
        ϵ = sqrt(norm(η)^2/length(η))
        println(t,", ϵ=$ϵ, α=$(ϵ/maximum(η))"); sleep(0.5)
        refinemesh!(g,p,e,t,select(η,ϵ),"regular")
    end
    return g,p,e,t
end
PoissonAdaptive(refinemesh(0:0.25:1)...,1,0,x->exp(-100abs2(x[2]-0.5)))

More general problems for finite element boundary value problems are also enabled by mesh representations imported from external sources. These methods can automatically generalize to higher dimensional manifolds and are compatible with discrete differential geometry.