Adaptive multistep numerical solvers with Grassmann element assembly
Author chakravala
Popularity
10 Stars
Updated Last
1 Year Ago
Started In
November 2019

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)
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)
A = assemblestiffness(t,ϵ,m,g)
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
ϵ = 1.0
while ϵ > 5e-5 && length(t) < 10000
m = volumes(t)
A,M,b = assemble(t,c,a,f,m,h)
ξ = solvedirichlet(A+M,b,e)
η = jumps(t,c,a,f,ξ,m,h)
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