Adapode.jl

Adaptive P/ODE numerics with Grassmann.jl element TensorField assembly

This project originally started as a FORTRAN 95 project called adapode and evolved with Grassmann.jl and Cartan.jl.

Adapode.jl is an extension to the Cartan.jl package with support for solving ordinary differential equations and partial differential equations. It has various types of time-stepping methods (including multistep and adaptive), spectral element methods, and finite element methods.

Cartan.jl introduces a pioneering unified numerical framework for comprehensive differential geometric algebra, purpose-built for the formulation and solution of partial differential equations on manifolds with non-trivial topological structure and Grassmann.jl algebra. Written in Julia, Cartan.jl unifies differential geometry, geometric algebra, and tensor calculus with support for fiber product topology; enabling directly executable generalized treatment of geometric PDEs over grids, meshes, and simplicial decompositions.

The system supports intrinsic formulations of differential operators (including the exterior derivative, codifferential, Lie derivative, interior product, and Hodge star) using a coordinate-free algebraic language grounded in Grassmann-Cartan multivector theory. Its core architecture accomodates numerical representations of fiber bundles, product manifolds, and submanifold immersion, providing native support for PDE models defined on structured or unstructured domains.

Cartan.jl integrates naturally with simplex-based finite element exterior calculus, allowing for geometrical discretizations of field theories and conservation laws. With its synthesis of symbolic abstraction and numerical execution, Cartan.jl empowers researchers to develop PDE models that are simultaneously founded in differential geometry, algebraically consistent, and computationally expressive, opening new directions for scientific computing at the interface of geometry, algebra, and analysis.

   _______     __                          __
  |   _   |.--|  |.---.-..-----..-----..--|  |.-----.
  |       ||  _  ||  _  ||  _  ||  _  ||  _  ||  -__|
  |___|___||_____||___._||   __||_____||_____||_____|
                         |__|

developed by chakravala with Grassmann.jl and Cartan.jl

For GridBundle initialization it is typical to invoke a combination of ProductSpace and QuotientTopology, while optional Julia packages extend SimplexBundle initialization, such as Meshes.jl, GeometryBasics.jl, Delaunay.jl, QHull.jl, MiniQhull.jl, Triangulate.jl, TetGen.jl, MATLAB.jl, FlowGeometry.jl.

Ordinary differential equations

Example (Lorenz). Observe vector fields by integrating streamlines

using Grassmann, Cartan, Adapode, Makie # GLMakie
Lorenz(s,r,b) = x -> Chain(
    s*(x[2]-x[1]), x[1]*(r-x[3])-x[2], x[1]*x[2]-b*x[3])
p = TensorField(ProductSpace(-40:0.2:40,-40:0.2:40,10:0.2:90))
vf = Lorenz(10.0,60.0,8/3).(p) # pick Lorenz parameters, apply
streamplot(vf,gridsize=(10,10)) # visualize vector field

ODE solvers in the Adapode package are built on Cartan, providing both Runge-Kutta and multistep methods with optional adaptive time stepping.

fn,x0 = Lorenz(10.0,28.0,8/3),Chain(10.0,10.0,10.0)
ic = InitialCondition(fn,x0,2pi) # tmax = 2pi
lines(odesolve(ic,MultistepIntegrator{4}(2^-15)))

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.

When there is a Levi-Civita connection with zero torsion related to a metrictensor, then there exist Christoffel symbols of the secondkind. In particular, these can be expressed in terms of the metrictensor g to express local geodesic differential equations for Riemannian geometry.

Example. using Grassmann, Cartan, Adapode, Makie # GLMakie

torus(x) = Chain(
    (2+0.5cos(x[1]))*cos(x[2]),
    (2+0.5cos(x[1]))*sin(x[2]),
    0.5sin(x[1]))
tor = torus.(TorusParameter(60,60))
tormet = surfacemetric(tor) # intrinsic metric
torcoef = secondkind(tormet) # Christoffel symbols
ic = geodesic(torcoef,Chain(1.0,1.0),Chain(1.0,sqrt(2)),10pi)
sol = geosolve(ic,ExplicitIntegrator{4}(2^-7)) # Runge-Kutta
lines(torus.(sol))

totalarclength(sol) # apparent length of parameter path
@basis MetricTensor([1 1; 1 1]) # abstract non-Euclidean V
solm = TensorField(tormet(sol),Chain{V}.(value.(fiber(sol))))
totalarclength(solm) # 2D estimate totalarclength(torus.(sol))
totalarclength(torus.(sol)) # compute in 3D Euclidean metric
lines(solm) # parametric curve can have non-Euclidean metric
lines(arclength(solm)); lines!(arclength(sol))

