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