Adapode.jl

Adaptive multistep numerical ODE solver with Grassmann element assembly
Author chakravala
Popularity
3 Stars
Updated Last
20 Days Ago
Started In
November 2019

DirectSum.jl

Adapode.jl

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

DOI Docs Stable Docs Dev Gitter Build Status Build status Coverage Status codecov.io

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

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 = detsimplex(t)
    b = assemblefunction(t,f,m)
    A = assemblestiffness(t,c,m)
    R,r = assemblerobin(e,κ,gD,gN)
    return (A+R)\(b+r)
end
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 = detsimplex(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)
        scene = mesh(t,color=ξ,shading=false); display(scene)
        if typeof(g)<:AbstractRange
            scatter!(p,ξ,markersize=0.01)
        else
            wireframe!(scene[end][1],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 is compatible with discrete differential geometry.

Used By Packages

No packages found.