Example (Klein geodesic). General ImmersedTopology are supported

klein(x) = klein(x[1],x[2]/2)
function klein(v,u)
    x = cos(u)*(-2/15)*(3cos(v)-30sin(u)+90sin(u)*cos(u)^4-
        60sin(u)*cos(u)^6+5cos(u)*cos(v)*sin(u))
    y = sin(u)*(-1/15)*(3cos(v)-3cos(v)*cos(u)^2-
        48cos(v)*cos(u)^4+48cos(v)*cos(u)^6-
        60sin(u)+5cos(u)*cos(v)*sin(u)-
        5cos(v)*sin(u)*cos(u)^3-80cos(v)*sin(u)*cos(u)^5+
        80cos(v)*sin(u)*cos(u)^7)
    z = sin(v)*(2/15)*(3+5cos(u)*sin(u))
    Chain(x,y,z)
end # too long paths over QuotientTopology can stack overflow
kle = klein.(KleinParameter(100,100))
klecoef = secondkind(surfacemetric(kle))
ic = geodesic(klecoef,Chain(1.0,1.0),Chain(1.0,2.0),2pi)
lines(geosolve(ic,ExplicitIntegrator{4}(2^-7)));wireframe(kle)

Example (Upper half plane). Intrinsic hyperbolic Lobachevsky metric

halfplane(x) = TensorOperator(Chain(
    Chain(Chain(0.0,inv(x[2])),Chain(-inv(x[2]),0.0)),
    Chain(Chain(-inv(x[2]),0.0),Chain(0.0,-inv(x[2])))))
z1 = geosolve(halfplane,Chain(1.0,1.0),Chain(1.0,2.0),10pi,7)
z2 = geosolve(halfplane,Chain(1.0,0.1),Chain(1.0,2.0),10pi,7)
z3 = geosolve(halfplane,Chain(1.0,0.5),Chain(1.0,2.0),10pi,7)
z4 = geosolve(halfplane,Chain(1.0,1.0),Chain(1.0,1.0),10pi,7)
z5 = geosolve(halfplane,Chain(1.0,1.0),Chain(1.0,1.5),10pi,7)
lines(z1); lines!(z2); lines!(z3); lines!(z4); lines!(z5)

Example (da Rios). The Cartan abstractions enable easily integrating

start(x) = Chain(cos(x),sin(x),cos(1.5x)*sin(1.5x)/5)
x1 = start.(TorusParameter(180));
darios(t,dt=tangent(fiber(t))) = hodge(wedge(dt,tangent(dt)))
sol = odesolve(darios,x1,1.0,2^-11)
mesh(sol,normalnorm)

Other kinds of solvers include the LeapIntegrator applied on a LeapCondition specification of first or second order.

fprime(v) = dirichlet!(laplacian_chebyshevfft(v))
function leapfrog(u0,N=25,tmax=1,tplot=1/3)
    dt = 6/(N-1)^2
    lc = LeapCondition(fprime,u0,u0,dt,tmax)
    odesolve(lc,LeapIntegrator{2}(Int(round(tplot/dt))))
end

Example with initial LeapCondition at rest:

x = Chebyshev(41)
XY = TensorField(ProductSpace{2}(x,x))
X,Y = split(XY)
ex = exp(-40((X-0.4)^2+Y^2))
lf = leapfrog(ex,41,4,1/30)
contour(lf,alpha=0.03) # previous page figure
surface(lf[:,:,10]) # 10 or 20, 30, 60, 70, 80

This same method generalizes to higher dimensions also:

x = Chebyshev(31)
XYZ = TensorField(ProductSpace{3}(x,x,x))
X,Y,Z = split(XYZ)
ex = exp(-40((X-0.4)^2+Y^2+(Z-0.7)^2))
lf = leapfrog(ex,31,2,1/30)
contour(lf[:,:,:,10],alpha=0.02,levels=5)

Spectral element methods

Based on the Cartan interoperability with multidimensional FFTW, there are several PDE solvers based on Dirichlet, Neumann, or periodic boundary conditions for the wave/heat/biharmonic/Riesz/Laplace differential operator. A comprehensive list of such methods has been fully implemented, here are a few fundamental examples:

x =

Adapode.jl

Install / Use

README

Adapode.jl

Ordinary differential equations

Spectral element methods